Merge branch 'set_source3D' into 'master'

A source term can be add in the horizontal momentum equation

See merge request !92

dgpy/scripts/slim3d.py

@@ -58,6 +58,8 @@ class Slim3d_equations :
         self._nuh = None
         self._add_shear = None
         self._verticalAdaptation = {}
+        self._forcing_x = None
+        self._forcing_y = None
 
     def set_implicit_vertical(self, flag):
         self._slimSolver.setImplicitVertical(flag)
@@ -207,6 +209,21 @@ class Slim3d_equations :
             dgpy.Msg.Fatal("Unknown mode for wind stress: "+mode)
         self._wind_stress = (mode, wind_x, wind_y, density_air)
 
+    def set_forcing(self, forcing_x, forcing_y):
+        """Add a momentum source term derived from a gravitational potential in the shallow water equation.
+        
+        keyword arguments: 
+        
+        * forcing_x 
+            netcdf or .msh file containing the forcing term along the x-axis in the local basis for the whole domain [in ms^-2].  
+        * forcing_y 
+            netcdf or .msh file containing the forcing term along the y-axis in the local basis for the whole domain [in ms^-2].
+        """
+        self._forcing_x = forcing_x
+        self._forcing_y = forcing_y
+        if forcing_x is None or forcing_y is None:
+            dgpy.Msg.Fatal("You must set a forcing along the x axis and along the y axis")
+            
     def set_boundary_coast(self, physical_tags):
         if slim_private._is_string(physical_tags):
             self._boundary_coast.append(physical_tags)

dgpy/scripts/slim3d_private.py

@@ -151,7 +151,7 @@ def slim3d_setup(loop):
             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.functorNumpy(lambda cm : mergeUVNumpy(cm,eq._wind_stress_x, eq._wind_stress_y))
-
+          
     if eq._initial_elevation:
         f.etaInitFunc = slim_private._load_function(eq._initial_elevation, slimSolver.groups2d)
     else:
@@ -257,7 +257,7 @@ def slim3d_setup(loop):
         eq._coriolis_PC = dgpy.functionPrecomputedExtrusion(slimSolver.extrusion(), 3)
         eq._coriolis_PC.compute(coriolisFunc)
         f.coriolisFunc = eq._coriolis_PC
-    
+        
     if eq._atmPress:    
         f.atmPressFunc = slim_private._load_function(eq._atmPress, slimSolver.groups2d)
 
@@ -305,7 +305,13 @@ def slim3d_setup(loop):
         horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
     else :
         horMomEq.setBoundary0Flux(topTags)
-
+        
+    if eq._forcing_x: 
+        eq._forcing_x = slim_private._load_function(eq._forcing_x, slimSolver.groups3d)
+        eq._forcing_y = slim_private._load_function(eq._forcing_y, slimSolver.groups3d)
+        eq._forcing = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._forcing_x, eq._forcing_y))
+        horMomEq.setSource(eq._forcing)
+        
     for openBnd in eq._boundary_open:
         horMomEq.setBoundaryNeumann(openBnd.tag, openBnd.horizontal_momentum_flux_f)