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

The user can now set different threshold velocities (v0) for settling and

erosion (u0)
parent ab08d470
Pipeline #1974 failed with stage
in 41 minutes and 20 seconds
...@@ -981,7 +981,7 @@ class ShallowWaterTracer2d: ...@@ -981,7 +981,7 @@ class ShallowWaterTracer2d:
self._mass = dgpy.fullMatrixDouble(1, 1) self._mass = dgpy.fullMatrixDouble(1, 1)
self._integrator_mass.compute(self._mass,"") self._integrator_mass.compute(self._mass,"")
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): def set_sediment(self, initial_bottom_concentration, windU, windV, factorW = 1e-2, wMax = 1e-4, u0 = 0.3, , v0 = 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. """Define wind velocity for sediment transport.
keyword argument: keyword argument:
...@@ -997,7 +997,9 @@ class ShallowWaterTracer2d: ...@@ -997,7 +997,9 @@ class ShallowWaterTracer2d:
* wMax * wMax
maximum settling viscosity (default: 1e-4 m/s) maximum settling viscosity (default: 1e-4 m/s)
* u0 * u0
threshold velocity for erosion and settling (default: 0.3 m/s, see Lambrechts et al., 2010, eq 1 and 6) threshold velocity for erosion (default: 0.3 m/s, see Lambrechts et al., 2010, eq 1 and 6)
* v0
threshold velocity for settling (default: 0.3 m/s, see Lambrechts et al., 2010, eq 1 and 6)
* w0 * w0
wind speed threshold (default: 10 m/s) wind speed threshold (default: 10 m/s)
* a01 * a01
...@@ -1024,6 +1026,7 @@ class ShallowWaterTracer2d: ...@@ -1024,6 +1026,7 @@ class ShallowWaterTracer2d:
self._factorW = factorW self._factorW = factorW
self._wMax = wMax self._wMax = wMax
self._u0 = u0 self._u0 = u0
self._v0 = v0
self._w0 = w0 self._w0 = w0
self._a01 = a01 self._a01 = a01
self._a02 = a02 self._a02 = a02
...@@ -1043,7 +1046,7 @@ class ShallowWaterTracer2d: ...@@ -1043,7 +1046,7 @@ class ShallowWaterTracer2d:
if self._wind is None or self._sediment_bottom is None: 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") 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_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, self._factorW, self._wMax, self._u0, self._w0, self._a01, self._a02, self._a1, self._n, self._omega2, self._eros0Fact) _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._v0, self._w0, self._a01, self._a02, self._a1, self._n, self._omega2, self._eros0Fact)
_flux_dof.interpolate(_flux) _flux_dof.interpolate(_flux)
_flux_dof.scale(dt) _flux_dof.scale(dt)
self._solution.axpy(_flux_dof,1) self._solution.axpy(_flux_dof,1)
......
...@@ -1335,7 +1335,7 @@ void getHC::operator()(functorCache &cm, fullMatrix<double> &val) const { ...@@ -1335,7 +1335,7 @@ 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, double factorW, double wMax, double u0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact) { 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 v0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact) {
_sed = sed; _sed = sed;
_sol = sol; _sol = sol;
_wind = wind; _wind = wind;
...@@ -1345,6 +1345,7 @@ sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *so ...@@ -1345,6 +1345,7 @@ sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *so
_factorW = factorW; _factorW = factorW;
_wMax = wMax; _wMax = wMax;
_u0 = u0; _u0 = u0;
_v0 = v0;
_w0 = w0; _w0 = w0;
_a01 = a01; _a01 = a01;
_a02 = a02; _a02 = a02;
...@@ -1370,7 +1371,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val) ...@@ -1370,7 +1371,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
wS = (uWind < 10) ? wS : 0; wS = (uWind < 10) ? wS : 0;
if (sed(i,0) < -1e-9) Msg::Warning("Wrong!! Negative sediment concentration"); if (sed(i,0) < -1e-9) Msg::Warning("Wrong!! Negative sediment concentration");
double sed0 = std::max(0., sed(i,0)); double sed0 = std::max(0., sed(i,0));
double settling = (uBot < _u0) ? - sed0 * wS * (1.- _pow2(uBot/_u0)) : 0; double settling = (uBot < _v0) ? - sed0 * wS * (1.- _pow2(uBot/_v0)) : 0;
if (sedGround(i,0) < 1e-6) if (sedGround(i,0) < 1e-6)
val(i,0) = settling/(bath(i,0)+sol(i,0)); val(i,0) = settling/(bath(i,0)+sol(i,0));
else{ else{
......
...@@ -388,10 +388,10 @@ class getHC: public functor { ...@@ -388,10 +388,10 @@ class getHC: public functor {
class sedimentBottomFlux2D : public functor { class sedimentBottomFlux2D : public functor {
const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity; const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
double _factorW, _wMax, _u0, _w0, _a01, _a02, _a1, _n, _omega2, _eros0Fact; double _factorW, _wMax, _u0, _w0, _a01, _a02, _a1, _n, _omega2, _eros0Fact, _v0;
public : public :
void operator()(functorCache &cm, fullMatrix<double> &val) const ; 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, double factorW, double wMax, double u0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact); 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 v0, 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