Commit 1bd9dd71 authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

remove stereo

parent e6dfcdf1
Pipeline #1335 passed with stage
in 26 minutes and 55 seconds
...@@ -89,20 +89,15 @@ class dgConservationLawShallowWater2d::clipToPhysics : public function { ...@@ -89,20 +89,15 @@ class dgConservationLawShallowWater2d::clipToPhysics : public function {
}; };
class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function { class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
fullMatrix<double> bath, sol, _xyz, _g; fullMatrix<double> bath, sol, _g;
bool _linear, _filterLinear; bool _linear, _filterLinear;
double _R;
public: public:
maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear, const functor *g):function(1){ maxConvectiveSpeed (const functor *bathymetry, bool linear, bool filterLinear, const functor *g):function(1){
setArgument(bath,bathymetry); setArgument(bath,bathymetry);
setArgument(sol,function::getSolution()); setArgument(sol,function::getSolution());
_linear = linear; _linear = linear;
_filterLinear = filterLinear; _filterLinear = filterLinear;
_R = R;
setArgument(_g, g); setArgument(_g, g);
if(_R>0) {
setArgument(_xyz, function::getCoordinates());
}
} }
void call(dataCacheMap *m, fullMatrix<double> &val) { void call(dataCacheMap *m, fullMatrix<double> &val) {
if (_linear){ if (_linear){
...@@ -121,15 +116,6 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function { ...@@ -121,15 +116,6 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
else{ else{
for (int i=0; i<val.size1(); i++) for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(sol(i,0)*_g(i,0))+hypot(sol(i,1),sol(i,2)); } val(i,0) = sqrt(sol(i,0)*_g(i,0))+hypot(sol(i,1),sol(i,2)); }
/*if (_R>0)
for (int i=0; i<val.size1(); i++) {
double x = _xyz(i,0);
double y = _xyz(i,1);
double J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
val(i,0) /= J;
}*/
} }
}; };
...@@ -138,21 +124,13 @@ void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<d ...@@ -138,21 +124,13 @@ void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<d
const fullMatrix<double> &sol = cm.get(*function::getSolution()); const fullMatrix<double> &sol = cm.get(*function::getSolution());
const fullMatrix<double> *diff = _diffusiveFlux ? &cm.get(*_diffusiveFlux) : NULL; const fullMatrix<double> *diff = _diffusiveFlux ? &cm.get(*_diffusiveFlux) : NULL;
const fullMatrix<double> &g = cm.get(*_g); const fullMatrix<double> &g = cm.get(*_g);
const fullMatrix<double> *xyz = _R ? &cm.get(*function::getCoordinates()) : NULL;
const fullMatrix<double> &movingBathFactor = cm.get(*_movingBathFactor); const fullMatrix<double> &movingBathFactor = cm.get(*_movingBathFactor);
const fullMatrix<double> &atmPress = cm.get(*_atmPress); const fullMatrix<double> &atmPress = cm.get(*_atmPress);
const fullMatrix<double> &rhoSurf = cm.get(*_rhoSurf); const fullMatrix<double> &rhoSurf = cm.get(*_rhoSurf);
const fullMatrix<double> *etaOld = _linearSolverFrom3D ? &cm.get(*_etaOld) : NULL; const fullMatrix<double> *etaOld = _linearSolverFrom3D ? &cm.get(*_etaOld) : NULL;
size_t nQP = cm.nEvaluationPoint(); size_t nQP = cm.nEvaluationPoint();
val.resize(nQP, 9, false); val.resize(nQP, 9, false);
double J = 1.0;
for(size_t i=0; i< nQP; i++) { for(size_t i=0; i< nQP; i++) {
if (_R > 0){
double x = (*xyz)(i,0);
double y = (*xyz)(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
}
double H = bath(i,0); double H = bath(i,0);
double eta = sol(i,0); double eta = sol(i,0);
double u = sol(i,1); double u = sol(i,1);
...@@ -163,18 +141,18 @@ void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<d ...@@ -163,18 +141,18 @@ void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<d
else if (_linearSolverFrom3D) else if (_linearSolverFrom3D)
H+=(*etaOld)(i,0); H+=(*etaOld)(i,0);
// flux_x // flux_x
val(i,0) = H*u * J/A; val(i,0) = H*u /A;
val(i,1) = ( g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density )*J; val(i,1) = g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density;
val(i,2) = 0; val(i,2) = 0;
// flux_y // flux_y
val(i,3) = H*v * J/A; val(i,3) = H*v /A;
val(i,4) = 0; val(i,4) = 0;
val(i,5) = ( g(i,0)*(1+rhoSurf(i,0)/_density) *eta + atmPress(i,0)/_density )*J; val(i,5) = g(i,0)*(1+rhoSurf(i,0)/_density) *eta + atmPress(i,0)/_density;
if (!_linear && !_from3D) { if (!_linear && !_from3D) {
val(i,1) += u*u*J; val(i,1) += u*u;
val(i,2) += u*v*J; val(i,2) += u*v;
val(i,4) += v*u*J; val(i,4) += v*u;
val(i,5) += v*v*J; val(i,5) += v*v;
} }
if (diff) { if (diff) {
val(i,1) += (*diff)(i,1); val(i,1) += (*diff)(i,1);
...@@ -203,31 +181,19 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double ...@@ -203,31 +181,19 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double
const fullMatrix<double> *solGradient = (_nu || !_linear) ? &cm.get(*function::getSolutionGradient()) : NULL; const fullMatrix<double> *solGradient = (_nu || !_linear) ? &cm.get(*function::getSolutionGradient()) : NULL;
const fullMatrix<double> *nu = _nu ? &cm.get(*_nu) : NULL; const fullMatrix<double> *nu = _nu ? &cm.get(*_nu) : NULL;
const fullMatrix<double> *bathymetryGradient = _nu ? &cm.get(*_bathymetryGradient) : NULL; const fullMatrix<double> *bathymetryGradient = _nu ? &cm.get(*_bathymetryGradient) : NULL;
const fullMatrix<double> *xyz = _R ? &cm.get(*function::getCoordinates()) : NULL;
const fullMatrix<double> *movingBathFactor = _useMovingBathWD ? &cm.get(*_movingBathFactor) : NULL; const fullMatrix<double> *movingBathFactor = _useMovingBathWD ? &cm.get(*_movingBathFactor) : NULL;
const fullMatrix<double> *movingBathFactorGradient = _useMovingBathWD ? &cm.get(*_movingBathFactorGradient) : NULL; const fullMatrix<double> *movingBathFactorGradient = _useMovingBathWD ? &cm.get(*_movingBathFactorGradient) : NULL;
const fullMatrix<double> &rhoSurfGrad = cm.get(*_rhoSurfGrad); const fullMatrix<double> &rhoSurfGrad = cm.get(*_rhoSurfGrad);
size_t nQP = cm.nEvaluationPoint(); size_t nQP = cm.nEvaluationPoint();
val.resize(nQP, 3, false); val.resize(nQP, 3, false);
double J = 1.0;
for(size_t i=0; i< nQP; i++) { for(size_t i=0; i< nQP; i++) {
if(_R > 0){
double x = (*xyz)(i,0);
double y = (*xyz)(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
}
double eta = sol(i,0); double eta = sol(i,0);
double u = sol(i,1); double u = sol(i,1);
double v = sol(i,2); double v = sol(i,2);
double H = eta + bathymetry(i,0); double H = eta + bathymetry(i,0);
val (i,0) = 0;//1/_R/_R*(_xyz(i,0)*u + _xyz(i,1)*v);//0.; val (i,0) = 0;
val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + source(i,0) + windStress(i,0)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,0)/_density; val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + source(i,0) + windStress(i,0)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,0)/_density;
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + source(i,1) + windStress(i,1)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,1)/_density; val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + source(i,1) + windStress(i,1)/(H*_density) + g(i,0)*eta*rhoSurfGrad(i,1)/_density;
if (_R>0) {
val (i,1) -= g(i,0)*(eta)*(*xyz)(i,0)/2.0/_R/_R;
val (i,2) -= g(i,0)*(eta)*(*xyz)(i,1)/2.0/_R/_R;
}
if(_nu){ if(_nu){
double H = bathymetry(i,0); double H = bathymetry(i,0);
double dHdx = (*bathymetryGradient)(i, 0); double dHdx = (*bathymetryGradient)(i, 0);
...@@ -243,12 +209,8 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double ...@@ -243,12 +209,8 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double
if(!_linear){ if(!_linear){
double div = (*solGradient)(i,3) + (*solGradient)(i,7); double div = (*solGradient)(i,3) + (*solGradient)(i,7);
double bottomFriction = hypot(u,v)*quadraticDissipation(i,0); double bottomFriction = hypot(u,v)*quadraticDissipation(i,0);
val(i,1) += u*(- bottomFriction + div/J); val(i,1) += u*(- bottomFriction + div);
val(i,2) += v*(- bottomFriction + div/J); val(i,2) += v*(- bottomFriction + div);
if (_R > 0) {
val(i,1)+= (-(*xyz)(i,0)*(u*u+v*v))/(2.0*_R*_R);
val(i,2)+= (-(*xyz)(i,1)*(u*u+v*v))/(2.0*_R*_R);
}
} }
if(_useMovingBathWD){ if(_useMovingBathWD){
double H = bathymetry(i,0) + eta; double H = bathymetry(i,0) + eta;
...@@ -257,9 +219,6 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double ...@@ -257,9 +219,6 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double
double dAdy = (*movingBathFactorGradient)(i,1); double dAdy = (*movingBathFactorGradient)(i,1);
val(i,0) += - H / ( A*A ) * ( dAdx * u + dAdy * v); val(i,0) += - H / ( A*A ) * ( dAdx * u + dAdy * v);
} }
val(i,0) *= J*J;
val(i,1) *= J*J;
val(i,2) *= J*J;
if (_nudgingCoeff) { if (_nudgingCoeff) {
double uref = (*nudgingVel)(i,1); double uref = (*nudgingVel)(i,1);
double vref = (*nudgingVel)(i,2); double vref = (*nudgingVel)(i,2);
...@@ -285,27 +244,6 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double ...@@ -285,27 +244,6 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double
} }
} }
class massFactor: public function {
fullMatrix<double> _xyz;
double _R;
public :
massFactor(double R, const functor *xyz) : function(3) {
_R = R;
setArgument(_xyz, xyz ? xyz : function::getCoordinates());
}
void call(dataCacheMap *m, fullMatrix<double> &val) {
int nQP=val.size1();
for (int i=0; i<nQP; i++) {
double x = _xyz(i,0);
double y = _xyz(i,1);
double J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
val(i,0) = J*J;
val(i,1) = J*J;
val(i,2) = J*J;
}
}
};
void dgConservationLawShallowWater2d::diffusiveFlux(functorCache &cm,fullMatrix<double> &val) const { void dgConservationLawShallowWater2d::diffusiveFlux(functorCache &cm,fullMatrix<double> &val) const {
const fullMatrix<double> &solGrad = cm.get(*function::getSolutionGradient()); const fullMatrix<double> &solGrad = cm.get(*function::getSolutionGradient());
const fullMatrix<double> &nu = cm.get(*_nu); const fullMatrix<double> &nu = cm.get(*_nu);
...@@ -353,25 +291,7 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -353,25 +291,7 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
bool _isDiffusive = _ipTerm; bool _isDiffusive = _ipTerm;
double _upwindFactor = _upwindFactorRiemann; double _upwindFactor = _upwindFactorRiemann;
const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL; const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL;
const fullMatrix<double> *xyz = _R>0 ? &cache.get(*_xyz) : NULL;
int rev = 1;
double J = 1.0;
const dgGroupOfFaces &iGroup = *cache.interfaceGroup();
int connectionR = iGroup.nConnection() == 2 ? 1 : 0;
const dgGroupOfElements *groupR = &iGroup.elementGroup(connectionR);
const dgGroupOfElements *groupL = &iGroup.elementGroup(0);
for(size_t i=0; i< nQP; i++) { for(size_t i=0; i< nQP; i++) {
rev = 1;
if(_R > 0)
{
double x = (*xyz)(i,0);
double y = (*xyz)(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
if(groupR->getPhysicalTag() != groupL->getPhysicalTag())
rev = -1;
}
double pa = 0.5*(atmPressL(i,0) + atmPressR(i,0)); double pa = 0.5*(atmPressL(i,0) + atmPressR(i,0));
double rhoSurf = 0.5*(rhoSurfL(i,0) + rhoSurfR(i,0)); double rhoSurf = 0.5*(rhoSurfL(i,0) + rhoSurfR(i,0));
double Fun, Fut, Feta; double Fun, Fut, Feta;
...@@ -387,20 +307,14 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -387,20 +307,14 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
nxR = -1.*normalsR(i,0); nxR = -1.*normalsR(i,0);
nyR = -1.*normalsR(i,1); nyR = -1.*normalsR(i,1);
} }
/* if(rev<0)
{
double RR = sqrt(_xyz.get(i,0)*_xyz.get(i,0) + _xyz.get(i,1)*_xyz.get(i,1));
printf("nx = %.16g nx2 = %.16g, ny = %.16g ny2 = %.16g, 2R = %.16g , RR = %.16g\n ", nx,_xyz.get(i,0)/RR,ny,_xyz.get(i,1)/RR, 2*_R, RR);
}*/
double uL = nxL*solL(i,1) + nyL*solL(i,2); double uL = nxL*solL(i,1) + nyL*solL(i,2);
double uR = nxR*solR(i,1) + nyR*solR(i,2); double uR = nxR*solR(i,1) + nyR*solR(i,2);
double vL = -nyL*solL(i,1) + nxL*solL(i,2); double vL = -nyL*solL(i,1) + nxL*solL(i,2);
double vR = -nyR*solR(i,1) + nxR*solR(i,2); double vR = -nyR*solR(i,1) + nxR*solR(i,2);
double hL = bathL(i,0); double hL = bathL(i,0);
double hR = bathR(i,0); double hR = bathR(i,0);
uR *= rev;
if (_laxFriedrichs) { if (_laxFriedrichs) {
//double Amb = (movingBathFactorL(i,0) + movingBathFactorR(i,0))/2.;
double etaR = solR(i,0); double etaR = solR(i,0);
double etaL = solL(i,0); double etaL = solL(i,0);
const double etaM = (etaR + etaL) /2; const double etaM = (etaR + etaL) /2;
...@@ -422,12 +336,13 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -422,12 +336,13 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
utM = (vR + vL) /2; utM = (vR + vL) /2;
c = sqrt(g(i,0) * H + fabs(unM)); c = sqrt(g(i,0) * H + fabs(unM));
} }
Fun = J * (- g(i,0) * (1+rhoSurf/_density) * etaM - pa/_density ) + c * (uR - uL) / 2; Fun = - g(i,0) * (1+rhoSurf/_density) * etaM - pa/_density + c * (uR - uL) / 2;
Feta = J * (- H * unM) + c * (etaR - etaL) / 2; Feta = - H * unM + c * (etaR - etaL) / 2;
Fut = 0.; Fut = 0.;
if (!_linear && !_from3D){ if (!_linear && !_from3D){
Fun += J * (- unM * unM); Fun += - unM * unM;
Fut += J * (- unM * utM) + c * (vR - vL) / 2; //Feta /= Amb; TODO verify if Wetting and drying and Lax Friedrichs !
Fut += - unM * utM + c * (vR - vL) / 2; // jump only if advection
} }
} }
else { else {
...@@ -441,9 +356,9 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -441,9 +356,9 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double sq_g_h = sqrt(g(i,0)/hMin); double sq_g_h = sqrt(g(i,0)/hMin);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2; double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2; double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -(g(i,0)*(1+rhoSurf/_density)*etaRiemann + pa/_density)*J; Fun = -g(i,0)*(1+rhoSurf/_density)*etaRiemann + pa/_density;
Fut = 0; Fut = 0;
Feta = -(hMin)*unRiemann*J; Feta = -hMin*unRiemann;
}else{ }else{
const double hMin = std::min(hL, hR); const double hMin = std::min(hL, hR);
double HR = solR(i,0) + hMin, HL = solL(i,0) + hMin; double HR = solR(i,0) + hMin, HL = solL(i,0) + hMin;
...@@ -453,9 +368,9 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -453,9 +368,9 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double u,v,H,Amb; double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor); roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
double eta = H-hMin; double eta = H-hMin;
Fun = -(g(i,0)*(1+rhoSurf/_density)*eta + pa/_density)*J - u*u*J; Fun = -g(i,0)*(1+rhoSurf/_density)*eta + pa/_density - u*u;
Fut = -u*v*J; Fut = -u*v;
Feta = -H*u*J/Amb; Feta = -H*u/Amb;
} }
} }
val(i,0) = Feta; val(i,0) = Feta;
...@@ -467,8 +382,8 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou ...@@ -467,8 +382,8 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
if (_isDiffusive) { if (_isDiffusive) {
val(i,1)+=(*ip)(i,1); val(i,1)+=(*ip)(i,1);
val(i,2)+=(*ip)(i,2); val(i,2)+=(*ip)(i,2);
val(i,4)+=(*ip)(i,4)*rev; val(i,4)+=(*ip)(i,4);
val(i,5)+=(*ip)(i,5)*rev; val(i,5)+=(*ip)(i,5);
} }
} }
} }
...@@ -498,11 +413,6 @@ void roeSolver(double uL, double uR, double vL, double vR, double HL, double HR, ...@@ -498,11 +413,6 @@ void roeSolver(double uL, double uR, double vL, double vR, double HL, double HR,
vStar = Hv / HStar; vStar = Hv / HStar;
AStar = (sqHL * mbFactL + sqHR * mbFactL)/(sqHL + sqHR); AStar = (sqHL * mbFactL + sqHR * mbFactL)/(sqHL + sqHR);
} }
void dgConservationLawShallowWater2d::setCoordinatesFunction(functor *xyz){
_xyz=xyz;
}
/*============================================================================== /*==============================================================================
* Conservation law : constructor - setup - destructor * Conservation law : constructor - setup - destructor
...@@ -525,7 +435,6 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat ...@@ -525,7 +435,6 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_windStress = _fzerov; _windStress = _fzerov;
_linear = false; _linear = false;
_constantJac = false; _constantJac = false;
_R = -1.0;
static const functor *gdefault = new functionConstant(9.80616); static const functor *gdefault = new functionConstant(9.80616);
_g = gdefault; _g = gdefault;
_density = 1025; _density = 1025;
...@@ -589,7 +498,7 @@ void dgConservationLawShallowWater2d::setup() { ...@@ -589,7 +498,7 @@ void dgConservationLawShallowWater2d::setup() {
_sourceTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::source); _sourceTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::source);
_gradPsiTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::gradPsiTerm); _gradPsiTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::gradPsiTerm);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann); _riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false, _g); _maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, false, _g);
_volumeTerm0[""] =_sourceTerm; _volumeTerm0[""] =_sourceTerm;
_volumeTerm1[""]=_gradPsiTerm; _volumeTerm1[""]=_gradPsiTerm;
_interfaceTerm0[""] = _riemannTerm; _interfaceTerm0[""] = _riemannTerm;
...@@ -597,11 +506,9 @@ void dgConservationLawShallowWater2d::setup() { ...@@ -597,11 +506,9 @@ void dgConservationLawShallowWater2d::setup() {
// _sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel, _nudgingVelIsTransport); // _sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel, _nudgingVelIsTransport);
// _gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g); // _gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g);
// _riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann); // _riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
// _maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g); // _maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, true, _g);
// setImexMode(_imexMode); // setImexMode(_imexMode);
// setLinearFilterMode(_linearFilter); // setLinearFilterMode(_linearFilter);
/*if (_R>0)
_massFactor[""] = std::make_pair(new massFactor(_R, _xyz), true);*/
//_interfaceTerm0[""] = _nu ? : new riemann(_bathymetry,NULL, _linear); //_interfaceTerm0[""] = _nu ? : new riemann(_bathymetry,NULL, _linear);
...@@ -645,75 +552,6 @@ dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() { ...@@ -645,75 +552,6 @@ dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
// BC : Wall // BC : Wall
/*class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition {
class term : public function {
fullMatrix<double> normals, sol, _xyz, bath;
bool _linear;
double _R;
public:
term(bool linear, double R, const functor *bathF) : function(3) {
setArgument(sol,function::getSolution());
setArgument(_xyz,function::getCoordinates());
setArgument(normals,function::getNormals());
setArgument(bath,bathF);
_linear = linear;
_R = R;
}
void call(dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = sol.size1();
for (size_t i=0; i<nQP; i++) {
double nx = normals(i,0);
double ny = normals(i,1);
double eta = sol(i,0);
double J = 1.0;
if (_R>0) {
double x = _xyz(i,0);
double y = _xyz(i,1);
J = 4*_R*_R/(4*_R*_R+x*x+y*y);
J = 1.0;
}
double u = sol(i, 1);
double v = sol(i, 2);
double un = u * nx + v * ny;
// val(i, 0) = 0;
// val(i, 1) = -_g*eta*nx/J;
// val(i, 2) = -_g*eta*ny/J;
// if(!_linear && un > 0){
// val(i, 1)-=un*u/J;
// val(i, 2)-=un*v/J;
// }
if ( _linear) {
val(i, 0) = 0;//bath(i, 0) * un * J;
val(i, 1) = -_g*eta*nx;
val(i, 2) = -_g*eta*ny;
} else {
double uL = un;
double vL = u * (-ny) + v * nx;
double h = bath(i,0);
double HL = eta + h;
double HStar,uStar,vStar,Amb;
//roeSolver(uL, 0, vL, 0, HL, HL, uStar, vStar, HStar, Amb, _g);
roeSolver(uL, -uL, vL, vL, HL, HL, uStar, vStar, HStar, Amb, _g);
val(i,0) = 0;//(-HStar * uStar + (HL * un)) * J;
double FunL = - _g * ((HStar -h)) * J - (uStar*uStar) * J;
double FutL = - (uStar*vStar) * J;
val(i,1) = FunL * nx - FutL * ny;
val(i,2) = FunL * ny + FutL * nx;
// if (uStar>0) {
// val(i, 1) = (uStar * uStar - uL*uL) / J;
// val(i, 2) = (uStar * vStar - uL*vL) / J;
//
}
}
}
};
term _term;
public:
boundaryWall(dgConservationLawShallowWater2d *claw) : _term(claw->isLinear(), claw->getRadius(), claw->getBathymetry()) {
_term0 = &_term;
}
};*/
void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue(functorCache &cm, fullMatrix<double> &val) const { void dgConservationLawShallowWater2d::dgConservationLawShallowWater2dWallExtValue(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0); const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0);
const fullMatrix<double> &solIn = cm.get(*function::getSolution(), 0); const fullMatrix<double> &solIn = cm.get(*function::getSolution(), 0);
...@@ -803,19 +641,17 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(bool slip) ...@@ -803,19 +641,17 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(bool slip)
class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBoundaryCondition { class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBoundaryCondition {
class term : public function { class term : public function {
fullMatrix<double> normals, sol, _xyz, bath, discharge, mbFact, _g; fullMatrix<double> normals, sol, bath, discharge, mbFact, _g;
bool _linear; bool _linear;
double _R, S0; double S0;
public: public:
term(bool linear, double R, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, const functor *g):function(3){ term(bool linear, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, const functor *g):function(3){
setArgument(sol,function::getSolution()); setArgument(sol,function::getSolution());
setArgument(_xyz,function::getCoordinates());
setArgument(normals,function::getNormals()); setArgument(normals,function::getNormals());
setArgument(bath,bathF); setArgument(bath,bathF);
setArgument(discharge,dischargeF); setArgument(discharge,dischargeF);
setArgument(mbFact,movingBathF); setArgument(mbFact,movingBathF);
_linear = linear; _linear = linear;
_R = R;
setArgument(_g, g); setArgument(_g, g);
dgFunctionIntegratorInterface integrator(solution->getGroups(), bathF, solution); dgFunctionIntegratorInterface integrator(solution->getGroups(), bathF, solution);
fullMatrix<double> Smat; fullMatrix<double> Smat;
...@@ -823,15 +659,8 @@ class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBounda ...@@ -823,15 +659,8 @@ class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBounda
S0 = Smat(0,0); S0 = Smat(0,0);
} }
void call(dataCacheMap *m, fullMatrix<double> &val) { void call(dataCacheMap *m, fullMatrix<double> &val) {
double J = 1;
double Fun, Fut, Feta; double Fun, Fut, Feta;
for (int i=0; i<sol.size1(); i++) { for (int i=0; i<sol.size1(); i++) {
if (_R>0) {
double x = _xyz(i,0);