Commit c0b7ceca authored by Philippe Delandmeter's avatar Philippe Delandmeter Committed by Valentin Vallaeys
Browse files

slim3d: relaxation at the surface for TEq

parent 4ac1ea96
Pipeline #1757 passed with stage
in 32 minutes and 53 seconds
......@@ -12,7 +12,7 @@ shutil.copyfile(data_dir+'mesh3d_show.msh', output_dir+'mesh3d_show.msh')
domain = slim3d.Domain(data_dir+'mesh3d.msh', data_dir+'periodicMesh.txt')
equations = slim3d.Slim3d_equations(domain, temperature=True, salinity=True)
equations.set_implicit_vertical_diffusion(True)
equations.set_implicit_vertical(True)
equations.set_vertical_viscosity('gotm')
equations.set_vertical_diffusivity('gotm')
equations.set_boundary_coast(['coastL', 'coastR'])
......
......@@ -176,6 +176,15 @@ def slim3d_setup(loop):
f.TInitFunc = dgpy.functorNumpy(lambda cm : eq._initial_temperature[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._initial_temperature[3])
else:
dgpy.Msg.Fatal('Unknown mode for initial temperature:' + eq._initial_temperature[0])
if eq._temperature_boundary_relaxation:
if eq._temperature_boundary_relaxation[0] == 'netcdf':
eq._TRelaxFunc = slim_private._load_function(eq._temperature_boundary_relaxation[1], slimSolver.groups3d)
f.TRelaxFunc = eq._TRelaxFunc
elif eq._temperature_boundary_relaxation[0] == 'vertical_gradient':
f.TRelaxFunc = dgpy.functorNumpy(lambda cm : eq._temperature_boundary_relaxation[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._temperature_boundary_relaxation[3])
else:
dgpy.Msg.Fatal('Unknown mode for temperature relaxation:' + eq._temperature_boundary_relaxation[0])
elif (not eq._linear_density) or (eq._linear_density[0] == 'temperature'):
dgpy.Msg.Fatal('Initial temperature must be set for rhoFunc')
......@@ -338,6 +347,8 @@ def slim3d_setup(loop):
e.TEq.setBoundary0Flux(bottomTags)
e.TEq.setBoundary0Flux(eq._boundary_coast)
e.TEq.setBoundary0Flux('vertical_bottom')
if f.TRelaxFunc:
e.TEq.setRelaxation(f.TRelaxFunc, eq._temperature_boundary_relaxation[4])
for openBnd in eq._boundary_open:
e.TEq.setBoundaryNeumann(openBnd.tag, openBnd.t_flux_f)
......
......@@ -294,9 +294,6 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double _upwindFactor = _upwindFactorRiemann;
const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL;
for(size_t i=0; i< nQP; i++) {
//if (cache.interfaceGroup()->physicalTag() == "coastL"){
// printf("sw2dEq: etaR %g, uR %g\n", solR(i,0), solR(i,1));
//}
double pa = 0.5*(atmPressL(i,0) + atmPressR(i,0));
double rhoSurf = 0.5*(rhoSurfL(i,0) + rhoSurfR(i,0));
double Fun, Fut, Feta;
......
......@@ -39,7 +39,6 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
const fullMatrix<double> &bathL = cache.get(*_bathymetry, 0);
const fullMatrix<double> &bathR = cache.get(*_bathymetry, 1);
double g = slim3dParameters::g;
// determine which element is below
for(int i=0; i< val.size1(); i++) {
double nx = normals(i,0);
double ny = normals(i,1);
......@@ -62,7 +61,6 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
dgConservationLawSW3dContinuity::dgConservationLawSW3dContinuity( int propagateFromTop ) :
dgConservationLawFunction(1), _propagateFromTop(propagateFromTop){
_laxFriedrichsFactor = 0.0;
_uv = NULL;
_eta = NULL;
}
......
......@@ -14,7 +14,6 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
// they are thus constants in z
const functor *_uv, *_eta, *_bathymetry;
bool _propagateFromTop;
double _laxFriedrichsFactor;
void riemann(functorCache &cache, fullMatrix<double> &val) const;
void advection(functorCache &cache, fullMatrix<double> &val) const;
void boundaryWall(functorCache &cache, fullMatrix<double> &val) const;
......@@ -29,9 +28,7 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
const functor *getEta()const {return _eta;}
const functor *getUV()const {return _uv;}
inline void setBathymetry(const functor *bathymetry) { _bathymetry = bathymetry;}
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
void finalize(const dgGroupOfElements &group, dgDofContainer &residual) const;
//inline void setW(const functor *f) { _w = f; };
};
/**Pressure equation for 3D shallow water equation.
......
......@@ -310,30 +310,14 @@ dgBoundaryCondition *dgConservationLawSW3dMomentum::newBoundaryBottom(const func
void dgConservationLawSW3dMomentum::boundarySurface(functorCache &cm, fullMatrix<double> &val, const functor *windF) const {
const fullMatrix<double> &normals = cm.get(*function::getNormals());
// const fullMatrix<double> &UV = cm.get(*function::getSolution());
const fullMatrix<double> &UVgrad = cm.get(*function::getSolutionGradient());
const fullMatrix<double> &wind = cm.get(*windF);
const fullMatrix<double> &nuH = cm.get(*_nuH);
const fullMatrix<double> &nuV = cm.get(*_nuV);
double rho0 = slim3dParameters::rho0;
val.resize(cm.nEvaluationPoint(), 2, false);
// bnd condition for linear momentum
for(int i=0;i<cm.nEvaluationPoint(); i++) {
double nx = normals(i,0);
double ny = normals(i,1);
double nz = normals(i,2);
val(i,0) = 0;
val(i,1) = 0;
if (_nuV) {
val(i,0) += wind(i,0)/rho0*nz;
val(i,1) += wind(i,1)/rho0*nz;
double dudn = UVgrad(i,0)*nx + UVgrad(i,1)*ny;// + UVgrad(i,2)*nz;
double dvdn = UVgrad(i,3)*nx + UVgrad(i,4)*ny;// + UVgrad(i,5)*nz;
double nGradUVn = dudn*nx + dvdn*ny;
double dudzn = UVgrad(i,2)*nx + UVgrad(i,5)*ny;
val(i,0) += -nuH(i,0)*nGradUVn*nx - nuV(i,0)*dudzn*nx*nz;// - jumpPenaltyX;
val(i,1) += -nuH(i,0)*nGradUVn*ny - nuV(i,0)*dudzn*ny*nz;// - jumpPenaltyY;
}
val(i,0) = wind(i,0)/rho0*nz;
val(i,1) = wind(i,1)/rho0*nz;
}
}
......
......@@ -29,8 +29,12 @@ void dgConservationLawSW3dTracer::maxConvectiveSpeed(functorCache &cm, fullMatri
void dgConservationLawSW3dTracer::source(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &sol = cm.get(cm.solution());
const fullMatrix<double> &cE = cm.get(_explicitSol ? *_explicitSol : cm.solution());
const fullMatrix<double> *sourceFunc = _sourceFunc ? &cm.get(*_sourceFunc) : NULL;
const fullMatrix<double> *ref = _ref ? &cm.get(*_ref) : NULL;
const fullMatrix<double> *relaxRef = _relaxRef ? &cm.get(*_relaxRef) : NULL;
const fullMatrix<double> &xyz = cm.get(*function::getCoordinates());
double zMin = 12;
int nQP = cm.nEvaluationPoint();
val.resize(nQP, 1, true); //already set to zero in resize
for(int i=0; i< nQP; i++) {
......@@ -38,6 +42,10 @@ void dgConservationLawSW3dTracer::source(functorCache &cm, fullMatrix<double> &v
val(i,0) += (*sourceFunc)(i,0);
if(ref)
val(i,0) += ((*ref)(i,0)-sol(i,0))*_dt/_T;
if(relaxRef){
if (xyz(i,2) > -zMin)
val(i,0) += ((*relaxRef)(i,0)-cE(i,0)) * (zMin+xyz(i,2))/zMin / _relaxTau;
}
}
}
......@@ -209,7 +217,7 @@ void dgConservationLawSW3dTracer::setup () {
Msg::Fatal("dgConservationLawSW3dTracer: Bath function must be defined");
_maximumConvectiveSpeed[""] = newFunctorMember(*this, &dgConservationLawSW3dTracer::maxConvectiveSpeed);
_diffusiveFlux = haveDiff ? newFunctorMember(*this, &dgConservationLawSW3dTracer::diffusiveFlux) : NULL;
_volumeTerm0[""] = (_sourceFunc || _ref) ? newFunctorMember(*this, &dgConservationLawSW3dTracer::source) : NULL;
_volumeTerm0[""] = (_sourceFunc || _ref || _relaxRef) ? newFunctorMember(*this, &dgConservationLawSW3dTracer::source) : NULL;
_volumeTerm1[""] = newFunctorMember(*this, &dgConservationLawSW3dTracer::advection);
_interfaceTerm0[""] = newFunctorMember(*this, &dgConservationLawSW3dTracer::interfaceTerm);
}
......@@ -267,6 +275,8 @@ dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawSW
_bathymetry = NULL;
_diffusiveFlux = NULL;
_ref=NULL;
_relaxRef = NULL;
_dt=0;
_T=1;
_relaxTau=0;
}
......@@ -24,6 +24,8 @@ class dgConservationLawSW3dTracer : public dgConservationLawSW3dIMEX {
const functor *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_sourceFunc;
int _dt, _T;
const functor *_ref;
double _relaxTau;
const functor *_relaxRef;
void source(functorCache &cm, fullMatrix<double> &val) const ;
void advection(functorCache &cm, fullMatrix<double> &val) const ;
void interfaceTerm(functorCache &cm, fullMatrix<double> &val) const ;
......@@ -38,6 +40,7 @@ class dgConservationLawSW3dTracer : public dgConservationLawSW3dIMEX {
/**A new advection-diffusion law */
dgConservationLawSW3dTracer();
inline void setSource(const functor *f) {_sourceFunc = f;};
inline void setRelaxation(const functor *f, double tau) {_relaxRef = f; _relaxTau = tau;};
inline void setCallRef(int dt, int T, const functor *ref){_dt=dt;_T=T; _ref=ref;};
inline void setKappaH(const functor *kap) {_kappaH = kap;};
inline void setKappaV(const functor *kap) {_kappaV = kap;};
......
......@@ -322,6 +322,7 @@ slim3dSolverFunctions::slim3dSolverFunctions(slim3dSolver *solver) : _solver(sol
bathFunc2d = bathCGMinFunc2d = NULL;
coriolisFunc = NULL;
TInitFunc = SInitFunc = TFunc = SFunc = sedInitFunc = sedGroundInitFunc = NULL;
TRelaxFunc = NULL;
uvInitFunc = etaInitFunc = wInitFunc = NULL;
newZFunc = wMeshSurfFunc = wMeshFunc = wMeshDzFunc = zBotDistFunc = NULL;
nuvFunc = nuh = kappavFunc = kappavS = kappavT = kappavJumpS = kappavJumpT = kappahTotal = kappahUser = kappahJumpS = kappahJumpT = kappahS = kappahT = NULL;
......@@ -673,6 +674,7 @@ slim3dSolverFunctions::~slim3dSolverFunctions()
if (kappahJumpT) delete kappahJumpT;
if (kappahS) delete kappahS;
if (kappahT) delete kappahT;
if (TRelaxFunc) delete TRelaxFunc;
}
......
......@@ -37,6 +37,7 @@ public:
functionPrecomputedExtrusion *bathFunc2d, *bathCGMinFunc2d;
const functor *coriolisFunc;
const functor *TInitFunc, *SInitFunc, *TInitGradientFunc, *SInitGradientFunc, *sedInitFunc, *sedGroundInitFunc;
const functor *TRelaxFunc;
const functor *uvInitFunc, *etaInitFunc, *wInitFunc;
const functor *SFunc, *TFunc, *SCenterGradientFunc, *TCenterGradientFunc;
// moving mesh
......
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