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

Change function to functor in slimMovingBathWettingDrying and allowing for...

Change function to functor in slimMovingBathWettingDrying and allowing for using Lax Friedrichs for tracer equation (needed for tracer consistency if Lax Friedrichs is used for SW2d)

lez saisir le message de validation pour vos modifications. Les lignes
parent 66ba8445
......@@ -6,6 +6,7 @@ import param
import dgpy
do = ["gravity", "bathy", "forcing", "coriolis","manning"]
do = ["manning"]
check_result = True
......
......@@ -157,10 +157,6 @@ class ShallowWater2d:
self._wetting_drying = wetting_drying
self._equation.setBathymetry(self._domain._bath_function)
self._equation.setBathymetryGradient(self._domain._bath_gradient_function)
if self._wetting_drying is not None:
self._equation.setMovingBathWettingDrying(self._wetting_drying)
self._domain._bath_function = self._equation.getBathymetry()
self._domain._bath_gradient_function = self._equation.getBathymetryGradient()
self._equation.setGravity(self._domain._g)
self._coriolis = None
self._solution = dgpy.dgDofContainer(self._domain._groups, self._equation.getNbFields())
......@@ -168,6 +164,10 @@ class ShallowWater2d:
self._solution.setFieldName(0, 'eta')
self._solution.setFieldName(1, 'u')
self._solution.setFieldName(2, 'v')
if self._wetting_drying is not None:
self._equation.setMovingBathWettingDrying(self._solution.getFunction(), self._solution.getFunctionGradient(), self._wetting_drying)
self._domain._bath_function = self._equation.getBathymetry()
self._domain._bath_gradient_function = self._equation.getBathymetryGradient()
self._name = "sw2d"
if temporal_scheme == "explicit" or temporal_scheme == "implicit" or temporal_scheme == "multirate":
self._temporal_scheme = temporal_scheme
......@@ -603,13 +603,10 @@ class ShallowWaterTracer2d:
self._equation.setGravity(self._domain._g)
self._wetting_drying = wetting_drying
self._equation.setBathymetry(self._domain._bath_function)
if self._wetting_drying is not None:
#dgpy.Msg.Fatal("Wetting Drying is not implemented yet for tracer, Sorry") #TODO faire le wetting drying pour tracer
self._equation.setMovingBathWettingDrying(self._wetting_drying)
self._domain._bath_function = self._equation.getBathymetry()
self._solution = dgpy.dgDofContainer(self._domain._groups, 1)
self._solution.setAll(0.0)
self._solution.setFieldName(0, name)
self._hydro_sol = None
self._name = name
if temporal_scheme == "explicit" or temporal_scheme == "implicit":
self._temporal_scheme = temporal_scheme
......@@ -627,6 +624,9 @@ class ShallowWaterTracer2d:
if not isinstance(hydro_sol, ShallowWater2d):
dgpy.Msg.Fatal("You have to provide the hydrodynamic equation to solve the tracer!")
self._hydro_sol = hydro_sol._solution.getFunction()
if self._wetting_drying is not None:
self._equation.setMovingBathWettingDrying(self._hydro_sol, None, self._wetting_drying)
self._domain._bath_function = self._equation.getBathymetry()
self._equation.setHydroSolution(self._hydro_sol)
self._export_every_sub_time_step = export_every_sub_time_step
self._evaluator = None
......
......@@ -577,6 +577,8 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_fzerov = new functionConstant(std::vector<double>(2,0));
_bathymetry = NULL;
_bathymetryGradient = NULL;
_hydro = NULL;
_hydroGradient = NULL;
_originalBathymetry = NULL;
_originalBathymetryGradient = NULL;
_linearDissipation = _fzero;
......@@ -658,15 +660,17 @@ void dgConservationLawShallowWater2d::setup() {
_clipToPhysics[""] = new clipToPhysics(_bathymetry, 0.01);
}
void dgConservationLawShallowWater2d::setMovingBathWettingDrying(double alphaMovingBathWD){
void dgConservationLawShallowWater2d::setMovingBathWettingDrying(const functor *hydro, const functor *hydroGradient, 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);
_hydro = hydro;
_hydroGradient = hydroGradient;
_bathymetry = new movingBath(_originalBathymetry, _hydro, _alphaMovingBathWD);
_bathymetryGradient = new movingBathGradient(_originalBathymetry, _originalBathymetryGradient, _hydro, _hydroGradient, _alphaMovingBathWD);
_movingBathFactor = new movingBathFactor(_originalBathymetry, _hydro, _alphaMovingBathWD);
_movingBathFactorGradient = new movingBathFactorGradient(_originalBathymetry,_originalBathymetryGradient, _hydro, _hydroGradient, _alphaMovingBathWD);
}
dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
......
......@@ -37,6 +37,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
double _density;
const functor *_bathymetry, *_bathymetryGradient;
const functor *_hydro, *_hydroGradient;
const functor *_originalBathymetry, *_originalBathymetryGradient;
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_windStress, *_coriolisFactor, *_coordinatesF, *_g;
const functor *_movingBathFactor, *_movingBathFactorGradient;
......@@ -94,7 +95,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
_nudgingVelIsTransport = transport;
}
/**use the moving bath WD (Karna et al, 2011) and define the bathymetry and the moving bathymetry */
void setMovingBathWettingDrying(double alphaMovingBathWD);
void setMovingBathWettingDrying(const functor *hydro, const functor *hydroGradient, double alphaMovingBathWD);
/** a factor to smooth the if in Hv term in the roe riemann solver */
inline void setUpwindFactorRiemann(double upwindFactorRiemann) {
_upwindFactorRiemann=upwindFactorRiemann;
......
......@@ -105,19 +105,38 @@ void dgConservationLawShallowWaterTracer2d::riemannRoeCB(functorCache &m, fullMa
nyR = -normalsR(i,1);
}
//printf("nl=(%g,%g), nR=(%g,%g)\n",nxL,nyL,nxR,nyR);
double cL = solL(i,0), cR = solR(i,0);
double uL = nxL*hydroSolL(i,1) + nyL*hydroSolL(i,2);
double vL = -nyL*hydroSolL(i,1) + nxL*hydroSolL(i,2);
double uR = nxR*hydroSolR(i,1) + nyR*hydroSolR(i,2);
double vR = - nyR*hydroSolR(i,1) + nxR*hydroSolR(i,2);
double h = (bathL(i,0)+bathR(i,0))/2.;
double HR = hydroSolR(i,0) + h;
double HL = hydroSolL(i,0) + h;
if(HR<0 || HL<0) printf("in sw2d tracer Riemann: HR = %e HL = %e\n", HR,HL);
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), _upwindFactorRiemann);
double Feta = -H*u/J/Amb;
val(i,0) = Feta * (uL > 0 ? cL : cR) ;
double Feta, u;
if (_laxFriedrichs) {
double hL = bathL(i,0);
double hR = bathR(i,0);
double etaR = hydroSolR(i,0);
double etaL = hydroSolL(i,0);
const double etaM = (etaR + etaL) /2;
const double hMin = std::min(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 c = sqrt(g(i,0) * (hMin + etaM) + fabs(unM));
Feta = J * (- (hMin + etaM) * unM) + c * (etaR - etaL) / 2;
u =unM;
}
else{
double h = (bathL(i,0)+bathR(i,0))/2.;
double HR = hydroSolR(i,0) + h;
double HL = hydroSolL(i,0) + h;
if(HR<0 || HL<0) printf("in sw2d tracer Riemann: HR = %e HL = %e\n", HR,HL);
double v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactorRiemann);
Feta = -H*u/J/Amb;
}
double cL = solL(i,0), cR = solR(i,0);
val(i,0) = Feta * (u > 0 ? cL : cR) ;
val(i,1) = -val(i,0);
if(ip){
val(i,0) += (*ip)(i,0);
......@@ -177,8 +196,8 @@ dgBoundaryCondition *dgConservationLawShallowWaterTracer2d::newForceTracer(const
//**********************************************
void dgConservationLawShallowWaterTracer2d::setup () {
if(_useMovingBathWD){
_bathymetry = new movingBath(_originalBathymetry, _alphaMovingBathWD);
_movingBathFactor = new movingBathFactor(_originalBathymetry, _alphaMovingBathWD);
_bathymetry = new movingBath(_originalBathymetry, _hydro, _alphaMovingBathWD);
_movingBathFactor = new movingBathFactor(_originalBathymetry, _hydro, _alphaMovingBathWD);
}
_diffusivity[""] = _kappa;
......@@ -208,6 +227,7 @@ dgConservationLawShallowWaterTracer2d::dgConservationLawShallowWaterTracer2d() :
_alphaMovingBathWD = 0.0;
_useMovingBathWD = false;
_bathymetry = NULL;
_hydro = NULL;
_originalBathymetry = NULL;
_movingBathFactor = new functionConstant(1.);
_upwindFactorRiemann = -1;
......
......@@ -6,11 +6,12 @@
class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
const functor *_kappa;
const functor *_diffusion;
const functor *_source;
const functor *_source;
const functor *_hydroSolution;
const functor *_diffusiveFlux;
const functor *_ipTerm;
const functor *_movingBathFactor;
const functor *_hydro;
const functor *fzero, *fzerov;
const functor *_g;
double _R;
......@@ -19,6 +20,7 @@ class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
bool _useMovingBathWD;
double _alphaMovingBathWD;
bool _isConstantJac, _isLinear;
bool _laxFriedrichs;
void diffusiveFluxCB(functorCache &m, fullMatrix<double> &val) const;
void psiTermCB(functorCache &m, fullMatrix<double> &val) const;
......@@ -33,9 +35,11 @@ class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
bool isConstantJac() const {return _isConstantJac;}
bool isLinear() const {return _isLinear;}
void setIsLinear(bool isLinear) {_isLinear = isLinear;}
void setIsConstantJac(bool isConstantJac) {_isConstantJac = isConstantJac;}
/* sw2D solved with Lax Friedrichs*/
void hydroWithLaxFriedrichs(bool laxFriedrichs) { _laxFriedrichs = laxFriedrichs;}
/*set the gravity as a functor */
inline void setGravity(const functor *gravityF) { _g = gravityF;}
/**set the function to evaluate the source in the tracer equation */
......@@ -49,10 +53,11 @@ class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
inline void setUpwindFactorRiemann(double upwindFactorRiemann){
_upwindFactorRiemann=upwindFactorRiemann;
}
inline void setMovingBathWettingDrying(double alphaMovingBathWD){
inline void setMovingBathWettingDrying(const functor *hydro, const functor *hydroGradient, double alphaMovingBathWD){
_useMovingBathWD=true;
_alphaMovingBathWD = alphaMovingBathWD;
_originalBathymetry = _bathymetry;
_hydro = hydro;
}
inline const functor *getBathymetry() {
return _bathymetry;
......
#include "function.h"
#include "functor.h"
#include "slimMovingBathWettingDrying.h"
movingBathFactor::movingBathFactor(const functor *originalBathF, double alphaF):function(1){
setArgument(sol, function::getSolution());
setArgument(bath, originalBathF);
alpha = alphaF;
}
void movingBathFactor::call (dataCacheMap *m, fullMatrix<double> &val) {
for (int k = 0 ; k < bath.size1(); k++ ){
double H = bath(k,0)+sol(k,0);
val(k,0) = 0.5 * H / sqrt( H*H + alpha*alpha ) + 0.5;
movingBathFactor::movingBathFactor(const functor *originalBathF, const functor *hydro, double alphaF){
_hydro = hydro;
_bath = originalBathF;
_alpha = alphaF;
}
void movingBathFactor::operator() (functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &hydro = cm.get(*_hydro);
const fullMatrix<double> &bath = cm.get(*_bath);
val.resize(cm.nEvaluationPoint(),1,false);
for (int k = 0 ; k < cm.nEvaluationPoint(); k++ ){
double H = bath(k,0)+hydro(k,0);
val(k,0) = 0.5 * H / sqrt( H*H + _alpha*_alpha ) + 0.5;
}
}
movingBathFactorGradient::movingBathFactorGradient(const functor *originalBathF, const functor *originalBathGradientF, double alphaF):function(2) {
setArgument(sol, function::getSolution());
setArgument(solGradient, function::getSolutionGradient());
setArgument(bath, originalBathF);
setArgument(bathGradient, originalBathGradientF);
alpha = alphaF;
}
void movingBathFactorGradient::call (dataCacheMap *m, fullMatrix<double> &val){
for (int k = 0 ; k < bath.size1(); k++ ){
double H = bath(k,0) + sol(k,0);
double dHdx = bathGradient(k,0) + solGradient(k,0);
double dHdy = bathGradient(k,1) + solGradient(k,1);
val(k,0) = 0.5 * dHdx / sqrt( H*H + alpha*alpha ) - 0.5 * H*H * dHdx / pow(H*H + alpha*alpha, 3.0/2);
val(k,1) = 0.5 * dHdy / sqrt( H*H + alpha*alpha ) - 0.5 * H*H * dHdy / pow(H*H + alpha*alpha, 3.0/2);
movingBathFactorGradient::movingBathFactorGradient(const functor *originalBathF, const functor *originalBathGradientF, const functor *hydro, const functor *hydroGradient, double alphaF){
_hydro = hydro;
_hydroGradient = hydroGradient;
_bath = originalBathF;
_bathGradient = originalBathGradientF;
_alpha = alphaF;
}
void movingBathFactorGradient::operator() (functorCache &cm, fullMatrix<double> &val) const{
const fullMatrix<double> &hydro = cm.get(*_hydro);
const fullMatrix<double> &hydroGradient = cm.get(*_hydroGradient);
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &bathGradient = cm.get(*_bathGradient);
val.resize(cm.nEvaluationPoint(),2,false);
for (int k = 0 ; k < cm.nEvaluationPoint(); k++ ){
double H = bath(k,0) + hydro(k,0);
double dHdx = bathGradient(k,0) + hydroGradient(k,0);
double dHdy = bathGradient(k,1) + hydroGradient(k,1);
val(k,0) = 0.5 * dHdx / sqrt( H*H + _alpha*_alpha ) - 0.5 * H*H * dHdx / pow(H*H + _alpha*_alpha, 3.0/2);
val(k,1) = 0.5 * dHdy / sqrt( H*H + _alpha*_alpha ) - 0.5 * H*H * dHdy / pow(H*H + _alpha*_alpha, 3.0/2);
}
}
movingBath::movingBath(const functor *originalBathF, double alphaF):function(1){
setArgument(sol, function::getSolution());
setArgument(bath, originalBathF);
alpha = alphaF;
}
void movingBath::call (dataCacheMap *m, fullMatrix<double> &val) {
for (int k = 0 ; k < bath.size1(); k++ ){
double h = bath(k,0) + sol(k,0);
val(k,0) = bath(k,0) + 0.5 * ( sqrt( h*h + alpha*alpha ) - h );
movingBath::movingBath(const functor *originalBathF, const functor *hydro, double alphaF){
_hydro = hydro;
_bath = originalBathF;
_alpha = alphaF;
}
void movingBath::operator() (functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &hydro = cm.get(*_hydro);
const fullMatrix<double> &bath = cm.get(*_bath);
val.resize(cm.nEvaluationPoint(),1,false);
for (int k = 0 ; k < cm.nEvaluationPoint(); k++ ){
double h = bath(k,0) + hydro(k,0);
val(k,0) = bath(k,0) + 0.5 * ( sqrt( h*h + _alpha*_alpha ) - h );
}
}
movingBathGradient::movingBathGradient(const functor *originalBathF, const functor *originalBathGradientF, double alphaF):function(2){
setArgument(sol, function::getSolution());
setArgument(solGradient, function::getSolutionGradient());
setArgument(bath, originalBathF);
setArgument(bathGradient, originalBathGradientF);
alpha = alphaF;
}
void movingBathGradient::call (dataCacheMap *m, fullMatrix<double> &val) {
for (int k = 0 ; k < bath.size1(); k++ ){
double H = bath(k,0) + sol(k,0);
double HGradient_x = bathGradient(k,0) + solGradient(k,0);
double HGradient_y = bathGradient(k,1) + solGradient(k,1);
val(k,0) = bathGradient(k,0) + 0.5 * ( H * HGradient_x / sqrt( H*H + alpha*alpha ) - HGradient_x );
val(k,1) = bathGradient(k,1) + 0.5 * ( H * HGradient_y / sqrt( H*H + alpha*alpha ) - HGradient_y );
movingBathGradient::movingBathGradient(const functor *originalBathF, const functor *originalBathGradientF, const functor *hydro, const functor *hydroGradient, double alphaF){
_hydro = hydro;
_hydroGradient = hydroGradient;
_bath = originalBathF;
_bathGradient = originalBathGradientF;
_alpha = alphaF;
}
void movingBathGradient::operator() (functorCache &cm, fullMatrix<double> &val) const{
const fullMatrix<double> &hydro = cm.get(*_hydro);
const fullMatrix<double> &hydroGradient = cm.get(*_hydroGradient);
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &bathGradient = cm.get(*_bathGradient);
val.resize(cm.nEvaluationPoint(),2,false);
for (int k = 0 ; k < cm.nEvaluationPoint(); k++ ){
double H = bath(k,0) + hydro(k,0);
double HGradient_x = bathGradient(k,0) + hydroGradient(k,0);
double HGradient_y = bathGradient(k,1) + hydroGradient(k,1);
val(k,0) = bathGradient(k,0) + 0.5 * ( H * HGradient_x / sqrt( H*H + _alpha*_alpha ) - HGradient_x );
val(k,1) = bathGradient(k,1) + 0.5 * ( H * HGradient_y / sqrt( H*H + _alpha*_alpha ) - HGradient_y );
}
}
#ifndef _SLIM_MOVING_BATH_WETTING_DRYING_H_
#define _SLIM_MOVING_BATH_WETTING_DRYING_H_
#include "function.h"
#include "functor.h"
#include "stdint.h"
#include "dgDofContainer.h"
class movingBathFactor : public function {
fullMatrix<double> sol, bath;
class movingBathFactor : public functor {
const functor *_bath, *_hydro;
double _alpha;
public:
double alpha;
movingBathFactor(const functor *originalBathF, double alphaF);
void call (dataCacheMap *m, fullMatrix<double> &val);
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
movingBathFactor(const functor *originalBathF, const functor *hydro, double alphaF);
};
class movingBathFactorGradient : public function {
fullMatrix<double> sol, bath, solGradient, bathGradient;
double alpha;
class movingBathFactorGradient : public functor {
const functor *_hydro, *_bath, *_hydroGradient, *_bathGradient;
double _alpha;
public:
movingBathFactorGradient(const functor *originalBathF, const functor *originalBathGradientF, double alphaF);
void call (dataCacheMap *m, fullMatrix<double> &val);
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
movingBathFactorGradient(const functor *originalBathF, const functor *originalBathGradientF, const functor *hydro, const functor *hydroGradient, double alphaF);
};
class movingBath : public function {
fullMatrix<double> bath, sol;
class movingBath : public functor {
const functor *_bath, *_hydro;
double _alpha;
public:
double alpha;
movingBath(const functor *originalBathF, double alphaF);
void call (dataCacheMap *m, fullMatrix<double> &val);
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
movingBath(const functor *originalBathF, const functor *hydro, double alphaF);
};
class movingBathGradient : public function {
fullMatrix<double> bath, bathGradient, sol, solGradient;
class movingBathGradient : public functor {
const functor *_bath, *_bathGradient, *_hydro, *_hydroGradient;
double _alpha;
public:
double alpha;
movingBathGradient(const functor *originalBathF, const functor *originalBathGradientF, double alphaF);
void call (dataCacheMap *m, fullMatrix<double> &val);
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
movingBathGradient(const functor *originalBathF, const functor *originalBathGradientF, const functor *hydro, const functor *hydroGradient, double alphaF);
};
#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