Commit 20fed63b authored by lambrechts's avatar lambrechts
Browse files

dg : shallow water 2d : add a lax-friedrichs option

git-svn-id: https://geuz.org/svn/gmsh/trunk@20021 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 22a28ae4
......@@ -439,26 +439,41 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double h = bath(i,0);
uR *= rev;
if (_linear) {
//linear equations
if (_laxFriedrichs) {
if (_linear) {
Msg::Fatal("no need for lax friedrichs in linear equation");
}
double etaR = solR(i,0);
double etaL = solL(i,0);
double sq_g_h = sqrt(_g/h);
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;
}else{
double HR = solR(i,0) + h, HL = solL(i,0) + h;
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;
Fun = -_g*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
const double etaM = (etaR + etaL) /2;
const double unM = (uR + uL) /2;
const double utM = (vR + vL) /2;
const double c = sqrt(_g * (h + etaM) + fabs(unM));
Fun = J * (- _g * etaM - unM * unM) + c * (uR - uL) / 2;
Feta = J * (- (h + etaM) * unM) + c * (etaR - etaL) / 2;
Fut = J * (-unM * utM) + c * (vR - vL) / 2;
}
else {
if (_linear) {
//linear equations
double etaR = solR(i,0);
double etaL = solL(i,0);
double sq_g_h = sqrt(_g/h);
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;
}else{
double HR = solR(i,0) + h, HL = solL(i,0) + h;
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;
Fun = -_g*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
}
}
val(i,0) = Feta;
......
......@@ -31,6 +31,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
bool _useMovingBathWD;
double _alphaMovingBathWD;
double _upwindFactorRiemann;
bool _laxFriedrichs;
double _R;
double _g;
......@@ -57,6 +58,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
inline void setIsLinear(bool linear) {_linear = linear; _constantJac = linear;}
virtual void setImexMode(imexMode mode);
void setLinearFilterMode(bool set);
inline void setLaxFriedrichs(bool flag) { _laxFriedrichs = flag;}
virtual bool isLinear() const {return (_linear || _linearIMEX);}
virtual bool isConstantJac() const {return (_constantJac || _constantJacIMEX);}
/**if this flag is true, a spherical version of the equations are solved*/
......
......@@ -105,7 +105,7 @@ const dgDofContainer &dgNewton::solve (const dgDofContainer *U0, double alpha, c
double newtonResidual = 0;
newtonResidual = _dof.getNormInfRHS();
if (newtonResidual != newtonResidual) {
if (_verb > 0)
if (_verb > 0 && Msg::GetCommRank() == 0)
Msg::Info("Newton failed");
_converged = 0;
if(_returnZeroIfNotConverged)
......@@ -119,7 +119,7 @@ const dgDofContainer &dgNewton::solve (const dgDofContainer *U0, double alpha, c
conv = log(1. / newtonResidual) / log(10);
else if (ii > 0. && newtonResidual != 0.0 )
conv = log(firstNewtonResidual / newtonResidual) / log(10);
if (_verb > 1)
if (_verb > 1 && Msg::GetCommRank() == 0)
Msg::Info("Newton residual (%i): %g (%.2f)", ii, newtonResidual, conv);
//todo : this check ii>0 are for boundary condition but we have to find another solution
if (ii > 0 && ((firstNewtonResidual != 0.0 && newtonResidual / firstNewtonResidual < _rtol) || (newtonResidual < _atol))) {
......@@ -156,7 +156,7 @@ const dgDofContainer &dgNewton::solve (const dgDofContainer *U0, double alpha, c
}
}
if (!hasConverged && _verb>0) {
if (!hasConverged && _verb>0 && Msg::GetCommRank() == 0) {
Msg::Info("Newton not converged after %i iters !!",_maxit);
_converged = 0;
if(_returnZeroIfNotConverged)
......
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