Commit 5fafcc74 authored by David Vincent's avatar David Vincent
Browse files

fix problem in sedimentBottomFlux2D:

* the settling term has also to be divided by H (= bathy+eta)
* the erosion term has to be divided by H instead of bath

use interpolate instead of L2Projection in compute_sediment function of
Tracer class
parent fa2c146d
Pipeline #1930 passed with stage
in 37 minutes and 21 seconds
......@@ -1044,7 +1044,7 @@ class ShallowWaterTracer2d:
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, 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.interpolate(_flux)
_flux_dof.scale(dt)
self._solution.axpy(_flux_dof,1)
_HC = dgpy.dgDofContainer(self._domain._groups, 1)
......
......@@ -1372,7 +1372,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
double sed0 = std::max(0., sed(i,0));
double settling = (uBot < _u0) ? - sed0 * wS * (1.- _pow2(uBot/_u0)) : 0;
if (sedGround(i,0) < 1e-6)
val(i,0) = settling;
val(i,0) = settling/(bath(i,0)+sol(i,0));
else{
double erosion0 = 0;
if (uWind > 1e-3){
......@@ -1385,7 +1385,7 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
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");
val(i,0) = (settling + erosionTotal)/bath(i,0);
val(i,0) = (settling + erosionTotal)/(bath(i,0)+sol(i,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