Commit c9ea0532 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

Wetting&Drying and dg bath

parent 7484d26c
......@@ -393,10 +393,11 @@ class dgConservationLawShallowWater2d::diffusivity : public function {
void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<double> &val) const {
const size_t nQP = cache.nEvaluationPoint();
val.resize(nQP,3);
val.resize(nQP,6);
const fullMatrix<double> &solL = cache.get(cache.solution(), 0);
const fullMatrix<double> &solR = cache.get(cache.solution(), 1);
const fullMatrix<double> &bath = cache.get(*_bathymetry,0);
const fullMatrix<double> &bathL = cache.get(*_bathymetry,0);
const fullMatrix<double> &bathR = cache.get(*_bathymetry,1);
const fullMatrix<double> &normals = cache.get(*function::getNormals(),0);
const fullMatrix<double> &movingBathFactorL = cache.get(*_movingBathFactor,0);
const fullMatrix<double> &movingBathFactorR = cache.get(*_movingBathFactor,1);
......@@ -435,7 +436,8 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double uL = nx*solL(i,1) + ny*solL(i,2), uR = nx*solR(i,1) + ny*solR(i,2);
double vL = -ny*solL(i,1) + nx*solL(i,2), vR = - ny*solR(i,1) + nx*solR(i,2);
double h = bath(i,0);
double hL = bathL(i,0);
double hR = bathR(i,0);
uR *= rev;
if (_laxFriedrichs) {
......@@ -443,21 +445,29 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double etaR = solR(i,0);
double etaL = solL(i,0);
const double etaM = (etaR + etaL) /2;
const double hMin = std::(hL, hR);
uR = uR * hR/hMin;
uL = uL * hL/hMin;
const double unM = (uR + uL) /2;
const double c = sqrt(_g * (h + etaM) + fabs(unM));
Fun = J * (- _g * etaM ) + c * (uR - uL) / 2;
Feta = J * (- h * unM) + c * (etaR - etaL) / 2;
Fut = c * (vR - vL) / 2;
const double c = sqrt(_g * hMin);
Fun = J * (- _g * etaM ) + c * (uR - uL) / 2;
Feta = J * (- hMin * unM) + c * (etaR - etaL) / 2;
Fut = 0;
}
else {
double etaR = solR(i,0);
double etaL = solL(i,0);
const double etaM = (etaR + etaL) /2;
const double hMin = std::(hL, hR);
uR = uR * hR/hMin;
uL = uL * hL/hMin;
vR = vR * hR/hMin;
vL = vL * hL/hMin;
const double unM = (uR + uL) /2;
const double utM = (vR + vL) /2;
const double c = sqrt(_g * (h + etaM) + fabs(unM));
const double c = sqrt(_g * (hMin + etaM) + fabs(unM));
Fun = J * (- _g * etaM - unM * unM) + c * (uR - uL) / 2;
Feta = J * (- (h + etaM) * unM) + c * (etaR - etaL) / 2;
Feta = J * (- (hMin + etaM) * unM) + c * (etaR - etaL) / 2;
Fut = J * (-unM * utM) + c * (vR - vL) / 2;
}
}
......@@ -466,18 +476,22 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
//linear equations
double etaR = solR(i,0);
double etaL = solL(i,0);
double sq_g_h = sqrt(_g/h);
const double hMin = std::(hL, hR);
double sq_g_h = sqrt(_g/hMin);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -_g*etaRiemann*J;
Fut = 0;
Feta = -(h)*unRiemann*J;
Feta = -(hMin)*unRiemann*J;
}else{
double HR = solR(i,0) + h, HL = solL(i,0) + h;
const double hMin = std::(hL, hR);
double HR = solR(i,0) + hMin, HL = solL(i,0) + hMin;
if ( fabs(hL-hR) > 1e-3 )
Msg::Error("Discontinuous Bathymetry is not implemented for non-linear equation with Roe Solver. Please use the linear equations and/or the LaxFriedrichs Flux.");
if(HR<0 || HL<0) printf(" HR = %e HL =%e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g, movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
double eta = H-h;
double eta = H-hMin;
Fun = -_g*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
......@@ -487,11 +501,10 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
val(i,0) = Feta;
val(i,1) = Fun * nx - Fut * ny;
val(i,2) = Fun * ny + Fut * nx;
/*val(i,3) = -val(i,0);
val(i,4) = - Fun * nx * rev + Fut * ny;
val(i,5) = - Fun * ny * rev - Fut * nx;*/
val(i,3) = -val(i,0);
val(i,4) = -val(i,1);
val(i,5) = -val(i,2);
//printf("%f\n", val(i,0));
if (_isDiffusive) {
val(i,1)+=(*ip)(i,1);
......@@ -595,14 +608,7 @@ void dgConservationLawShallowWater2d::setLinearFilterMode(bool set) {
}
void dgConservationLawShallowWater2d::setup() {
if (_useMovingBathWD) {
_bathymetry = new movingBath(_originalBathymetry, _alphaMovingBathWD);
_bathymetryGradient = new movingBathGradient(_originalBathymetry, _originalBathymetryGradient, _alphaMovingBathWD);
_movingBathFactor = new movingBathFactor(_originalBathymetry, _alphaMovingBathWD);
_movingBathFactorGradient = new movingBathFactorGradient(_originalBathymetry,_originalBathymetryGradient, _alphaMovingBathWD);
}
//_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear);
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz, _g);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
......@@ -626,6 +632,17 @@ void dgConservationLawShallowWater2d::setup() {
_clipToPhysics[""] = new clipToPhysics(_bathymetry, 0.01);
}
void dgConservationLawShallowWater2d::setMovingBathWettingDrying(double alphaMovingBathWD){
_useMovingBathWD = true;
_alphaMovingBathWD = alphaMovingBathWD;
_originalBathymetry = _bathymetry;
_originalBathymetryGradient= _bathymetryGradient;
_bathymetry = new movingBath(_originalBathymetry, _alphaMovingBathWD);
_bathymetryGradient = new movingBathGradient(_originalBathymetry, _originalBathymetryGradient, _alphaMovingBathWD);
_movingBathFactor = new movingBathFactor(_originalBathymetry, _alphaMovingBathWD);
_movingBathFactorGradient = new movingBathFactorGradient(_originalBathymetry,_originalBathymetryGradient, _alphaMovingBathWD);
}
dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
delete _fzero;
delete _fzerov;
......
......@@ -84,12 +84,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
/**set an absorbing layer. Param: TAG, AbsCoef, ExtFields */
void setAbsLayer(const std::string, const functor *, const functor *);
/**use the moving bath WD (Karna et al, 2011) and define the bathymetry and the moving bathymetry */
inline void setMovingBathWettingDrying(double alphaMovingBathWD) {
_useMovingBathWD = true;
_alphaMovingBathWD = alphaMovingBathWD;
_originalBathymetry = _bathymetry;
_originalBathymetryGradient= _bathymetryGradient;
}
void setMovingBathWettingDrying(double alphaMovingBathWD);
/** a factor to smooth the if in Hv term in the roe riemann solver */
inline void setUpwindFactorRiemann(double upwindFactorRiemann) {
_upwindFactorRiemann=upwindFactorRiemann;
......@@ -99,10 +94,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
return functionSumNew(_bathymetry,functionExtractCompNew(function::getSolution(),0));
}
inline const functor *getBathymetry() {
if(_useMovingBathWD)
return new movingBath(_originalBathymetry,_alphaMovingBathWD);
else
return _bathymetry;
return _bathymetry;
}
inline const functor *getElevation() {
return functionExtractCompNew(function::getSolution(),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