shallowWater2DC : conservative shallow water with new wetting-drying

algorithm

benchmarks/wetting-drying/balzano1/balzano.geo

+lc = 600;
+Point(1) = {0, -3600, 0, lc};
+Point(2) = {0, 3600, 0, lc};
+Point(3) = {13800, 3600, 0, lc};
+Point(4) = {13800, -3600, 0, lc};
+Point(5) = {3600, -3600, 0, lc};
+Point(6) = {3600, 3600, 0, lc};
+Point(7) = {4800, -3600, 0, lc};
+Point(8) = {4800, 3600, 0, lc};
+Point(9) = {6000, -3600, 0, lc};
+Point(10) = {6000, 3600, 0, lc};
+Line(4) = {2, 6};
+Line(5) = {6, 8};
+Line(6) = {8, 10};
+Line(7) = {10, 3};
+Line(8) = {3, 4};
+Line(9) = {4, 9};
+Line(10) = {9, 7};
+Line(11) = {7, 5};
+Line(12) = {5, 1};
+Line(13) = {1, 2};
+Line Loop(14) = {4, 5, 6, 7, 8, 9, 10, 11, 12, 13};
+Plane Surface(15) = {14};
+Line(16) = {6, 5};
+Line(17) = {7, 8};
+Line(18) = {10, 9};
+Line {16,17,18} In Surface {15};
+Physical Line("Wall") = {4,5,6,7,9,10,11,12,13};
+Physical Line("OpenSea") = {8};
+Physical Surface("Domain") = {15};

benchmarks/wetting-drying/balzano1/balzano1.py

+import slimPre
+import slim
+import numpy
+
+#########################
+### global parameters ###
+#########################
+
+meshfile = "balzano.msh"
+tend = 4*12*3600
+
+dt = 900
+
+def ssebnd(t) :
+    t = t+5*3600
+    return numpy.sin(numpy.pi*2*t/(12*3600))*2
+
+def bathf(x,y) :
+    return 5*(x/13800)
+
+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")
+
+### open boundary ###
+time_vector = numpy.arange(0,tend+2*dt,dt,numpy.float64)
+time = slimPre.Time(time_vector)
+eta = numpy.zeros_like(time_vector)
+for i in range(eta.shape[0]) :
+    eta[i] = ssebnd(time_vector[i])
+slimPre.write_file('eta.nc', region=None, time=time, data=[('eta', eta)])
+
+### initial condition ###
+zeros = numpy.zeros_like(xyz)
+etaInit = numpy.ones((xyz.shape[0],1))*ssebnd(0)
+etaInit = numpy.maximum(-bath,etaInit)
+slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
+
+###########
+### run ###
+###########
+
+domain = slim.Domain(meshfile,('bath.nc','bath'))
+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_open('OpenSea', sse=('eta.nc','eta'))
+eq.set_boundary_coast('Wall')
+loop = slim.Loop(maximum_time_step=dt, export_time=3600)
+loop.add_equation(eq)
+loop.run()

benchmarks/wetting-drying/balzano1/view.idx

+Merge "balzano.msh";
+Merge "output/bath/bath.idx";
+Merge "output/sw2d/eta/eta.idx";
+Merge "output/sw2d/uv/uv.idx";
+Merge "balzano1.opt";

benchmarks/wetting-drying/balzano3/balzano3.py

