Commit 087a703b authored by Valentin Vallaeys's avatar Valentin Vallaeys Committed by Jonathan Lambrechts
Browse files

new mode splitting

parent 4d2a6feb
import slim3d
from slimPre import make_directory
output_directory = 'output'
make_directory(output_directory)
domain = slim3d.Domain('data/mesh3d.msh')
......
......@@ -4,7 +4,7 @@ from dgpy.scripts import slim_private
class Domain :
"""Create the slim3d Domain"""
def __init__(self, mesh_file_name, periodic_map_file='') :
def __init__(self, mesh_file_name, periodic_map_file='', reference_density=1027, gravity=9.81) :
"""keyword arguments:
* mesh_file_name
......@@ -14,6 +14,8 @@ class Domain :
self._topTags = topTags
self._bottomTags = bottomTags
self._slimSolver = dgpy.slim3dSolver(mesh_file_name, bottomTags, topTags, periodic_map_file)
self._density = reference_density
self._g = gravity
class Slim3d_equations :
"""Create the slim3d Equation"""
......@@ -40,7 +42,8 @@ class Slim3d_equations :
self._coriolis = None
self._atmPress = None
self._linear_density = None
self._reference_density = 1027
self._reference_density = self._domain._density
self._gravity = self._domain._g
self._initial_elevation = None
self._initial_salinity = None
self._initial_temperature = None
......@@ -49,6 +52,8 @@ class Slim3d_equations :
self._boundary_open = []
self._gotm_option_file = None
self._z0 = [0.005, 0.02]
self._linear2d = False
self._linear_solver_2D = True
def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver
......@@ -61,6 +66,12 @@ class Slim3d_equations :
if slimSolver.getSolveT() and self._vertical_diffusivity:
slimSolver.setSolveTImplicitVerticalDiffusion(True)
def set_linear_2d_equations(self, flag):
self._linear2d = flag
def set_linear_solver_2d(self, flag):
self._linear_solver_2D = flag
def set_limiter(self, flag):
slimSolver = self._slimSolver
slimSolver.setFlagUVLimiter(flag)
......@@ -125,7 +136,7 @@ class Slim3d_equations :
if z0S:
self._z0[1] = z0S
def set_linear_density(self, salinity=False, temperature=False, constant_coefficient=0, linear_coefficient=0):
def set_linear_density(self, mode, constant_coefficient=0, linear_coefficient=0):
"""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:
......@@ -142,9 +153,6 @@ class Slim3d_equations :
dgpy.Msg.Fatal('density is a linear function of either salinity, temperature or z coordinate')
self._linear_density = (mode, constant_coefficient, linear_coefficient)
def set_reference_density(self, density):
self._reference_density = density
def set_initial_elevation(self, elevation):
self._initial_elevation = elevation
......@@ -202,11 +210,10 @@ class Slim3d_equations :
class Loop:
"""Temporal solver"""
def __init__(self, equations, maximum_2d_time_step=3600, ratio_3dvs2d_time_step=1, export_time_step=3600, ratio_full_export=-1, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output"):
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"):
self._slimSolver = equations._slimSolver
self._equations = equations
self._max_2d_dt = maximum_2d_time_step
self._dt2d3d_ratio = ratio_3dvs2d_time_step
self._dt = time_step
self._export_dt = export_time_step
self._ratio_full_export = ratio_full_export
fmt = '%Y-%m-%d %H:%M:%S'
......@@ -237,6 +244,7 @@ class Loop:
self._ncExport = []
self._restart_dir = None
self._restart_ind = -1
self._timer = False
def export_uv(self, vect=True):
self._export_uv = True
......@@ -259,6 +267,9 @@ class Loop:
def restart(self, directory, index):
self._restart_dir = directory
self._restart_ind = index
def set_timer(self):
self._timer = True
def setup(self):
slim3d_private.slim3d_setup(self)
......@@ -352,6 +363,7 @@ class Loop:
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))
return uvNorm
def terminate(self):
self._timeIntegrator.terminate()
......@@ -360,6 +372,9 @@ class Loop:
def loop(self):
if not self._timeIntegrator:
self.setup()
if self._timer:
timer = dgpy.dgTimer.root()
timer.start()
while self.get_time() < self.final_time:
self.advance_one_time_step()
self.check_sanity()
......@@ -369,6 +384,8 @@ class Loop:
self.export_full()
self.check_tracer_consistency()
self.print_iter_info()
if self._timer:
timer.printFull()
self.terminate()
......
......@@ -15,15 +15,13 @@ def slim3d_setup(loop):
slimSolver.setSolveTurbulence(True)
slimSolver.setAdvectTurbulentEnergy(True)
if eq._gotm_option_file:
slimSolver.turbulenceSetupFile = 'eq._gotm_option_file'
slimSolver.turbulenceSetupFile = eq._gotm_option_file
else:
gotmFileStr = gotmOptionFile
f = open('gotmturb.nml', 'w')
f.write(gotmFileStr)
f.close()
# TODO download automatically or in SLIM
slimSolver.turbulenceSetupFile = 'gotmturb.nml'
slimSolver.turbulenceSetupFile = 'gotmturb.nml'
if (slimSolver.getSolveUVImplicitVerticalDiffusion() or slimSolver.getSolveTurbulence()) and not slimSolver.getComputeBottomFriction():
dgpy.Msg.Fatal('If Vertical implicit diffusion or using gotm, bottom friction must be true')
if eq._horizontal_viscosity == 'smagorinsky' :
......@@ -62,9 +60,9 @@ def slim3d_setup(loop):
val[:] = eq._initial_salinity[2] + xyz[:,2] * eq._initial_salinity[3]
f.SInitFunc = dgpy.functionNumpy(1, zFunc, [d.xyzOrigDof.getFunction()])
else:
dgpy.Msg.Fatal('this mode does not exist')
elif (not eq._linear_density) or (eq._linear_density[0]):
dgpy.Msg.Fatal('Initial salinity must be set if rhoFunc is not set')
dgpy.Msg.Fatal('Unknown mode for initial salinity:' + eq._initial_salinity[0])
elif (not eq._linear_density) or (eq._linear_density[0] == 'salinity'):
dgpy.Msg.Fatal('Initial salinity must be set for rhoFunc')
if eq._initial_temperature:
if eq._initial_temperature[0] == 'netcdf':
......@@ -74,33 +72,38 @@ def slim3d_setup(loop):
val[:] = eq._initial_temperature[2] + xyz[:,2] * eq._initial_temperature[3]
f.TInitFunc = dgpy.functionNumpy(1, zFunc, [d.xyzOrigDof.getFunction()])
else:
dgpy.Msg.Fatal('this mode does not exist')
elif (not eq._linear_density) or (not eq._linear_density[0]):
dgpy.Msg.Fatal('Initial temperature must be set if rhoFunc is not set')
dgpy.Msg.Fatal('Unknown mode for initial temperature:' + eq._initial_temperature[0])
elif (not eq._linear_density) or (eq._linear_density[0] == 'temperature'):
dgpy.Msg.Fatal('Initial temperature must be set for rhoFunc')
if eq._linear_density :
if eq._linear_density[0] == 'salinity':
f.rhoFunc = dgpy.linearStateEquation(d.SDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
if slimSolver.getSolveS():
f.rhoFunc = dgpy.linearStateEquation(d.SDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
else:
f.rhoFunc = dgpy.linearStateEquation(f.SInitFunc, eq._linear_density[2], eq._linear_density[1])
elif eq._linear_density[0] == 'temperature' :
f.rhoFunc = dgpy.linearStateEquation(d.TDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
if slimSolver.getSolveT():
f.rhoFunc = dgpy.linearStateEquation(d.TDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
else:
f.rhoFunc = dgpy.linearStateEquation(f.TInitFunc, eq._linear_density[2], eq._linear_density[1])
else:
f.rhoFunc = dgpy.linearStateEquation(d.zDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
if eq._vertical_diffusivity == 'gotm':
dgpy.Msg.Fatal('GOTM does not work with linear density. S and T must be set, and the automatic Jackett-McDougall state equation will be computed')
dgpy.slim3dParameters.setRho0(eq._reference_density)
dgpy.slim3dParameters.setG(eq._gravity)
eq._gravityFunc = dgpy.functionConstant(eq._gravity)
eq._areaFunc = None
if eq._horizontal_viscosity == 'constant':
f.nuh = dgpy.functionConstant(eq._hor_visc_const)
f.nuh2d = dgpy.functionConstant(eq._hor_visc_const)
elif eq._horizontal_viscosity == 'smagorinsky':
area2d = dgpy.dgDofContainer(slimSolver.groups2d, 1)
slimSolver.copy2d3d.computeArea2d( area2d)
area2d.scatter()
eq._areaFunc = dgpy.functionPrecomputedExtrusion(slimSolver.extrusion(), 3, area2d.getFunction())
f.nuh = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof.getFunction(), eq._areaFunc)
f.nuh2d = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof2d.getFunction(), eq._areaFunc)
if eq._horizontal_diffusivity == 'constant':
f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const)
......@@ -147,7 +150,7 @@ def slim3d_setup(loop):
bottomTags = eq._domain._bottomTags
topTags = eq._domain._topTags
g = dgpy.slim3dParameters.getG()
g = eq._gravity
# Open Boundaries
def uvOpenGivenEtaNumpy(cmap, val, uv, normals) :
un = uv[:,0] * normals[:,0] + uv[:,1] * normals[:,1]
......@@ -161,64 +164,64 @@ def slim3d_setup(loop):
val[:] = numpy.sqrt( (bath[:]+eta[:]) / g) * un
def transport2UVNumpy(cmap, val, uv, bath, eta) :
val[:,0] = uv[:,0] / ( bath[:] + eta[:] )
val[:,1] = uv[:,1] / ( bath[:] + eta[:] )
val[:,1] = uv[:,1] / ( bath[:] + eta[:] )
def mergeEtaUV(cmap, val, eta, uv):
val[:,0] = eta[:]
val[:,1] = uv[:,0]
val[:,2] = uv[:,1]
eta = f.eta2dFunc
uv = d.uvDof.getFunction()
uvAv2d = f.uvAv2dFunc
uvInt2d = f.uvInt2dFunc
eq._uv_open = {}
eq._uvInt2d_open = {}
eq._uvDev_open = {}
eq._uvAv2d_open = {}
eq._uvTransport_open = {}
eq._uvAv2dTransport_open = {}
eq._u_open = {}
eq._uAv2d_open = {}
eq._v_open = {}
eq._uAv2d_open = {}
eq._vAv2d_open = {}
eq._flux_section = {}
eq._flux_section = {}
eq._flux = {}
eq._flux2d = {}
eq._eta_open = {}
eq._S_open = {}
eq._T_open = {}
eq._sw2D_open = {}
for openBnd in eq._boundary_open:
eq._uvDev_open[openBnd] = dgpy.functionConstant([0.,0.])
if openBnd.u and openBnd.eta:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvInt2d_open[openBnd] = eq._uvAv2d_open[openBnd]
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = eq._uvAv2d_open[openBnd]
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.u:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvInt2d_open[openBnd] = eq._uvAv2d_open[openBnd]
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = eq._uvAv2d_open[openBnd]
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._eta_open[openBnd] = dgpy.functionNumpy(1, etaOpenGivenUVNumpy, [uvAv2d, f.bathFunc2d, eta, dgpy.function.getNormals()])
elif openBnd.eta:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uv, dgpy.function.getNormals()])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uvInt2d, dgpy.function.getNormals()])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uv, dgpy.function.getNormals()])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uvAv2d, dgpy.function.getNormals()])
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.flux:
......@@ -228,13 +231,11 @@ def slim3d_setup(loop):
eq._flux_section[openBnd] = sectionFM(0,0)
eq._flux[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups3d)
eq._flux2d[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups2d)
eq._uv_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups3d, eq._flux_section[openBnd], eq._flux[openBnd], f.bathFunc2d, eta)
eq._uvInt2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._uv_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups3d, eq._flux_section[openBnd], eq._flux[openBnd], f.bathFunc2d, eta)
eq._uvAv2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._eta_open[openBnd] = eta
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._eta_open[openBnd] = eta
if openBnd.salinity:
......@@ -247,7 +248,18 @@ def slim3d_setup(loop):
eq._T_open[openBnd] = f.TInitFunc
else:
eq._T_open[openBnd] = slim_private._load_function(openBnd.temperature, slimSolver.groups3d)
eq._sw2D_open[openBnd] = dgpy.functionNumpy(3, mergeEtaUV, [eq._eta_open[openBnd], eq._uvAv2d_open[openBnd]])
sw2DEq = e.sw2DEq
sw2DEq.setFrom3D(True)
sw2DEq.setIsLinear(eq._linear2d)
sw2DEq.setLinearSolverFrom3D(eq._linear_solver_2D)
sw2DEq.setGravity(eq._gravityFunc)
sw2DEq.setDensity(eq._reference_density)
sw2DBndWall = sw2DEq.newBoundaryWall(True)
sw2DEq.addBoundaryCondition(eq._boundary_coast, sw2DBndWall)
for openBnd in eq._boundary_open:
sw2DEq.addBoundaryCondition(openBnd.tag, sw2DEq.newOutsideValueBoundary(eq._sw2D_open[openBnd]))
horMomEq = e.horMomEq
horMomEq.setLaxFriedrichsFactor(0.0)
......@@ -257,8 +269,8 @@ def slim3d_setup(loop):
else :
horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
for openBnd in eq._boundary_open:
toReplaceHorMomOpen = dgpy.VectorFunctorConst([uv, uvInt2d, eta])
replaceByHorMomOpen = dgpy.VectorFunctorConst([eq._uv_open[openBnd], eq._uvInt2d_open[openBnd], eq._eta_open[openBnd]])
toReplaceHorMomOpen = dgpy.VectorFunctorConst([dgpy.function.getSolution(), uvAv2d, eta])
replaceByHorMomOpen = dgpy.VectorFunctorConst([eq._uvDev_open[openBnd], eq._uvAv2d_open[openBnd], eq._eta_open[openBnd]])
horMomBndOpen = horMomEq.newOutsideValueBoundaryGeneric("", toReplaceHorMomOpen, replaceByHorMomOpen)
horMomEq.addBoundaryCondition(openBnd.tag, horMomBndOpen)
horMomBndWall = horMomEq.newBoundaryWall()
......@@ -330,42 +342,25 @@ def slim3d_setup(loop):
turbEq.setBoundary0Flux('vertical_bottom')
for openBnd in eq._boundary_open:
turbEq.addBoundaryCondition(openBnd.tag, turbEq.newInFluxBoundary(f.tinyFunc))
eta2dEq = e.eta2dEq
eta2dEq.setBoundary0Flux(eq._boundary_coast)
for openBnd in eq._boundary_open:
toReplaceEtaOpen = dgpy.VectorFunctorConst([uvAv2d, dgpy.function.getSolution()])
replaceByEtaOpen = dgpy.VectorFunctorConst([eq._uvAv2d_open[openBnd], eq._eta_open[openBnd]])
etaBndOpen = eta2dEq.newOutsideValueBoundaryGeneric("", toReplaceEtaOpen, replaceByEtaOpen)
eta2dEq.addBoundaryCondition(openBnd.tag, etaBndOpen)
uv2dEq = e.uv2dEq
uv2dBndWall = uv2dEq.newBoundaryWall()
uv2dEq.addBoundaryCondition(eq._boundary_coast,uv2dBndWall)
for openBnd in eq._boundary_open:
toReplaceUV2dOpen = dgpy.VectorFunctorConst([dgpy.function.getSolution(), eta])
replaceByUV2dOpen = dgpy.VectorFunctorConst([eq._uvAv2d_open[openBnd], eq._eta_open[openBnd]])
uv2dBndOpen = uv2dEq.newOutsideValueBoundaryGeneric("", toReplaceUV2dOpen, replaceByUV2dOpen)
uv2dEq.addBoundaryCondition(openBnd.tag, uv2dBndOpen)
# .-------------------------------.
# | setup time integrator |
# '-------------------------------'
timeIntegrator = dgpy.slim3dTimeIntegratorPC(slimSolver)
rk = dgpy.dgERK(e.eta2dEq, None, dgpy.DG_ERK_EULER)
dtRK = rk.computeInvSpectralRadius(d.etaDof2d)
CFL = 0.7236
dt = CFL*dtRK
maxUV = 1.5 # max advective velocity
c2d = slim_private.math.sqrt(9.81*20.)
dt1 = c2d/(c2d+maxUV) * dt
dt2 = min(dt1, loop._max_2d_dt)
n_iter_per_export = numpy.ceil(loop._export_dt/(dt2*loop._dt2d3d_ratio))
dt = loop._export_dt/(n_iter_per_export*loop._dt2d3d_ratio)
# rk = dgpy.dgERK(e.eta2dEq, None, dgpy.DG_ERK_EULER)
# dtRK = rk.computeInvSpectralRadius(d.etaDof2d)
# CFL = 0.7236
# dt = CFL*dtRK
# maxUV = 1.5 # max advective velocity
# c2d = slim_private.math.sqrt(9.81*20.)
# dt1 = c2d/(c2d+maxUV) * dt
# dt2 = min(dt1, loop._max_2d_dt)
# n_iter_per_export = numpy.ceil(loop._export_dt/(dt2*loop._dt2d3d_ratio))
# dt = loop._export_dt/(n_iter_per_export*loop._dt2d3d_ratio)
n_iter_per_export = numpy.ceil(loop._export_dt/loop._dt)
dt = loop._export_dt/n_iter_per_export
if dgpy.Msg.GetCommRank() == 0 :
dgpy.Msg.Info('Runge-Kutta dt :' + str(dtRK) + ' CFL:' + str(CFL) + ' dt:' + str(dt))
dgpy.Msg.Info('dt_user :' + str(loop._dt) + ' dt:' + str(dt))
# .-------------------------------.
# | setup export options |
......@@ -375,7 +370,6 @@ def slim3d_setup(loop):
#timeIntegrator.setFullExportAtEveryTimeStep(False)
timeIntegrator.setFullExportIntervalRatio(loop._ratio_full_export )
timeIntegrator.setTimeStep(dt)
timeIntegrator.set3dTimeStepRatio(loop._dt2d3d_ratio)
timeIntegrator.setInitTime(loop.initial_time)
timeIntegrator.setEndTime(loop.final_time)
loop._timeIntegrator = timeIntegrator
......@@ -468,7 +462,7 @@ gotmOptionFile = \
/
!-------------------------------------------------------------------------------
! What turbulence parameters have been described?
! What turbulence parameters have been described?
!
! cm0_fix -> value of cm0 for turb_method=2
! Prandtl0_fix -> value of the turbulent Prandtl-number for stab_method=1-4
......
......@@ -229,7 +229,7 @@ def _new_default_linear_system(domain, equation) :
except :
dgpy.Msg.Fatal("No linear system available")
return sys
def _findTopAndBottomTags(mesh_file_name):
bottomTags = []
topTags = []
......
set(SRC
dgConservationLawShallowWater1d.cpp
dgConservationLawShallowWater2d.cpp
dgConservationLawShallowWaterAbsLayer2d.cpp
dgConservationLawShallowWaterTracer1d.cpp
dgConservationLawShallowWaterTracer2d.cpp
dgConservationLawShallowWaterEta.cpp
dgConservationLawShallowWater2d.cpp
slimMovingBathWettingDrying.cpp
)
......
......@@ -88,19 +88,6 @@ class dgConservationLawShallowWater2d::clipToPhysics : public function {
}*/
};
class dgConservationLawShallowWater2d::maxDiffusiveSpeed: public function {
fullMatrix<double> nu;
public:
/* maxConvectiveSpeed (const functor *nuF):function(3){
setArgument(nu, nuF);
}*/
void call(dataCacheMap *m, fullMatrix<double> &val) {
/*for(int i=0; i< val.size1(); i++) {
val(i,0) = sqrt(bath(i,0)*_g);
}*/
}
};
class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
fullMatrix<double> bath, sol, _xyz, _g;
bool _linear, _filterLinear;
......@@ -146,36 +133,25 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
}
};
class dgConservationLawShallowWater2d::gradPsiTerm: public function {
fullMatrix<double> bath, sol, diff, _xyz, movingBathFactor, _g;
bool _linear, _isDiffusive;
double _R;
public:
gradPsiTerm (const functor *bathymetryF, const functor *diffusiveFluxF, bool linear, double R, const functor *movingBathFactorF, const functor *xyz, const functor *g):function(9){
setArgument(bath,bathymetryF);
setArgument(sol,function::getSolution());
_isDiffusive = diffusiveFluxF;
if (_isDiffusive) {
setArgument(diff,diffusiveFluxF);
}
_linear = linear;
_R = R;
setArgument(_g, g);
if(_R>0) {
setArgument(_xyz, xyz ? xyz : function::getCoordinates());
}
setArgument(movingBathFactor, movingBathFactorF);
}
void call(dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = val.size1();
void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &bath = cm.get(*_bathymetry);
const fullMatrix<double> &sol = cm.get(*function::getSolution());
const fullMatrix<double> *diff = _diffusiveFlux ? &cm.get(*_diffusiveFlux) : NULL;
const fullMatrix<double> &g = cm.get(*_g);
const fullMatrix<double> *xyz = _R ? &cm.get(*function::getCoordinates()) : NULL;
const fullMatrix<double> &movingBathFactor = cm.get(*_movingBathFactor);
const fullMatrix<double> &atmPress = cm.get(*_atmPress);
const fullMatrix<double> &rhoSurf = cm.get(*_rhoSurf);
const fullMatrix<double> *etaOld = _linearSolverFrom3D ? &cm.get(*_etaOld) : NULL;
size_t nQP = cm.nEvaluationPoint();
val.resize(nQP, 9, false);
double J = 1.0;
for(size_t i=0; i< nQP; i++) {
if (_R > 0){
double x = _xyz(i,0);
double y = _xyz(i,1);
double x = (*xyz)(i,0);
double y = (*xyz)(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
J = 1.0;
}
double H = bath(i,0);
double eta = sol(i,0);
......@@ -184,146 +160,115 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
double A = movingBathFactor(i,0);//equals 1 if no WD
if (!_linear)
H+=eta;
else if (_linearSolverFrom3D)
H+=(*etaOld)(i,0);
// flux_x
val(i,0) = H*u * J/A;
val(i,1) = _g(i,0)*eta*J;
val(i,1) = ( g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density )*J;
val(i,2) = 0;
// flux_y
val(i,3) = H*v * J/A;
val(i,4) = 0;
val(i,5) = _g(i,0)*eta*J;
if (!_linear) {
val(i,5) = ( g(i,0)*(1+rhoSurf(i,0)/_density) *eta + atmPress(i,0)/_density )*J;
if (!_linear && !_from3D) {
val(i,1) += u*u*J;
val(i,2) += u*v*J;
val(i,4) += v*u*J;
val(i,5) += v*v*J;
}
if (_isDiffusive) {
val(i,1) += diff(i,1);
val(i,2) += diff(i,2);
val(i,4) += diff(i,4);
val(i,5) += diff(i,5);
if (diff) {
val(i,1) += (*diff)(i,1);
val(i,2) += (*diff)(i,2);
val(i,4) += (*diff)(i,4);
val(i,5) += (*diff)(i,5);
}
// flux_z
val(i,6) = 0;
val(i,7) = 0;
val(i,8) = 0;
}
}
};
class dgConservationLawShallowWater2d::source: public function {
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _windStress, _nu, _xyz;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _nudgingCoeff, _nudgingVel, _g;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer, _nudgingVelIsTransport;
double _R, _density;
public :
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF, const functor *windStressF,
const functor *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R, const functor *gravityF, double density,
const functor *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *nudgingCoeffF, const functor *nudgingVelF, bool nudgingVelIsTransport) :function (3){
_linear = linear;
_R = R;
setArgument(_g, gravityF);