slim3d interface is updated

dgpy/scripts/slim3d.py

@@ -26,6 +26,7 @@ class Slim3d_equations :
         slimSolver.setSolveUVImplicitVerticalDiffusion(False)
         slimSolver.setSolveSImplicitVerticalDiffusion(False)
         slimSolver.setSolveTImplicitVerticalDiffusion(False)
+        self._verticalDiffImplicit = False
         slimSolver.setFlagUVLimiter(True)
         if salinity:
             slimSolver.setFlagSLimiter(True)
@@ -51,11 +52,14 @@ class Slim3d_equations :
 
     def set_implicit_vertical_diffusion(self, flag):
         slimSolver = self._slimSolver
-        slimSolver.setSolveUVImplicitVerticalDiffusion(flag)
-        if slimSolver.getSolveS():
-            slimSolver.setSolveSImplicitVerticalDiffusion(flag)
-        if slimSolver.getSolveT():
-            slimSolver.setSolveTImplicitVerticalDiffusion(flag)
+        self._verticalDiffImplicit = flag
+        if flag:
+            if self._vertical_viscosity:
+                slimSolver.setSolveUVImplicitVerticalDiffusion(True)
+            if slimSolver.getSolveS() and self._vertical_diffusivity:
+                slimSolver.setSolveSImplicitVerticalDiffusion(True)
+            if slimSolver.getSolveT() and self._vertical_diffusivity:
+                slimSolver.setSolveTImplicitVerticalDiffusion(True)
 
     def set_limiter(self, flag):
         slimSolver = self._slimSolver
@@ -81,17 +85,26 @@ class Slim3d_equations :
         self._hor_diff_fact = factor
         self._hor_diff_max = maximum
 
-    def set_vertical_viscosity(self, mode, constant_value=1e-5):
-        if (mode != "gotm") and (mode != "constant") and (mode != "PP"):
+    def set_vertical_viscosity(self, mode, constant_value=1e-5, paca_values=[1e-2, 1e-5, 10, 2]):
+        if (mode != "gotm") and (mode != "constant") and (mode != "pacanowski_philander"):
             dgpy.Msg.Fatal("Unknown viscosity for vertical_viscosity : "+mode)
         self._vertical_viscosity = mode
-        self._vertical_viscosity_value = constant_value
+        self._vertical_viscosity_constant_value = constant_value
+        self._vertical_viscosity_paca_values = paca_values
+        if self._verticalDiffImplicit:
+            self._slimSolver.setSolveUVImplicitVerticalDiffusion(True)
+
         
-    def set_vertical_diffusivity(self, mode, constant_value=1e-5):
-        if (mode != "gotm") and (mode != "constant") and (mode != "PP"):
+    def set_vertical_diffusivity(self, mode, constant_value=1e-5, paca_values=[1e-5, 10]):
+        if (mode != "gotm") and (mode != "constant") and (mode != "pacanowski_philander"):
             dgpy.Msg.Fatal("Unknown diffusivity for vertical_diffusivity : "+mode)
         self._vertical_diffusivity = mode
-        self._vertical_diffusivity_value = constant_value
+        self._vertical_diffusivity_constant_value = constant_value
+        self._vertical_diffusivity_paca_values = paca_values
+        if self._slimSolver.getSolveS() and self._vertical_diffusivity:
+            self._slimSolver.setSolveSImplicitVerticalDiffusion(True)
+        if self._slimSolver.getSolveT() and self._vertical_diffusivity:
+            self._slimSolver.setSolveTImplicitVerticalDiffusion(True)
 
     def set_gotm_option_file(self,f):
         self._gotm_option_file = f
@@ -110,9 +123,21 @@ class Slim3d_equations :
             self._z0[1] = z0S
 
     def set_linear_density(self, salinity=False, temperature=False, constant_coefficient=0, linear_coefficient=0):
