Commit 9f4477ad authored by David Vincent's avatar David Vincent Committed by Jonathan Lambrechts
Browse files

shallowWater : ip on local basis

parent e44754bb
Pipeline #452 failed with stage
in 25 minutes and 27 seconds
......@@ -8,6 +8,8 @@
#include "slimMovingBathWettingDrying.h"
#include "functorMember.h"
#include "functor.h"
......@@ -635,7 +637,6 @@ void dgConservationLawShallowWater2d::setup() {
_maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g);
setImexMode(_imexMode);
setLinearFilterMode(_linearFilter);
/*if (_R>0)
_massFactor[""] = std::make_pair(new massFactor(_R, _xyz), true);*/
//_interfaceTerm0[""] = _nu ? : new riemann(_bathymetry,NULL, _linear);
......@@ -929,3 +930,102 @@ public:
dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryCouplingSW1D(const functor *solExtF){
return new boundaryCouplingSW1D(this, solExtF);
}
/*==============================================================================
* IP Term
*============================================================================*/
class dgIPTermOnSphere : public functor {
const functor *_diffusiveFlux, *_diffusivity;
bool _isotropic;
public:
void operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &diffusiveFluxL = cm.get(*_diffusiveFlux, 0);
const fullMatrix<double> &diffusiveFluxR = cm.get(*_diffusiveFlux, 1);
const fullMatrix<double> &diffusivityL = cm.get(*_diffusivity, 0);
const fullMatrix<double> &diffusivityR = cm.get(*_diffusivity, 1);
const fullMatrix<double> &solutionL = cm.get(cm.solution(), 0);
const fullMatrix<double> &solutionR = cm.get(cm.solution(), 1);
const fullMatrix<double> &normalsL = cm.get(*function::getNormals(),0);
const fullMatrix<double> &normalsR = cm.get(*function::getNormals(),1);
int nbFields = solutionL.size2();
fullMatrix<double> ipf = dgConservationLawFunction::muFactor(&cm);
val.resize(cm.nEvaluationPoint(), nbFields * 2);
double nxR, nyR, nxL, nyL;
for (int iPt = 0; iPt < cm.nEvaluationPoint(); ++iPt) {
// the sign of the right normals changes because the basis are different (except for element on the boundary)
if (cm.interfaceGroup()->nConnection() ==1){
nxR = normalsR(iPt,0);
nyR = normalsR(iPt,1);
}
else {
nxR = -1.*normalsR(iPt,0);
nyR = -1.*normalsR(iPt,1);
}
nxL = normalsL(iPt,0);
nyL = normalsL(iPt,1);
double mufactor = ipf(iPt, 0);
// viscosities for u and v
double nuU = std::max(diffusivityR(iPt, 1), diffusivityL(iPt, 1));
double nuV = std::max(diffusivityR(iPt, 2), diffusivityL(iPt, 2));
// left velocities
double unL = nuU * solutionL(iPt,1)*nxL + nuV * solutionL(iPt,2)*nyL;
double utL = nuV * solutionL(iPt,2)*nxL - nuU * solutionL(iPt,1)*nyL;
double unxL = diffusiveFluxL(iPt,2)*nxL + diffusiveFluxL(iPt,4)*nyL;
double unyL = diffusiveFluxL(iPt,3)*nxL + diffusiveFluxL(iPt,5)*nyL;
double utxL = diffusiveFluxL(iPt,4)*nxL - diffusiveFluxL(iPt,2)*nyL;
double utyL = diffusiveFluxL(iPt,5)*nxL - diffusiveFluxL(iPt,3)*nyL;
double unnL = unxL*nxL + unyL*nyL;
double utnL = utxL*nxL + utyL*nyL;
// right velocities
double unR = nuU * solutionR(iPt,1)*nxR + nuV * solutionR(iPt,2)*nyR;
double utR = nuV * solutionR(iPt,2)*nxR - nuU * solutionR(iPt,1)*nyR;
double unxR = diffusiveFluxR(iPt,2)*nxR + diffusiveFluxR(iPt,4)*nyR;
double unyR = diffusiveFluxR(iPt,3)*nxR + diffusiveFluxR(iPt,5)*nyR;
double utxR = diffusiveFluxR(iPt,4)*nxR - diffusiveFluxR(iPt,2)*nyR;
double utyR = diffusiveFluxR(iPt,5)*nxR - diffusiveFluxR(iPt,3)*nyR;
double unnR = unxR*nxR + unyR*nyR;
double utnR = utxR*nxR + utyR*nyR;
//mean velocities and jump
double meanNormalFlux = 0.5*(unnR+unnL);
double meanTangentFlux = 0.5*(utnR+utnL);
double jumpNormalVel = 0.5*(unL-unR)* mufactor;// nu is in un and ut
double jumpTangentVel = 0.5*(utL-utR)* mufactor;
double Fun = meanNormalFlux + jumpNormalVel;
double Fut = meanTangentFlux + jumpTangentVel;
double FuL = Fun * nxL + Fut * nyL;
double FvL = -Fun * nyL + Fut * nxL;
double FuR = Fun * nxR + Fut * nyR;
double FvR = -Fun * nyR + Fut * nxR;
//left values
val(iPt,0) = 0.;
val(iPt,1) = FuL;
val(iPt,2) = FvL;
//right values
val(iPt,3) = 0.;
val(iPt,4) = -FuR;
val(iPt,5) = -FvR;
}
}
dgIPTermOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity, bool isotropic){
_diffusiveFlux = diffusiveFlux;
_diffusivity = diffusivity;
_isotropic = isotropic;
}
};
functor *dgNewIpTermIsotropicOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity) {
return new dgIPTermOnSphere (nbFields, diffusiveFlux, diffusivity, true);
}
......@@ -112,4 +112,9 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
dgBoundaryCondition *newForcedDischarge(dgDofContainer *solution, const functor *discharge, std::string tag);
dgBoundaryCondition *newBoundaryCouplingSW1D(const functor *solExtF);
};
/*==============================================================================
* IP Term
*============================================================================*/
functor *dgNewIpTermIsotropicOnSphere (int nbFields, const functor *diffusiveFlux, const functor *diffusivity);
#endif
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