+import slimPre
+import slim
+import numpy
+
+#########################
+### global parameters ###
+#########################
+
+meshfile = "balzano.msh"
+tend = 10*12*3600
+
+dt = 900
+
+def ssebnd(t) :
+    if t < 6*3600 : return 3*numpy.sin(numpy.pi*(2*t/(12*3600)+0.5))
+    else : return -3 
+
+def bathf(x,y) :
+    if x <= 3600 or x >= 6000 : return x/2760
+    if x <= 4800 : return -x/2760 + 60/23
+    else : return x/920 - 100/23
+
+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")
+
+### open boundary ###
+time_vector = numpy.arange(0,tend+2*dt,dt,numpy.float64)
+time = slimPre.Time(time_vector)
+eta = numpy.zeros_like(time_vector)
+for i in range(eta.shape[0]) :
+    eta[i] = ssebnd(time_vector[i])
+slimPre.write_file('eta.nc', region=None, time=time, data=[('eta', eta)])
+
+### initial condition ###
+zeros = numpy.zeros_like(xyz)
+etaInit = numpy.ones((xyz.shape[0],1))*ssebnd(0)
+etaInit = numpy.maximum(-bath,etaInit)
+slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
+
+###########
+### run ###
+###########
+
+domain = slim.Domain(meshfile,('bath.nc','bath'))
+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_open('OpenSea', sse=('eta.nc','eta'))
+eq.set_boundary_coast('Wall')
+loop = slim.Loop(maximum_time_step=dt, export_time=3600)
+loop.add_equation(eq)
+loop.run()

benchmarks/wetting-drying/balzano3/balzanoN.geo

+lc = 600;
+Point(1) = {0, -3600, 0, lc};
+Point(2) = {0, 3600, 0, lc};
+Point(3) = {13800, 3600, 0, lc};
+Point(4) = {13800, -3600, 0, lc};
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+Line Loop(1) = {1,2,3,4};
+Plane Surface(1) = {1};
+Physical Line("Wall") = {1,2,4};
+Physical Line("OpenSea") = {3};
+Physical Surface("Domain") = {1};

benchmarks/wetting-drying/balzano3/view.idx

+Merge "balzano.msh";
+Merge "output/bath/bath.idx";
+Merge "output/sw2d/eta/eta.idx";
+Merge "output/sw2d/uv/uv.idx";
+Merge "balzano3.opt";

benchmarks/wetting-drying/rough/balzano.geo

+lc = 600;
+Point(1) = {0, -3600, 0, lc};
+Point(2) = {0, 3600, 0, lc};
+Point(3) = {13800, 3600, 0, lc};
+Point(4) = {13800, -3600, 0, lc};
+Point(5) = {3600, -3600, 0, lc};
+Point(6) = {3600, 3600, 0, lc};
+Point(7) = {4800, -3600, 0, lc};
+Point(8) = {4800, 3600, 0, lc};
+Point(9) = {6000, -3600, 0, lc};
+Point(10) = {6000, 3600, 0, lc};
+Line(4) = {2, 6};
+Line(5) = {6, 8};
+Line(6) = {8, 10};
+Line(7) = {10, 3};
+Line(8) = {3, 4};
+Line(9) = {4, 9};
+Line(10) = {9, 7};
+Line(11) = {7, 5};
+Line(12) = {5, 1};
+Line(13) = {1, 2};
+Line Loop(14) = {4, 5, 6, 7, 8, 9, 10, 11, 12, 13};
+Plane Surface(15) = {14};
+Line(16) = {6, 5};
+Line(17) = {7, 8};
+Line(18) = {10, 9};
+Line {16,17,18} In Surface {15};
+Physical Line("Wall") = {4,5,6,7,9,10,11,12,13};
+Physical Line("OpenSea") = {8};
+Physical Surface("Domain") = {15};

benchmarks/wetting-drying/rough/rough.opt

+General.BoundingBoxSize = 16970.56274847714;
+General.ColorScheme = 3;
+General.MaxX = 6000;
+General.MaxY = 6000;
+General.MinX = -6000;
+General.MinY = -6000;
+General.MouseSelection = 0;
+General.RotationX = 286.4702775693401;
+General.RotationY = 1.229464019360118;
+General.RotationZ = 358.7450137902272;
+General.TrackballQuaternion0 = 0.5985561937506934;
+General.TrackballQuaternion1 = -0.002039873344194671;
+General.TrackballQuaternion2 = 0.01519401212968505;
+General.TrackballQuaternion3 = 0.8009341195347022;
+General.TranslationX = -42.57475924987375;
+General.TranslationY = -2043.58844399392;
+View[0].RaiseZ = -250;
+View[1].RaiseZ = 250;