-        if (temperature and salinity) or (not temperature and not salinity):
-            dgpy.Msg.Fatal('density is a linear function of either salinity or temperature, not both, not none') 
-        self._linear_density = (salinity, constant_coefficient, linear_coefficient)
+        """Set the density as a linear function of either salinity, temperature or z coordinate. If this function is not set, density follows Jackett-McDougall
+        
+        keyword arguments:
+        
+        * mode                  
+            "salinity", "temperature" or "z_coordinate"
+        * constant_coefficient
+            value of the salinity, temperature or z coord for which rho = rho_0
+        * linear_coefficient
+            linear_coefficient of the function
+            
+        """
+        if mode != "salinity" and mode != "temperature" and mode != "z_coordinate":
+            dgpy.Msg.Fatal('density is a linear function of either salinity, temperature or z coordinate') 
+        self._linear_density = (mode, constant_coefficient, linear_coefficient)
 
     def set_reference_density(self, density):
         self._reference_density = density

dgpy/scripts/slim3d_private.py

@@ -79,10 +79,12 @@ def slim3d_setup(loop):
         dgpy.Msg.Fatal('Initial temperature must be set if rhoFunc is not set')
 
     if eq._linear_density :
-        if eq._linear_density[0]:
+        if eq._linear_density[0] == 'salinity':
             f.rhoFunc = dgpy.linearStateEquation(d.SDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
-        else :
+        elif eq._linear_density[0] == 'temperature' :
             f.rhoFunc = dgpy.linearStateEquation(d.TDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
+        else:
+            f.rhoFunc = dgpy.linearStateEquation(d.zDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
         if eq._vertical_diffusivity == 'gotm':
             dgpy.Msg.Fatal('GOTM does not work with linear density. S and T must be set, and the automatic Jackett-McDougall state equation will be computed')
     dgpy.slim3dParameters.setRho0(eq._reference_density)
@@ -111,13 +113,17 @@ def slim3d_setup(loop):
         f.kappahTotal = dgpy.okubo(eq._hor_diff_fact, eq._hor_diff_max, eq._areaFunc)
 
     if eq._vertical_viscosity == 'constant':
-        f.nuvFunc = dgpy.functionConstant(eq._vertical_viscosity_value)
-    elif eq._vertical_viscosity == 'PP':
-        f.nuvFunc = dgpy.nuPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), 1e-2, 1e-5, 10, 2)
+        f.nuvFunc = dgpy.functionConstant(eq._vertical_viscosity_constant_value)
+    elif eq._vertical_viscosity == 'pacanowski_philander':
+        paca_values = eq._vertical_viscosity_paca_values
+        f.nuvFunc = dgpy.nuPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), paca_values[0], paca_values[1], paca_values[2], paca_values[3])
     if eq._vertical_diffusivity == 'constant':
-        f.kappavFunc = dgpy.functionConstant(eq._vertical_diffusivity_value)
-    elif eq._vertical_diffusivity == 'PP':
-        f.kappavFunc = dgpy.kappaPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), f.nuvFunc, 1e-5, 10)
+        f.kappavFunc = dgpy.functionConstant(eq._vertical_diffusivity_constant_value)
+    elif eq._vertical_diffusivity == 'pacanowski_philander':
+        paca_values = eq._vertical_diffusivity_paca_values
+        if not f.nuvFunc:
+            dgpy.Msg.Fatal('To give a Pacanowski-Philander vertical diffusivity, a vertical viscosity is required')
+        f.kappavFunc = dgpy.kappaPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), f.nuvFunc, paca_values[0], paca_values[1])
     
     if slimSolver.getComputeBottomFriction():
         f.z0BFunc = dgpy.functionConstant(eq._z0[0])
@@ -180,9 +186,9 @@ def slim3d_setup(loop):
     for openBnd in eq._boundary_open:
         if openBnd.u and openBnd.eta:
             eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d) 
-            eq._v_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d) 
+            eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d) 
             eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d) 
-            eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d) 
+            eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d) 
             if openBnd.transport:
                 eq._uvTransport_open[openBnd]     = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])  
                 eq._uv_open[openBnd]              = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta]) 
@@ -196,9 +202,9 @@ def slim3d_setup(loop):
             eq._eta_open[openBnd]     = slim_private._load_function(openBnd.eta, slimSolver.groups2d)  
         elif openBnd.u:
             eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d) 
-            eq._v_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d) 
+            eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d) 
             eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d) 
-            eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d) 
+            eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d) 
             if openBnd.transport:
                 eq._uvTransport_open[openBnd]     = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])  
                 eq._uv_open[openBnd]              = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])