slim.py 77.3 KB
Newer Older
1
2
3
"""Second-generation Louvain-la-Neuve Ice-ocean model - User interface

    Contact: jonathan.lambrechts@uclouvain.be
4
    webpage: www.climate.be/slim
5

6
    SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.
7

8
9
10
11
    To simulate the model, run the script my_run_script.py:
        - On a single core: rundgpy my_run_script.py
        - On multiple cores: mpirun -np N rundgpy my_run_script.py (replace N by the number of cores)

12
13
    Equations: 
        - 2D:
14
            * Shallow water equations (see :class:`~ShallowWater2d`)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
15
            * Tracer equation (see :class:`~slim.ShallowWaterTracer2d`)
16
17
            * Particle Transport (see class TODO)
        - 3D:
18
            * Shallow water equations (see class TODO)
19
20
21
            * Tracer equation (see class TODO)
            
    Domain: 
22
        The domain must be represented by a mesh using the .msh format (generated by gmsh or unref)
23
        Various forms of physical domain can be used:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
24
25
26
        
        - 2D domain:
            * domain which are projected by mean of cartographic transformation (such as utm projection, stereographic projection)
27
            * sphere. In this case, the local poblems are written in a local orthonormal curvilinear system. All the local problems are assembled in the global discrete algebraic system (for further details, see Comblen et al., 2009, A finite element method for solving the shallow water equations on the sphere).
Valentin Vallaeys's avatar
Valentin Vallaeys committed
28
29
30
        - 3D domain:
            * domain which are projected by mean of cartographic transformation (such as utm projection, stereographic projection)

31
32
    Data: 
        Various file format can be used: 
33
            * netcdf file (data.nc) 
34
                Those files have a special structure: data are given to each nodes at one or several time(s). The file 'slimPre' helps generating those files.
35
            * msh file (data_COMP_0.msh)
36
                 
37
            * idx file (data.idx or data_COMP_0.idx) 
38
                 
39
        Plane domain: In this case, data such as initial condition, wind stress, forcing ... must be expressed in the local basis (x,y) which is defined by the projection.
40
        
41
        Sphere: In this case, data such as initial condition, wind stress, forcing ... must be expressed in the global orthonormal basis (x,y,z) with z pointing upwards.
42
    
43
44
45
46
47
    partition_id:
        This function returns the number of the current partition as a string
    
    partition_nb:
        This function returns the number of partitions used as a string
48
"""
49

50
from dgpy.scripts import slim_private
Valentin Vallaeys's avatar
Valentin Vallaeys committed
51
import dgpy
52
 
David Vincent's avatar
David Vincent committed
53
class Domain:
54
55
    """Create the numerical domain"""
    
56
    def __init__(self, mesh, bathy, g=None, density=1025, solve_on_sphere=False): 
57
        """keyword arguments:
58
        
59
        * mesh            
60
            path to the mesh file (.msh format). The partitioned mesh will automatically be loaded in multiprocessing
61
62
63
        * bathy           
            bathymetric file [in meters, positive] (.msh, .idx or .nc format)
        * g               
David Vincent's avatar
David Vincent committed
64
            map of the mean gravitational acceleration (.msh, .idx or .nc format) (default: 9.81 m/s^2)
65
66
67
68
        * density         
            density of the liquid (default: 1025 kg/m^3)
        * solve_on_sphere 
            boolean stating if you want to solve the shallow water equation on a sphere or on any 2d manifold instead of solving it on a plane (default: False)
69
        """
70
        self._solve_on_manifold = solve_on_sphere
71
72
        if (dgpy.Msg.GetCommSize() > 1) and (not slim_private._check_if_partionned(mesh)):
            dgpy.dgMeshPartition(mesh, dgpy.Msg.GetCommSize())
