Commit e3ee8d75 authored by David Vincent's avatar David Vincent
Browse files

create new testcase: check if tracer offline is ok (solve hydro and tracer,...

create new testcase: check if tracer offline is ok (solve hydro and tracer, then tracer offline and compare tracer and tracer offline)
It also uses importIdx (on a tracer)
modification in slim.py:
	* hydro_sol in tracer has the same initial conition as the hydro (it was 0)
	* the file for offline are imported by means of importMsh instead of importIdx
parent fef24e83
Pipeline #1194 passed with stage
in 28 minutes
......@@ -617,9 +617,9 @@ class ShallowWaterTracer2d:
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!")
self._datafile = hydro_sol+'/offline_sw2d/'
self._hydro_sol_dof = dgpy.dgDofContainer(self._domain._groups, 3)
self._hydro_sol_dof.setAll(0)
self._datafile = hydro_sol+'/'+slim_private.path.basename(hydro_sol)
self._hydro_sol_dof.importMsh(self._datafile+"offline_sw2d-000000_000000/offline_sw2d-000000_000000")
self._hydro_sol = self._hydro_sol_dof.getFunction()
else:
if not isinstance(hydro_sol, ShallowWater2d):
......@@ -887,6 +887,7 @@ class Loop:
self._periodic = periodic
self._index_start = index_start
self._n_index_per_period = n_index_per_period
self._offline = True
if periodic and n_index_per_period < 0:
dgpy.Msg.Fatal("in set_time_step_offline: number of file per period not defined")
......@@ -1016,7 +1017,7 @@ class Loop:
index = self._index_start+self._n_iter_offline%self._n_index_per_period
else:
index = self._n_iter_offline
i._hydro_sol_dof.importIdx(i._datafile+"-%06d_%06d.idx" %(index,j))
i._hydro_sol_dof.importMsh(i._datafile+"offline_sw2d-%06d_%06d/offline_sw2d-%06d_%06d" %(index,j,index,j))
i._hydro_sol_dof.scatter()
newton_converged = i._temporal_solver.subiterate(i._solution, self._dt, self._time)
if not newton_converged:
......@@ -1026,7 +1027,7 @@ class Loop:
index = self._index_start+self._n_iter_offline%self._n_index_per_period
else:
index = self._n_iter_offline
i._hydro_sol_dof.importIdx(i._datafile+"-%06d_%06d.idx" %(index,j))
i._hydro_sol_dof.importMsh(i._datafile+"offline_sw2d-%06d_%06d/offline_sw2d-%06d_%06d" %(index,j,index,j))
i._hydro_sol_dof.scatter()
for e in self._export_dofs_full:
e[0].exportIdxSubTimeStep(self._index_exporter_full, j, self._time - e[1])
......
d = 5e5;
Point(1) = {-d, -d, 0};
Point(2) = {d, -d, 0};
Point(3) = {d, d, 0};
Point(4) = {-d, d, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "8e4";
//Field[1].F = "1e5";
Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.LcIntegrationPrecision = 1e-3;
Mesh.CharacteristicLengthExtendFromBoundary=0;
Mesh.CharacteristicLengthFromPoints=0;
#################################################
########## PREPROCESSING ##########
#################################################
import slimPre
import numpy as np
from dgpy.scripts.termcolor import colored
from dgpy.scripts import Common
import dgpy
#mesh
Common.genMesh("square.geo", 2)
mesh = slimPre.Mesh("square.msh")
#bathymetry
slimPre.write_file('bathymetry.nc', region=None, time=None, data=[('bathymetry', 10)])
#gravity
slimPre.write_file('gravity.nc', region=None, time=None, data=[('gravity', 9.80616)])
#manning
slimPre.write_file('manning.nc', region=None, time=None, data=[('manning', 0.)])
#coriolis
slimPre.write_file('coriolis.nc', region=None, time=None, data=[('coriolis', 0.)])
#initial condition
region_global = slimPre.Region(mesh, '')
xyz = region_global.coordinates
nPt = len(xyz)
eta = np.empty(nPt)
u = np.empty(nPt)
v = np.empty(nPt)
tracer = np.empty(nPt)
for i in range(nPt):
eta[i]=(xyz[i,0]/500000.)
u[i]=0
v[i]=0
tracer[i]=(1.+(xyz[i,0]/500000.))/2.
slimPre.write_file('sse.nc', region=region_global, time=None, data=[('sse', eta)])
slimPre.write_file('u.nc', region=region_global, time=None, data=[('u', u)])
slimPre.write_file('v.nc', region=region_global, time=None, data=[('v', v)])
slimPre.write_file('tracer_CI.nc', region=region_global, time=None, data=[('CI', tracer)])
#parameter
dt = 15*60*10
nbSteps = 50
n_print = 10
#################################################
########## SIMULATION ##########
#################################################
import slim
domain = slim.Domain("square.msh",("bathymetry.nc","bathymetry"), g= ("gravity.nc","gravity"))
eq = slim.ShallowWater2d(domain, "implicit", initial_time=0, final_time=nbSteps*dt, export_every_sub_time_step = True)
eq.set_initial_condition(sse=("sse.nc","sse"), ux=("u.nc","u"), uy=("v.nc","v"))
eq.set_viscosity("constant", constant_viscosity=0.00)
eq.set_boundary_coast("Wall", slip=True)
eq.compute_mass("mass.dat")
eq2 = slim.ShallowWaterTracer2d(domain, "implicit",eq, name="tracer", offline = False, initial_time=0, final_time=nbSteps*dt)
eq2.set_initial_condition(("tracer_CI.nc","CI"))
eq2.set_diffusivity("okubo",okubo_coefficient=0.0)
eq2.set_boundary_coast("Wall")
eq2.compute_mass("tracer.dat")
loop=slim.Loop(maximum_time_step = dt, export_time = dt,path="./output")
loop.add_equation(eq)
loop.add_equation(eq2)
loop.run()
eq3 = slim.ShallowWaterTracer2d(domain, "implicit",'./output', name="tracer_offline", offline = True, initial_time=0, final_time=nbSteps*dt)
eq3.set_initial_condition(("tracer_CI.nc","CI"))
eq3.set_diffusivity("okubo",okubo_coefficient=0.0)
eq3.set_boundary_coast("Wall")
eq3.compute_mass("tracer.dat")
loop2=slim.Loop(maximum_time_step = dt, export_time = dt,path="./output")
loop2.set_time_step_offline(dt, periodic = False, index_start=1, n_index_per_period=-1)
loop2.add_equation(eq3)
loop2.run()
tracer = dgpy.dgDofContainer(mesh._groups, 1)
tracer.importIdx("./output/tracer/tracer-%06d.idx"%nbSteps)
tracer_offline = dgpy.dgDofContainer(mesh._groups, 1)
tracer_offline.importIdx("./output/tracer_offline/tracer_offline-%06d.idx"%nbSteps)
def check_error(cm):
t1 = cm.get(tracer.getFunction())
t2 = cm.get(tracer_offline.getFunction())
return (t1[:] - t2[:])**2
errorF = dgpy.functorNumpy(check_error)
integrator = dgpy.dgFunctionIntegrator(mesh._groups, errorF)
integral = dgpy.fullMatrixDouble(3,1)
integrator.compute(integral)
error = np.sqrt(integral(0,0))/5e5/5e5
print("L2 error on tracer: %e\n" % error)
if (error > 1e-12):
dgpy.Msg.Error("Error is too high");
dgpy.Msg.Exit(1)
dgpy.Msg.Exit(0)
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