slim3d.py 19.6 KB
Newer Older
1
import dgpy
2
3
from dgpy.scripts import slim3d_private
from dgpy.scripts import slim_private
4
5
6

class Domain :
    """Create the slim3d Domain"""
Valentin Vallaeys's avatar
Valentin Vallaeys committed
7
    def __init__(self, mesh_file_name, periodic_map_file='', reference_density=1027, gravity=9.81) :
8
9
10
11
12
13
14
15
16
        """keyword arguments:
        
        * mesh_file_name            
            path to the mesh file (.msh format). The partitioned mesh will automatically be loaded in multiprocessing
        """
        (topTags, bottomTags) = slim_private._findTopAndBottomTags(mesh_file_name)
        self._topTags = topTags
        self._bottomTags = bottomTags
        self._slimSolver = dgpy.slim3dSolver(mesh_file_name, bottomTags, topTags, periodic_map_file)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
17
18
        self._density = reference_density
        self._g = gravity
19
20
21
22
23
24
25
26
27
28
29
30

class Slim3d_equations :
    """Create the slim3d Equation"""
    def __init__(self, Domain, salinity=False, temperature=False) :
        self._domain = Domain
        self._slimSolver = Domain._slimSolver
        slimSolver = self._slimSolver
        slimSolver.setSolveS(salinity)
        slimSolver.setSolveT(temperature)
        slimSolver.setSolveUVImplicitVerticalDiffusion(False)
        slimSolver.setSolveSImplicitVerticalDiffusion(False)
        slimSolver.setSolveTImplicitVerticalDiffusion(False)
31
        self._verticalDiffImplicit = False
32
33
34
35
36
37
38
39
40
41
42
        slimSolver.setFlagUVLimiter(True)
        if salinity:
            slimSolver.setFlagSLimiter(True)
        if temperature:
            slimSolver.setFlagTLimiter(True)
        slimSolver.setComputeBottomFriction(True)
        self._horizontal_viscosity = None
        self._horizontal_diffusivity = None
        self._vertical_viscosity = None
        self._vertical_diffusivity = None
        self._coriolis = None
43
        self._atmPress = None
44
        self._linear_density = None
Valentin Vallaeys's avatar
Valentin Vallaeys committed
45
46
        self._reference_density = self._domain._density
        self._gravity = self._domain._g
47
48
49
        self._initial_elevation = None
        self._initial_salinity = None
        self._initial_temperature = None
50
        self._wind_stress = None
51
52
        self._boundary_coast = []
        self._boundary_open = []
53
        self._gotm_option_file = None
54
        self._z0 = [0.005, 0.02]
Valentin Vallaeys's avatar
Valentin Vallaeys committed
55
56
        self._linear2d = False
        self._linear_solver_2D = True
57
58
59
60
61
        self._LFF = 0.
        self._additional_shear = None
        self._areaFunc = None
        self._nuh = None
        self._add_shear = None
62
63
64

    def set_implicit_vertical_diffusion(self, flag):
        slimSolver = self._slimSolver
65
66
67
68
69
70
71
72
        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)
73

Valentin Vallaeys's avatar
Valentin Vallaeys committed
74
75
76
77
78
79
    def set_linear_2d_equations(self, flag):
        self._linear2d = flag

    def set_linear_solver_2d(self, flag):
        self._linear_solver_2D = flag

80
81
82
83
84
85
86
87
    def set_limiter(self, flag):
        slimSolver = self._slimSolver
        slimSolver.setFlagUVLimiter(flag)
        if slimSolver.getSolveS():
            slimSolver.setFlagSLimiter(flag)
        if slimSolver.getSolveT():
            slimSolver.setFlagTLimiter(flag)

88
89
90
    def set_lax_friedrichs(self, factor=1.):
        self._LFF = factor

91
92
93
    def set_conservative_ale(self, flag):
        self._slimSolver.setUseConservativeALE(flag)

94
    def set_horizontal_viscosity(self, mode, constant=1, factor=0.02, maximum=1e9):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
95
96
        if (mode != "smagorinsky") and (mode != "constant"):
            dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode)
97
98
99
100
101
        self._horizontal_viscosity = mode
        self._hor_visc_const = constant
        self._hor_visc_fact = factor
        self._hor_visc_max = maximum

102
103
104
105

    def set_additional_artificial_horizontal_shear(self, additional_shear):
        self._additional_shear = additional_shear

106
    def set_horizontal_diffusivity(self, mode, constant=1, factor=0.03, maximum=1e9):
Valentin Vallaeys's avatar
Valentin Vallaeys committed
107
108
        if (mode != "okubo") and (mode != "constant"):
            dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode)