73
            self._mesh = dgpy.dgMesh(mesh[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh')
74
        else:
75
            self._mesh = dgpy.dgMesh(mesh)
76
77
        #TODO add option to run in high order
        order = 1
78
        self._groups = dgpy.dgGroupCollection(self._mesh, order)
79
        if self._solve_on_manifold:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
80
            self._space_transform = dgpy.dgSpaceTransformFlatTriangle(self._mesh); 
81
            self._mesh.setSpaceTransform(self._space_transform, self._groups)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
82
        self._bath = dgpy.dgDofContainer(self._groups, 1)
83
        slim_private._load(self._bath, bathy)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
84
        self._bath_PC = dgpy.functionPrecomputed(self._groups, 3)
85
        self._bath_PC.compute(self._bath.getFunction())
David Vincent's avatar
David Vincent committed
86
        self._bath_PC.compute(self._bath.getFunction(),dgpy.dataCacheMap.NODE_GROUP_MODE)
87
        self._bath_gradient = self._bath.getFunctionGradient()
Valentin Vallaeys's avatar
Valentin Vallaeys committed
88
        self._bath_gradient_PC = dgpy.functionPrecomputed(self._groups, 3)
89
        self._bath_gradient_PC.compute(self._bath.getFunctionGradient())
90
91
        self._bath_function = self._bath_PC
        self._bath_gradient_function = self._bath_gradient_PC
David Vincent's avatar
David Vincent committed
92
        if g is None:
93
            self._g = dgpy.functionConstant(9.81)
David Vincent's avatar
David Vincent committed
94
        else:
95
            self._g = slim_private._load_function(g,self._groups)
96
        self._density = density
97
 
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
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:
        
        * physical_tag  
            tag of the part of the boundary on which this boundary condition is applied.
        * discharge    
            netcdf or .msh file containing the discharge which will be applied at the boundary [in m^3s^-1] (default: None)
        * sse            
            netcdf or .msh file containing the sea surface elevation which will be applied at the boundary [in m] (default: None, i.e. the sea surface elevation computed by the model is set to this boundary)
        * ux             
            netcdf or .msh file containing the velocity along the x axis of the local basis  or the velocity along the x axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
        * uy             
            netcdf or .msh file containing the velocity along the y axis of the local basis  or the velocity along the y axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
        * uz             
            netcdf or .msh file containing the velocity along the z axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
        * transport_flux 
            set whether ux and uy (and uz) are velocities (False) or transport (True) (Default: False) 
        * ramp_period  
            period for the linear ramp before applying the nudging [in s] (Default: 0s, no ramp) 
        """
        #TODO add option for external gradient
        _uz_= _sse = _ux = _uy = _velocity_local = _discharge = _depth_integral = None        
        if (ux is None and uy is not None) or (ux is not None and uy is None):
            dgpy.Msg.Fatal("Both ux and uy should be given for boundary condition for tag: "+physical_tag)
        if ux is not None and discharge is not None:
            dgpy.Msg.Fatal("Impossible to impose the velocity and the discharge for boundary condition for tag: "+physical_tag)
        if sse is not None and discharge is not None:
            dgpy.Msg.Fatal("Impossible to impose the sea surface elevation and the discharge for boundary condition for tag: "+physical_tag)
            #TODO add option to prescribe etaExt
        if sse is None and ux is None and discharge is None:
            _external_value = dgpy.outsideValueFree(self._solution.getFunction(), self._domain._bath_function, self._domain._g)
        elif discharge is not None:
            _discharge = slim_private._load_function(discharge, self._domain._groups)
            _depth_integral = dgpy.fullMatrixDouble(1, 1)
            dgpy.dgFunctionIntegratorInterface(self._domain._groups, self._domain._bath_PC).compute(physical_tag, _depth_integral)
            _external_value = dgpy.dischargeToBoundaryVelocity(self._solution.getFunction(), _discharge, _depth_integral(0, 0), self._domain._bath_function)
            if (ramp_period > 0.):
                _external_value_ramp = dgpy.ramp3d(_external_value, dgpy.function.getTime(), self._initial_time, ramp_period)
            else:
                _external_value_ramp = _external_value
        elif ux is None:
            _sse = slim_private._load_function(sse, self._domain._groups)
            _external_value = dgpy.outsideValueEta(self._solution.getFunction(), _sse)
        else:
            _ux_= slim_private._load_function(ux, self._domain._groups)
            _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)    
            _ux = slim_private._get_x(_velocity_local).functor
            _uy = slim_private._get_y(_velocity_local).functor
            if transport_flux:
                if sse is not None:
                    _sse = slim_private._load_function(sse, self._domain._groups)    
                    _external_value = dgpy.transportToBoundaryVelocity(_sse, _ux, _uy, self._domain._bath_function)
                else:
                    _external_value = dgpy.transportToBoundaryVelocity(self._solution.getFunction(), _ux, _uy, self._domain._bath_function)
            else:
                if sse is not None:
                    _sse = slim_private._load_function(sse, self._domain._groups)    
                _external_value = dgpy.mergeVector(_sse, _ux, _uy, self._solution.getFunction())
        if (ramp_period > 0.):
            _external_value_ramp = dgpy.ramp3d(_external_value, dgpy.function.getTime(), self._initial_time, ramp_period)
        else:
            _external_value_ramp = _external_value
        def nonConservativeToConservative(cm, solF, bathF) :
            sol = cm.get(solF)
            n = cm.get(dgpy.function.getNormals())
            solIn = cm.get(cm.solution())
            H = cm.get(bathF) + sol[:,[0]]
            solExt = slim_private.np.hstack([H,H*sol[:,[1,2]]])
            solExt[:,[1,2]] = solIn[:,[1,2]]- n[:,[0,1]] * (9.81/H)**0.5*(solIn[:,[0]]-solExt[:,[0]])/2
            return solExt
        conservative = dgpy.functorNumpy(lambda cm : nonConservativeToConservative(cm,_external_value_ramp, self._domain._bath_function))
        self._equation.addBoundaryCondition(physical_tag, self._equation.newOutsideValueBoundary(conservative))
        self._ref += [_sse, _ux, _uy, _velocity_local, _external_value, _external_value_ramp, ramp_period, _discharge, _depth_integral, conservative] 
 
    def export_value_at_point(self, output_name, x, y, z=None):
        """set the coordinates where a time serie of the solution is exported

        keyword arguments:
        
        * output_name 
            name of the output file
        * coord_x     
            vector containing the x coordinates where the solution is evaluated
        * coord_y     
            vector containing the y coordinates where the solution is evaluated
        * coord_z     
            vector containing the z coordinates where the solution is evaluated (Default: None, replaced by vector of 0)
        """
        if self._evaluator is None:
            self._evaluator = dgpy.dgFunctionEvaluator(self._domain._groups, self._solution.getFunction())
        self._evaluator_points_x = x
        self._evaluator_points_y = y
        if z is None:
            z = [0.]*len(x)
        self._evaluator_points_z = z
        self._evaluator_output_name = output_name
        f = open(output_name, 'w')
        sx = "x,"
        sy = "y,"
        sz = "z,"
        for i in range(len(x)):
            sx += str(x[i])+","
            sy += str(y[i])+","
            sz += str(z[i])+","
        f.write(sx[:-1]+'\n')
        f.write(sy[:-1]+'\n')
        f.write(sz[:-1]+'\n')
        f.close()
    
    def compute_mass(self,output_name):
        """compute the total mass of the domain
        
        keyword arguments:
        
        * output_name 
            name of the output file
        """
        self._output_mass = output_name
        self._compute_mass = True
        self._height = dgpy.functorNumpy(lambda cm : (cm.get(self._solution.getFunction())[:,[0]])+cm.get(self._domain._bath_function))
        self._integrator_mass = dgpy.dgFunctionIntegrator(self._domain._groups, self._height)
        self._integrator_mass.setSolution(self._solution.getFunction(), self._solution.getFunctionGradient())
        self._mass = dgpy.fullMatrixDouble(1, 1)
        self._integrator_mass.compute(self._mass,"")

David Vincent's avatar
David Vincent committed
357
class ShallowWater2d:
358
    """Create the shallow water equation with various options"""
David Vincent's avatar
David Vincent committed
359

360
    def __init__(self, domain, temporal_scheme, linear_equation=False, wetting_drying=None, export_every_sub_time_step=False, initial_time=None, final_time=None): 
361
        """keyword arguments:
362
        
363
364
        * domain                     
            object of class Domain
365
366
367
368
369
370
371
372
373
374
375
        * 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
376
377
378
379
380
381
382
383
384
385
386
387
        * linear_equation            
            boolean stating if you want tos solve the linear version of the  shallow water equation (Default: False)
        * wetting_drying             
            value for the wetting & drying algorithm (see Karna et al. 2011) (Default: None)
        * 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]
388
            - date string in format "year-month-day hours:minutes:seconds" (e.g. "2000-02-15 10:05:00")
389
390
391
392
393
394
        * final_time
            final time of the simulation (default: 0)
            
            Format:
                
            - float [in seconds]
395
396
397
            - 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'   
398
399
        if initial_time is None:
            self._initial_time = 0.
400
        elif slim_private._is_string(initial_time):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
401
            date = slim_private.datetime.datetime.strptime(initial_time, fmt)
402
            self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
403
404
405
406
        else:
            self._initial_time = float(initial_time)
        if final_time is None:
            self._final_time = 0.
407
        elif slim_private._is_string(final_time):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
408
            date = slim_private.datetime.datetime.strptime(final_time, fmt)
409
            self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
410
411
        else:
            self._final_time = float(final_time)
David Vincent's avatar
David Vincent committed
412
        self._domain = domain
Valentin Vallaeys's avatar
Valentin Vallaeys committed
413
        self._equation = dgpy.dgConservationLawShallowWater2d()
David Vincent's avatar
David Vincent committed
414
        self._equation.setIsLinear(linear_equation)
David Vincent's avatar
David Vincent committed
415
        self._equation.setLaxFriedrichs(True)
416
        self._wetting_drying = wetting_drying
417
418
        self._equation.setBathymetry(self._domain._bath_function)
        self._equation.setBathymetryGradient(self._domain._bath_gradient_function)
419
        self._equation.setGravity(self._domain._g)
420
        self._coriolis = None
Valentin Vallaeys's avatar
Valentin Vallaeys committed
421
        self._solution = dgpy.dgDofContainer(self._domain._groups, self._equation.getNbFields())
David Vincent's avatar
David Vincent committed
422
423
424
425
        self._solution.setAll(0.0)
        self._solution.setFieldName(0, 'eta')
        self._solution.setFieldName(1, 'u')
        self._solution.setFieldName(2, 'v')
426
427
428
429
        if self._wetting_drying is not None:
            self._equation.setMovingBathWettingDrying(self._solution.getFunction(), self._solution.getFunctionGradient(), self._wetting_drying)
            self._domain._bath_function = self._equation.getBathymetry()
            self._domain._bath_gradient_function = self._equation.getBathymetryGradient()
430
        self._name = "sw2d"
431
432
433
434
        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 !")
435
        self._export_every_sub_time_step = export_every_sub_time_step
436
        self._evaluator = None
437
        self._ref = []
438
        self._compute_mass = False
David Vincent's avatar
David Vincent committed
439
        self._export_on_structured_grid =False
David Vincent's avatar
David Vincent committed
440
  
441
    def set_initial_condition(self, sse, ux, uy, uz=None):
442
        """
443
        Initial conditions
444

445
        keyword arguments:
446
447
        
        * sse 
448
            sea surface elevation [in m] (.msh or .nc format) 
449
        * ux  
450
            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) 
451
        * uy  
452
            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) 
453
        * uz  
454
            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)
455
        """
456
        _sse = _ux = _uy = _uz = _ux_ = _uy_ = None
457
        if sse is not None:
458
            _sse = slim_private._load_function(sse, self._domain._groups)
459
        if ux is not None:
460
            _ux = slim_private._load_function(ux, self._domain._groups)
461
        if uy is not None:
462
            _uy = slim_private._load_function(uy, self._domain._groups)
463
        if uz is not None:
464
            _uz = slim_private._load_function(uz, self._domain._groups)
465
        _velocity_local = slim_private._check_vector(_ux,_uy,_uz,self._domain)
466
467
        _ux_ = slim_private._get_x(_velocity_local).functor
        _uy_ = slim_private._get_y(_velocity_local).functor
Valentin Vallaeys's avatar
Valentin Vallaeys committed
468
469
        _f0 = dgpy.functionConstant([0.,0.,0.])
        _initial_condition = dgpy.mergeVector(_sse, _ux_, _uy_, _f0)
470
471
        self._solution.interpolate(_initial_condition)

472
    def set_viscosity(self, mode, smagorinsky_coefficient=0.1, constant_viscosity=1e-6):
473
        """Turbulent viscosity : Smagorinsky scheme or constant value
474
475
        
        keyword arguments:
476
477
478
479
480
481
482
483
484
485
486
487
        
        * mode       
            type of turbulent viscosity
            
            * "smagorinsky" 
                Smagorinsky scheme (see Smagorinsky 1962)
            * "constant"    
                constant value 
        * smagorinsky_coefficient 
            coefficient value for Smagorinsky scheme (default: 0.1)
        * constant_viscosity      
            constant value [in m^2/s] (default: 1e-6 m^2/s)
488
        """
489
        if mode == "smagorinsky":
490
491
            self._smagorinsky_coefficient = smagorinsky_coefficient
            self._nu = dgpy.smagorinsky2d(self._smagorinsky_coefficient, self._solution.getFunctionGradient(), dgpy.function.getElementVolume())
492
            self._equation.setDiffusivity(self._nu)
493
        elif mode == "constant":
494
495
            self._constant_viscosity = constant_viscosity
            self._nu = dgpy.functionConstant(self._constant_viscosity)
496
            self._equation.setDiffusivity(self._nu)
497
        else: 
Valentin Vallaeys's avatar
Valentin Vallaeys committed
498
            dgpy.Msg.Fatal("Viscosity scheme not implemented, use smagorinsky or a constant viscosity scheme")
499

500
    def set_dissipation(self, mode, coefficient=None):
501
        """Dissipation due to bottom friction: linear, basic quadratic scheme or Chezy-Manning scheme
502
503
        
        keyword arguments:
504
505
506
507
508
        
        * mode
            scheme used for the bottom friction
                                      
            * "linear"    
509
                linear dissipation tau_b/(rho*H) = gamma*u (default: gamma = 1e-6 s^-1)
510
            * "quadratic" 
511
                basic quadratic dissipation tau_b/(rho*H) = (C_d/H)*|u|*u (default: C_d = 2.5e-3)
512
            * "manning"   
513
514
515
                Chezy-Manning-Strickler formulation tau_b/(rho*H) = (n*n*g/(H^4/3))*|u|*u (default: n = 0.03 sm^(-1/3))
        * coefficient
            netcdf or .msh file containing the coefficient value (gamma, C_d or n depending on the mode) over the whole domain (default: None, the default value of the mode will be used)
516
        """
517
518
        if coefficient is not None:
            self._coefficient = slim_private._load_function(coefficient, self._domain._groups)
519
        if mode == "linear":
520
521
522
            if coefficient is None:
                self._coefficient = dgpy.functionConstant(1e-6)
            self._equation.setLinearDissipation(self._coefficient)
523
        elif mode == "quadratic":
524
525
526
            if coefficient is None:
                self._coefficient = dgpy.functionConstant(2.5e-3)
            self._cd = dgpy.quadraticBottomFriction(self._coefficient, self._domain._bath_function, self._solution.getFunction())
527
            self._equation.setQuadraticDissipation(self._cd)
528
        elif mode == "manning":
529
530
531
            if coefficient is None:
                self._coefficient = dgpy.functionConstant(0.03)
            self._cd = dgpy.manningBottomFriction(self._coefficient, self._domain._g, self._domain._bath_function, self._solution.getFunction())
532
            self._equation.setQuadraticDissipation(self._cd)
533
        else: 
Valentin Vallaeys's avatar
Valentin Vallaeys committed
534
            dgpy.Msg.Fatal("Dissipation scheme not implemented!")
535
536
    
    def set_coriolis(self, coriolis):
537
        """Add a coriolis term (2*Omega*cos(latitude)) in the shallow water equation 
538
        
539
        keyword argument: 
540
541
542
        
       * coriolis 
            netcdf or .msh file containing the coriolis term for the whole domain. 
543
        """
544
        self._coriolis = slim_private._load_function(coriolis, self._domain._groups)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
545
        self._coriolis_PC = dgpy.functionPrecomputed(self._domain._groups, 3)
546
547
        self._coriolis_PC.compute(self._coriolis)
        self._equation.setCoriolisFactor(self._coriolis_PC)
David Vincent's avatar
David Vincent committed
548
    
549
    def set_wind_stress(self, mode, wind_x, wind_y, density_air=-1, wind_z= None):
550
        """Add a wind stress term in the shallow water equation.
551
        
552
553
        keyword arguments:
        
554
555
        * mode                  
            type of wind forcing
556
            
557
558
559
560
561
562
563
564
565
566
567
568
            * "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 
569
            
570
            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain.  
571
572
        * density_air  
            density of the ambiant air (only necessary for "speed" mode) (default: -1)
573
        * wind_z 
574
            netcdf or .msh file containing the wind term along the z-axis in the global basis (x,y,z) for the whole domain (only for equation on a sphere). (default: None)
575
        """
576
577
578
579
580
        if mode == 'speed':
            if density_air < 0:
                dgpy.Msg.Fatal("Setting wind stress. in 'speed' mode, density_air must be set")
            _wind_speed_x = slim_private._load_function(wind_x, self._domain._groups)
            _wind_speed_y = slim_private._load_function(wind_y, self._domain._groups)
581
            if wind_z is not None:
582
583
584
585
                _wind_speed_z = slim_private._load_function(wind_z, self._domain._groups)
            else:
                _wind_speed_z = None
            self._wind_speed_local = slim_private._check_vector(_wind_speed_x, _wind_speed_y, _wind_speed_z, self._domain)
586
            self._wind_stress_local = dgpy.windCubicStress(self._wind_speed_local, density_air)
587
588
589
590
591
            self._ref += [density_air]    
            self._equation.setWindStress(self._wind_stress_local)   
        elif mode == 'stress':
            _wind_stress_x = slim_private._load_function(wind_x, self._domain._groups)
            _wind_stress_y = slim_private._load_function(wind_y, self._domain._groups)
592
            if wind_z is not None:
593
594
595
596
597
                _wind_stress_z = slim_private._load_function(wind_z, self._domain._groups)
            else:
                _wind_stress_z = None
            self._wind_stress_local = slim_private._check_vector(_wind_stress_x, _wind_stress_y, _wind_stress_z, self._domain)
            self._equation.setWindStress(self._wind_stress_local)
598
        
599
    def set_forcing(self, forcing_x, forcing_y, forcing_z = None):
600
        """Add a momentum source term derived from a gravitational potential in the shallow water equation.
601
        
602
        keyword arguments: 
603
604
605
606
        
        * forcing_x 
            netcdf or .msh file containing the forcing term along the x-axis in the local basis 
            
607
            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in ms^-2].  
608
609
610
        * forcing_y 
            netcdf or .msh file containing the forcing term along the y-axis in the local basis 
            
611
            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in ms^-2].
612
613
614
        * forcing_z 
            netcdf or .msh file containing the forcing term along the z-axis in the local basis 
            
615
            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in ms^-2].
616
        """
617
618
        _forcing_x = slim_private._load_function(forcing_x, self._domain._groups)
        _forcing_y = slim_private._load_function(forcing_y, self._domain._groups)
David Vincent's avatar
David Vincent committed
619
        _forcing_z = None
620
        if forcing_z is not None:
621
            _forcing_z = slim_private._load_function(forcing_z, self._domain._groups)
622
        self._forcing_local = slim_private._check_vector(_forcing_x, _forcing_y, _forcing_z, self._domain)
623
        self._equation.setSource(self._forcing_local)
624

Valentin Vallaeys's avatar
Valentin Vallaeys committed
625
    def set_nudging(self, coefficient, external_sse, external_ux, external_uy, external_uz=None, ramp_period=0, transport_flux=False):
626
        """Add a nudging term as a source term in the shallow water equation.
627
        
628
        keyword arguments: 
629
630
631
632
633
634
        
        * coefficient    
            netcdf or .msh file containing the nudging coefficient for the whole domain [in s^-1]
        * external_sse   
            netcdf or .msh file containing the external sea surface elevation for the whole domain [in m]
        * external_ux    
635
            netcdf or .msh file containing the external velocity along the x-axis expressed 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 [in m/s]
636
        * external_uy    
637
            netcdf or .msh file containing the external velocity along the y-axis expressed 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 [in m/s]
638
639
640
641
642
643
        * external_uz    
            netcdf or .msh file containing the external velocity along the z-axis expressed in the global basis (x,y,z) for the whole domain [in m/s]
        * ramp_period    
            period for the linear ramp before applying the nudging [in s] (Default: 0s, no ramp) 
        * transport_flux 
            set whether ux and uy (and uz) are velocities (False) or transport (True) (Default: False) 
644
        """
645
646
647
648
        _coeff_nudging = slim_private._load_function(coefficient, self._domain._groups)
        _external_sse = slim_private._load_function(external_sse, self._domain._groups)
        _external_ux = slim_private._load_function(external_ux, self._domain._groups)
        _external_uy = slim_private._load_function(external_uy, self._domain._groups)
649
650
        _external_uz = None
        if external_uz is not None:
651
            _external_uz = slim_private._load_function(external_uz, self._domain._groups)
652
653
654
655
656
        _external_vel_local_ = slim_private._check_vector(_external_ux, _external_uy, _external_uz, self._domain)
        _external_vel_local_x = slim_private._get_x(_external_vel_local_).functor
        _external_vel_local_y = slim_private._get_y(_external_vel_local_).functor
        _external_sol_local = slim_private._merger_3d(_external_sse, _external_vel_local_x, _external_vel_local_y).functor
        self._ref += [_coeff_nudging, ramp_period, transport_flux]    
657
        self._equation.setNudgingCoefficient(_coeff_nudging)
658
        if (ramp_period > 0.):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
659
            self._external_sol_local = dgpy.ramp3d(_external_sol_local, dgpy.function.getTime(), self._initial_time, ramp_period)
660
661
        else:
            self._external_sol_local = _external_sol_local
Valentin Vallaeys's avatar
Valentin Vallaeys committed
662
        self._equation.setNudgingVelocity(self._external_sol_local, transport_flux)
663
664

    def set_boundary_coast(self, physical_tag, slip=False):
665
        """Boundary condition: normal speed set to zero and tangential speed set to zero or free (no slip or free slip condition)
666
        
667
        keyword arguments:
668
669
670
671
672
        
        * physical_tag 
            tag of the part of the boundary on which this boundary condition is applied.
        * slip         
            flag wheter applying slip (True) or no slip (False) condition at the boundary (Default: False)
673
        """
674
675
        self._equation.addBoundaryCondition(physical_tag, self._equation.newBoundaryWall(slip))

676
677
    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).
