Commit fa2c146d authored by David Vincent's avatar David Vincent

The user is now able to modify the parameter of the sediment settling

and erosion parametrisation
parent 976d64a7
Pipeline #1911 passed with stage
in 39 minutes and 36 seconds
......@@ -981,7 +981,7 @@ class ShallowWaterTracer2d:
self._mass = dgpy.fullMatrixDouble(1, 1)
self._integrator_mass.compute(self._mass,"")
def set_sediment(self, initial_bottom_concentration, windU, windV):
def set_sediment(self, initial_bottom_concentration, windU, windV, factorW = 1e-2, wMax = 1e-4, u0 = 0.3, w0 = 10, a01 = 2.8e-2, a02 = 1.44e-1, a1 = 1e-6, n = 4., omega2 = 2.4538, eros0Fact = 0.245):
"""Define wind velocity for sediment transport.
keyword argument:
......@@ -992,6 +992,28 @@ class ShallowWaterTracer2d:
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].
* factorW
dependence factor between the settling viscosity and the suspended sediment concentration (default: 0.01 m^4/kg/s)
* wMax
maximum settling viscosity (default: 1e-4 m/s)
* u0
threshold velocity for erosion and settling (default: 0.3 m/s, see Lambrechts et al., 2010, eq 1 and 6)
* w0
wind speed threshold (default: 10 m/s)
* a01
empirical constant that depend only on the characteristics of the fine sediment on the seafloor (default: 2.8e-2).
It is used to compute A1 in eq 8 of Lambrechts et al., 2010: A1= a01 + (a02-a01) * (0.5 - atan(10*alpha)/(2*atan(10))) where alpha is the wind orientation factor
* a02
empirical constant that depend only on the characteristics of the fine sediment on the seafloor (default: 1.44e-1)
It is used to compute A1 in eq 8 of Lambrechts et al., 2010: A1= a01 + (a02-a01) * (0.5 - atan(10*alpha)/(2*atan(10))) where alpha is the wind orientation factor
* a1
empirical constant that depend only on the characteristics of the fine sediment on the seafloor (default: 1 e-6, A2 in eq 9 of Lambrechts et al., 2010)
* n
constant in erosion flux parametrisation (default: 4, see Lambrechts et al., 2010, eq 1)
* omega2
square of the wave frequency (default: 2.4538 Hz, see Lambrechts et al., 2010, eq 8-10)
* eros0Fact
erosion factor multiplying F (defualt: 0.245, see Lambrechts et al., 2010, eq 10)
"""
self._is_sediment = True
self._sediment_bottom = dgpy.dgDofContainer(self._domain._groups, 1)
......@@ -999,6 +1021,16 @@ class ShallowWaterTracer2d:
_windU=slim_private._load_function(windU, self._domain._groups)
_windV=slim_private._load_function(windV, self._domain._groups)
self._wind = slim_private._check_vector(_windU, _windV, None, self._domain)
self._factorW = factorW
self._wMax = wMax
self._u0 = u0
self._w0 = w0
self._a01 = a01
self._a02 = a02
self._a1 = a1
self._n = n
self._omega2 = omega2
self._eros0Fact = eros0Fact
def compute_sediment(self,dt):
"""Compute the sediment flux
......@@ -1011,7 +1043,7 @@ class ShallowWaterTracer2d:
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 = dgpy.sedimentBottomFlux2D(self._solution.getFunction(), self._hydro_sol, self._wind, self._domain._bath_function, self._sediment_bottom.getFunction(), self._domain._g, self._factorW, self._wMax, self._u0, self._w0, self._a01, self._a02, self._a1, self._n, self._omega2, self._eros0Fact)
_flux_dof.L2Projection(_flux)
_flux_dof.scale(dt)
self._solution.axpy(_flux_dof,1)
......
......@@ -1335,13 +1335,23 @@ void getHC::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *sol, 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, double factorW, double wMax, double u0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact) {
_sed = sed;
_sol = sol;
_wind = wind;
_bath = bath;
_sedGround = sedGround;
_gravity = gravity;
_factorW = factorW;
_wMax = wMax;
_u0 = u0;
_w0 = w0;
_a01 = a01;
_a02 = a02;
_a1 = a1;
_n = n;
_omega2 = omega2;
_eros0Fact = eros0Fact;
}
void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val) const {
......@@ -1351,37 +1361,27 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &sedGround = cm.get(*_sedGround);
const fullMatrix<double> &gravity = cm.get(*_gravity);
double factorW = 1e-2;
double wMax = 1e-4;
double u0 = 0.3;
double w0 = 10;
double a01 = 2.8e-2;
double a02 = 1.44e-1;
double a1 = 1e-6;
double n = 4.;
double omega2 = 2.4538;
double eros0Fact = 0.245;
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(sol(i,1), sol(i,2));
double wS = std::max(0., std::min(wMax, factorW * sed(i,0)));
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;
if (sed(i,0) < -1e-9) Msg::Warning("Wrong!! Negative sediment concentration");
double sed0 = std::max(0., sed(i,0));
double settling = (uBot < u0) ? - sed0 * wS * (1.- _pow2(uBot/u0)) : 0;
double settling = (uBot < _u0) ? - sed0 * wS * (1.- _pow2(uBot/_u0)) : 0;
if (sedGround(i,0) < 1e-6)
val(i,0) = settling;
else{
double erosion0 = 0;
if (uWind > 1e-3){
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;
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;
}
double erosion1 = (uBot > u0) ? a1*(pow(uBot/u0, n) - 1) : 0;
double erosion1 = (uBot > _u0) ? _a1*(pow(uBot/_u0, _n) - 1) : 0;
double erosionTotal = std::min(0.01*sedGround(i,0), erosion0+erosion1); // supposing dt3d < 100, there will be no more erosion than available
if (erosionTotal < 0)
Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
......
......@@ -388,9 +388,10 @@ class getHC: public functor {
class sedimentBottomFlux2D : public functor {
const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
double _factorW, _wMax, _u0, _w0, _a01, _a02, _a1, _n, _omega2, _eros0Fact;
public :
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
sedimentBottomFlux2D(const functor *sed, const functor *sol, 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, double factorW, double wMax, double u0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact);
};
......
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