109
110
111
112
113
        self._horizontal_diffusivity = mode
        self._hor_diff_const = constant
        self._hor_diff_fact = factor
        self._hor_diff_max = maximum

114
115
    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"):
116
117
            dgpy.Msg.Fatal("Unknown viscosity for vertical_viscosity : "+mode)
        self._vertical_viscosity = mode
118
119
120
121
122
        self._vertical_viscosity_constant_value = constant_value
        self._vertical_viscosity_paca_values = paca_values
        if self._verticalDiffImplicit:
            self._slimSolver.setSolveUVImplicitVerticalDiffusion(True)

123
        
124
125
    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"):
126
127
            dgpy.Msg.Fatal("Unknown diffusivity for vertical_diffusivity : "+mode)
        self._vertical_diffusivity = mode
128
129
130
131
132
133
        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)
134

Valentin Vallaeys's avatar
Valentin Vallaeys committed
135
    def set_gotm_option_file(self,f):
136
137
        self._gotm_option_file = f

138
139
    def set_coriolis(self, coriolis):
        self._coriolis = coriolis
140
141
142
    
    def set_atmospheric_pressure(self, atmPress):
        self._atmPress = atmPress
143

144
    def set_bottom_friction(self, flag, z0B=None, z0S=None):
145
        self._slimSolver.setComputeBottomFriction(flag)
146
147
148
149
        if z0B:
            self._z0[0] = z0B
        if z0S:
            self._z0[1] = z0S
150