678
        
679
        keyword arguments:
680
681
682
        
        * physical_tag  
            tag of the part of the boundary on which this boundary condition is applied.
683
684
        * discharge    
            netcdf or .msh file containing the discharge which will be applied at the boundary [in m^3s^-1] (default: None)
685
686
687
        * sse            
            netcdf or .msh file containing the sea surface elevation which will be applied at the boundary [in m] (default: None, i.e. the sea surface elevation computed by the model is set to this boundary)
        * ux             
688
            netcdf or .msh file containing the velocity along the x axis of the local basis  or the velocity along the x axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
689
        * uy             
690
            netcdf or .msh file containing the velocity along the y axis of the local basis  or the velocity along the y axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
691
        * uz             
692
            netcdf or .msh file containing the velocity along the z axis of the global basis if you solve the equation on a sphere. This velocity will be applied at the boundary [in m/s] (default: None, i.e. the velocity computed by the model is set to this boundary)
693
694
695
696
        * transport_flux 
            set whether ux and uy (and uz) are velocities (False) or transport (True) (Default: False) 
        * ramp_period  
            period for the linear ramp before applying the nudging [in s] (Default: 0s, no ramp) 
697
698
        """
        #TODO add option for external gradient
699
        _uz_= _sse = _ux = _uy = _velocity_local = _discharge = _depth_integral = None        
700
701
702
703
704
705
706
707
708
709
710
711
        if (ux is None and uy is not None) or (ux is not None and uy is None):
            dgpy.Msg.Fatal("Both ux and uy should be given for boundary condition for tag: "+physical_tag)
        if ux is not None and discharge is not None:
            dgpy.Msg.Fatal("Impossible to impose the velocity and the discharge for boundary condition for tag: "+physical_tag)
        if sse is not None and discharge is not None:
            dgpy.Msg.Fatal("Impossible to impose the sea surface elevation and the discharge for boundary condition for tag: "+physical_tag)
            #TODO add option to prescribe etaExt
        if sse is None and ux is None and discharge is None:
            _external_value = dgpy.outsideValueFree(self._solution.getFunction(), self._domain._bath_function, self._domain._g)
        elif discharge is not None:
            _discharge = slim_private._load_function(discharge, self._domain._groups)
            _depth_integral = dgpy.fullMatrixDouble(1, 1)
712
            dgpy.dgFunctionIntegratorInterface(self._domain._groups, self._domain._bath_PC).compute(physical_tag, _depth_integral)
713
            _external_value = dgpy.dischargeToBoundaryVelocity(self._solution.getFunction(), _discharge, _depth_integral(0, 0), self._domain._bath_function)
714
715
716
717
            if (ramp_period > 0.):
                _external_value_ramp = dgpy.ramp3d(_external_value, dgpy.function.getTime(), self._initial_time, ramp_period)
            else:
                _external_value_ramp = _external_value
718
        elif ux is None:
719
            _sse = slim_private._load_function(sse, self._domain._groups)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
720
            _external_value = dgpy.outsideValueEta(self._solution.getFunction(), _sse)
721
        else:
722
723
724
725
726
            _ux_= slim_private._load_function(ux, self._domain._groups)
            _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)    
David Vincent's avatar
David Vincent committed
727
728
            _ux = slim_private._get_x(_velocity_local).functor
            _uy = slim_private._get_y(_velocity_local).functor
729
730
            if transport_flux:
                if sse is not None:
731
                    _sse = slim_private._load_function(sse, self._domain._groups)    
732
                    _external_value = dgpy.transportToBoundaryVelocity(_sse, _ux, _uy, self._domain._bath_function)
733
                else:
734
                    _external_value = dgpy.transportToBoundaryVelocity(self._solution.getFunction(), _ux, _uy, self._domain._bath_function)
735
736
            else:
                if sse is not None:
737
                    _sse = slim_private._load_function(sse, self._domain._groups)    
Valentin Vallaeys's avatar
Valentin Vallaeys committed
738
                _external_value = dgpy.mergeVector(_sse, _ux, _uy, self._solution.getFunction())
Valentin Vallaeys's avatar
Valentin Vallaeys committed
739
        if (ramp_period > 0.):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
740
            _external_value_ramp = dgpy.ramp3d(_external_value, dgpy.function.getTime(), self._initial_time, ramp_period)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
741
742
743
        else:
            _external_value_ramp = _external_value
        self._equation.addBoundaryCondition(physical_tag, self._equation.newOutsideValueBoundary(_external_value_ramp))
744
        self._ref += [_sse, _ux, _uy, _velocity_local, _external_value, _external_value_ramp, ramp_period, _discharge, _depth_integral] 
745
 
746
    def export_value_at_point(self, output_name, x, y, z=None):
747
        """set the coordinates where a time serie of the solution is exported
