sw2dC: add option for manning coefficient

benchmarks/wetting-drying/thacker/thacker.py

-import slimPre
-import slim
-import numpy
-
-#########################
-### global parameters ###
-#########################
-
-meshfile = "thacker.msh"
-
-tend = 4*12*3600
-
-dt = 50
-
-h0 = 50
-eta0 = 2
-R = 430.62e3
-g = 9.81
-
-def ssebnd(x,y) :
-    r = numpy.sqrt(x**2 + y**2)
-    A = ((h0+eta0)**2 - h0**2)/((h0+eta0)**2 + h0**2)
-    o = numpy.sqrt(8*g*h0/(R**2)) 
-    t = numpy.pi/o
-    eta = h0*(numpy.sqrt(1-A**2)/(1-A*numpy.cos(o*t)) - 1 - (r**2/R**2)*((1-A**2)/(1-A*numpy.cos(o*t))**2 - 1) )
-    return eta
-
-def bathf(x,y) :
-    r2 = x**2 + y**2
-    return h0 * (R**2 - r2)/(R**2)
-
-hlimt = 1
-hlimd = 0.1
-linDrag = 0.01
-
-###################
-### pre process ###
-###################
-
-### mesh ###
-mesh = slimPre.Mesh(meshfile)
-region = slimPre.Region(mesh)
-xyz = region.coordinates
-
-### bathymetry ###
-bath = numpy.zeros((xyz.shape[0]))
-for i in range(xyz.shape[0]) :
-    bath[i] = bathf(xyz[i,0], xyz[i,1])
-slimPre.write_file('bath.nc', region=region, time=None, data=[('bath', bath)])
-slimPre.netcdf_to_msh(meshfile, "bath.nc", "bath", "output/bath")
-
-### init elev ###
-etaInit = numpy.zeros((xyz.shape[0]))
-for i in range(xyz.shape[0]) :
-    etaInit[i] = numpy.maximum(ssebnd(xyz[i,0],xyz[i,1]), -bath[i]+.05)
-zeros = numpy.zeros_like(etaInit)
-slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
-
-slimPre.write_file('g.nc', region=None, time=None, data=[('g', g)])
-
-###########
-### run ###
-###########
-
-domain = slim.Domain(meshfile,('bath.nc','bath'), g=("g.nc","g"))
-eq = slim.ShallowWater2dWD(domain, "implicit", final_time=tend, hlimt=hlimt, hlimd=hlimd, linDrag=linDrag)
-eq.set_initial_condition(('init.nc','eta'),('init.nc','u'),('init.nc','v'))   
-eq.set_boundary_coast('Wall')
-loop = slim.Loop(maximum_time_step=dt, export_time=600)
-loop.add_equation(eq)
-loop.run()
+import slimPre
+import slim
+import numpy
+
+#########################
+### global parameters ###
+#########################
+
+meshfile = "thacker.msh"
+
+tend = 4*12*3600
+
+dt = 50
+
+h0 = 50
+eta0 = 2
+R = 430.62e3
+g = 9.81
+
+def ssebnd(x,y) :
+    r = numpy.sqrt(x**2 + y**2)
+    A = ((h0+eta0)**2 - h0**2)/((h0+eta0)**2 + h0**2)
+    o = numpy.sqrt(8*g*h0/(R**2)) 
+    t = numpy.pi/o
+    eta = h0*(numpy.sqrt(1-A**2)/(1-A*numpy.cos(o*t)) - 1 - (r**2/R**2)*((1-A**2)/(1-A*numpy.cos(o*t))**2 - 1) )
+    return eta
+
+def bathf(x,y) :
+    r2 = x**2 + y**2
+    return h0 * (R**2 - r2)/(R**2)
+
+hlimt = 1
+hlimd = 0.1
+linDrag = 0.
+manning = 0.
+
+###################
+### pre process ###
+###################
+
+### mesh ###
+mesh = slimPre.Mesh(meshfile)
+region = slimPre.Region(mesh)
+xyz = region.coordinates
+
+### bathymetry ###
+bath = numpy.zeros((xyz.shape[0]))
+for i in range(xyz.shape[0]) :
+    bath[i] = bathf(xyz[i,0], xyz[i,1])
+slimPre.write_file('bath.nc', region=region, time=None, data=[('bath', bath)])
+slimPre.netcdf_to_msh(meshfile, "bath.nc", "bath", "output/bath")
+
+### init elev ###
+etaInit = numpy.zeros((xyz.shape[0]))
+for i in range(xyz.shape[0]) :
+    etaInit[i] = numpy.maximum(ssebnd(xyz[i,0],xyz[i,1]), -bath[i]+.05)
+zeros = numpy.zeros_like(etaInit)
+slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
+
+slimPre.write_file('g.nc', region=None, time=None, data=[('g', g)])
+
+###########
+### run ###
+###########
+
+domain = slim.Domain(meshfile,('bath.nc','bath'), g=("g.nc","g"))
+eq = slim.ShallowWater2dWD(domain, "implicit", final_time=tend, hlimt=hlimt, hlimd=hlimd, linDrag=linDrag, manning=manning)
+eq.set_initial_condition(('init.nc','eta'),('init.nc','u'),('init.nc','v'))   
+eq.set_boundary_coast('Wall')
+loop = slim.Loop(maximum_time_step=dt, export_time=600)
+loop.add_equation(eq)
+loop.run()