Valentin Vallaeys's avatar
Valentin Vallaeys committed
151
    def set_linear_density(self, mode, constant_coefficient=0, linear_coefficient=0):
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        """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)
167
168
169
170

    def set_initial_elevation(self, elevation):
        self._initial_elevation = elevation

171
    def set_initial_salinity(self, mode, salinity=None, surface_salinity=0, vertical_gradient=0):
172
173
        self._initial_salinity = (mode, salinity, surface_salinity, vertical_gradient)

174
    def set_initial_temperature(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0):
175
176
        self._initial_temperature = (mode, temperature, surface_temperature, vertical_gradient)

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
    # TODO check value of air_density
    def set_wind_stress(self, mode, wind_x, wind_y, density_air=1.25):
        """Add a wind stress term in the shallow water equation.
        
        keyword arguments:
        
        * mode                  
            type of wind forcing
            
            * "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 
            
            (except for equation on a sphere for which it has to be expressed in the global basis (x,y,z)) for the whole domain.  
        * density_air  
            density of the ambiant air (only necessary for "speed" mode) (default: 1.25 [kg m^-3])
        """
        if mode != "speed" and mode != "stress":
            dgpy.Msg.Fatal("Unknown mode for wind stress: "+mode)
        self._wind_stress = (mode, wind_x, wind_y, density_air)

207
208
209
210
211
212
213
    def set_boundary_coast(self, physical_tags):
        if slim_private._is_string(physical_tags):
            self._boundary_coast.append(physical_tags)
        else:
          for tag in physical_tags:
              self._boundary_coast.append(tag)

214
    def set_boundary_open(self, physical_tags, eta=None, u=None, v=None, salinity=None, temperature=None, flux=None, transport=False):
215
216
217
218
219
220
221
222
223
224
225
226
        groups2 = self._slimSolver.groups2d
        groups3 = self._slimSolver.groups3d
        if (u and not v) or (v and not u):
            dgpy.Msg.Fatal('Error. Either u and v or none of them must be set in open boundary conditions')
        uvF = functorMerge([slim_private._load_function(u, groups3), slim_private._load_function(v, groups3)]) if u else None
        uv2F = functorMerge([slim_private._load_function(u, groups2), slim_private._load_function(u, groups2)]) if u else None
        etaF = slim_private._load_function(eta, groups3) if eta else None
        eta2F = slim_private._load_function(eta,groups2) if eta else None
        flowF  = slim_private._load_function(flux, groups3) if flux else None
        flow2F = slim_private._load_function(flux, groups2) if flux else None
        SF = slim_private._load_function(salinity, groups3) if salinity else None
        TF = slim_private._load_function(temperature, groups3) if temperature else None
227
        if slim_private._is_string(physical_tags):
228
229
230
            physical_tags = [physical_tags]
        for tag in physical_tags:
            openBnd = slim3d_private.OpenBoundary(self._slimSolver, tag, etaF, eta2F, uvF, uv2F, SF, TF, flowF, flow2F, transport)
231
232
233
234
            self._boundary_open.append(openBnd)

class Loop:
    """Temporal solver"""
Valentin Vallaeys's avatar
Valentin Vallaeys committed
235
    def __init__(self, equations, time_step=3600, export_time_step=3600, ratio_full_export=-1, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output"): 
236
237
        self._slimSolver = equations._slimSolver
        self._equations = equations
Valentin Vallaeys's avatar
Valentin Vallaeys committed
238
        self._dt = time_step
239
        self._export_dt = export_time_step
240
        self._ratio_full_export = ratio_full_export
241
242
        fmt = '%Y-%m-%d %H:%M:%S'   
        date = slim_private.datetime.datetime.strptime(initial_time, fmt)
243
        self.initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
244
        self._initial_date = date
245
246
        if slim_private._is_string(final_time):
            date = slim_private.datetime.datetime.strptime(final_time, fmt)
247
            self.final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
248
        else:
249
            self.final_time = self.initial_time + final_time
250
251
252
        self._odir = output_directory
        self._export_uv       = False
        self._export_w        = False
253
        self._export_z        = True
254
255
256
257
258
259
        self._export_S        = False
        self._export_T        = False
        self._export_uvAv2d   = False
        self._export_eta      = False
        self._export_kappav   = False
        self._export_nuv      = False
260
261
262
263
        self._export_uv_nc    = False
        self._export_S_nc     = False
        self._export_eta_nc   = False
        self._export_z_nc     = False
264
        self._export_rho      = False
265
        self._export_nuh      = False
266
267
        self._timeIntegrator = None
        self._export = []
268
269
270
        self._ncExport = []
        self._restart_dir = None
        self._restart_ind = -1
Valentin Vallaeys's avatar
Valentin Vallaeys committed
271
        self._timer = False
272
    
273
274
275
    def export_uv(self, vect=True):
        self._export_uv = True
        self._export_uv_vect = vect
276
    def export_w(self): self._export_w = True
277
    def export_z_coordinate(self): self._export_z = True
278
279
    def export_salinity(self): self._export_S = True
    def export_temperature(self): self._export_T = True
280
281
282
    def export_uv2d(self, vect=True):
        self._export_uvAv2d = True
        self._export_uvAv2d_vect = vect
283
284
285
    def export_elevation(self): self._export_eta = True
    def export_vertical_diffusivity(self): self._export_kappav = True
    def export_vertical_viscosity(self): self._export_nuv = True
286
    def export_rho(self): self._export_rho = True    
287
    def export_nuh(self): self._export_nuh = True    
288
289
290
291
    def export_uv_nc(self): self._export_uv_nc = True
    def export_salinity_nc(self): self._export_S_nc = True
    def export_elevation_nc(self): self._export_eta_nc = True
    def export_z_nc(self): self._export_z_nc = True
292
    
293
294
295
    def restart(self, directory, index):
        self._restart_dir = directory
        self._restart_ind = index
Valentin Vallaeys's avatar
Valentin Vallaeys committed
296
297
298
        
    def set_timer(self):
        self._timer = True
299

300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
    def export_fields(self):
        timeIntegrator = self._timeIntegrator
        if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
        d = self._slimSolver.getDofs()
        if self._export_nuh:
            d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
        if self._export_uvAv2d:
            timeIntegrator.exporter_computeUVMean()
        if self._export_w:
            timeIntegrator.exporter_computeW()
        if self._export_rho:
            timeIntegrator.exporter_computeRho()
        if self._export_eta:
            timeIntegrator.exporter_computeEta()
        if self._export_z:
            timeIntegrator.exporter_computeZ()
        for e in self._export :
            e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
        if self._ncWriter3d:
            for dof,fname,field in self._ncExport :
                filename = fname+'_{0:05d}'.format(timeIntegrator.getExportCount())
                self._ncWriter3d.exportUGRID(dof, filename, timeIntegrator.getTime() - self.initial_time, timeIntegrator.getIter(),field)

323
    def setup(self):
324
        slim3d_private.slim3d_setup(self)
325
        slimSolver = self._slimSolver
326
        if dgpy.Msg.GetCommRank() == 0:
327
328
            if not slim_private.path.exists(self._odir):
                slim_private.makedirs(self._odir)
329
        dgpy.Msg.Barrier()
330
        d = slimSolver.getDofs()
331
        if self._export_uv:     self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", self._export_uv_vect))
332
333
334
        if self._export_w:      self._export.append(dgpy.dgIdxExporter(d.wDof3dExport, self._odir + "/w"))
        if self._export_z:      self._export.append(dgpy.dgIdxExporter(d.zDofExport, self._odir + "/z"))
        if self._export_rho:    self._export.append(dgpy.dgIdxExporter(d.rhoDof3dExport, self._odir + "/rho"))
335
        if self._export_nuh:    self._export.append(dgpy.dgIdxExporter(d.nuhDof3d, self._odir + "/nuh"))
336
337
        if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2dExport, self._odir + "/uvAv2d", self._export_uvAv2d_vect))
        if self._export_eta:    self._export.append(dgpy.dgIdxExporter(d.etaDof2dExport, self._odir + "/eta"))
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
        if self._export_S: 
            if slimSolver.getSolveS():
                self._export.append(dgpy.dgIdxExporter(d.SDof, self._odir + "/salinity"))
            else: 
                dgpy.Msg.Warning("Salinity is not solved. It won't be exported")
        if self._export_T: 
            if slimSolver.getSolveT():
                self._export.append(dgpy.dgIdxExporter(d.TDof, self._odir + "/temperature"))
            else: 
                dgpy.Msg.Warning("Temperature is not solved. It won't be exported")
        if self._export_kappav: 
            if slimSolver.getSolveTurbulence():
                self._export.append(dgpy.dgIdxExporter(d.kappavDof, self._odir + "/kappav"))
            else: 
                dgpy.Msg.Warning("Turbulence is not solved with gotm. kappa_v won't be exported")
        if self._export_nuv: 
            if slimSolver.getSolveTurbulence():
                self._export.append(dgpy.dgIdxExporter(d.nuvDof, self._odir + "/nuv"))
            else: 
                dgpy.Msg.Warning("Turbulence is not solved with gotm. nu_v won't be exported")
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374


        if slimSolver.havePNETCDF:
            self._ncWriter2d = dgpy.slimNetCDFIO(slimSolver.groups2d)
            self._ncWriter2d.setInitTime(self._initial_date.year,self._initial_date.month,self._initial_date.day,self._initial_date.hour,self._initial_date.minute,self._initial_date.second,0,0)
            self._ncWriter3d = dgpy.slim3dNetCDFIO(slimSolver.columnInfo,self._ncWriter2d)
            if self._export_eta_nc: self._ncExport.append([d.etaDof2d, self._odir + '/eta','eta'])
            if self._export_uv_nc: self._ncExport.append([d.uvDof, self._odir + '/uv','uv'])
            if self._export_z_nc: self._ncExport.append([d.zDof, self._odir + '/z','z'])
            if self._export_S_nc: 
                if slimSolver.getSolveS():
                    self._ncExport.append([d.SDof, self._odir + '/S','S']) 
                else: 
                    dgpy.Msg.Warning("Salinity is not solved. It won't be exported")       
        else:
            self._ncWriter2d = None 
            self._ncWriter3d = None
375
376

        timeIntegrator = self._timeIntegrator
377
378
        if self._restart_ind > -1:
            timeIntegrator.reStart(self._restart_dir, self._restart_ind)
379
        self.export_fields()
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
        self._tic = slim_private.time.time()

    def get_time(self):
        return self._timeIntegrator.getTime()
    def advance_one_time_step(self):
        self._timeIntegrator.advanceOneTimeStep()
    def check_sanity(self):
        self._timeIntegrator.checkSanity()
    def check_export(self):
        return self._timeIntegrator.checkExport()
    def check_full_export(self):
        return self._timeIntegrator.checkFullExport()
    def export_full(self):
        self._timeIntegrator.exportAll()
    def check_tracer_consistency(self):
        return self._timeIntegrator.checkTracerConsistency()
    def check_mass_conservation(self):
        return self._timeIntegrator.checkMassConservation()
    def print_iter_info(self):
        d = self._slimSolver.getDofs()
        uvNorm = d.uvDof.norm()
        if dgpy.Msg.GetCommRank() == 0 : dgpy.Msg.Info("%5i [ %s ] normuv %6.12e clock %.1fs" % (self._timeIntegrator.getIter(),slim_private.time.strftime("%d %b %Y %H:%M:%S", slim_private.time.gmtime(self._timeIntegrator.getTime())), uvNorm, slim_private.time.time() - self._tic))
Valentin Vallaeys's avatar
Valentin Vallaeys committed
402
        return uvNorm
403
404
405
406
407
408
    def terminate(self):
        self._timeIntegrator.terminate()

    def loop(self):
        if not self._timeIntegrator:
            self.setup()
Valentin Vallaeys's avatar
Valentin Vallaeys committed
409
410
411
        if self._timer:
            timer = dgpy.dgTimer.root()
            timer.start()
412
413
414
415
416
417
418
419
420
        while self.get_time() < self.final_time:
            self.advance_one_time_step()
            self.check_sanity()
            if self.check_export(): # Currently, we cannot export without using check_export! Since it is this func that updates the export id's
                self.export_fields()
            if self.check_full_export():
                self.export_full()
            self.check_tracer_consistency()
            self.print_iter_info()
Valentin Vallaeys's avatar
Valentin Vallaeys committed
421
422
        if self._timer:
            timer.printFull()
423
424
425
        self.terminate()


426