interface3D: add wind + fix turbulence

dgpy/scripts/slim3d.py

@@ -41,6 +41,7 @@ class Slim3d_equations :
         self._initial_elevation = None
         self._initial_salinity = None
         self._initial_temperature = None
+        self._wind_stress = None
         self._boundary_coast = []
         self._boundary_open = []
 
@@ -109,6 +110,36 @@ class Slim3d_equations :
     def set_initial_temperature(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0):
         self._initial_temperature = (mode, temperature, surface_temperature, vertical_gradient)
 
+    # TODO check value of air_density
+    def set_wind_stress(self, mode, wind_x, wind_y, density_air=1.25):
+        """Add a wind stress term in the shallow water equation.
+        
+        keyword arguments:
+        
+        * mode                  
+            type of wind forcing
+            
+            * "speed"   
+                wind speed given in m s^-1. wind_stress will be computed with Smith-Banke formula.
+                  stress = C_D(speed)*density_air*speed
+                Setting the air density is mandatory 
+            * "stress" 
+                wind stress given in kg m^-1 s^-1. (air density will not be used)
+        * wind_x 
+            netcdf or .msh file containing the wind term along the x-axis in the local basis 
+            
+            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain.  
+        * wind_y 
+            netcdf or .msh file containing the wind term along the y-axis in the local basis 
+            
+            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain.  
+        * density_air  
+            density of the ambiant air (only necessary for "speed" mode) (default: 1.25 [kg m^-3])
+        """
+        if mode != "speed" and mode != "stress":
+            dgpy.Msg.Fatal("Unknown mode for wind stress: "+mode)
+        self._wind_stress = (mode, wind_x, wind_y, density_air)
+
     def set_boundary_coast(self, physical_tags):
         if slim_private._is_string(physical_tags):
             self._boundary_coast.append(physical_tags)

dgpy/scripts/slim3d_setup.py

@@ -11,7 +11,8 @@ def slim3d_setup(loop):
     if eq._vertical_viscosity == 'gotm' :
         slimSolver.setSolveTurbulence(True)
         slimSolver.setAdvectTurbulentEnergy(True)
-        slimSolver.turbulenceSetupFile = 'gotmTurb.nml'
+        # TODO download automatically or in SLIM
+        slimSolver.turbulenceSetupFile = 'gotmturb.nml'
     if (slimSolver.getSolveUVImplicitVerticalDiffusion() or slimSolver.getSolveTurbulence()) and not slimSolver.getComputeBottomFriction():
         dgpy.Msg.Fatal('If Vertical implicit diffusion or using gotm, bottom friction must be true')
     if eq._horizontal_viscosity == 'smagorinsky' :
@@ -20,10 +21,22 @@ def slim3d_setup(loop):
     d = slimSolver.getDofs(True)
     f = slimSolver.functions
     f.initializeBath()
-
-    ### TODO
-    #f.windFunc = 
-    #f.windStressFunc =
+    
+    def mergeUVNumpy(cmap, val, u, v) :
+        val[:,0] = u[:] 
+        val[:,1] = v[:]  
+    
+    if eq._wind_stress:
+        if eq._wind_stress[0] == 'speed':
+            eq._wind_speed_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
+            eq._wind_speed_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d) 
+            f.windFunc = dgpy.functionNumpy(2, mergeUVNumpy, [eq._wind_speed_x, eq._wind_speed_y])
+            f.windStressFunc = dgpy.windCubicStress(f.windFunc, eq._wind_stress[3])   
+        # TODO Fatal if no wind speed and solving sediment
+        elif eq._wind_stress[0] == 'stress':
+            eq._wind_stress_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
+            eq._wind_stress_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d)
+            f.windStressFunc = dgpy.functionNumpy(2, mergeUVNumpy, [eq._wind_stress_x, eq._wind_stress_y])
 
     if eq._initial_elevation:
         f.etaInitFunc = slim_private._load_function(eq._initial_elevation, slimSolver.groups2d)
@@ -119,10 +132,7 @@ def slim3d_setup(loop):
         val[:] = numpy.sqrt( (bath[:]+eta[:]) / g) * un
     def transport2UVNumpy(cmap, val, uv, bath, eta) :
         val[:,0] = uv[:,0] / ( bath[:] + eta[:] )  
-        val[:,1] = uv[:,1] / ( bath[:] + eta[:] )  
-    def mergeUVNumpy(cmap, val, u, v) :
-        val[:,0] = u[:] 
-        val[:,1] = v[:] 
+        val[:,1] = uv[:,1] / ( bath[:] + eta[:] )
 
     eta = f.eta2dFunc
     uv = d.uvDof.getFunction()
@@ -213,7 +223,7 @@ def slim3d_setup(loop):
     horMomEq = e.horMomEq
     horMomEq.setLaxFriedrichsFactor(0.0) 
     if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc :
-        #horMomEq.setBoundary0Flux(topTags) # CHECK WHICH TODO     !!!
+        #horMomEq.setBoundary0Flux(topTags) # TODO CHECK WHICH      !!!
         horMomEq.setBoundarySymmetry(topTags)
     else :
         horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
@@ -253,7 +263,7 @@ def slim3d_setup(loop):
     e.newRGradEq.setBoundarySymmetry(bottomTags)
     e.newRGradEq.setBoundarySymmetry(eq._boundary_coast)
     for openBnd in eq._boundary_open:
-        e.newRGradEq.setBoundarySymmetry(openBnd.tag) # MUST BE CHANGED !!!
+        e.newRGradEq.setBoundarySymmetry(openBnd.tag) # TODO MUST BE CHANGED !!!
     #rGradBndOut = e.newRGradEq.newOutsideValueBoundary("", f.rhoFunc)
     #e.newRGradEq.addBoundaryCondition(['open', 'openSouth'], rGradBndOut)
     
@@ -282,14 +292,14 @@ def slim3d_setup(loop):
             e.TEq.addBoundaryCondition(openBnd.tag, TBndOpen)
     
     if slimSolver.getAdvectTurbulentEnergy() :
-        for eq in [e.tkeAdvEq, e.epsAdvEq] :
-            eq.setLaxFriedrichsFactor(0.0)
-            eq.setBoundary0Flux(topTags)
-            eq.setBoundary0Flux(bottomTags)
-            eq.setBoundary0Flux(eq._boundary_coast)
-            eq.setBoundary0Flux('vertical_bottom')
+        for turbEq in [e.tkeAdvEq, e.epsAdvEq] :
+            turbEq.setLaxFriedrichsFactor(0.0)
+            turbEq.setBoundary0Flux(topTags)
+            turbEq.setBoundary0Flux(bottomTags)
+            turbEq.setBoundary0Flux(eq._boundary_coast)
+            turbEq.setBoundary0Flux('vertical_bottom')
             for openBnd in eq._boundary_open:
-                eq.addBoundaryCondition(openBnd.tag, eq.newInFluxBoundary(f.tinyFunc))
+                turbEq.addBoundaryCondition(openBnd.tag, turbEq.newInFluxBoundary(f.tinyFunc))
     
     eta2dEq = e.eta2dEq
     eta2dEq.setBoundary0Flux(eq._boundary_coast)