748
749

        keyword arguments:
750
751
752
753
754
755
756
757
758
        
        * output_name 
            name of the output file
        * coord_x     
            vector containing the x coordinates where the solution is evaluated
        * coord_y     
            vector containing the y coordinates where the solution is evaluated
        * coord_z     
            vector containing the z coordinates where the solution is evaluated (Default: None, replaced by vector of 0)
759
760
        """
        if self._evaluator is None:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
761
            self._evaluator = dgpy.dgFunctionEvaluator(self._domain._groups, self._solution.getFunction())
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
        self._evaluator_points_x = x
        self._evaluator_points_y = y
        if z is None:
            z = [0.]*len(x)
        self._evaluator_points_z = z
        self._evaluator_output_name = output_name
        f = open(output_name, 'w')
        sx = "x,"
        sy = "y,"
        sz = "z,"
        for i in range(len(x)):
            sx += str(x[i])+","
            sy += str(y[i])+","
            sz += str(z[i])+","
        f.write(sx[:-1]+'\n')
        f.write(sy[:-1]+'\n')
        f.write(sz[:-1]+'\n')
        f.close()
780
    
781
    def compute_mass(self,output_name):
782
        """compute the total mass of the domain
783
784
785
786
787
        
        keyword arguments:
        
        * output_name 
            name of the output file
