Commit c7210b0a authored by Jonathan Lambrechts's avatar Jonathan Lambrechts Committed by Valentin Vallaeys
Browse files

list of tracer equations

parent 7dd562c2
......@@ -346,11 +346,6 @@ def slim3d_setup(loop):
for openBnd in eq._boundary_open:
e.SEq.setBoundaryNeumann(openBnd.tag, openBnd.s_flux_f)
if slimSolver.getSolveSImplicitVerticalDiffusion() or slimSolver.getSolveSImplicitVerticalAdvection() :
e.SDiffEq.setBoundary0Flux(topTags)
e.SDiffEq.setBoundary0Flux(bottomTags)
e.SDiffEq.setBoundary0Flux(eq._boundary_coast)
if slimSolver.getSolveT() :
e.TEq.setBoundary0Flux(topTags)
e.TEq.setBoundary0Flux(bottomTags)
......@@ -359,25 +354,6 @@ def slim3d_setup(loop):
for openBnd in eq._boundary_open:
e.TEq.setBoundaryNeumann(openBnd.tag, openBnd.t_flux_f)
if slimSolver.getSolveTImplicitVerticalDiffusion() :
if eq._temperature_boundary_relaxation:
if eq._temperature_boundary_relaxation[0] == 'netcdf':
eq._temp_bnd_relax = slim_private._load_function(eq._temperature_boundary_relaxation[1], slimSolver.groups3d)
elif eq._temperature_boundary_relaxation[0] == 'vertical_gradient':
def zFunc(cmap, val, xyz) :
val[:] = eq._temperature_boundary_relaxation[2] + xyz[:,2] * eq._temperature_boundary_relaxation[3]
eq._temp_bnd_relax = dgpy.functionNumpy(1, zFunc, [d.xyzOrigDof.getFunction()])
else:
dgpy.Msg.Fatal('Unknown mode for initial temperature:' + eq._temperature_boundary_relaxation[0])
TDiffBndSurf = e.TDiffEq.newBoundarySurfaceRelaxation(eq._temp_bnd_relax, eq._temperature_boundary_relaxation[4])
e.TDiffEq.addBoundaryCondition(topTags, TDiffBndSurf)
TDiffBndBottom = e.TDiffEq.newBoundaryBottomRelaxation(eq._temp_bnd_relax, eq._temperature_boundary_relaxation[4], eq._temperature_boundary_relaxation[5])
e.TDiffEq.addBoundaryCondition(bottomTags, TDiffBndBottom)
else:
e.TDiffEq.setBoundary0Flux(topTags)
e.TDiffEq.setBoundary0Flux(bottomTags)
e.TDiffEq.setBoundary0Flux(eq._boundary_coast)
if slimSolver.getAdvectTurbulentEnergy() :
for turbEq in [e.tkeAdvEq, e.epsAdvEq] :
turbEq.setBoundary0Flux(topTags)
......
......@@ -2316,12 +2316,10 @@ bool dgSlopeLimiterOriginal::boundaryValues(const dgGroupOfFaces * faces, dgDofC
}
}
dgSlopeLimiterOriginal::dgSlopeLimiterOriginal(dgConservationLaw* claw, dgGroupCollection* groups): dgLimiter(claw)
dgSlopeLimiterOriginal::dgSlopeLimiterOriginal(dgConservationLaw* claw, const dgGroupCollection* groups): dgLimiter(claw)
{
_groups = groups;
_nodalVolume = new dgDofContainer(*_groups, 1);
dgNodalVolume::compute(_nodalVolume);
}
......
......@@ -60,12 +60,12 @@ public :
/** original limiter weighted by the volume */
class dgSlopeLimiterOriginal : public dgLimiter{
dgGroupCollection *_groups;
const dgGroupCollection *_groups;
dgDofContainer* _nodalVolume;
public :
virtual bool boundaryValues(const dgGroupOfFaces * faces, dgDofContainer * solution,
fullMatrix<double> & valueOnFaces, const functor * secondary=NULL);
dgSlopeLimiterOriginal (dgConservationLaw *claw, dgGroupCollection *groups);
dgSlopeLimiterOriginal (dgConservationLaw *claw, const dgGroupCollection *groups);
virtual ~dgSlopeLimiterOriginal();
virtual void apply (dgDofContainer *solution);
};
......
#ifndef _DG_CONSERVATION_LAW_SW_3D_TRACER_H
#define _DG_CONSERVATION_LAW_SW_3D_TRACER_H
#include "dgConservationLawFunction.h"
class function;
/**Advection and diffusion of a scalar, the advection and diffusion are provided by functions.
* \f[\frac{dc}{dt} + v. \nabla c - nu \Delta c = 0\f]
*/
class dgConservationLawSW3dTracer : public dgConservationLawFunction {
const functor *_uv, *_w, *_wOld, *_wM, *_eta, *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_ipTerm, *_slope, *_Kedd, *_isopycnalDiffusion, *_interfaceIso, *_ipTermIso, *_sourceFunc, *_uGM;
const functor *_fzero, *_fzeroVec;
double _laxFriedrichsFactor;
int _dt,_T;
const functor *_uv, *_w, *_eta, *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_sourceFunc, *_explicitC, *_explicitCGrad;
int _dt, _T;
const functor *_ref;
void source(functorCache &cm, fullMatrix<double> &val) const ;
void advection(functorCache &cm, fullMatrix<double> &val) const ;
void interfaceTerm(functorCache &cm, fullMatrix<double> &val) const ;
void maxConvectiveSpeed(functorCache &cm, fullMatrix<double> &val) const ;
void diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const ;
//class boundaryRelaxation;
void isopycnalDiffusion(functorCache &cm, fullMatrix<double> &val) const ;
void interfaceIso(functorCache &cm, fullMatrix<double> &val) const ;
void diffusivityTensor(functorCache &cm, fullMatrix<double> &val) const ;
//class boundarySurfaceT;
//class boundarySurfaceS;
protected:
dgGroupCollection *_groups;
void addIPTerm(functorCache &cm, fullMatrix<double> &val)const;
public:
void setup();
/**A new advection-diffusion law.
* The horizontal advection velocity is given by 'uv' vector function while the vertical advection velocity is w.
* Horizontal and vertical diffusivity are given by the functions kappaH and kappaV. */
dgConservationLawSW3dTracer(const functor *uv,const functor *w);
/**A new advection-diffusion law */
dgConservationLawSW3dTracer();
/** set factor for scaling Lax-Friedrichs penalty term [default LFF = 0]. */
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;}
inline void setSource(const functor *f) {_sourceFunc = f;};
inline void setCallRef(int dt, int T, const functor *ref){_dt=dt;_T=T; _ref=ref;};
inline void setUV(const functor *f) {_uv = f;};
inline void setW(const functor *f) {_w = f;};
inline void setWOld(const functor *f) {_wOld = f;};
inline void setKappaH(const functor *kap) {_kappaH = kap;};
inline void setKappaV(const functor *kap) {_kappaV = kap;};
inline void setWMesh(const functor *wM) {_wM = wM;};
inline void setEta(const functor *eta) {_eta = eta;};
inline void setBathymetry(const functor *bath) {_bathymetry = bath;};
inline void setIsoDiff(const functor *S, const functor *Kedd) {_slope = S; _Kedd = Kedd;};
inline void setGMVel(const functor *uGM) {_uGM = uGM; };
void setFullExplicit();
void setImplicitVertical(const functor *explicitC, const functor *explicitGradC);
const functor *getUV()const {return _uv;}
const functor *getEta()const {return _eta;}
~dgConservationLawSW3dTracer();
......@@ -61,40 +47,4 @@ class dgConservationLawSW3dTracer : public dgConservationLawFunction {
//dgBoundaryCondition *newRelaxationBoundary(const functor* outValue, double TRelax);
};
/**Diffusion of a scalar, the diffusivity provided by a function.
* \f[\frac{dc}{dt} - nu \Delta c = 0\f]
* This is lightweight version contains only vertical diffusion
* to speed up the implict solver!
*/
class dgConservationLawSW3dTracerVDiff : public dgConservationLawFunction {
const functor *_kappaV, *_advectiveDiffusiveFlux, *_diffusiveFlux, *_sourceTerm, *_interfaceTerm, *_ipTerm;
const functor *_w, *_wM, *_dwMdz;
const functor *_fzero;
void diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const ;
void advectiveDiffusiveFlux(functorCache &cm, fullMatrix<double> &val) const ;
void diffusivityTensor(functorCache &cm, fullMatrix<double> &val) const ;
void source(functorCache &cm, fullMatrix<double> &val) const ;
void interfaceTerm(functorCache &cm, fullMatrix<double> &val) const ;
double _laxFriedrichsFactor;
bool _useConservativeALE;
void boundarySurfaceRelaxation(functorCache &cm, fullMatrix<double> &val, const functor *surfaceVal, const double tau) const;
void boundaryBottomRelaxation(functorCache &cm, fullMatrix<double> &val, const functor *bottomVal, const double tau, const double belowIsBottom) const;
public:
void setup();
/**A new advection-diffusion law.
* Vertical diffusivity is given by the functions kappaV. */
dgConservationLawSW3dTracerVDiff();
virtual bool isLinear() const {return true;}
inline void setKappaV(const functor *kap) {_kappaV = kap;};
inline void setW(const functor *f) {_w = f;};
inline void setWMesh(const functor *wM) {_wM = wM;};
inline void setWMeshDz(const functor *dwMdz) {_dwMdz = dwMdz;};
~dgConservationLawSW3dTracerVDiff();
dgFaceTermAssembler *getFaceTerm(const dgGroupOfFaces &group) const;
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;}
inline void setUseConservativeALE(bool b) {_useConservativeALE = b;};
dgBoundaryCondition *newBoundarySurfaceRelaxation(const functor* surfaceVal, const double tau);
dgBoundaryCondition *newBoundaryBottomRelaxation(const functor* surfaceVal, const double tau, const double belowIsBottom);
};
#endif
......@@ -703,9 +703,6 @@ slim3dSolverEquations::~slim3dSolverEquations()
delete TEq;
delete tkeAdvEq;
delete epsAdvEq;
delete SDiffEq;
delete TDiffEq;
delete sedDiffEq;
delete depthIntegratorRho;
delete sw2DEq;
delete uvHorDivEq;
......@@ -735,9 +732,6 @@ slim3dSolverEquations::slim3dSolverEquations(slim3dSolver* solver) : _solver(sol
TEq = new dgConservationLawSW3dTracer();
tkeAdvEq = new dgConservationLawSW3dTracer();
epsAdvEq = new dgConservationLawSW3dTracer();
SDiffEq = new dgConservationLawSW3dTracerVDiff();
TDiffEq = new dgConservationLawSW3dTracerVDiff();
sedDiffEq = new dgConservationLawSW3dTracerVDiff();
//depthIntegratorRho = new dgSW3dDepthIntegrate(1,INTEGRATE_FROM_TOP);
sw2DEq = new dgConservationLawShallowWater2d();
//uvHorDivEq = new dgUVHorDivergence();
......@@ -857,7 +851,6 @@ void slim3dSolverEquations::initialize()
}
}
if ( _solver->_solveSDiff ) {
SDiffEq->setKappaV(funcs->kappavS);
//SDiffEq->setBoundary0Flux(_solver->physicalTags2d);
}
......@@ -871,7 +864,6 @@ void slim3dSolverEquations::initialize()
}
}
if ( _solver->_solveTDiff ) {
TDiffEq->setKappaV(funcs->kappavT);
//TDiffEq->setBoundary0Flux(_solver->physicalTags2d);
}
......
......@@ -127,12 +127,6 @@ public:
dgConservationLawSW3dTracer *tkeAdvEq;
/** 3d transport equation for turbulence dissipation rate */
dgConservationLawSW3dTracer *epsAdvEq;
/** Vertical diffusion for salinity */
dgConservationLawSW3dTracerVDiff *SDiffEq;
/** Vertical diffusion for temperature */
dgConservationLawSW3dTracerVDiff *TDiffEq;
/** Vertical diffusion for sediment */
dgConservationLawSW3dTracerVDiff *sedDiffEq;
/** Vertical integrator for integrating r (baroclinic head) */
dgSW3dDepthIntegrate *depthIntegratorRho;
/** new way to compute rGRad */
......
......@@ -48,7 +48,7 @@ slim3dSolver::slim3dSolver(const std::string meshFile3d, std::vector<std::string
prismConverter = NULL;
movingMesh = NULL;
//zDiff1dEq = NULL;
uvLimiter = SLimiter = TLimiter = tkeLimiter = epsLimiter = sedLimiter= NULL;
uvLimiter = tkeLimiter = epsLimiter = sedLimiter= NULL;
jumpDiffTracerS = jumpDiffTracerT = jumpDiffTracerVertS = jumpDiffTracerVertT = NULL;
dgCG3d = dgCG2d = NULL;
copyField2d = NULL;
......@@ -122,8 +122,6 @@ slim3dSolver::~slim3dSolver()
if ( _verticalAdaptMesh ) delete _verticalAdaptMesh;
// if ( zDiff1dEq ) delete zDiff1dEq;
if ( uvLimiter ) delete uvLimiter;
if ( SLimiter ) delete SLimiter;
if ( TLimiter ) delete TLimiter;
if ( tkeLimiter ) delete tkeLimiter;
if ( epsLimiter ) delete epsLimiter;
if ( sedLimiter ) delete sedLimiter;
......@@ -169,26 +167,9 @@ slim3dSolverEquations* slim3dSolver::createEquations()
uvLimiter = new dgSlopeLimiterUV (_equations->horMomEq,groups3d);
//uvLimiter = new dgSlopeLimiterOriginal (_equations->horMomEq,groups3d);
}
if ( _solveSDiff && !_solveS && _useSLimiter){
//SLimiter = new dgSlopeLimiter4 (_equations->SDiffEq,groups3d, _dofs->SDof);
SLimiter = new dgSlopeLimiterOriginal (_equations->SDiffEq,groups3d);
}
else if (_useSLimiter && _solveS){
//SLimiter = new dgSlopeLimiter4 (_equations->SEq,groups3d, _dofs->SDof);
SLimiter = new dgSlopeLimiterOriginal (_equations->SEq,groups3d);
//SLimiter = new dgSlopeLimiter2 (_equations->SEq,groups3d, nodalVolume3d->get(), true);
}
if ( _solveTDiff && !_solveT && _useTLimiter){
//TLimiter = new dgSlopeLimiter4 (_equations->TDiffEq,groups3d, _dofs->TDof);
TLimiter = new dgSlopeLimiterOriginal (_equations->TDiffEq,groups3d);
}
else if (_useTLimiter && _solveT){
//TLimiter = new dgSlopeLimiter4 (_equations->TEq,groups3d, _dofs->TDof);
TLimiter = new dgSlopeLimiterOriginal (_equations->TEq,groups3d);
}
if ( _solveSedDiff && !_solveSed && _useSedLimiter){
//sedLimiter = new dgSlopeLimiter4 (_equations->sedDiffEq,groups3d, _dofs->sedDof);
sedLimiter = new dgSlopeLimiterOriginal (_equations->sedDiffEq,groups3d);
sedLimiter = new dgSlopeLimiterOriginal (_equations->sedEq[0],groups3d);
}
else if (_useSedLimiter && _solveSed){
//sedLimiter = new dgSlopeLimiter4 (_equations->sedEq,groups3d, _dofs->sedDof);
......
......@@ -72,7 +72,7 @@ public:
dgSW3dMovingMesh *movingMesh;
//slim3dZDiff1dEq *zDiff1dEq;
dgSW3dVerticalBottomRemover *verticalBottomRemover;
dgLimiter *uvLimiter, *SLimiter, *TLimiter, *tkeLimiter, *epsLimiter, *sedLimiter, *sedLimiterGround2d;
dgLimiter *uvLimiter, *tkeLimiter, *epsLimiter, *sedLimiter, *sedLimiterGround2d;
dgJumpDiffusionTracer *jumpDiffTracerS, *jumpDiffTracerT, *jumpDiffTracerVertS, *jumpDiffTracerVertT;
dgCGMeanFilter *dgCG3d, *dgCG2d;
// dgCGStructure *dgCG3d, *dgCG2d;
......
This diff is collapsed.
......@@ -8,8 +8,10 @@
#include "dgDIRK.h"
/** This is the base class for all time integration instants */
class slim3dTracerSolver;
class slim3dTimeIntegrator {
protected:
std::vector<slim3dTracerSolver*> _tracerSolvers;
void _initializeFields();
clock_t _startClock;
......@@ -22,15 +24,9 @@ protected:
int _iter;
bool _exportEveryIter;
bool _exploded;
linearSystem<fullMatrix<double> > *sw2DSys, *uSys, *SDiffSys, *TDiffSys, *sedDiffSys;
linearSystem <double> *streamFuncSys;
dgDofManager *dofManager2D, *dofManagerU, *dofManagerSDiff, *dofManagerTDiff, *dofManagerSedDiff;
dgDofManager *dofManagerStreamFunc;
linearSystem<fullMatrix<double> > *sw2DSys, *uSys;
dgDofManager *dofManager2D, *dofManagerU;
dgDIRK *rk;
//for GMVel
dgDofContainer *rhoGradDDof3d, *streamFuncDDof3d;
// TODO remove dgDIRK, create newton directly
dgSteady *GMStreamFuncSolver;
// compute mass conservation and tracer deviation
bool _checkMassConservation, _checkDeviation;
class massConservationChecker;
......@@ -41,10 +37,8 @@ protected:
/* for adaptivity. vJump advection */
dgDofContainer *vJumpDof_K;
public:
dgNewton *momUNewton, *momVNewton, *SDiffNewton, *TDiffNewton;
std::vector<dgNewton *> sedDiffNewton;
dgDofContainer *sw2DSol, *SDiffNewtonSol, *TDiffNewtonSol, *momUNewtonSol, *momVNewtonSol;
std::vector<dgDofContainer *> sedDiffNewtonSol;
dgNewton *momUNewton, *momVNewton;
dgDofContainer *sw2DSol, *momUNewtonSol, *momVNewtonSol;
/** Constructor */
slim3dTimeIntegrator( slim3dSolver* solver );
virtual ~slim3dTimeIntegrator();
......@@ -96,37 +90,13 @@ public:
void terminate(int exitFlag=0);
};
/** Simple predictor-corrector time integration */
class slim3dTimeIntegratorPC : public slim3dTimeIntegrator {
bool _started;
std::vector<dgIdxExporter*> _fullExporter;
public:
dgDofContainer *uvDof_K;
dgDofContainer *SDof_K;
dgDofContainer *TDof_K;
dgDofContainer *sedDof_K;
const dgDofContainer *uKim, *vKim, *SKim, *TKim;
std::vector<const dgDofContainer *>sedKim;
dgDofContainer *tkeDof_K, *tkeDof_K1, *tkeDof_K2, *tkeDof_K3;
dgDofContainer *epsDof_K, *epsDof_K1, *epsDof_K2, *epsDof_K3;
// for predictor corrector
dgDofContainer *etaDof2d_tav_new;
dgDofContainer *uvAvDof2d_tav_new, *uvAvDof2d_tav_newHalf;
dgDofContainer *uvDof_new, *uvDof_old;
dgDofContainer *SDof_new, *SDof_old;
dgDofContainer *TDof_new, *TDof_old;
const dgDofContainer *uKim, *vKim;
dgDofContainer *uvCorr;
std::vector<dgDofContainer *> sedDof_new, sedDof_old;
std::vector<dgDofContainer *> sedGroundDof_new, sedGroundDof_old;
/** Adams-Bashfort coefficients */
double betaAB3;
/** Leapfrog-AM3 coefficient (Shchepetkin and McWilliams 2005) **/
double gammaLFAM3;
/** implicity of semi-implicit vertical diffusion */
double CrNiTheta;
slim3dTimeIntegratorPC( slim3dSolver* solver );
~slim3dTimeIntegratorPC();
virtual void iterate();
......
......@@ -53,7 +53,7 @@ import slim3d
domain = slim3d.Domain('box.msh')
equations = slim3d.Slim3d_equations(domain, temperature=True, salinity=True)
equations.set_implicit_vertical_diffusion(False)
equations.set_implicit_vertical_advection(False)
equations.set_implicit_vertical_advection(True)
equations.set_lax_friedrichs_factor(1.)
equations.set_bottom_friction(False)
......
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