benchmarks/wetting-drying/rough/rough.py

+import slimPre
+import slim
+import numpy
+import scipy.interpolate
+
+#########################
+### global parameters ###
+#########################
+
+meshfile = "balzano.msh"
+tend = 10*12*3600
+
+dt = 900
+
+def ssebnd(t) :
+    h = 21-t/(3600)
+    return max(h,3)
+
+N=10
+pts = numpy.ndarray((N+4,2))
+vals = numpy.zeros((N+4,))
+L = 13800
+l = 3600
+pts[:4,:] = [[0,-l],[0,l],[L,-l],[L,l]]
+vals[0:2] = -20
+vals[2:4]=0
+numpy.random.seed(0)
+pts[4:,0] = numpy.random.rand(N)*0.79*L
+pts[4:,1] = numpy.random.rand(N)*2*l-l
+vals[4:] = numpy.random.rand(N)*-20+1
+interp = scipy.interpolate.LinearNDInterpolator(pts,vals,0)
+def bathf(x,y) :
+    return interp(x,y)
+
+hlimt = 1
+hlimd = 0.5
+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")
+
+### open boundary ###
+time_vector = numpy.arange(0,tend+2*dt,dt,numpy.float64)
+time = slimPre.Time(time_vector)
+eta = numpy.zeros_like(time_vector)
+for i in range(eta.shape[0]) :
+    eta[i] = ssebnd(time_vector[i])
+slimPre.write_file('eta.nc', region=None, time=time, data=[('eta', eta)])
+
+### initial condition ###
+zeros = numpy.zeros_like(xyz)
+etaInit = numpy.ones((xyz.shape[0],1))*ssebnd(0)
+etaInit = numpy.maximum(-bath,etaInit)
+slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
+
+###########
+### run ###
+###########
+
+domain = slim.Domain(meshfile,('bath.nc','bath'))
+eq = slim.ShallowWater2dWD(domain, "implicit", final_time=tend, hlimt=hlimt, hlimd=hlimd)
+eq.set_initial_condition(('init.nc','eta'),('init.nc','u'),('init.nc','v'))    
+eq.set_boundary_open('OpenSea', sse=('eta.nc','eta'))
+eq.set_boundary_coast('Wall')
+loop = slim.Loop(maximum_time_step=dt, export_time=3600)
+loop.add_equation(eq)
+loop.run()

benchmarks/wetting-drying/rough/view.idx

+Merge "balzano.msh";
+Merge "output/bath/bath.idx";
+Merge "output/sw2d/eta/eta.idx";
+Merge "output/sw2d/uv/uv.idx";
+Merge "rough.opt";

benchmarks/wetting-drying/volcano/view.idx

+Merge "volcano.msh";
+Merge "output/bath/bath.idx";
+Merge "output/sw2d/eta/eta.idx";
+Merge "output/sw2d/uv/uv.idx";
+Merge "view_volcano.opt";

benchmarks/wetting-drying/volcano/view_volcano.opt

+General.BoundingBoxSize = 16970.56274847714;
+General.ColorScheme = 3;
+General.MaxX = 6000;
+General.MaxY = 6000;
+General.MinX = -6000;
+General.MinY = -6000;
+General.MouseSelection = 0;
+General.RotationX = 286.4702775693401;
+General.RotationY = 1.229464019360118;
+General.RotationZ = 358.7450137902272;
+General.TrackballQuaternion0 = 0.5985561937506934;
+General.TrackballQuaternion1 = -0.002039873344194671;
+General.TrackballQuaternion2 = 0.01519401212968505;
+General.TrackballQuaternion3 = 0.8009341195347022;
+General.TranslationX = -42.57475924987375;
+General.TranslationY = -2043.58844399392;
+View[0].RaiseZ = -250;
+View[1].RaiseZ = 250;