788
        """
789
        self._output_mass = output_name
790
        self._compute_mass = True
791
        self._height = dgpy.functorNumpy(lambda cm : (cm.get(self._solution.getFunction())[:,[0]])+cm.get(self._domain._bath_function))
792
        self._integrator_mass = dgpy.dgFunctionIntegrator(self._domain._groups, self._height)
793
        self._integrator_mass.setSolution(self._solution.getFunction(), self._solution.getFunctionGradient())
Valentin Vallaeys's avatar
Valentin Vallaeys committed
794
        self._mass = dgpy.fullMatrixDouble(1, 1)
795
        self._integrator_mass.compute(self._mass,"")
796

797

798
class ShallowWaterTracer2d:
799
    """Create the shallow water equation for a tracer"""
800
    def __init__(self, domain, temporal_scheme, hydro_sol, linear_equation=False, name="Tracer", offline=False, wetting_drying=None, export_every_sub_time_step=False, initial_time=None, final_time=None): 
801
        """keyword arguments:
802
        
803
804
        * domain                     
            object of class Domain
805
806
807
808
809
810
811
812
813
814
815
816
817
        * temporal_scheme 
            scheme for the temporal integration
                               
            * "explicit" 
                explicit Euler scheme
            * "implicit" 
                implicit order 2 Runge-Kutta
        * hydro_sol
            hydrodynamics solution
            
            Format:
            
            - equation : ShallowWater2d object (online mode)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
818
            - string : path to the folder containing the hydrodynamics folder 'offline_sw2d' (offline mode).
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
        * linear_equation            
            boolean stating linear shallow water tracer equation (default: False)
        * name                       
            name of the Tracer (default: "Tracer")
        * offline                    
            resolve shallow water equation (False) or use file to provide hydrodynamics (True) (default: False)
        * wetting_drying             
            value for the wetting & drying algorithm. It has to be the same as for the shallow water equation (see Karna et al. 2011) (default: None)
        * export_every_sub_time_step 
            boolean stating if you want to export the hydrodynamics every sub time step for further use of offline tracer (default: False)
        * initial_time               
            initial time of the simulation (default: 0)
            
            Format:
            
            - float [in seconds]
835
            - date string in format "year-month-day hours:minutes:seconds" (e.g. "2000-06-20 10:05:00")
836
837
838
839
840
841
        * final_time                 
            final time of the simulation (default: 0)
            
            Format:
            
            - float [in seconds]
842
            - date string in format "year-month-day hours:minutes:seconds" (e.g. "2000-06-20 22:42:17")
843
        """
844
        fmt = '%Y-%m-%d %H:%M:%S'   
845
846
        if initial_time is None:
            self._initial_time = 0.
847
        elif slim_private._is_string(initial_time):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
848
            date = slim_private.datetime.datetime.strptime(initial_time, fmt)
849
            self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
850
851
852
853
        else:
            self._initial_time = float(initial_time)
        if final_time is None:
            self._final_time = 0.
854
        elif slim_private._is_string(final_time):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
855
            date = slim_private.datetime.datetime.strptime(final_time, fmt)
856
            self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
857
858
        else:
            self._final_time = float(final_time)
859
        self._domain = domain
Valentin Vallaeys's avatar
Valentin Vallaeys committed
860
        self._equation = dgpy.dgConservationLawShallowWaterTracer2d()
861
        self._equation.setIsLinear(linear_equation)
David Vincent's avatar
David Vincent committed
862
        self._equation.hydroWithLaxFriedrichs(True)
863
        self._equation.setGravity(self._domain._g)
864
        self._wetting_drying = wetting_drying
865
        self._equation.setBathymetry(self._domain._bath_function)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
866
        self._solution = dgpy.dgDofContainer(self._domain._groups, 1)
867
868
        self._solution.setAll(0.0)
        self._solution.setFieldName(0, name)
869
        self._hydro_sol = None
870
        self._name = name
871
872
873
874
        if temporal_scheme == "explicit" or temporal_scheme == "implicit":
            self._temporal_scheme = temporal_scheme
        else:
            dgpy.Msg.Fatal("Unknown temporal scheme "+temporal_scheme+" for shallow water tracer equation !")
875
        self._offline = offline
876
877
878
        if self._offline:
            if not slim_private._is_string(hydro_sol):
                dgpy.Msg.Fatal("You have to provide the sea surface elevation and the velocity to solve the tracer in offline mode!")
879
            self._datafile = hydro_sol+'/offline_sw2d/'
880
            self._hydro_sol_dof = dgpy.dgDofContainer(self._domain._groups, 3)
881
            self._hydro_sol_dof.importMsh(self._datafile+"offline_sw2d-000000_000000/offline_sw2d-000000_000000")
882
883
884
885
886
            self._hydro_sol = self._hydro_sol_dof.getFunction()
        else:
            if not isinstance(hydro_sol, ShallowWater2d):
                dgpy.Msg.Fatal("You have to provide the hydrodynamic equation to solve the tracer!")
            self._hydro_sol = hydro_sol._solution.getFunction()
887
888
889
        if self._wetting_drying is not None:
            self._equation.setMovingBathWettingDrying(self._hydro_sol, None, self._wetting_drying)
            self._domain._bath_function = self._equation.getBathymetry()
890
        self._equation.setHydroSolution(self._hydro_sol) 
891
        self._export_every_sub_time_step = export_every_sub_time_step
Valentin Vallaeys's avatar
Valentin Vallaeys committed
892
        self._evaluator = None
893
        self._ref = []
894
        self._compute_mass = False
David Vincent's avatar
David Vincent committed
895
        self._export_on_structured_grid =False
896
897
           
    def set_initial_condition(self, initial_condition):
898
       """Initial conditions
