Commit fef24e83 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

testcase swe3dMovinMesh works now with the interface

parent 6ca3412d
Pipeline #1189 passed with stage
in 33 minutes and 36 seconds
......@@ -69,6 +69,9 @@ class Slim3d_equations :
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_conservative_ale(self, flag):
self._slimSolver.setUseConservativeALE(flag)
def set_horizontal_viscosity(self, mode, constant=1, factor=0.02, maximum=1e9):
if (mode != "smagorinsky") and (mode != "constant"):
dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode)
......@@ -208,13 +211,13 @@ class Loop:
self._ratio_full_export = ratio_full_export
fmt = '%Y-%m-%d %H:%M:%S'
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])
self.initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
self._initial_date = date
if 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])
self.final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
else:
self._final_time = self._initial_time + final_time
self.final_time = self.initial_time + final_time
self._odir = output_directory
self._export_uv = False
self._export_w = False
......@@ -307,39 +310,66 @@ class Loop:
else:
self._ncWriter2d = None
self._ncWriter3d = None
return self._timeIntegrator
def loop(self):
if not self._timeIntegrator:
self._timeIntegrator = self.setup()
timeIntegrator = self._timeIntegrator
slimSolver = self._slimSolver
d = slimSolver.getDofs()
if self._restart_ind > -1:
timeIntegrator.reStart(self._restart_dir, self._restart_ind)
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
for e in self._export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self._initial_time)
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)
tic = slim_private.time.time()
while timeIntegrator.getTime() < timeIntegrator.getEndTime() :
timeIntegrator.advanceOneTimeStep()
if timeIntegrator.checkExport() :
timeIntegrator.checkSanity()
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
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)
if timeIntegrator.checkFullExport() :
timeIntegrator.exportAll()
timeIntegrator.checkTracerConsistency()
uvNorm = d.uvDof.norm()
if dgpy.Msg.GetCommRank() == 0 : dgpy.Msg.Info("%5i [ %s ] normuv %6.12e clock %.1fs" % (timeIntegrator.getIter(),slim_private.time.strftime("%d %b %Y %H:%M:%S", slim_private.time.gmtime(timeIntegrator.getTime())), uvNorm, slim_private.time.time() - tic))
timeIntegrator.terminate(0)
self._ncWriter3d.exportUGRID(dof, filename, timeIntegrator.getTime() - self.initial_time, timeIntegrator.getIter(),field)
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 export_fields(self):
timeIntegrator = self._timeIntegrator
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
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)
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))
def terminate(self):
self._timeIntegrator.terminate()
def loop(self):
if not self._timeIntegrator:
self.setup()
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()
self.terminate()
......@@ -376,8 +376,8 @@ def slim3d_setup(loop):
timeIntegrator.setFullExportIntervalRatio(loop._ratio_full_export )
timeIntegrator.setTimeStep(dt)
timeIntegrator.set3dTimeStepRatio(loop._dt2d3d_ratio)
timeIntegrator.setInitTime(loop._initial_time)
timeIntegrator.setEndTime(loop._final_time)
timeIntegrator.setInitTime(loop.initial_time)
timeIntegrator.setEndTime(loop.final_time)
loop._timeIntegrator = timeIntegrator
gotmOptionFile = \
......
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from dgpy import *
from shutil import copyfile
import os, math
from dgpy.scripts import Common
# .-------------------------------.
# | Description |
# '-------------------------------'
......@@ -20,227 +13,91 @@ relVolTolerance = 2e-14
relTracerMassTolerance = 1e-14
TracerDevTolerance = 1.2e-4 #arbitrary
# .-------------------------------.
# | generate mesh |
# '-------------------------------'
Common.genMesh("box2d.geo", 2)
nbLayers = 4
def getLayersSigma (e,v) :
h = 50 *(1-v[0]/5000.)/2. + 50 *(1+v[0]/5000.)/2.
#h = 20 *(1-x/5000.)/2. + 5 *(1+x/5000.)/2.
#h = 50 * (math.sin(6*3.14*(1-x/5000.)/2.)+1) + 80
return [z * h /nbLayers for z in range(nbLayers + 1)]
from dgpy.scripts.extrude import *
extrude(mesh('box2d.msh'), getLayersSigma).write('box.msh')
# .-------------------------------.
# | mesh and output directory |
# '-------------------------------'
meshFile = 'box.msh'
odir = 'movingMesh/'
# create outputdir if !exists
if Msg.GetCommRank()==0 and not os.path.exists(odir) :
os.mkdir(odir)
Msg.Barrier()
copyfile(meshFile,odir + '/mesh.msh')
# .-------------------------------.
# | create solver |
# '-------------------------------'
s = slim3dSolver(meshFile, ["bottom"], ["top"])
# options
s.setSolveS(True)
s.setSolveT(True)
s.setSolveSImplicitVerticalDiffusion(False)
s.setSolveTImplicitVerticalDiffusion(False)
s.setSolveUVImplicitVerticalDiffusion(False)
s.setSolveUVImplicitVerticalAdvection(False)
s.setSolveSImplicitVerticalAdvection(True)
s.setSolveTImplicitVerticalAdvection(False)
s.setSolveTurbulence(False)
s.setAdvectTurbulentEnergy(False)
s.setComputeBottomFriction(False)
s.exportDirectory = odir
s.turbulenceSetupFile = "gotmEstuaryTurb.nml"
s.setUseConservativeALE(True)
# .-------------------------------.
# | set functions |
# '-------------------------------'
f = s.functions
def SF(val,xyz) :
for i in range(val.size1()) :
#val.set(i,0, 5*(-xyz(i,2)/100) + 3 *(1+xyz(i,2)/100) )
#val.set(i,0, -42 *(1-xyz(i,0)/5000)/2. + 4 *(1+xyz(i,0)/5000)/2. )
val.set(i,0, 4.00 )
f.SInitFunc = functionPython(1,SF,[f.xyzFunc3d])
def TF(val,xyz) :
for i in range(val.size1()) :
val.set(i,0, 5*(-xyz(i,2)/50) + 3 *(1+xyz(i,2)/50) )
f.TInitFunc = functionPython(1,TF,[f.xyzFunc3d])
slim3dParameters.setRho0(1002.83351)
etaAmp = 0.1
def etaF(val,xyz) :
for i in range(val.size1()) :
val.set(i,0, etaAmp*math.exp(- (xyz(i,0)/2000)**2) )
f.etaInitFunc = functionPython(1,etaF,[f.xyzFunc2d])
f.z0BFunc = functionConstant(0.005)
f.z0SFunc = functionConstant(0.02)
# .-------------------------------.
# | create equations |
# '-------------------------------'
e = s.createEquations()
d = s.getDofs()
eta2d = d.etaDof2d.getFunction()
uvAv2d = d.uvAvDof2d.getFunction()
eta2dCG = d.etaDof2dCG.getFunction()
horMomEq = e.horMomEq
horMomBndZero = horMomEq.new0FluxBoundary()
horMomBndSymm = horMomEq.newSymmetryBoundary('')
horMomBndWall = horMomEq.newBoundaryWall()
horMomEq.addBoundaryCondition('bottom',horMomBndWall) # must be wall for bath!
horMomEq.addBoundaryCondition('sideL',horMomBndWall)
horMomEq.addBoundaryCondition('river',horMomBndWall)
horMomEq.addBoundaryCondition('sideR',horMomBndWall)
horMomEq.addBoundaryCondition('sea',horMomBndWall)
horMomEq.addBoundaryCondition('top',horMomBndSymm) # zero for nonconst tracers!
if s.getSolveUVImplicitVerticalDiffusion() or s.getSolveUVImplicitVerticalAdvection():
vertMomUEq = e.vertMomUEq
vertMomUBndZero = vertMomUEq.new0FluxBoundary()
vertMomUBndSymm = vertMomUEq.newSymmetryBoundary('')
if s.getSolveUVImplicitVerticalDiffusion():
vertMomUBndBottom = vertMomUEq.newBoundaryBottom(f.bottomFriction2dPrecompFunc,f.zBotDist2dPrecompFunc)
vertMomUEq.addBoundaryCondition('bottom',vertMomUBndBottom)
else :
vertMomUEq.addBoundaryCondition('bottom',vertMomUBndZero)
vertMomUEq.addBoundaryCondition('top',vertMomUBndZero)
vertMomUEq.addBoundaryCondition('sideL',vertMomUBndZero)
vertMomUEq.addBoundaryCondition('river',vertMomUBndZero)
vertMomUEq.addBoundaryCondition('sideR',vertMomUBndZero)
vertMomUEq.addBoundaryCondition('sea',vertMomUBndZero)
wEq = e.wEq
wBndZero = wEq.new0FluxBoundary()
wBndSymm = wEq.newSymmetryBoundary('')
wEq.addBoundaryCondition('sideL',wBndZero)
wEq.addBoundaryCondition('river',wBndZero)
wEq.addBoundaryCondition('sideR',wBndZero)
wEq.addBoundaryCondition('sea',wBndZero)
wEq.addBoundaryCondition('top',wBndSymm) # must be symm!!
wEq.addBoundaryCondition('bottom',wBndZero)
if s.getSolveS() :
SEq = e.SEq
SBndZero = SEq.new0FluxBoundary()
SEq.addBoundaryCondition(['sideL','river','sideR','sea','top','bottom'],SBndZero)
if s.getSolveSImplicitVerticalDiffusion() or s.getSolveSImplicitVerticalAdvection():
SDiffEq = e.SDiffEq
SDiffBndZero = SDiffEq.new0FluxBoundary()
SDiffEq.addBoundaryCondition('bottom',SDiffBndZero)
SDiffEq.addBoundaryCondition('top',SDiffBndZero)
#SDiffEq.addBoundaryCondition('sideL',SDiffBndZero)
#SDiffEq.addBoundaryCondition('river',SDiffBndZero)
#SDiffEq.addBoundaryCondition('sideR',SDiffBndZero)
#SDiffEq.addBoundaryCondition('sea',SDiffBndZero)
e.newRGradEq.setBoundarySymmetry(['sideL','river','sideR','sea','top','bottom'])
if s.getSolveT() :
TEq = e.TEq
TBndZero = TEq.new0FluxBoundary()
TEq.addBoundaryCondition(['sideL','river','sideR','sea','top','bottom'],TBndZero)
if s.getSolveTImplicitVerticalDiffusion() or s.getSolveTImplicitVerticalAdvection():
TDiffEq = e.TDiffEq
TDiffBndZero = TDiffEq.new0FluxBoundary()
TDiffEq.addBoundaryCondition('bottom',TDiffBndZero)
TDiffEq.addBoundaryCondition('top',TDiffBndZero)
#TDiffEq.addBoundaryCondition('sideL',TDiffBndZero)
#TDiffEq.addBoundaryCondition('river',TDiffBndZero)
#TDiffEq.addBoundaryCondition('sideR',TDiffBndZero)
#TDiffEq.addBoundaryCondition('sea',TDiffBndZero)
if s.getAdvectTurbulentEnergy() :
tkeAdvEq = e.tkeAdvEq
tkeAdvBndZero = tkeAdvEq.new0FluxBoundary()
tkeAdvEq.addBoundaryCondition(['sideL','river','sideR','sea','top','bottom'],tkeAdvBndZero)
epsAdvEq = e.epsAdvEq
epsAdvBndZero = epsAdvEq.new0FluxBoundary()
epsAdvEq.addBoundaryCondition(['sideL','river','sideR','sea','top','bottom'],epsAdvBndZero)
eta2dEq = e.eta2dEq
etaBndZero = eta2dEq.new0FluxBoundary()
eta2dEq.addBoundaryCondition(['sideL','river','sideR','sea'],etaBndZero)
uv2dEq = e.uv2dEq
uv2dBndWall = uv2dEq.newBoundaryWall()
uv2dEq.addBoundaryCondition(['sideL','river','sideR','sea'],uv2dBndWall)
# .-------------------------------.
# | setup time integrator |
# '-------------------------------'
t = slim3dTimeIntegratorPC(s)
rk = dgERK(e.eta2dEq, None, DG_ERK_EULER)
dtRK = rk.computeInvSpectralRadius(d.etaDof2d)
CFL = 0.7263
dt = CFL*dtRK
#dt = 0.2 # for convergence
if Msg.GetCommRank() == 0 :
print('Runge-Kutta dt :' + str(dtRK) + ' CFL:' + str(CFL) + ' dt:' + str(dt))
M = 10
t.setTimeStep(dt)
t.set3dTimeStepRatio(M)
t.setExportInterval(200)
t.setEndTime(8e3)
#t.setExportAtEveryTimeStep(True)
## .-------------------------------.
# | setup export options |
# '-------------------------------'
### Generate Mesh ###
from dgpy.scripts import Common
Common.genMesh("box2d.geo", 2)
# list of fields to export
export = []
export.append(dgIdxExporter(d.SDof, odir + "/salinity"))
export.append(dgIdxExporter(d.uvDof, odir + "/uv", True))
export.append(dgIdxExporter(d.etaDof2d, odir + "/eta"))
export.append(dgIdxExporter(d.etaDof2dCG, odir + "/etaCG"))
export.append(dgIdxExporter(d.zDof, odir + "/z"))
#initial export
#for e in export :
# e.exportIdx(t.getExportCount(), t.getTime())
tic = time.time()
while t.getTime() < t.getEndTime() :
t.advanceOneTimeStep()
t.checkSanity()
if t.checkExport() :
Msg.Info("export %i" % t.getExportCount())
for e in export :
e.exportIdx(t.getExportCount(), t.getTime())
Msg.Info("%5i time : %g normuv %6.12e clock %.1fs" % (t.getIter(), t.getTime(), d.uvDof.norm(), time.time() - tic))
[relVolErr, relSErr, relTErr] = t.checkMassConservation()
[absSDev, absTDev] = t.checkTracerConsistency()
if abs(relVolErr) > relVolTolerance :
print ('Volume conservation is broken!!')
Msg.Exit(-1)
if abs(absSDev) > TracerDevTolerance :
print ('Tracer consistency is broken!!')
Msg.Exit(-1)
if max(abs(relSErr), abs(relTErr)) > relTracerMassTolerance :
print ('Tracer mass conservation is broken!!')
Msg.Exit(-1)
if t.checkFullExport() :
t.exportAll()
t.terminate(0)
print('mesh done')
### Prepro ###
import slimPre
mesh_file = 'box2d.msh'
nPart = slimPre.dgpy.Msg.GetCommSize()
if nPart > 1:
if slimPre.dgpy.Msg.GetCommRank() == 0:
slimPre.dgpy.dgMeshPartition(mesh_file, nPart)
slimPre.dgpy.Msg.Barrier()
mesh_file = mesh_file[:-4] + '_' + str(nPart) + '.msh'
mesh2d = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh2d)
xyz = region_global.coordinates
bath = 50 *(1-xyz[:,0]/5000.)/2. + 50 *(1+xyz[:,0]/5000.)/2.
slimPre.write_file('bath_2d.nc', region=region_global, time=None, data=[('bath',bath)])
slimPre.extrude(mesh_file, ('bath_2d.nc','bath'), nb_layers = 4, mesh_file_name_out='box.msh')
mesh_file = 'box.msh'
mesh = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh, 'top')
xyz = region_global.coordinates
import numpy
etaInit = 0.1 * numpy.exp(-(xyz[:,0]/2000)**2)
slimPre.write_file('etaInit.nc', region=region_global, time=None, data=[('eta',etaInit)])
print('prepro done')
### Prepro ###
import slim3d
domain = slim3d.Domain('box.msh')
equations = slim3d.Slim3d_equations(domain, temperature=True, salinity=True)
equations.set_implicit_vertical_diffusion(False)
equations.set_conservative_ale(True)
equations.set_boundary_coast(['sideL', 'sideR', 'river', 'sea'])
equations.set_limiter(True)
equations.set_initial_temperature('vertical_gradient', None, 3, -2/50.) #temporary. Change size after
equations.set_initial_salinity('vertical_gradient', None, 4, 0.)
equations.set_initial_elevation(('etaInit.nc', 'eta'))
time_loop = slim3d.Loop(equations,
maximum_2d_time_step=1,
ratio_3dvs2d_time_step=10,
export_time_step=200,
final_time=8e3,
output_directory='output')
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv(False)
time_loop.export_w()
time_loop.export_uv2d(False)
#time_loop.loop()
time_loop.setup()
while time_loop.get_time() < time_loop.final_time:
time_loop.advance_one_time_step()
time_loop.check_sanity()
#print('******************')
#print(time_loop.equations._slimSolver.getUseConservativeALE())
if time_loop.check_export():
time_loop.print_iter_info()
[relVolErr, relSErr, relTErr] = time_loop.check_mass_conservation()
[absSDev, absTDev] = time_loop.check_tracer_consistency()
time_loop.export_fields()
if abs(relVolErr) > relVolTolerance :
print ('Volume conservation is broken!!')
slim3d.dgpy.Msg.Exit(-1)
if abs(absSDev) > TracerDevTolerance :
print ('Tracer consistency is broken!!')
slim3d.dgpy.Msg.Exit(-1)
if max(abs(relSErr), abs(relTErr)) > relTracerMassTolerance :
print ('Tracer mass conservation is broken!!')
slim3d.dgpy.Msg.Exit(-1)
time_loop.terminate()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment