Commit fcb54ee8 authored by David Vincent's avatar David Vincent Committed by Valentin Vallaeys
Browse files

modify dgConservationShallowWater2D to add a wind stress source term...

modify dgConservationShallowWater2D to add a wind stress source term (additionnally to setSource) in order to allow adding wind AND gravitational potential as source term
modification in slimFunction are added (smagorinsky, bottom friction ....) It is function created in this file to allow setting such term via the interface.
Modifying CMakeList to add slim.py among the dgpy script
parent a3052201
......@@ -6,6 +6,7 @@ set(SCRIPTS
Richards.py
SWGW.py
SWNI.py
slim.py
dgftp.py
extrude.py
extrudePeriodic.py
......
......@@ -166,7 +166,7 @@ class ShallowWater2d:
"""
_wind_stress = dgDofContainer(self._groups, 2)
_load(_wind_stress, wind_stress)
self._solution.setSource(_wind_stress.getFunction())
self._solution.setWindStress(_wind_stress.getFunction())
def set_gravitational_potential(self, gravitational_potential):
"""
......
......@@ -215,12 +215,12 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
class dgConservationLawShallowWater2d::source: public function {
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _nu, _xyz;
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _windStress, _nu, _xyz;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _absLayerCoef, _absLayerExtForc;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer;
double _R, _g;
public :
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF,
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF, const functor *windStressF,
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){
......@@ -245,6 +245,7 @@ class dgConservationLawShallowWater2d::source: public function {
Msg::Info("AbsLayer is used.");
}
setArgument(_source,sourceF);
setArgument(_windStress,windStressF);
setArgument(_bathymetry,bathymetryF);
if(_isDiffusive || !_linear){
setArgument(solGradient,function::getSolutionGradient());
......@@ -276,8 +277,8 @@ class dgConservationLawShallowWater2d::source: public function {
double u = sol(i,1);
double v = sol(i,2);
val (i,0) = 0;//1/_R/_R*(_xyz(i,0)*u + _xyz(i,1)*v);//0.;
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);
val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + _source(i,0) + _windStress(i,0);
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + _source(i,1) + _windStress(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;
......@@ -577,6 +578,7 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_coriolisFactor = _fzero;
_quadraticDissipation = _fzero;
_source = _fzerov;
_windStress = _fzerov;
_linear = false;
_constantJac = false;
_coordinatesF = NULL;
......@@ -628,11 +630,11 @@ void dgConservationLawShallowWater2d::setup() {
_diffusion = _nu ? new diffusivity(_nu) : NULL;
_ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL;
//_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear);
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz, _g);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false, _g);
_sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g);
_riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g);
......@@ -680,7 +682,7 @@ void dgConservationLawShallowWater2d::setAbsLayer(const std::string tag, const f
if (!absCoef || !extFunc)
Msg::Fatal("Dg/ShallowWater/AbsLayer: Bad 'absCoef' and/or 'extFunc'.");
if (_volumeTerm0[tag]) delete _volumeTerm0[tag];
_volumeTerm0[tag] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient,
_volumeTerm0[tag] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient,
_nu, false || _linear, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, absCoef, extFunc);
}
......
......@@ -37,7 +37,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
const functor *_bathymetry, *_bathymetryGradient;
const functor *_originalBathymetry, *_originalBathymetryGradient;
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_coordinatesF;
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_windStress, *_coriolisFactor, *_coordinatesF;
const functor *_movingBathFactor, *_movingBathFactorGradient;
const functor *_fzero, *_fzerov, *_nu, *_diffusiveFlux, *_diffusion, *_ipTerm, *_xyz;
......@@ -75,6 +75,8 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
inline void setQuadraticDissipation(const functor *quadraticDissipation) { _quadraticDissipation = quadraticDissipation;}
/**set the function to evaluate the source term \f{eqnarray*}\frac{du}{dt} &=& s(0)\\ \frac{dv}{dt} &=& s(1)\f} */
inline void setSource(const functor *source) { _source = source;}
/**set the function to evaluate the wind induced surface stress term \f{eqnarray*}\frac{du}{dt} &=& s(0)\\ \frac{dv}{dt} &=& s(1)\f} */
inline void setWindStress(const functor *windStress) { _windStress = windStress;}
/**set the function to evaluate the bathymetry h (H = h+eta) */
inline void setBathymetry(functor *bathymetry) { _bathymetry= bathymetry; }
/**set the function to evaluate the gradient of the bathymetry: grad(h) */
......
......@@ -1131,8 +1131,9 @@ void mergeVector::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
dischargeToBoundaryVelocity::dischargeToBoundaryVelocity(const functor *solution, const functor *discharge, const functor *section) {
dischargeToBoundaryVelocity::dischargeToBoundaryVelocity(const functor *solution, const functor *discharge, const functor *section, const functor *bathy) {
_section = section;
_bathy = bathy;
_solution = solution;
_discharge = discharge;
_normals = NULL;
......@@ -1140,12 +1141,13 @@ dischargeToBoundaryVelocity::dischargeToBoundaryVelocity(const functor *solution
void dischargeToBoundaryVelocity::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &solution = cm.get(*_solution);
const fullMatrix<double> &bathy = cm.get(*_bathy);
const fullMatrix<double> &discharge = cm.get(*_discharge);
const fullMatrix<double> &section = cm.get(*_section);
const fullMatrix<double> &normals = cm.get(*function::getNormals());
val.resize(cm.nEvaluationPoint(),3,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
double c = discharge(i,0) / section(i,0);
double c = discharge(i,0) / section(i,0)*bathy(i,0)/(bathy(i,0)+solution(i,0));
val(i,0) = solution(i,0);
val(i,1) = - normals(i,0) * c;
val(i,2) = - normals(i,1) * c;
......
......@@ -316,10 +316,10 @@ class mergeVector: public functor {
};
class dischargeToBoundaryVelocity: public functor {
const functor *_solution, *_normals, *_discharge, *_section;
const functor *_solution, *_normals, *_discharge, *_section, *_bathy;
public:
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
dischargeToBoundaryVelocity(const functor *solution, const functor *discharge, const functor *section);
dischargeToBoundaryVelocity(const functor *solution, const functor *discharge, const functor *section, const functor *bathy);
};
#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