Commit 0d2ab455 authored by David Vincent's avatar David Vincent
Browse files

adapt slim2d interface in order to solve the sediment transport

parent 28681ed8
Pipeline #1805 passed with stage
in 32 minutes and 43 seconds
......@@ -805,7 +805,7 @@ class ShallowWaterTracer2d:
boolean stating if you want to export the hydrodynamics every sub time step for further use of offline tracer (default: False)
* initial_time
initial time of the simulation (default: 0)
Format:
- float [in seconds]
......@@ -817,6 +817,8 @@ class ShallowWaterTracer2d:
- float [in seconds]
- date string in format "year-month-day hours:minutes:seconds" (e.g. "2000-06-20 22:42:17")
* is_sediment
use shallow water tracer law to predict the sediment transport (default: False)
"""
fmt = '%Y-%m-%d %H:%M:%S'
if initial_time is None:
......@@ -834,6 +836,9 @@ class ShallowWaterTracer2d:
else:
self._final_time = float(final_time)
self._domain = domain
self._is_sediment = False
self._wind = None
self._sediment_bottom = None
self._equation = dgpy.dgConservationLawShallowWaterTracer2d()
self._equation.setIsLinear(linear_equation)
self._equation.hydroWithLaxFriedrichs(True)
......@@ -978,6 +983,43 @@ class ShallowWaterTracer2d:
self._mass = dgpy.fullMatrixDouble(1, 1)
self._integrator_mass.compute(self._mass,"")
def set_sediment(self, initial_bottom_concentration, windU, windV):
"""Define wind velocity for sediment transport.
keyword argument:
*initial_bottom_concentration
netcdf or .msh file containing the initial sediment concentration at the bottom of the domain for the whole domain [in c/m^2].
* windU
netcdf or .msh file containing the surface wind velocity along the x-axis in the local basis [in m*s^-1].
* windV
netcdf or .msh file containing the surface wind velocity along the y-axis in the local basis [in m*s^-1].
"""
self._sediment_bottom = dgpy.dgDofContainer(self._domain._groups, 1)
slim_private._load(self._sediment_bottom, initial_bottom_concentration)
self._wind = slim_private._check_vector(windU, windV, None, self._domain)
def compute_sediment(self,dt):
"""Compute the sediment flux
keyword arguments:
* dt
time step
"""
if self._wind is None or self._sediment_bottom is None:
dgpy.Msg.Fatal("You have to define the surface wind velocity by using the function set_wind of class ShallowWaterTracer2d if you want to predict sediment transport")
_flux_dof = dgpy.dgDofContainer(self._domain._groups, 1)
_flux = dgpy.sedimentBottomFlux2D(self._solution.getFunction(), self._hydro_sol, self._wind, self._domain._bath_function, self._sediment_bottom.getFunction(), self._domain._g)
_flux_dof.L2Projection(_flux)
_flux_dof.scale(dt)
self._solution.axpy(_flux_dof,1)
_HC = dgpy.dgDofContainer(self._domain._groups, 1)
_H = slim_private._sum(self._solution.getFunction(), self._domain._bath_function).functor
_H_dof= dgpy.dgDofContainer(self._domain._groups, 1)
_H_dof.interpolate(_H)
_HC.multiply(_flux_dof,_H_dof)
self._sediment_bottom.axpy(_HC,1)
class Loop:
"""Temporal solver"""
......@@ -1106,12 +1148,14 @@ class Loop:
equation._sys_tracer2 = slim_private._new_default_linear_system(equation._domain, equation._equation)
equation._temporal_solver = dgpy.dgMDIRK(equation._equation, equation._sys_tracer1.dof_manager, equation._sys_tracer2.dof_manager, 2)
equation._temporal_solver.getNewton().setVerb(10)
equation._temporal_solver.getNewton().setAtol(1.e-8)
equation._temporal_solver.getNewton().setAtol(1e-8)
equation._temporal_solver.getNewton().setRtol(1e-8)
else:
dgpy.Msg.Fatal("Temporal scheme '"+equation._temporal_scheme+"' not implemented!")
self._number_of_tracer += 1
self._export_dofs.append([dgpy.dgIdxExporter(equation._solution, self._output_directory +'/'+ equation._name), equation._initial_time, equation._final_time])
if equation._is_sediment:
self._export_dofs.append([dgpy.dgIdxExporter(equation._sediment_bottom, self._output_directory +'/'+ equation._name + "_sediment_bottom"), equation._initial_time, equation._final_time])
if equation._initial_time < self._min_initial_time or self._number_of_equations == 0 :
self._min_initial_time = equation._initial_time
self._time = self._min_initial_time
......@@ -1261,6 +1305,9 @@ class Loop:
self._equations[0]._temporal_solver.iterate(self._equations[0]._solution, self._dt, self._time)
else:
for i in self._equations:
if isinstance(i,ShallowWaterTracer2d):
if i._is_sediment:
i.compute_sediment(self._dt)
if self._time > i._initial_time-self._eps and self._time < i._final_time+self._eps:
i._temporal_solver.initialize(i._solution, self._dt, self._time)
for j in range(self._equations[0]._temporal_solver.nstep()):
......
......@@ -24,7 +24,20 @@ class _velocity_extractor :
val[:,0] = sol[:,1]
val[:,1] = sol[:,2]
val[:,2] = 0
return val
return val
class _sum:
def __init__(self,sol,bath) :
self._sol = sol
self._bath = bath
self.functor = dgpy.functorNumpy(self.call)
def call(self, cm) :
sol = cm.get(self._sol)
bath = cm.get(self._bath)
nPt = sol.shape
val = np.zeros([nPt[0], 1])
val[:,0] = sol[:,0] + bath[:]
return val
class _get_x:
def __init__(self,sol) :
......
......@@ -102,7 +102,7 @@ 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
error = np.sqrt(integral(0,0))/1e6/1e6
print("L2 error on tracer: %e\n" % error)
if (error > 1e-12):
dgpy.Msg.Error("Error is too high");
......
......@@ -1335,9 +1335,9 @@ void getHC::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *uvBot, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity) {
sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity) {
_sed = sed;
_uvBot = uvBot;
_sol = sol;
_wind = wind;
_bath = bath;
_sedGround = sedGround;
......@@ -1346,7 +1346,7 @@ sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *uv
void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &sed = cm.get(*_sed);
const fullMatrix<double> &uvBot = cm.get(*_uvBot);
const fullMatrix<double> &sol = cm.get(*_sol);
const fullMatrix<double> &wind = cm.get(*_wind);
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &sedGround = cm.get(*_sedGround);
......@@ -1364,7 +1364,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
val.resize(cm.nEvaluationPoint(),1,true);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
if (sed(i,0) > 3) Msg::Warning("vertical Velocity for sediment parametrization not valid anymore if sediment concentration greater than 3kg m^-3");
double uBot = hypot(uvBot(i,0), uvBot(i,1));
double uBot = hypot(sol(i,1), sol(i,2));
double wS = std::max(0., std::min(wMax, factorW * sed(i,0)));
double uWind = hypot(wind(i,0), wind(i,1));
wS = (uWind < 10) ? wS : 0;
......@@ -1376,7 +1376,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
else{
double erosion0 = 0;
if (uWind > 1e-3){
double alpha = (wind(i,0)*uvBot(i,0) + wind(i,1)*uvBot(i,1))/(uWind*uBot);
double alpha = (wind(i,0)*sol(i,1) + wind(i,1)*sol(i,2))/(uWind*uBot);
double a0 = a01 + (a02 - a01) * (0.5 - atan(10*alpha)/(2*atan(10)));
double F = eros0Fact*exp(-0.905 * omega2 * bath(i,0) / gravity(i,0) - 0.0207);
erosion0 = a0 * pow(uWind/w0, 3)*pow(uBot/u0, n)*F;
......
......@@ -387,10 +387,10 @@ class getHC: public functor {
class sedimentBottomFlux2D : public functor {
const functor *_sed, *_uvBot, *_wind, *_bath, *_sedGround, *_gravity;
const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
public :
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
sedimentBottomFlux2D(const functor *sed, const functor *uvBot, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity);
sedimentBottomFlux2D(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity);
};
......
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