Commit 113a8aa0 authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

update test cases

parent 5a65f5e8
...@@ -363,10 +363,10 @@ def slim3d_setup(loop): ...@@ -363,10 +363,10 @@ def slim3d_setup(loop):
# dt2 = min(dt1, loop._max_2d_dt) # dt2 = min(dt1, loop._max_2d_dt)
# n_iter_per_export = numpy.ceil(loop._export_dt/(dt2*loop._dt2d3d_ratio)) # n_iter_per_export = numpy.ceil(loop._export_dt/(dt2*loop._dt2d3d_ratio))
# dt = loop._export_dt/(n_iter_per_export*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) n_iter_per_export = int(numpy.ceil(loop._export_dt/loop._dt))
dt = loop._export_dt/n_iter_per_export dt = loop._export_dt/n_iter_per_export
if dgpy.Msg.GetCommRank() == 0 : if dgpy.Msg.GetCommRank() == 0 :
dgpy.Msg.Info('dt_user :' + str(loop._dt) + ' dt:' + str(dt)) dgpy.Msg.Info('time step: ' + str(dt) + ', number of iterations between export: ' + str(n_iter_per_export))
# .-------------------------------. # .-------------------------------.
# | setup export options | # | setup export options |
......
...@@ -168,7 +168,6 @@ horMomEq = e.horMomEq ...@@ -168,7 +168,6 @@ horMomEq = e.horMomEq
#horMomEq.setIsSpherical #horMomEq.setIsSpherical
horMomEq.setLaxFriedrichsFactor(0.0) horMomEq.setLaxFriedrichsFactor(0.0)
#horMomEq.setNuV(f.nuVFunc) #horMomEq.setNuV(f.nuVFunc)
horMomBndSides = horMomEq.newBoundaryWallLinear(0.0)
horMomBndZero = horMomEq.new0FluxBoundary() horMomBndZero = horMomEq.new0FluxBoundary()
horMomBndSymm = horMomEq.newSymmetryBoundary('') horMomBndSymm = horMomEq.newSymmetryBoundary('')
#horMomBndTop = horMomEq.newBoundarySurface(fWind) #horMomBndTop = horMomEq.newBoundarySurface(fWind)
...@@ -224,43 +223,14 @@ if s.getSolveT() : ...@@ -224,43 +223,14 @@ if s.getSolveT() :
TEq.addBoundaryCondition('top_Surface',TBndZero) TEq.addBoundaryCondition('top_Surface',TBndZero)
TEq.addBoundaryCondition('bottom_Surface',TBndZero) TEq.addBoundaryCondition('bottom_Surface',TBndZero)
eta2dEq = e.eta2dEq sw2DWall = e.sw2DEq.newBoundaryWall()
#eta2dEq.setIsSpherical(R) e.sw2DEq.addBoundaryCondition("Wall", sw2DWall)
etaBndZero = eta2dEq.new0FluxBoundary()
eta2dEq.addBoundaryCondition('Wall',etaBndZero)
#eta2dEq.addBoundaryCondition('cut',etaBndZero)
#eta2dEq.addBoundaryCondition('paste',etaBndZero)
uv2dEq = e.uv2dEq
#uv2dEq.setIsSpherical(R)
#uv2dEq.setWindStress(fWind)
uv2dBndWall = uv2dEq.newBoundaryWall()
uv2dEq.addBoundaryCondition('Wall',uv2dBndWall)
#uv2dEq.addBoundaryCondition('cut',uv2dBndWall)
#uv2dEq.addBoundaryCondition('paste',uv2dBndWall)
# .-------------------------------. # .-------------------------------.
# | setup time integrator | # | setup time integrator |
# '-------------------------------' # '-------------------------------'
timeIntegrator = slim3dTimeIntegratorPC(s) timeIntegrator = slim3dTimeIntegratorPC(s)
rk = dgERK(e.eta2dEq, None, DG_ERK_EULER)
dtRK = rk.computeInvSpectralRadius(d.etaDof2d)
CFL = 0.7236
dt = CFL*dtRK
maxUV = 1.5 # max advective velocity
c2d = math.sqrt(9.81*20.)
dt = c2d/(c2d+maxUV) * dt
if Msg.GetCommRank() == 0 :
Msg.Info('Runge-Kutta dt :' + str(dtRK) + ' CFL:' + str(CFL) + ' dt:' + str(dt))
#dt3d = dtRK * math.sqrt(9.81*20) / maxUV
#M = int(0.72*dt3d/dt)
M = 20
if Msg.GetCommRank() == 0 :
Msg.Info('Mode split ratio M :' + str(M))
# .-------------------------------. # .-------------------------------.
# | setup export options | # | setup export options |
...@@ -268,7 +238,6 @@ if Msg.GetCommRank() == 0 : ...@@ -268,7 +238,6 @@ if Msg.GetCommRank() == 0 :
timeIntegrator.setExportInterval(1e4) timeIntegrator.setExportInterval(1e4)
timeIntegrator.setExportAtEveryTimeStep(True) timeIntegrator.setExportAtEveryTimeStep(True)
timeIntegrator.setFullExportAtEveryTimeStep(False) timeIntegrator.setFullExportAtEveryTimeStep(False)
timeIntegrator.setFullExportIntervalRatio(10)
# list of fields to export # list of fields to export
export = [] export = []
...@@ -279,13 +248,14 @@ export.append(dgIdxExporter(d.etaDof2d, odir + "/eta")) ...@@ -279,13 +248,14 @@ export.append(dgIdxExporter(d.etaDof2d, odir + "/eta"))
for e in export : for e in export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime()) e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime())
dt = 10000
# .--------------------. # .--------------------.
# | main loop | # | main loop |
# '--------------------' # '--------------------'
timeIntegrator.setTimeStep(dt) timeIntegrator.setTimeStep(dt)
timeIntegrator.set3dTimeStepRatio(M)
timeIntegrator.setInitTime(False) timeIntegrator.setInitTime(False)
timeIntegrator.setEndTime(dt*M) timeIntegrator.setEndTime(dt)
tic = time.time() tic = time.time()
timeIntegrator.advanceOneTimeStep() timeIntegrator.advanceOneTimeStep()
timeIntegrator.checkSanity() timeIntegrator.checkSanity()
......
...@@ -64,8 +64,7 @@ equations.set_boundary_coast(['Wall']) ...@@ -64,8 +64,7 @@ equations.set_boundary_coast(['Wall'])
equations.set_linear_density("z_coordinate", 0., 0.) equations.set_linear_density("z_coordinate", 0., 0.)
time_loop = slim3d.Loop(equations, time_loop = slim3d.Loop(equations,
maximum_2d_time_step=100, time_step=400,
ratio_3dvs2d_time_step=20,
export_time_step=2000, export_time_step=2000,
final_time=4000, final_time=4000,
output_directory='output') output_directory='output')
......
...@@ -13,7 +13,6 @@ oprefix = 'mixed_layer_depth' ...@@ -13,7 +13,6 @@ oprefix = 'mixed_layer_depth'
plotFigures = False plotFigures = False
saveFigures = False saveFigures = False
### Generate Mesh ### ### Generate Mesh ###
from dgpy.scripts import Common from dgpy.scripts import Common
meshStem = 'katoPhillips' meshStem = 'katoPhillips'
...@@ -51,7 +50,7 @@ print('prepro done') ...@@ -51,7 +50,7 @@ print('prepro done')
### Prepro ### ### Prepro ###
import slim3d import slim3d
domain = slim3d.Domain(mesh_file_3d, mapFilename) domain = slim3d.Domain(mesh_file_3d, mapFilename, reference_density=1027)
equations = slim3d.Slim3d_equations(domain, temperature=False, salinity=True) equations = slim3d.Slim3d_equations(domain, temperature=False, salinity=True)
equations.set_implicit_vertical_diffusion(True) equations.set_implicit_vertical_diffusion(True)
equations.set_vertical_viscosity('gotm') equations.set_vertical_viscosity('gotm')
...@@ -64,18 +63,14 @@ equations.set_limiter(False) ...@@ -64,18 +63,14 @@ equations.set_limiter(False)
slim3d.dgpy.slim3dParameters.setNuv0( 1.3e-6) slim3d.dgpy.slim3dParameters.setNuv0( 1.3e-6)
slim3d.dgpy.slim3dParameters.setKapv0T( 1.4e-7) slim3d.dgpy.slim3dParameters.setKapv0T( 1.4e-7)
slim3d.dgpy.slim3dParameters.setKapv0S( 1.1e-9) slim3d.dgpy.slim3dParameters.setKapv0S( 1.1e-9)
slim3d.dgpy.slim3dParameters.setRho0(1027)
equations.set_wind_stress('stress', ('windStress.nc','u'), ('windStress.nc','v') ) equations.set_wind_stress('stress', ('windStress.nc','u'), ('windStress.nc','v') )
equations.set_boundary_coast(['side1','side3']) equations.set_boundary_coast(['side1','side3'])
equations.set_initial_salinity('vertical_gradient', None, 0, -13.3e-3) equations.set_initial_salinity('vertical_gradient', None, 0, -13.3e-3)
equations.set_initial_temperature('vertical_gradient', None, 10) equations.set_initial_temperature('vertical_gradient', None, 10)
time_loop = slim3d.Loop(equations, time_loop = slim3d.Loop(equations,
maximum_2d_time_step=12, time_step=60,
ratio_3dvs2d_time_step=5,
export_time_step=100*60, export_time_step=100*60,
final_time=30*3600, final_time=30*3600,
output_directory='output') output_directory='output')
...@@ -103,7 +98,6 @@ import dgpy ...@@ -103,7 +98,6 @@ import dgpy
s = equations._slimSolver s = equations._slimSolver
d = s.getDofs(True) d = s.getDofs(True)
maxExportCount = int(time_loop.final_time/time_loop._export_dt)+1 maxExportCount = int(time_loop.final_time/time_loop._export_dt)+1
export_counts = range(1,maxExportCount) export_counts = range(1,maxExportCount)
nbTimes = len(export_counts) nbTimes = len(export_counts)
......
...@@ -2,14 +2,13 @@ ...@@ -2,14 +2,13 @@
# | Description | # | Description |
# '-------------------------------' # '-------------------------------'
# Lake at Rest on a smooth non-constant bathymetry # Lake at Rest on a smooth non-constant bathymetry
# 20 days, uv slowly grows (!) but remains small
# #
# .-------------------------------. # .-------------------------------.
# | test case tolerances | # | test case tolerances |
# '-------------------------------' # '-------------------------------'
normUTolerance = 5e-10 normUTolerance = 5e-10
### Generate Mesh ### ### Generate Mesh ###
from dgpy.scripts import Common from dgpy.scripts import Common
mesh_file = Common.genMesh("square.geo", 2) mesh_file = Common.genMesh("square.geo", 2)
...@@ -48,16 +47,13 @@ equations = slim3d.Slim3d_equations(domain, temperature=True, salinity=False) ...@@ -48,16 +47,13 @@ equations = slim3d.Slim3d_equations(domain, temperature=True, salinity=False)
equations.set_boundary_coast(['Wall']) equations.set_boundary_coast(['Wall'])
equations.set_limiter(False) equations.set_limiter(False)
equations.set_initial_temperature('vertical_gradient', None, 25, 8.2e-3) equations.set_initial_temperature('vertical_gradient', None, 25, 8.2e-3)
equations.set_linear_density('temperature', 20, -.1) equations.set_linear_density('temperature', 20, -.1)
time_loop = slim3d.Loop(equations, time_loop = slim3d.Loop(equations,
maximum_2d_time_step=10000, #let the model find it time_step=3600,
ratio_3dvs2d_time_step=10, export_time_step=86400,
export_time_step=86400*2, final_time=86400*20,
final_time=86400*2,
output_directory='output') output_directory='output')
time_loop.export_elevation() time_loop.export_elevation()
...@@ -66,14 +62,13 @@ time_loop.export_uv(False) ...@@ -66,14 +62,13 @@ time_loop.export_uv(False)
time_loop.export_w() time_loop.export_w()
time_loop.export_uv2d(False) time_loop.export_uv2d(False)
#time_loop.loop()
time_loop.setup() time_loop.setup()
it = 0 it = 0
while time_loop.get_time() < time_loop.final_time: while time_loop.get_time() < time_loop.final_time:
time_loop.advance_one_time_step() time_loop.advance_one_time_step()
it += 1 it += 1
time_loop.check_sanity() time_loop.check_sanity()
if it % 20 == 0: if it % 24 == 0:
normUV = time_loop.print_iter_info() normUV = time_loop.print_iter_info()
if normUV > normUTolerance : if normUV > normUTolerance :
print ('rest conservation is broken!!') print ('rest conservation is broken!!')
......
...@@ -62,7 +62,6 @@ equations.set_initial_temperature('vertical_gradient', None, 3, -2/50.) #tempora ...@@ -62,7 +62,6 @@ equations.set_initial_temperature('vertical_gradient', None, 3, -2/50.) #tempora
equations.set_initial_salinity('vertical_gradient', None, 4, 0.) equations.set_initial_salinity('vertical_gradient', None, 4, 0.)
equations.set_initial_elevation((pre_data_dir+'etaInit.nc', 'eta')) equations.set_initial_elevation((pre_data_dir+'etaInit.nc', 'eta'))
time_loop = slim3d.Loop(equations, time_loop = slim3d.Loop(equations,
time_step=5, time_step=5,
export_time_step=200, export_time_step=200,
...@@ -75,7 +74,6 @@ time_loop.export_uv(False) ...@@ -75,7 +74,6 @@ time_loop.export_uv(False)
time_loop.export_w() time_loop.export_w()
time_loop.export_uv2d(False) time_loop.export_uv2d(False)
#time_loop.loop()
time_loop.setup() time_loop.setup()
while time_loop.get_time() < time_loop.final_time: while time_loop.get_time() < time_loop.final_time:
time_loop.advance_one_time_step() time_loop.advance_one_time_step()
......
Supports Markdown
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