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

modifying slim.py to allow solving the sw2D on a manifold

modify shallowWater to multiply wind_stress by 1/(rho*H) --> add setDensity to the sw2d equation
parent 7e7a36e8
import slim
#slim.set_output_directory("output")
domain = slim.Domain("mesh.msh","bath_COMP_0.msh", g=9.81, order=1)
import dgpy
domain = slim.Domain("mesh.msh","bath_COMP_0.msh", g=9.81, order=1, solve_on_manifold=True)
#PREPRO domain.set_projection_UTM(55, a= 6370000, b= 6370000)
# PREPRO domain.set_projection("+proj=stere +ellps=WGS84 +lat_ts=69 +lat_0=90 +lon_0=85E +a=6370000 +b=6370000")
eq = slim.ShallowWater2d(domain, linear_equation = False, wetting_drying = None)
eq = slim.ShallowWater2d(domain, linear_equation = False, wetting_drying = None, export_every_sub_time_step = False)
#eq.set_initial_condition("sse.nc", "ux.nc", "uy.nc")
eq.set_viscosity("constant", constant_viscosity = 1e-6)
eq.set_dissipation("quadratic", quadratic_coefficient = 2.5e-3)
......@@ -19,16 +19,20 @@ eq.set_boundary_coast("Wall", slip=False)
#eq.set_boundary_river("physical_tag", "discharge.nc")
eq.set_temporal_scheme("explicit")
#eq = slim.ShallowWaterTracer2d(domain, linear_equation=False, temporal_scheme="explicit", offline=False)
#eq.set_initial_condition("concentration.nc")
#eq.set_diffusivity("okubo OU constant", coeff=0.01)
#eq.set_hydro_solution(
#eq.set_boundary_coast("physical_tag")
#eq.set_boundary_sea("physical_tag", "concentration.nc", "concentration_gradient.nc")
eq2 = slim.ShallowWaterTracer2d(domain, linear_equation=False, name="Temperature", offline=False, wetting_drying = None)
#eq2.set_initial_condition("concentration.nc")
eq2.set_diffusivity("okubo", okubo_coefficient=0.03)
eq2.set_hydro_solution(equation=eq)
eq2.set_boundary_coast("Wall")
#eq2.set_boundary_sea("physical_tag", "concentration.nc", "concentration_gradient.nc")
eq2.set_temporal_scheme("explicit")
loop=slim.Loop(initial_time="10 06 1991 22:00:00", final_time= "10 06 1991 23:00:00", maximum_time_step = 900, export_time = 3000)
loop.add_equation(eq)
#loop.add_equation(eq2)
loop.set_output_directory("./output")
loop.run()
#while not loop.end():
......
This diff is collapsed.
......@@ -126,16 +126,21 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
if (_linear){
if (_filterLinear)
val.setAll(0);
else
else{
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(bath(i,0)*_g);
}
}else
if (!_filterLinear)
for (int i=0; i<val.size1(); i++)
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));
else
}
}
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);
......@@ -218,15 +223,16 @@ class dgConservationLawShallowWater2d::source: public function {
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;
double _R, _g, _density;
public :
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 *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R, double g, double density,
const functor *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *absLayerCoefF, const functor *absLayerExtForcF) :function (3){
_linear = linear;
_R = R;
_g = g;
_density = density;
_isDiffusive = nuF;
_useMovingBathWD = movingBathFactorF;
......@@ -277,8 +283,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) + _windStress(i,0);
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + _source(i,1) + _windStress(i,1);
val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + _source(i,0) + _windStress(i,0)/((eta+_bathymetry(i,0))*_density);
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + _source(i,1) + _windStress(i,1)/((eta+_bathymetry(i,0))*_density);
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;
......@@ -584,6 +590,7 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_coordinatesF = NULL;
_R = -1.0;
_g = 9.80616;
_density = 1025;
_nu = NULL;
_useMovingBathWD = false;
_movingBathFactor = new functionConstant(1.);
......@@ -630,11 +637,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, _windStress, _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, _density, _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, _windStress, _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, _density, _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);
......@@ -683,7 +690,7 @@ void dgConservationLawShallowWater2d::setAbsLayer(const std::string tag, const f
Msg::Fatal("Dg/ShallowWater/AbsLayer: Bad 'absCoef' and/or 'extFunc'.");
if (_volumeTerm0[tag]) delete _volumeTerm0[tag];
_volumeTerm0[tag] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient,
_nu, false || _linear, _R, _g, _movingBathFactor, _movingBathFactorGradient, _xyz, absCoef, extFunc);
_nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, absCoef, extFunc);
}
......
......@@ -34,6 +34,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
bool _laxFriedrichs;
double _R;
double _g;
double _density;
const functor *_bathymetry, *_bathymetryGradient;
const functor *_originalBathymetry, *_originalBathymetryGradient;
......@@ -64,6 +65,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
/**if this flag is true, a spherical version of the equations are solved*/
inline void setIsSpherical(double R) {_R = R;}
void setGravity(double g) {_g = g;}
void setDensity(double density) {_density = density;}
inline double getRadius() {return _R;}
void setCoordinatesFunction(functor *xyz);
......
......@@ -1166,18 +1166,20 @@ void okubo2d::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
mergeConcentrationGradient::mergeConcentrationGradient(const functor *f1, const functor *f2) {
mergeConcentrationGradient::mergeConcentrationGradient(const functor *f1, const functor *f2, , const functor *f3) {
_f1 = f1;
_f2 = f2;
_f3 = f3;
}
void mergeConcentrationGradient::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &f1 = cm.get(*_f1);
const fullMatrix<double> &f2 = cm.get(*_f2);
const fullMatrix<double> &f3 = cm.get(*_f3);
val.resize(cm.nEvaluationPoint(),3,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
val(i,0) = f1(i,0);
val(i,1) = f2(i,0);
val(i,2) = 0;
val(i,2) = f3(i,0);
}
}
......@@ -331,10 +331,10 @@ class okubo2d : public functor {
};
class mergeConcentrationGradient: public functor {
const functor *_f1, *_f2;
const functor *_f1, *_f2, *_f3;
public:
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
mergeConcentrationGradient(const functor *f1, const functor *f2);
mergeConcentrationGradient(const functor *f1, const functor *f2, const functor *f3);
};
#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