899
900
        
        keyword argument:
901
902
903
        
        * tracer 
            tracer initial condition [-] (.msh or .nc format) (default: 0 )
904
        """
905
       slim_private._load(self._solution, initial_condition)
906
    
907
    def set_diffusivity(self, mode, okubo_coefficient=0.03, constant_diffusivity=1e-6):
908
        """Diffusivity
909
910
911
        
        keyword arguments:

912
        * mode                  
913
914
915
916
917
918
919
920
921
922
            type of diffusivity
            
            * "okubo"   
                Okubo scheme 
            * "constant" 
                constant value 
        * okubo_coefficient   
            coefficient for Okubo scheme [in m^0.85/s] (default: 0.03 m^0.85/s)
        * constant_diffusivity 
            constant value [in m^2/s] (default: 1e-6 m^2/s)
923
        """
924
        if mode == "okubo":
Valentin Vallaeys's avatar
Valentin Vallaeys committed
925
            self._kappa = dgpy.okubo2d(okubo_coefficient, dgpy.function.getElementVolume())
926
            self._equation.setDiffusivity(self._kappa)
927
        elif mode == "constant":
Valentin Vallaeys's avatar
Valentin Vallaeys committed
928
            self._kappa = dgpy.functionConstant(constant_diffusivity)
929
            self._equation.setDiffusivity(self._kappa)
930
        else: 
931
932
            dgpy.Msg.Fatal("Diffusivity scheme not implemented, use okubo or a constant diffusivity scheme") 

933
    def set_source(self, tracer_source):
934
        """Add a tracer source term in the tracer equation.