benchmarks/wetting-drying/volcano/volcano.geo

+L = 12000;
+lc = 1000;
+Point(1) = {-L/2,-L/2,0,lc}; 
+Point(2) = {-L/2,L/2,0,lc}; 
+Point(3) = {L/2+L,L/2,0,lc}; 
+Point(4) = {L/2+L,-L/2,0,lc}; 
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Line Loop(1) = {1,2,3,4};
+Plane Surface(1) = {1};
+
+Physical Line("Wall") = {1,2,4};
+Physical Line("OpenSea") = {3};
+
+Physical Surface("Domain") = {1};

benchmarks/wetting-drying/volcano/volcano.py

+import slimPre
+import slim
+import numpy
+
+#########################
+### global parameters ###
+#########################
+
+meshfile = "volcano.msh"
+dt = 900
+tend = 7*24*3600
+
+def ssebnd(t) :
+    ti = min(t,4.5*24*3600)
+    return 25-10*numpy.sin(numpy.pi*(ti/(24*3600)))
+
+def bathf(x,y) :
+    r = (x*x+y*y)**0.5
+    if r < 2500 : return -r/100
+    if r < 5000 : return -(5000-r)/100
+    return 0
+
+hlimt = 1
+hlimd = 1
+linDrag = 0.05
+
+###################
+### 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")
+
+### open boundary ###
+time_vector = numpy.arange(0,tend+dt,dt,numpy.float64)
+time = slimPre.Time(time_vector)
+eta = numpy.zeros_like(time_vector)
+for i in range(eta.shape[0]) :
+    eta[i] = ssebnd(time_vector[i])
+slimPre.write_file('eta.nc', region=None, time=time, data=[('eta', eta)])
+
+### initial condition ###
+zeros = numpy.zeros_like(xyz)
+etaInit = numpy.ones((xyz.shape[0],1))*ssebnd(0)
+etaInit = numpy.maximum(-bath+0.01,etaInit)
+slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
+
+###########
+### run ###
+###########
+
+domain = slim.Domain(meshfile,('bath.nc','bath'))
+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_open('OpenSea', sse=('eta.nc','eta'))
+eq.set_boundary_coast('Wall')
+loop = slim.Loop(maximum_time_step=dt, export_time=3600)
+loop.add_equation(eq)
+loop.run()

dgpy/scripts/slim.py

@@ -95,6 +95,265 @@ class Domain:
             self._g = slim_private._load_function(g,self._groups)
         self._density = density
  