dgpy/scripts/slim.py

@@ -98,7 +98,7 @@ class Domain:
 class ShallowWater2dWD:
     """Create the shallow water equation with various options"""
 
-    def __init__(self, domain, temporal_scheme, export_every_sub_time_step=False, initial_time=None, final_time=None, hlimt=1, hlimd=0.1, linDrag=0.1): 
+    def __init__(self, domain, temporal_scheme, export_every_sub_time_step=False, initial_time=None, final_time=None, hlimt=1, hlimd=0.1, linDrag=0.1, manning=0.02): 
         """keyword arguments:
         
         * domain                     
@@ -148,6 +148,7 @@ class ShallowWater2dWD:
             self._final_time = float(final_time)
         self._domain = domain
         self._equation = dgpy.dgConservationLawShallowWater2dC(self._domain._bath_function, self._domain._bath_gradient_function, hlimt, hlimd, linDrag)
+        self._equation.setManningCoefficient(manning)        
         self._wetting_drying = None
         self._solution = dgpy.dgDofContainer(self._domain._groups, self._equation.getNbFields())
         self._solution.setAll(0.0)
@@ -213,7 +214,6 @@ class ShallowWater2dWD:
         self._coriolis_PC.compute(self._coriolis)
         self._equation.setCoriolisFactor(self._coriolis_PC)
     
-        
     def set_boundary_coast(self, physical_tag):
         """Boundary condition: normal speed set to zero and tangential speed set to zero or free (no slip or free slip condition)
         

modules/shallowWater/dgConservationLawShallowWater2dC.cpp

 #include "dgConservationLawShallowWater2dC.h"
 #include "functorMember.h"
 
-static double tdrag(double g, double H, double U, double V, double HLIM, double linDrag) {
+static double tdrag(double g, double H, double U, double V, double HLIM, double linDrag, double manningCoeff) {
   const double Hc = std::max(HLIM, H);
   const double a = std::max(0., (HLIM-std::max(0.,H))/HLIM);
-  return 0.0004*g*(hypot(U,V))*pow(Hc, -7./3)+a*linDrag;
+  return manningCoeff*manningCoeff*g*(hypot(U,V))*pow(Hc, -7./3)+a*linDrag;
 }
 
 static void tflux(double g, double H, double U, double V, double &fluxU, double &fluxV, double &c, double &ct, double _HLIM) {
@@ -73,7 +73,7 @@ void dgConservationLawShallowWater2dC::psiterm(functorCache &cm, fullMatrix<doub
     double H = sol(i, 0);
     double U = sol(i, 1);
     double V = sol(i, 2);
-    double drag = tdrag(_g,H,U,V,_hlimd, _linDrag);
+    double drag = tdrag(_g,H,U,V,_hlimd, _linDrag, _manningCoefficient);
     r(i, 0) = 0;
     r(i, 1) = -drag * U +_g*H*dh(i,0);
     r(i, 2) = -drag * V +_g*H*dh(i,1);
@@ -122,6 +122,7 @@ dgConservationLawShallowWater2dC::dgConservationLawShallowWater2dC(const functor
   _linDrag = linDrag;
   _g = 9.81;
   _ws = NULL;
+  _manningCoefficient = 0.02;
   _bath = &bath;
   _hlimd = hlimd;
   _hlimt = hlimt;

modules/shallowWater/dgConservationLawShallowWater2dC.h

@@ -5,7 +5,7 @@
 class dgDofContainer;
 
 class dgConservationLawShallowWater2dC : public dgConservationLawFunction {
-  double _hlimt, _hlimd, _g, _linDrag;
+  double _hlimt, _hlimd, _g, _linDrag, _manningCoefficient;
   const functor *_ws, *_coriolis, *_bath, *_bathGradient;
   public :
   void gradpsiterm(functorCache &cm, fullMatrix<double> &r) const;
@@ -13,6 +13,7 @@ class dgConservationLawShallowWater2dC : public dgConservationLawFunction {
   void interfaceterm(functorCache &cm, fullMatrix<double> &r) const;
   void setCoriolisFactor(const functor *f) {_coriolis = f;}
   void setWindStress(const functor *f) {_ws = f;}
+  void setManningCoefficient(const double d) {_manningCoefficient = d;}
   dgConservationLawShallowWater2dC(const functor &bath, const functor &bathGradient, double hlimt, double hlimd, double linDrag);
   dgBoundaryCondition *newBoundaryWall();
 };