Commit 1c295516 authored by vallaeys's avatar vallaeys
Browse files

dg : shallowWater : add a possibility to change the gravitational acceleration for E.T. simulations

git-svn-id: https://geuz.org/svn/gmsh/trunk@19623 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent fa1ecec2
......@@ -3,7 +3,6 @@
#include "dgConservationLawShallowWater2d.h"
#include "MElement.h"
static double g = 9.81;
dgConservationLawShallowWater1d::dgConservationLawShallowWater1d() : dgConservationLawFunction(2)
{
_hydroSol=NULL;
......@@ -19,6 +18,7 @@ dgConservationLawShallowWater1d::dgConservationLawShallowWater1d() : dgConservat
_diffusivity = new functionConstant(0.);
_clipSectionValue = 0;
_useClipToPhysics = false;
_g = 9.80616;
}
//***************************************
// Depending field (elevation and width)
......@@ -102,12 +102,14 @@ public :
//*************************************
class dgConservationLawShallowWater1d::advection: public function {
fullMatrix<double> sol, solGradient, elevation, diffusivity;
double _g;
public:
advection(const function *hydroSolF, const function *elevationF, const function *diffusivityF) : function(6) {
advection(const function *hydroSolF, const function *elevationF, const function *diffusivityF, double g) : function(6) {
setArgument(sol, hydroSolF);
setArgument(solGradient, function::getSolutionGradient());
setArgument(elevation, elevationF);
setArgument(diffusivity, diffusivityF);
_g = g;
};
void call (dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = val.size1();
......@@ -119,7 +121,7 @@ public:
double dudx = solGradient(i,3);
// flux_x
val(i,0) = section*u;
val(i,1) = .5 * u*u + g*eta - nu*dudx ;
val(i,1) = .5 * u*u + _g*eta - nu*dudx ;
// flux_y
val(i,2) = 0;
val(i,3) = 0;
......@@ -177,8 +179,9 @@ class dgConservationLawShallowWater1d::riemannLAX:public function {
fullMatrix<double> solL_eta, solR_eta, etaL_w, etaR_w;
fullMatrix<double> etaL, etaR, wL, wR;
functionReplace replaceElevationL, replaceElevationR, replaceWidthL, replaceWidthR;
double _g;
public:
riemannLAX(const function *hydroSol, const function *width, const function *elevation, const function *widthData, const function *diffusivity): function(4) {
riemannLAX(const function *hydroSol, const function *width, const function *elevation, const function *widthData, const function *diffusivity, double g): function(4) {
setArgument(normalsL,function::getNormals(), 0);
setArgument(normalsR,function::getNormals(), 1);
setArgument(solL, hydroSol, 0);
......@@ -191,6 +194,7 @@ public:
setArgument(widthR, width, 1);
setArgument(diffusivityL, diffusivity, 0);
setArgument(diffusivityR, diffusivity, 1);
_g = g;
addFunctionReplace(replaceElevationL);//The elevation depends on the solution
replaceElevationL.replace(solL_eta, hydroSol, 0);
......@@ -210,7 +214,7 @@ public:
for (int iPt = 0; iPt < val.size1(); ++iPt) {
double sL = solL(iPt,0), sR = solR(iPt,0);
double uL = solL(iPt,1), uR = solR(iPt,1) ;
double pL = g*elevationL(iPt,0), pR = g*elevationR(iPt,0);
double pL = _g*elevationL(iPt,0), pR = _g*elevationR(iPt,0);
double nL = normalsL(iPt,0);
const dgGroupOfFaces &iGroup = *m->interfaceGroup();
int connectionR = iGroup.nConnection() == 2 ? 1 : 0;
......@@ -218,7 +222,7 @@ public:
const dgGroupOfElements *groupL = &iGroup.elementGroup(0);
if(groupR->getPhysicalTag() == groupL->getPhysicalTag()){
double fs, fu;
dgConservationLawShallowWater1dLaxSolver(sL, sR, uL, uR, widthL(iPt,0), widthR(iPt,0), pL, pR, nL, fs, fu);
dgConservationLawShallowWater1dLaxSolver(sL, sR, uL, uR, widthL(iPt,0), widthR(iPt,0), pL, pR, nL, fs, fu, _g);
val(iPt,0) = -fs;
val(iPt,1) = -fu;
val(iPt,2) = -val(iPt,0);
......@@ -251,8 +255,8 @@ public:
solR_eta(iPt,0)=u(1); replaceElevationR.compute(); eta(1)=etaR(iPt,0);
etaL_w(iPt,0)=eta(0); replaceWidthL.compute(); w(0)=wL(iPt,0);
etaR_w(iPt,0)=eta(1); replaceWidthR.compute(); w(1)=wR(iPt,0);
caracs(0) = 0.5*(u(0) + n(0)*u(2)*sqrt(u(0)*w(0)/g));
caracs(1) = 0.5*(u(1) + n(1)*u(3)*sqrt(u(1)*w(1)/g));
caracs(0) = 0.5*(u(0) + n(0)*u(2)*sqrt(u(0)*w(0)/_g));
caracs(1) = 0.5*(u(1) + n(1)*u(3)*sqrt(u(1)*w(1)/_g));
double maxNorm=1e-8;
double eps=1e-3;
double delta[4]={u(0)*eps,u(1)*eps,eps,eps};
......@@ -281,31 +285,33 @@ public:
}
if(nr==10) Msg::Fatal("Newton for discontinuity solver diverges | norm:%e u=%.3f %.3f eta=%.3f %.3f %s-%s\n",du.norm(),u(2),u(3),eta(0),eta(1),groupL->getPhysicalTag().c_str(),groupR->getPhysicalTag().c_str());
val(iPt,0)=-u(2)*u(0)*n(0);
val(iPt,1)=-(.5*u(2)*u(2) + g*eta(0))*n(0);
val(iPt,1)=-(.5*u(2)*u(2) + _g*eta(0))*n(0);
val(iPt,2)=-u(3)*u(1)*n(1);
val(iPt,3)=-(.5*u(3)*u(3) + g*eta(1))*n(1);
val(iPt,3)=-(.5*u(3)*u(3) + _g*eta(1))*n(1);
}
}
}
void systemRiemannDiscontinuity(fullVector<double> &u, fullVector<double> &eta, fullVector<double> &w, fullVector<double> &carac, fullVector<double> &n, fullVector<double> &f){
f(0) = .5*(u(0) + n(0)*u(2)*sqrt(u(0)*w(0)/g)) - carac(0);
f(1) = .5*(u(1) + n(1)*u(3)*sqrt(u(1)*w(1)/g)) - carac(1);
f(0) = .5*(u(0) + n(0)*u(2)*sqrt(u(0)*w(0)/_g)) - carac(0);
f(1) = .5*(u(1) + n(1)*u(3)*sqrt(u(1)*w(1)/_g)) - carac(1);
f(2) = n(0)*u(0)*u(2) + n(1)*u(1)*u(2);
f(3) = .5*u(0)*u(2) + g*eta(0) - .5*u(1)*u(3) - g*eta(1);
f(3) = .5*u(0)*u(2) + _g*eta(0) - .5*u(1)*u(3) - _g*eta(1);
}
};
void dgConservationLawShallowWater1dLaxSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu){
double cL = fabs(uL)+sqrt(g*sL/wL);
double cR = fabs(uR)+sqrt(g*sR/wR);
void dgConservationLawShallowWater1dLaxSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu, double g){
double _g = g;
double cL = fabs(uL)+sqrt(_g*sL/wL);
double cR = fabs(uR)+sqrt(_g*sR/wR);
double aMax = std::max(cL,cR);
fs = .5*( sL*uL + sR*uR )*n0 + aMax*(sL-sR) ;
fu = .5*( .5*uL*uL+pL + .5*uR*uR+pR )*n0 + aMax*(uL-uR);
}
void dgConservationLawShallowWater1dRoeSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu){
void dgConservationLawShallowWater1dRoeSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu, double g){
double _g = g;
double fLs = sL*uL*n0, fRs = sR*uR*n0;
double fLu = (.5*uL*uL+pL)*n0, fRu = (.5*uR*uR+pR)*n0;
double uroe = ( uL * sqrt(sL) + uR * sqrt(sR) ) / ( sqrt(sL) + sqrt(sR) );
double croe = sqrt( g * (sL/wL + sR/wR));
double croe = sqrt( _g * (sL/wL + sR/wR));
double sj = sR-sL;
double uj = uR-uL;
double wj = wR-wL;
......@@ -349,21 +355,23 @@ public:
class dgConservationLawShallowWater1d::boundaryWall : public dgBoundaryCondition {
class term : public function {
fullMatrix<double> normals, elevation;
double _g;
public:
term (const function *hydroSol, const function *elevationF) : function (2) {
term (const function *hydroSol, const function *elevationF, double g) : function (2) {
//setArgument (sol, hydroSol, 0);
setArgument (elevation, elevationF, 0);
setArgument (normals, function::getNormals(), 0);
_g = g;
}
void call (dataCacheMap *map, fullMatrix<double> &val) {
val(0,0) = 0;
val(0,1) = - g*elevation(0,0)*normals(0,0);
val(0,1) = - _g*elevation(0,0)*normals(0,0);
}
};
public:
term _term;
boundaryWall(dgConservationLawShallowWater1d *claw) : _term (claw->getHydroSolution(), claw->getElevation()) {
boundaryWall(dgConservationLawShallowWater1d *claw) : _term (claw->getHydroSolution(), claw->getElevation(), claw->_g) {
_term0 = &_term;
}
};
......@@ -378,15 +386,17 @@ class dgConservationLawShallowWater1d::forceSol:public dgBoundaryCondition {
fullMatrix<double> sectionExt,velocityExt;
fullMatrix<double> elevationExt, elevationExt_w, elevation;
functionReplace elevationExtReplace, widthExtReplace;
double _g;
public:
term(const function *sectionExtF, const function *velocityExtF, const function *solF, const function *elevationF, const function *widthF):function(2){
term(const function *sectionExtF, const function *velocityExtF, const function *solF, const function *elevationF, const function *widthF, double g):function(2){
setArgument(sectionExt, sectionExtF, 0);
setArgument(velocityExt, velocityExtF, 0);
setArgument(sol, solF, 0);
setArgument(normals, function::getNormals(), 0);
setArgument(elevation, elevationF, 0);
setArgument(width, widthF, 0);
_g = g;
addFunctionReplace(elevationExtReplace);
elevationExtReplace.replace(solExt_eta,solF);
elevationExtReplace.get(elevationExt, elevationF);
......@@ -402,21 +412,21 @@ class dgConservationLawShallowWater1d::forceSol:public dgBoundaryCondition {
solExt_eta(0,0) = sExt;
elevationExtReplace.compute();
double pIn = g*elevation(0,0), pExt = g*elevationExt(0,0);
double pIn = _g*elevation(0,0), pExt = _g*elevationExt(0,0);
elevationExt_w(0,0)=elevationExt(0,0);
widthExtReplace.compute();
double wIn = width(0,0), wExt = widthExt(0,0);
double F1,F2;
dgConservationLawShallowWater1dLaxSolver(sIn, sExt, uIn, uExt, wIn, wExt, pIn, pExt, n0, F1, F2);
dgConservationLawShallowWater1dLaxSolver(sIn, sExt, uIn, uExt, wIn, wExt, pIn, pExt, n0, F1, F2, _g);
val(0,0) = - F1;
val(0,1) = - F2;
}
};
public:
term _term;
forceSol (dgConservationLawShallowWater1d *claw, const function *forcedSection, const function *forcedVelocity) : _term (forcedSection,forcedVelocity,claw->getHydroSolution(), claw->getElevation(), claw->getWidth()) {
forceSol (dgConservationLawShallowWater1d *claw, const function *forcedSection, const function *forcedVelocity) : _term (forcedSection,forcedVelocity,claw->getHydroSolution(), claw->getElevation(), claw->getWidth(), claw->_g) {
_term0 = &_term;
}
};
......@@ -428,13 +438,15 @@ class dgConservationLawShallowWater1d::boundaryCouplingSW2D:public dgBoundaryCon
class term:public function {
fullMatrix<double> normals, solExt, sol;
double _width, _section0;
double _g;
public:
term(const function *solExtF, const function *solF, const double width, const double section0):function(2){
term(const function *solExtF, const function *solF, const double width, const double section0, double g):function(2){
setArgument(solExt, solExtF);
setArgument(sol, solF);
setArgument(normals, function::getNormals(), 0);
_width = width;
_section0 = section0;
_g = g;
}
void call (dataCacheMap *map, fullMatrix<double> &val) {
double uIn = sol(0,1), uExt = solExt(0,1);
......@@ -443,17 +455,17 @@ class dgConservationLawShallowWater1d::boundaryCouplingSW2D:public dgBoundaryCon
double HIn = sol(0,0)/_width, HExt = solExt(0,0)/_width;
double n0 = normals(0,0);
double uStar,vStar,HStar, AStar;
roeSolver(uExt, uIn, vExt, vIn, HExt, HIn, uStar, vStar, HStar, AStar);
roeSolver(uExt, uIn, vExt, vIn, HExt, HIn, uStar, vStar, HStar, AStar, _g);
double etaStar = HStar - _bath;
double ffs = n0 * uStar * HStar * _width;
double ffu = n0 * ( uStar * uStar / 2 + g * etaStar );
double ffu = n0 * ( uStar * uStar / 2 + _g * etaStar );
val(0,0) = -ffs;
val(0,1) = -ffu;
}
};
public:
term _term;
boundaryCouplingSW2D (dgConservationLawShallowWater1d *claw, const function *solExtF, double width, double section0): _term(solExtF, claw->getHydroSolution(), width, section0) {
boundaryCouplingSW2D (dgConservationLawShallowWater1d *claw, const function *solExtF, double width, double section0): _term(solExtF, claw->getHydroSolution(), width, section0, claw->_g) {
_term0 = &_term;
}
};
......@@ -466,29 +478,31 @@ class dgConservationLawShallowWater1d::forceUpstreamDischarge:public dgBoundaryC
fullMatrix<double> normals, solExt, sol, width;
fullMatrix<double> dischargeExt;
fullMatrix<double> elevation;
double _g;
public:
term(const function *dischargeExtF, const function *solF, const function *elevationF, const function *widthF):function(2){
term(const function *dischargeExtF, const function *solF, const function *elevationF, const function *widthF, double g):function(2){
setArgument(dischargeExt, dischargeExtF);
setArgument(sol, solF);
setArgument(normals, function::getNormals(), 0);
setArgument(elevation, elevationF);
setArgument(width, widthF);
_g = g;
}
void call (dataCacheMap *map, fullMatrix<double> &val) {
double sIn = sol(0,0), sExt = sIn;
double uIn = sol(0,1), uExt = -dischargeExt(0,0)/sIn;
double n0 = normals(0,0);
double pIn = g*elevation(0,0), pExt = pIn;
double pIn = _g*elevation(0,0), pExt = pIn;
double wIn = width(0,0), wExt = wIn;
double F1,F2;
dgConservationLawShallowWater1dLaxSolver(sIn, sExt, uIn, uExt, wIn, wExt, pIn, pExt, n0, F1, F2);
dgConservationLawShallowWater1dLaxSolver(sIn, sExt, uIn, uExt, wIn, wExt, pIn, pExt, n0, F1, F2, _g);
val(0,0) = - F1;
val(0,1) = - F2;
}
};
public:
term _term;
forceUpstreamDischarge (dgConservationLawShallowWater1d *claw, const function *forcedDischarge) : _term (forcedDischarge,claw->getHydroSolution(), claw->getElevation(), claw->getWidth()) {
forceUpstreamDischarge (dgConservationLawShallowWater1d *claw, const function *forcedDischarge) : _term (forcedDischarge,claw->getHydroSolution(), claw->getElevation(), claw->getWidth(), claw->_g) {
_term0 = &_term;
}
};
......@@ -507,9 +521,10 @@ class dgConservationLawShallowWater1d::riemannBifurcation:public function {
fullMatrix<double> elevation1,elevation2,elevation3;
fullMatrix<double> width1,width2,width3;
functionReplace replaceElevation1, replaceElevation2, replaceElevation3, replaceWidth1, replaceWidth2, replaceWidth3;
double _g;
public:
riemannBifurcation(const function *hydroSol, const function *elevation, const function *width, const function *widthData):function(6) {
riemannBifurcation(const function *hydroSol, const function *elevation, const function *width, const function *widthData, double g):function(6) {
setArgument (normals1, function::getNormals(), 0);
setArgument (normals2, function::getNormals(), 1);
setArgument (normals3, function::getNormals(), 2);
......@@ -522,6 +537,7 @@ public:
setArgument (width1, width, 0);
setArgument (width2, width, 1);
setArgument (width3, width, 2);
double _g = g;
addFunctionReplace(replaceElevation1);//The elevation depends on the solution
replaceElevation1.replace(sol1_eta, hydroSol, 0);
......@@ -569,9 +585,9 @@ public:
eta2_w(0,0)=eta(1); replaceWidth2.compute(); w(1)=w2(0,0);
eta3_w(0,0)=eta(2); replaceWidth3.compute(); w(2)=w3(0,0);
caracs(0) = 0.5*(u(0) + n(0)*u(3)*sqrt(u(0)*w(0)/g));
caracs(1) = 0.5*(u(1) + n(1)*u(4)*sqrt(u(1)*w(1)/g));
caracs(2) = 0.5*(u(2) + n(2)*u(5)*sqrt(u(2)*w(2)/g));
caracs(0) = 0.5*(u(0) + n(0)*u(3)*sqrt(u(0)*w(0)/_g));
caracs(1) = 0.5*(u(1) + n(1)*u(4)*sqrt(u(1)*w(1)/_g));
caracs(2) = 0.5*(u(2) + n(2)*u(5)*sqrt(u(2)*w(2)/_g));
double maxNorm=1e-8;
double eps=1e-3;
......@@ -579,7 +595,7 @@ public:
int nr;
for (nr=0; nr < 10; nr++){
systemRiemannBifurcation(u, eta, w, caracs, n, fu);
systemRiemannBifurcation(u, eta, w, caracs, n, fu, _g);
for(int i=0;i<6;i++){
u(i)+=delta[i];
sol1_eta(0,0)=u(0); replaceElevation1.compute(); eta(0)=eta1(0,0);
......@@ -588,7 +604,7 @@ public:
eta1_w(0,0)=eta(0); replaceWidth1.compute(); w(0)=w1(0,0);
eta2_w(0,0)=eta(1); replaceWidth2.compute(); w(1)=w2(0,0);
eta3_w(0,0)=eta(2); replaceWidth3.compute(); w(2)=w3(0,0);
systemRiemannBifurcation(u, eta, w, caracs, n, fub);
systemRiemannBifurcation(u, eta, w, caracs, n, fub, _g);
for(int j=0;j<6;j++)
dfdu.set(j,i,(fub(j)-fu(j))/delta[i]);
u(i)-=delta[i];
......@@ -607,20 +623,21 @@ public:
// if(nr==10) Msg::Fatal("Newton for bifurcation solver diverges | norm:%e u=%.3f %.3f %.3f eta=%.3f %.3f %.3f %s-%s-%s\n",du.norm(),u(3),u(4),u(5),eta(0),eta(1),eta(2),group1->getPhysicalTag().c_str(),group2->getPhysicalTag().c_str(),group3->getPhysicalTag().c_str());
val(0,0)=-u(3)*u(0)*n(0);
val(0,1)=-(.5*u(3)*u(3) + g*eta(0))*n(0);
val(0,1)=-(.5*u(3)*u(3) + _g*eta(0))*n(0);
val(0,2)=-u(4)*u(1)*n(1);
val(0,3)=-(.5*u(4)*u(4) + g*eta(1))*n(1);
val(0,3)=-(.5*u(4)*u(4) + _g*eta(1))*n(1);
val(0,4)=-u(5)*u(2)*n(2);
val(0,5)=-(.5*u(5)*u(5) + g*eta(2))*n(2);
val(0,5)=-(.5*u(5)*u(5) + _g*eta(2))*n(2);
}
void systemRiemannBifurcation(fullVector<double> &u, fullVector<double> &eta, fullVector<double> &w, fullVector<double> &carac, fullVector<double> &n, fullVector<double> &f){
f(0) = .5*(u(0) + n(0)*u(3)*sqrt(u(0)*w(0)/g)) - carac(0);
f(1) = .5*(u(1) + n(1)*u(4)*sqrt(u(1)*w(1)/g)) - carac(1);
f(2) = .5*(u(2) + n(2)*u(5)*sqrt(u(2)*w(2)/g)) - carac(2);
void systemRiemannBifurcation(fullVector<double> &u, fullVector<double> &eta, fullVector<double> &w, fullVector<double> &carac, fullVector<double> &n, fullVector<double> &f, double g){
double _g = g;
f(0) = .5*(u(0) + n(0)*u(3)*sqrt(u(0)*w(0)/_g)) - carac(0);
f(1) = .5*(u(1) + n(1)*u(4)*sqrt(u(1)*w(1)/_g)) - carac(1);
f(2) = .5*(u(2) + n(2)*u(5)*sqrt(u(2)*w(2)/_g)) - carac(2);
f(3) = n(0)*u(0)*u(3) + n(1)*u(1)*u(4) + n(2)*u(2)*u(5);
f(4) = .5*u(3)*u(3) + g*eta(0) - .5*u(4)*u(4) - g*eta(1);
f(5) = .5*u(3)*u(3) + g*eta(0) - .5*u(5)*u(5) - g*eta(2);
f(4) = .5*u(3)*u(3) + _g*eta(0) - .5*u(4)*u(4) - _g*eta(1);
f(5) = .5*u(3)*u(3) + _g*eta(0) - .5*u(5)*u(5) - _g*eta(2);
}
};
//**********************************************
......@@ -635,9 +652,9 @@ void dgConservationLawShallowWater1d::setWidth(){
void dgConservationLawShallowWater1d::setup () {
_volumeTerm0[""] = new source (_hydroSol, _linearDissipation, _quadraticDissipation,_sourceMass, _sourceMomentum, _diffusivity);
_volumeTerm1[""] = new advection (_hydroSol, _elevation,_diffusivity);
_interfaceTerm0[""] = new riemannLAX (_hydroSol,_width, _elevation, _widthData, _diffusivity);
_bifurcationTerm0[""] = new riemannBifurcation (_hydroSol,_elevation, _width, _widthData);
_volumeTerm1[""] = new advection (_hydroSol, _elevation,_diffusivity, _g);
_interfaceTerm0[""] = new riemannLAX (_hydroSol,_width, _elevation, _widthData, _diffusivity, _g);
_bifurcationTerm0[""] = new riemannBifurcation (_hydroSol,_elevation, _width, _widthData, _g);
if(_useClipToPhysics)
_clipToPhysics[""] = new clipToPhysics(_hydroSol, _clipSectionValue);
}
......
......@@ -3,26 +3,26 @@
#include "dgConservationLawFunction.h"
#include "functionGeneric.h"
/**function to compute the flux*/
void dgConservationLawShallowWater1dLaxSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu);
void dgConservationLawShallowWater1dRoeSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu);
void dgConservationLawShallowWater1dLaxSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu, double g);
void dgConservationLawShallowWater1dRoeSolver(double sL, double sR,double uL, double uR, double wL, double wR, double pL, double pR, double n0, double &fs, double &fu, double g);
/**The non-conservative shallow water conservation law. (section, u) */
class dgConservationLawShallowWater1d : public dgConservationLawFunction {
class advection;
class source;
class riemannLAX;
class riemannBifurcation;
class riemannBifurcation;
class clipToPhysics;
class boundaryWall;
class forceSol;
class forceSol;
class forceUpstreamDischarge;
class boundaryCouplingSW2D;
class elevationShallow;
class widthShallow;
class widthDataShallow;
class sectionDataShallow;
class widthShallow;
class widthDataShallow;
class sectionDataShallow;
const function *_hydroSol;
const function *_linearDissipation, *_quadraticDissipation;
const function *_linearDissipation, *_quadraticDissipation;
const function *_sourceMass, *_sourceMomentum;
const function *_diffusivity;
const function *_width, *_widthData, *_sectionData, *_elevation, *_elevationInterp;
......@@ -30,6 +30,7 @@ class dgConservationLawShallowWater1d : public dgConservationLawFunction {
std::string _typeRiemann;
bool _zeroFluxBifurcations;
double _clipSectionValue;
double _g;
bool _useClipToPhysics;
public:
void setup();
......@@ -47,6 +48,7 @@ class dgConservationLawShallowWater1d : public dgConservationLawFunction {
inline const function* getElevationInterp() {return _elevationInterp;};
void setElevation();
void setWidth();
void setGravity(double g){_g = g;}
/**set the function to evaluate the linear dissipation Kr u */
inline void setLinearDissipation(const function *linearDissipation){_linearDissipation = linearDissipation;}
/**set the function to evaluate the quadratic dissipation Kr u |u| */
......
......@@ -8,7 +8,6 @@
#include "slimMovingBathWettingDrying.h"
#include "functorMember.h"
static double g = 9.80616;
......@@ -99,7 +98,7 @@ class dgConservationLawShallowWater2d::maxDiffusiveSpeed: public function {
}*/
void call(dataCacheMap *m, fullMatrix<double> &val) {
/*for(int i=0; i< val.size1(); i++) {
val(i,0) = sqrt(bath(i,0)*g);
val(i,0) = sqrt(bath(i,0)*_g);
}*/
}
};
......@@ -108,13 +107,15 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
fullMatrix<double> bath, sol, _xyz;
bool _linear, _filterLinear;
double _R;
double _g;
public:
maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear):function(1){
maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear, double g):function(1){
setArgument(bath,bathymetry);
setArgument(sol,function::getSolution());
_linear = linear;
_filterLinear = filterLinear;
_R = R;
_g = g;
if(_R>0) {
setArgument(_xyz, function::getCoordinates());
}
......@@ -125,14 +126,14 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
val.setAll(0);
else
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(bath(i,0)*g);
val(i,0) = sqrt(bath(i,0)*_g);
}else
if (!_filterLinear)
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt((bath(i,0)+sol(i,0))*g)+hypot(sol(i,1),sol(i,2));
val(i,0) = sqrt((bath(i,0)+sol(i,0))*_g)+hypot(sol(i,1),sol(i,2));
else
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(sol(i,0)*g)+hypot(sol(i,1),sol(i,2));
val(i,0) = sqrt(sol(i,0)*_g)+hypot(sol(i,1),sol(i,2));
/*if (_R>0)
for (int i=0; i<val.size1(); i++) {
double x = _xyz(i,0);
......@@ -146,9 +147,9 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
class dgConservationLawShallowWater2d::gradPsiTerm: public function {
fullMatrix<double> bath, sol, diff, _xyz, movingBathFactor;
bool _linear, _isDiffusive;
double _R;
double _R, _g;
public:
gradPsiTerm (const functor *bathymetryF, const functor *diffusiveFluxF, bool linear, double R, const functor *movingBathFactorF, const functor *xyz):function(9){
gradPsiTerm (const functor *bathymetryF, const functor *diffusiveFluxF, bool linear, double R, const functor *movingBathFactorF, const functor *xyz, double g):function(9){
setArgument(bath,bathymetryF);
setArgument(sol,function::getSolution());
_isDiffusive = diffusiveFluxF;
......@@ -157,6 +158,7 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
}
_linear = linear;
_R = R;
_g = g;
if(_R>0) {
setArgument(_xyz, xyz ? xyz : function::getCoordinates());
}
......@@ -165,7 +167,7 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
void call(dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = val.size1();
double wd_g = g;
double wd_g = _g;
double J = 1.0;
for(size_t i=0; i< nQP; i++) {
if (_R > 0){
......@@ -214,14 +216,15 @@ class dgConservationLawShallowWater2d::source: public function {
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _nu, _xyz;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _absLayerCoef, _absLayerExtForc;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer;
double _R;
double _R, _g;
public :
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF,
const functor *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R,
const functor *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R, double g,
const functor *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *absLayerCoefF, const functor *absLayerExtForcF) :function (3){
_linear = linear;
_R = R;
_g = g;
_isDiffusive = nuF;
_useMovingBathWD = movingBathFactorF;
......@@ -274,8 +277,8 @@ class dgConservationLawShallowWater2d::source: public function {
val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + _source(i,0);
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + _source(i,1);
if (_R>0) {
val (i,1) -= g*(eta)*_xyz(i,0)/2.0/_R/_R;
val (i,2) -= g*(eta)*_xyz(i,1)/2.0/_R/_R;
val (i,1) -= _g*(eta)*_xyz(i,0)/2.0/_R/_R;
val (i,2) -= _g*(eta)*_xyz(i,1)/2.0/_R/_R;
}
if(_isDiffusive){
double H = _bathymetry(i,0);
......@@ -441,19 +444,19 @@ 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);
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;
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, movingBathFactorL(i,0),