+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): 
+        """keyword arguments:
+        
+        * domain                     
+            object of class Domain
+        * temporal_scheme 
+            scheme for the temporal integration
+            
+            * "explicit" 
+                explicit Euler scheme
+            * "implicit" 
+                implicit order 2 Runge-Kutta
+            * "multirate"
+                multirate scheme (see Seny et al. 2012, 2014)
+                
+                Warning: This temporal scheme is not yet implemented for tracers
+        * export_every_sub_time_step 
+            boolean stating if you want to export the hydrodynamics every sub time step for further use of offline equation (Default: False)
+        * initial_time               
+            initial time of the simulation (specific to each equation) (Default: 0)
+            
+            Format: 
+
+            - float [in seconds]
+            - date string in format "year-month-day hours:minutes:seconds" (e.g. "2000-02-15 10:05:00")
+        * final_time
+            final time of the simulation (default: 0)
+            
+            Format:
+                
+            - float [in seconds]
+            - date string in format "year-month-day hours:minutes:seconds" (e.g. "2007-06-25 22:42:17")                              
+        """     
+        fmt = '%Y-%m-%d %H:%M:%S'   
+        if initial_time is None:
+            self._initial_time = 0.
+        elif slim_private._is_string(initial_time):
+            date = slim_private.datetime.datetime.strptime(initial_time, fmt)
+            self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
+        else:
+            self._initial_time = float(initial_time)
+        if final_time is None:
+            self._final_time = 0.
+        elif slim_private._is_string(final_time):
+            date = slim_private.datetime.datetime.strptime(final_time, fmt)
+            self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
+        else:
+            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._wetting_drying = None
+        self._solution = dgpy.dgDofContainer(self._domain._groups, self._equation.getNbFields())
+        self._solution.setAll(0.0)
+        self._solution.setFieldName(0, 'H')
+        self._solution.setFieldName(1, 'Hu')
+        self._solution.setFieldName(2, 'Hv')
+        self._name = "sw2dC"
+        if temporal_scheme == "explicit" or temporal_scheme == "implicit" or temporal_scheme == "multirate":
+            self._temporal_scheme = temporal_scheme
+        else:
+            dgpy.Msg.Fatal("Unknown temporal scheme "+temporal_scheme+" for shallow water equation !")
+        self._export_every_sub_time_step = export_every_sub_time_step
+        self._evaluator = None
+        self._ref = []
+        self._compute_mass = False
+        self._export_on_structured_grid =False
+  
+    def set_initial_condition(self, sse, ux, uy, uz=None):
+        """
+        Initial conditions
+
+        keyword arguments:
+        
+        * sse 
+            sea surface elevation [in m] (.msh or .nc format) 
+        * ux  
+            velocity along the x-axis in the local basis or in the global basis (x,y,z) if you want to solve on a sphere [in m/s] (.idx, .msh or .nc format) 
+        * uy  
+            velocity along the y-axis in the local basis or in the global basis (x,y,z) if you want to solve on a sphere [in m/s] (.idx, .msh or .nc format) 
+        * uz  
+            velocity along the z-axis in the global basis (x,y,z) if you want to solve on a sphere [in m/s] (.idx, .msh or .nc format) (default: None)
+        """
+        _sse = _ux = _uy = _uz = None
+        if sse is not None:
+            _sse = slim_private._load_function(sse, self._domain._groups)
+        if ux is not None:
+            _ux = slim_private._load_function(ux, self._domain._groups)
+        if uy is not None:
+            _uy = slim_private._load_function(uy, self._domain._groups)
+        if uz is not None:
+            _uz = slim_private._load_function(uz, self._domain._groups)
+        velocity_local = slim_private._check_vector(_ux,_uy,_uz,self._domain)
+        def nonConservativeToConservative(cm, etaF, uvF, bathF) :
+            H = cm.get(bathF)
+            if etaF :
+                H = H+cm.get(etaF)
+            if uvF :
+                return slim_private.np.hstack([H,H*cm.get(uvF)])
+            else :
+                return slim_private.np.hstack([H,H*0,H*0])
+        self._solution.interpolate(dgpy.functorNumpy(lambda cm : nonConservativeToConservative(cm,_sse,velocity_local,self._domain._bath_function)))
+
+    def set_coriolis(self, coriolis):
+        """Add a coriolis term (2*Omega*cos(latitude)) in the shallow water equation 
+        
+        keyword argument: 
+        
+       * coriolis 
+            netcdf or .msh file containing the coriolis term for the whole domain. 
+        """
+        self._coriolis = slim_private._load_function(coriolis, self._domain._groups)
+        self._coriolis_PC = dgpy.functionPrecomputed(self._domain._groups, 3)
+        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)
+        
+        keyword arguments:
+        
+        * physical_tag 
+            tag of the part of the boundary on which this boundary condition is applied.
+        """
+        self._equation.addBoundaryCondition(physical_tag, self._equation.newBoundaryWall())
+
+    def set_boundary_open(self, physical_tag, discharge=None, sse=None, ux=None, uy=None, uz=None, transport_flux=False, ramp_period=0):
+        """Open boundary condition: one can impose the discharge or the sea surface elevation and/or the velocity (ux and uy have to be defined at the same time but they can be zero).
+        
+        keyword arguments:
+