935
936
        
        keyword argument: 
937
938
939
        
        * tracer_source 
            netcdf or .msh file containing the source term for the whole domain [in c*s^-1]. 
940
        """
941
        self._source = slim_private._load_function(tracer_source, self._domain._groups)
942
        self._equation.setSource(self._source)    
943

944
    def set_boundary_coast(self, physical_tag):
945
        """Boundary condition: no flux condition at the boundary
946
947
        
        keyword argument:
948
949
950
        
        * physical_tag 
            tag of the part of the boundary on which this boundary condition is applied.
951
952
953
        """
        self._equation.addBoundaryCondition(physical_tag, self._equation.newBoundaryWall())

954
955
    def set_boundary_open(self, physical_tag, concentration):
        """Open boundary condition: one can impose the external concentration
956
957
        
        keyword arguments:
958
959
960
961
962
        
        * physical_tag             
            tag of the part of the boundary on which this boundary condition is applied.
        * concentration            
            netcdf or .msh file containing the concentration which will be applied at the boundary [in c]
963
964
        """
        #TODO add option for external gradient
965
        _external_value = slim_private._load_function(concentration, self._domain._groups)
966
967
        self._equation.addBoundaryCondition(physical_tag, self._equation.newOutsideValueBoundary(_external_value))
        self._ref += [_external_value]
968
       
969
970
971
972
973
974
975
976
977
978
979
980
981
982
    def set_boundary_flux(self, physical_tag, flux):
        """Boundary condition for a tracer flux: one can impose the flux throw the boundary condition, in [c s^-1], (with the tracer units [c]).
        
        keyword arguments:
        
        * physical_tag             
            tag of the part of the boundary on which this boundary condition is applied.
        * flux            
            netcdf or .msh file containing the flux which will be applied at the boundary in [c s^-1]
        """       
        _flux = slim_private._load_function(flux, self._domain._groups)
        _length_interface = dgpy.fullMatrixDouble(1, 1)
        one = slim_private.dgpy.functionConstant(1)
        dgpy.dgFunctionIntegratorInterface(self._domain._groups, one).compute(physical_tag, _length_interface)
983
        _flux_HC = dgpy.tracerFluxToNeumann(_flux, _length_interface(0,0))
984
        self._equation.setBoundaryNeumann(physical_tag, _flux_HC)
985
        self._ref += [_flux, one, _length_interface, _flux_HC]