Skip to content
GitLab
Commit e259abb5 authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

first attempt with new mode splitting

parent 85b6b79f
Pipeline #1084 failed with stage
in 31 minutes and 59 seconds
Expand all Hide whitespace changes
Inline Side-by-side
dgpy/scripts/slim3d_setup.py
View file @ e259abb5
......@@ -98,7 +98,7 @@ def slim3d_setup(loop):
area2d.scatter()
eq._areaFunc = dgpy.functionPrecomputedExtrusion(slimSolver.extrusion(), 3, area2d.getFunction())
f.nuh = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof.getFunction(), eq._areaFunc)
f.nuh2d = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof2d.getFunction(), eq._areaFunc)
# f.nuh2d = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof2d.getFunction(), eq._areaFunc)
if eq._horizontal_diffusivity == 'constant':
f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const)
......
modules/slim3d/dgConservationLawSW3dMomentum.h
View file @ e259abb5
......@@ -12,7 +12,7 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
void maxConvectiveSpeed(functorCache &cm, fullMatrix<double> &val) const;
void boundaryWallLinear(functorCache &cm, fullMatrix<double> &val, double uvPenaltyFactor) const ;
void boundarySurface(functorCache &cm, fullMatrix<double> &val, const functor *windF) const;
const functor *_uv2d, *_eta, *_rGrad, *_rhoSurf, *_w, *_wM, *_dwMdz, *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry, *_bathymetryMinCG, *_atmPress, *_rhoSurfGrad;
const functor *_uv2d, *_uv2dGrad, *_eta, *_rGrad, *_rhoSurf, *_w, *_wM, *_dwMdz, *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry, *_bathymetryMinCG, *_rhoSurfGrad, *_correction;
const functor *_fzero, *_fzerov;
const functor *_nuH, *_nuV, *_diffusiveFlux, *_ipTerm, *_maxDiffusivity;
double _laxFriedrichsFactor;
......@@ -39,11 +39,11 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
inline void setRGrad(const functor *f){ _rGrad = f;}
inline void setRhoSurf(const functor *f){ _rhoSurf = f;}
inline void setRhoSurfGrad(const functor *f){ _rhoSurfGrad = f;}
/**set function for atmospheric pressure */
inline void setAtmPress(const functor *f){ _atmPress = f;}
/**set depth averaged uv function [uInt, vInt] */
inline void setUV2d(const functor *uv2d){ _uv2d = uv2d;}
inline const functor* getUV2d(){ return _uv2d;}
/**set depth averaged uvGrad function [dudx, dudy, dudz, dvdx, dvdy, dvdz] */
inline void setUV2dGrad(const functor *uv2dGrad){ _uv2dGrad = uv2dGrad;}
/**set vertical velocity function [w] */
inline void setW(const functor *w){ _w = w;}
inline const functor* getW(){ return _w;}
......@@ -73,6 +73,8 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
inline void setBathymetryMinCG(const functor *bathymetryMinCG) { _bathymetryMinCG = bathymetryMinCG;}
/**set no to use the advection of momentum */
inline void setWithoutAdvection(bool b){ _withoutAdvection = b;}
/** set the correction source term*/
inline void setCorrectionSourceTerm(const functor *f ) {_correction = f;}
dgConservationLawSW3dMomentum();
~dgConservationLawSW3dMomentum();
/** slip wall boundary */
......@@ -168,10 +170,8 @@ class dgConservationLawSW2dU : public dgConservationLawFunction {
void source(functorCache &cm, fullMatrix<double> &val) const ;
void riemann(functorCache &cm, fullMatrix<double> &val) const ;
//class boundaryWall;
void diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const ;
bool _linear, _laxFriedrichs;
const functor *_source, *_coriolisFactor, *_windStress, *_bathymetry, *_fzero, *_fzerov, *_nuH, *_diffusiveFlux, *_eta, *_ipTerm, *_bathymetryGradient, *_rGradInt, *_rhoSurf, *_atmPress, *_rhoSurfGrad;
const functor *_bfCoeff, *_uvBottom, *_uvDevAdvection;
const functor *_source, *_coriolisFactor, *_correction, *_bathymetry, *_fzero, *_fzerov, *_eta, *_rGradInt, *_rhoSurf, *_atmPress, *_rhoSurfGrad;
void boundaryOpenFree(functorCache &cm, fullMatrix<double> &val) const;
void wallExtValue(functorCache &cm, fullMatrix<double> &val) const;
public:
......@@ -182,29 +182,17 @@ class dgConservationLawSW2dU : public dgConservationLawFunction {
/**set the function to evaluate the coriolis factor.
* \f{eqnarray*}\frac{du}{dt} &=& -f v \\ \frac{dv}{dt} &=& f u \f} */
inline void setCoriolisFactor(const functor *f){ _coriolisFactor = f;}
/**set the function to evaluate the wind stress*/
inline void setWindStress(const functor *f){ _windStress = f;}
/**set the function to evaluate the bottom friction drag coefficient.
* The coefficient is computed in the 3D model and is needed for mode splitting.
* \f{eqnarray*}\frac{du}{dt} &=& -c_d u_b||u_b||\\ \frac{dv}{dt} &=& -c_d v_b||u_b|| \f}*/
inline void setBottomFrictionCoefficient(const functor *f){ _bfCoeff = f;}
/**set the function to evaluate the 3d velocity at the bottom.
* The bottom velocity is computed in the 3D model and is needed for mode splitting.
* \f{eqnarray*}\frac{du}{dt} &=& -c_d u_b||u_b||\\ \frac{dv}{dt} &=& -c_d v_b||u_b|| \f}*/
inline void setBottomUV(const functor *f){ _uvBottom = f;}
/**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 *f){ _source = f;}
/**Set vertical viscosity function (scalar) */
inline void setNuH(const functor *f){ _nuH = f;}
/**set the function to evaluate the bathymetry h (H = h+eta) */
inline void setBathymetry(const functor *f) { _bathymetry = f;}
/**set the function to evaluate the gradient of the bathymetry h (dhdx dhdy [dhdz]) */
inline void setBathymetryGradient(const functor *f) { _bathymetryGradient = f;}
/**set the function to evaluate the free surface elevation (eta) */
inline void setEta(const functor *f) { _eta = f;}
/** set if use Lax-Friedrichs or not */
inline void setLaxFriedrichs(bool lf) {_laxFriedrichs = lf; }
/** set the correction source term*/
inline void setCorrectionSourceTerm(const functor *f ) {_correction = f;}
/**set the function to evaluate depth integrated gradient of baroclinic head ( intdrdx intdrdy) */
inline void setRGradInt(const functor *f) { _rGradInt = f;}
......@@ -212,8 +200,6 @@ class dgConservationLawSW2dU : public dgConservationLawFunction {
inline void setRhoSurfGrad(const functor *f) { _rhoSurfGrad = f;}
/**set the function to evaluate the atmospheric pressure */
inline void setAtmPress(const functor *f) { _atmPress = f;}
/**set the function to take into account the advection of ũ (deviation of u) (A_H)*/
inline void setUVDevAdvection(const functor *f){ _uvDevAdvection = f;}
dgConservationLawSW2dU();
~dgConservationLawSW2dU();
/**slip wall boundary */
......
modules/slim3d/slim3dEquations.cpp
View file @ e259abb5
......@@ -7,26 +7,24 @@
slim3dSolverDofs::slim3dSolverDofs(slim3dSolver *solver): _solver(solver) {
// init all dofs to NULL
uDof = vDof = uvDof = uvDofCenter = uvSurfDof2d = uvAvDof2d = NULL;
uDof = vDof = uvDof = uvDofCenter = uvSurfDof2d = uvAvDof2d = uvDevDof = NULL;
uvHorDivDof = NULL;
uvIntDof2d = uvCorrDof3d = uvBotDof2d = NULL;
uvIntDof2d = uvCorrDof2d = uvCorrDof3d = uvBotDof2d = NULL;
windStressDof2dCenter = NULL;
uvTauBDof2d = uvTauBDof2dCenter = NULL;
wDof3d = wSurfDof2d = NULL;
//wMeshDof3d = wMeshSurfDof2d = NULL;
//dwMeshDzDof3d = NULL;
etaDof2d = etaDof2dCG = etaDof3dCenter = NULL;
rDof3d = rIntDof2d = rBotDof2d = NULL;
//rDof3d = rIntDof2d = rBotDof2d = NULL;
rhoDof3d = rGradDof3d = rhoGradDof3d = rhoSurfDof2d = rGradIntDof2d = NULL;
SDof = SDofCenter = TDof = TDofCenter = SBottomDof2d= NULL;
//sedDof = sedBottomDof2d = sedGroundDof2d = sedBottomFluxDof2d = NULL;
sedDof.resize(0); sedBottomDof2d.resize(0); sedGroundDof2d.resize(0); sedBottomFluxDof2d.resize(0);
zDof = deltaZDof3d = deltaZDofSurf2d = zDofCenter = zBotDof2d = NULL;
uvDeviationProdIntDof2d = NULL;
uvDeviationProdDof3d = NULL;
xyzOrigDof = nnDof = ssDof = nuvDof = nuvDofCenter = kappavDof = kappavDofCenter = NULL;
tkeDof = tkeDofCenter = epsDof = epsDofCenter = lDof = NULL;
uvGradDof = uvGradDof2d= NULL;
uvGradDof = NULL;
nbDof3d = rhoGradGMVelDof3d = GMVelDof3d = streamFuncDof3d = NULL;
bathDof2d = bathCGMinDof2d = NULL;
arbitraryMovingMeshDof3d = kappaZEqDof3d = NULL;
......@@ -43,14 +41,14 @@ void slim3dSolverDofs::initialize()
dgGroupCollection *groups2d = _solver->groups2d;
// TODO allocate according to solver options
uvDof = new dgDofContainer(*groups3d, 2);
uvDevDof = new dgDofContainer(*groups3d, 2);
uvHorDivDof = new dgDofContainer(*groups3d, 1);
wDof3d = new dgDofContainer(*groups3d, 1);
//wDof3dCopy = new dgDofContainer(*groups3d, 1);
uvAvDof2d = new dgDofContainer(*groups2d, 2);
uvIntDof2d = new dgDofContainer(*groups2d, 2);
uvCorrDof2d = new dgDofContainer(*groups2d, 2);
uvCorrDof3d = new dgDofContainer(*groups3d, 2);
uvDeviationProdDof3d = new dgDofContainer(*groups3d, 3);
uvDeviationProdIntDof2d = new dgDofContainer(*groups2d, 3); // this is needed
etaDof2d = new dgDofContainer(*groups2d, 1);
etaDof2dCG = new dgDofContainer(*groups2d, 1);
// pressure related
......@@ -64,11 +62,11 @@ void slim3dSolverDofs::initialize()
intDRhodzDof3d = new dgDofContainer(*groups3d, 1);
intDRhodzDofSurf2d = new dgDofContainer(*groups2d, 1);
}
else {
rDof3d = new dgDofContainer(*groups3d, 1);
rIntDof2d = new dgDofContainer(*groups2d, 1);
rBotDof2d = new dgDofContainer(*groups2d, 1);
}
//else {
// rDof3d = new dgDofContainer(*groups3d, 1);
// rIntDof2d = new dgDofContainer(*groups2d, 1);
// rBotDof2d = new dgDofContainer(*groups2d, 1);
//}
// bottom friction related
uvBotDof2d = new dgDofContainer(*groups2d, 2);
uvTauBDof2d = new dgDofContainer(*groups2d, 2);
......@@ -132,7 +130,6 @@ void slim3dSolverDofs::initialize()
}
if ( _solver->getSolveUVGrad()){
uvGradDof = new dgDofContainer(*groups3d, 6);
uvGradDof2d = new dgDofContainer(*groups2d, 6);
}
if ( _solver->getSolveGMVel() ) {
rhoGradGMVelDof3d = new dgDofContainer(*groups3d, 3);
......@@ -160,10 +157,12 @@ void slim3dSolverDofs::initialize()
slim3dSolverDofs::~slim3dSolverDofs()
{
if ( uvDof ) delete uvDof;
if ( uvDevDof ) delete uvDevDof;
if ( uvHorDivDof ) delete uvHorDivDof;
if ( wDof3d ) delete wDof3d;
if ( uvAvDof2d ) delete uvAvDof2d;
if ( uvIntDof2d ) delete uvIntDof2d;
if ( uvCorrDof2d ) delete uvCorrDof2d;
if ( uvCorrDof3d ) delete uvCorrDof3d;
if ( etaDof2d ) delete etaDof2d;
if ( etaDof2dCG ) delete etaDof2dCG;
......@@ -180,9 +179,9 @@ slim3dSolverDofs::~slim3dSolverDofs()
if ( zDof ) delete zDof;
if ( deltaZDof3d ) delete deltaZDof3d;
if ( deltaZDofSurf2d ) delete deltaZDofSurf2d;
if ( rDof3d ) delete rDof3d;
if ( rIntDof2d ) delete rIntDof2d;
if ( rBotDof2d ) delete rBotDof2d;
//if ( rDof3d ) delete rDof3d;
//if ( rIntDof2d ) delete rIntDof2d;
//if ( rBotDof2d ) delete rBotDof2d;
if ( rGradDof3d) delete rGradDof3d;
if ( rhoGradDof3d) delete rhoGradDof3d;
if ( rhoDof3d ) delete rhoDof3d;
......@@ -213,11 +212,8 @@ slim3dSolverDofs::~slim3dSolverDofs()
if ( epsDof ) delete epsDof;
if ( epsDofCenter ) delete epsDofCenter;
if ( lDof ) delete lDof;
if (uvDeviationProdDof3d) delete uvDeviationProdDof3d;
if (uvDeviationProdIntDof2d) delete uvDeviationProdIntDof2d;
if ( nbDof3d ) delete nbDof3d;
if ( uvGradDof ) delete uvGradDof;
if ( uvGradDof2d ) delete uvGradDof2d;
if ( rhoGradGMVelDof3d) delete rhoGradGMVelDof3d;
if ( streamFuncDof3d) delete streamFuncDof3d;
if ( GMVelDof3d) delete GMVelDof3d;
......@@ -252,16 +248,14 @@ slim3dSolverFunctions::slim3dSolverFunctions(slim3dSolver *solver) : _solver(sol
timeFunc = function::getTime();
// functions depending on user input
uvCorrectorFunc = NULL;
pHydroStatFunc = rhoFunc = rhoFuncCenter = eosAlphaFunc = eosBetaFunc = NULL;
SCenterGradientFunc = TCenterGradientFunc = SInitGradientFunc = TInitGradientFunc = NULL;
bathFunc2d = bathCGMinFunc2d = NULL;
coriolisFunc = NULL;
TInitFunc = SInitFunc = TFunc = SFunc = sedInitFunc = sedGroundInitFunc = NULL;
uvInitFunc = uvAvInitFunc = etaInitFunc = wInitFunc = NULL;
uvDevFunc = uvDevProdFunc = uvDevAdvectionFunc = NULL;
newZFunc = wMeshSurfFunc = wMeshFunc = wMeshDzFunc = zBotDistFunc = NULL;
nuvFunc = nuh = nuh2d = kappavFunc = kappavS = kappavT = kappavJumpS = kappavJumpT = kappahTotal = kappahUser = kappahJumpS = kappahJumpT = kappahS = kappahT = NULL;
nuvFunc = nuh = kappavFunc = kappavS = kappavT = kappavJumpS = kappavJumpT = kappahTotal = kappahUser = kappahJumpS = kappahJumpT = kappahS = kappahT = NULL;
maxDiffS = 35.;
maxDiffT = 20.;
jumpPercentS = 0.05;
......@@ -272,13 +266,13 @@ slim3dSolverFunctions::slim3dSolverFunctions(slim3dSolver *solver) : _solver(sol
_nuv0Func = _kapv0Func = NULL;
z0BFunc = z0SFunc = uvTauBFunc2d = windStressFunc = windFunc = NULL;
bottomFrictionFunc = NULL;
bottomFrictionDragCoeff2d = bottomUVDeviation = NULL;
bottomFrictionDragCoeff2d = NULL;
SInitGradientFunc = TInitGradientFunc = zeroFunc3;
slopeFunc = slopeTap = NULL;
_bathDof3d = NULL;
eadySpeedFunc = GMIndepTermFunc = GMVelFunc = NULL;
wBottomFunc = NULL;
eta2dFunc = uvInt2dFunc = uvAv2dFunc = rhoSurf2dFunc = rhoSurfGrad2dFunc = NULL;
eta2dFunc = uvInt2dFunc = u2dFunc = uvAv2dFunc = uvAvGrad2dFunc = rhoSurf2dFunc = rhoSurfGrad2dFunc = NULL;
atmPressFunc = NULL;
uvIntTarget2dFunc = bottomFriction2dPrecompFunc = zBotDist2dPrecompFunc = NULL;
eta2dCGFunc = eta2dCGFuncOld = eta2dCGGradFunc = NULL;
......@@ -474,9 +468,7 @@ void slim3dSolverFunctions::initialize()
// TODO rename
const functor *zBottomFunc2d = dofs->zBotDof2d->getFunction();
const functor *uvBottomFunc2d = dofs->uvBotDof2d->getFunction();
const functor *uvIntFunc2d = dofs->uvIntDof2d->getFunction();
// bottom friction velocity in 2d and 3d geometries
zBotDistFunc = new zBotDist(bathFunc2d, zBottomFunc2d);
zBotDist2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
if (!bottomFrictionDragCoeff2d ){
......@@ -490,8 +482,6 @@ void slim3dSolverFunctions::initialize()
bottomFrictionFunc = new dragNormUb(bottomFrictionDragCoeff2d, uvBottomFunc2d);
bottomFriction2dPrecompFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
bottomUVDeviation = new bottomVelocityDeviation(bathFunc2d,eta2dCGFunc,uvBottomFunc2d,
uvIntFunc2d);
}
if(_solver->getSolveGMVel()){
eadySpeedFunc = new eadySpeed(dofs->nbDof3d->getFunction(), bathFunc2d, dofs->etaDof2d->getFunction());
......@@ -529,18 +519,6 @@ void slim3dSolverFunctions::initialize()
wMeshPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1);
wMeshDzPrecompFunc = new functionPrecomputed(*_solver->groups3d, _solver->groups3d->getElementGroup(0)->getOrder() * 2 + 1);
}
/*uvCorrectorFunc = new uvCorrector( dofs->uvIntDof2d->getFunction(),
dofs->uvIntTargetDof3d->getFunction(),
eta2dCGFunc,
bathFunc2d );*/
uvDevFunc = new uvDeviation( dofs->uvDof->getFunction(),
uvInt2dFunc,
eta2dCGFunc,
bathFunc2d );
uvDevProdFunc = new uvDeviationProducts(uvDevFunc);
uvDevAdvectionFunc = new uvDevAdvectionTerm(dofs->uvDeviationProdIntDof2d->getFunctionGradient(),eta2dCGFunc,bathFunc2d);
_initialized= true;
}
......@@ -561,6 +539,8 @@ void slim3dSolverFunctions::initializeBath()
eta2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2d->getFunction());
uvInt2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->uvIntDof2d->getFunction());
uvAv2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->uvAvDof2d->getFunction());
u2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->uvAvDof2d->getFunction());
uvAvGrad2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
uvIntTarget2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);
rhoSurf2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->rhoSurfDof2d->getFunction());
rhoSurfGrad2dFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->rhoSurfDof2d->getFunction());
......@@ -570,7 +550,7 @@ void slim3dSolverFunctions::initializeBath()
eta2dCGFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunction());
eta2dCGFuncOld = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunction());
eta2dCGGradFunc = new functionPrecomputedExtrusion(_solver->extrusion(), _solver->groups2d->getElementGroup(0)->getOrder() * 2 + 1);//*dofs->etaDof2dCG->getFunctionGradient());
_solver->verticalBottomRemover->computeHeightAboveVertBottom(bathFunc2d);
//_solver->verticalBottomRemover->computeHeightAboveVertBottom(bathFunc2d);
}
slim3dSolverFunctions::~slim3dSolverFunctions()
......@@ -579,7 +559,6 @@ slim3dSolverFunctions::~slim3dSolverFunctions()
if (zeroFunc3) delete zeroFunc3;
if (tinyFunc) delete tinyFunc;
if (oneFunc) delete oneFunc;
if (uvCorrectorFunc) delete uvCorrectorFunc;
if (rhoFunc) delete rhoFunc;
if (rhoFuncCenter) delete rhoFuncCenter;
if (eosAlphaFunc) delete eosAlphaFunc;
......@@ -612,6 +591,8 @@ slim3dSolverFunctions::~slim3dSolverFunctions()
if (eta2dCGGradFunc) delete eta2dCGGradFunc;
if (uvInt2dFunc) delete uvInt2dFunc;
if (uvAv2dFunc) delete uvAv2dFunc;
if (u2dFunc) delete u2dFunc;
if (uvAvGrad2dFunc) delete uvAvGrad2dFunc;
if (rhoSurf2dFunc) delete rhoSurf2dFunc;
if (rhoSurfGrad2dFunc) delete rhoSurfGrad2dFunc;
if (atmPressFunc) delete atmPressFunc;
......@@ -629,7 +610,6 @@ slim3dSolverFunctions::~slim3dSolverFunctions()
if (jumpDiffTracerVertS) delete jumpDiffTracerVertS;
if (jumpDiffTracerVertT) delete jumpDiffTracerVertT;
if (nuh) delete nuh;
if (nuh2d) delete nuh2d;
if (kappavS) delete kappavS;
if (kappavT) delete kappavT;
if (kappavJumpS) delete kappavJumpS;
......@@ -767,10 +747,9 @@ void slim3dSolverEquations::initialize()
//const functor *eta2d = dofs->etaDof2d->getFunction();
const functor *eta2d = funcs->eta2dFunc;
const functor *uv3d = dofs->uvDof->getFunction();
//const functor *uvInt2d = dofs->uvIntDof2d->getFunction();
const functor *uvInt2d = funcs->uvInt2dFunc;
//const functor *uvAv2d = dofs->uvAvDof2d->getFunction();
const functor *uvAv2d = funcs->uvAv2dFunc;
const functor *uvCorr2d = dofs->uvCorrDof2d->getFunction();
const functor *uv2d = funcs->uvAv2dFunc;
const functor *uv2dGrad = funcs->uvAvGrad2dFunc;
const functor *w3d = dofs->wDof3d->getFunction();
//const functor *w3dCopy = dofs->wDof3dCopy->getFunction();
//const functor *wMesh = funcs->wMeshFunc;
......@@ -788,12 +767,11 @@ void slim3dSolverEquations::initialize()
const functor *rhoSurf2d = funcs->rhoSurf2dFunc;
const functor *rhoSurfGrad2d = funcs->rhoSurfGrad2dFunc;
//const functor *rBot2d = dofs->rBotDof2d->getFunction();
const functor* bottomFrictionDragCoeff2d = funcs->bottomFrictionDragCoeff2d;
const functor* bottomUVDeviation = funcs->bottomUVDeviation;
horMomEq->setBathymetry(bath2d);
horMomEq->setBathymetryMinCG(bathCGMin2d);
horMomEq->setEta(eta2d);
horMomEq->setUV2d(uvInt2d);
horMomEq->setUV2d(uv2d);
horMomEq->setUV2dGrad(uv2dGrad);
if (_solver->_solveUVAdvVert){
horMomEq->setW(funcs->zeroFunc);
horMomEq->setWMesh(funcs->zeroFunc);
......@@ -807,8 +785,6 @@ void slim3dSolverEquations::initialize()
horMomEq->setRGrad(rGrad);
horMomEq->setRhoSurf(rhoSurf2d);
horMomEq->setRhoSurfGrad(rhoSurfGrad2d);
if (funcs->atmPressFunc)
horMomEq->setAtmPress(funcs->atmPressFunc);
if (funcs->coriolisFunc)
horMomEq->setCoriolisFactor(funcs->coriolisFunc);
if (funcs->nuh)
......@@ -997,14 +973,10 @@ void slim3dSolverEquations::initialize()
}
eta2dEq->setBathymetry(bath2d);
if ( _solver->_useModeSplit )
eta2dEq->setUV(uvAv2d);
else {
eta2dEq->setUsesDepthIntegratedUV(true);
eta2dEq->setUV(uvInt2d);
}
eta2dEq->setUV(uv2d);
uv2dEq->setBathymetry(bath2d);
uv2dEq->setEta(eta2d);
uv2dEq->setCorrectionSourceTerm(uvCorr2d);
//uv2dEq->setRInt(rInt2d);
//uv2dEq->setRInt(funcs->zeroFunc);
//uv2dEq->setRBot(rBot2d);
......@@ -1013,19 +985,8 @@ void slim3dSolverEquations::initialize()
uv2dEq->setRhoSurfGrad(rhoSurfGrad2d);
if (funcs->atmPressFunc)
uv2dEq->setAtmPress(funcs->atmPressFunc);
if (funcs->windStressFunc)
uv2dEq->setWindStress(funcs->windStressFunc);
if ( _solver->getComputeBottomFriction() ) {
uv2dEq->setBottomFrictionCoefficient(bottomFrictionDragCoeff2d);
uv2dEq->setBottomUV(bottomUVDeviation);
}
if (funcs->coriolisFunc)
uv2dEq->setCoriolisFactor(funcs->coriolisFunc);
if (funcs->nuh2d)
uv2dEq->setNuH(funcs->nuh2d);
if (funcs->uvDevAdvectionFunc)
uv2dEq->setUVDevAdvection(funcs->uvDevAdvectionFunc);
bool linear = false;
eta2dEq->setIsLinear(linear);
......
modules/slim3d/slim3dEquations.h
View file @ e259abb5
......@@ -30,7 +30,6 @@ private:
public:
const functor *zeroFunc, *zeroFunc2, *zeroFunc3, *tinyFunc, *oneFunc, *xyzFunc3d, *xyzFunc2d, *zFunc, *timeFunc;
const functor *uvCorrectorFunc;
const functor *pHydroStatFunc;
const functor *rhoFunc, *rhoFuncCenter, *eosAlphaFunc, *eosBetaFunc;
const functor *rhoGradFunc;
......@@ -53,7 +52,6 @@ public:
const functor *windStressFunc, *windFunc; //windStressFunc is the stress function used for horMomEq. windFunc is the actual wind.
const functor *atmPressFunc;
//const functor *rGradFunc;
const functor *uvDevFunc, *uvDevProdFunc, *uvDevAdvectionFunc;
const functor *eadySpeedFunc, *GMIndepTermFunc, *GMVelFunc;
const functor *slopeFunc, *slopeTap;
const functor *arbitraryMovingMeshFunc;
......@@ -67,6 +65,8 @@ public:
functionPrecomputedExtrusion *eta2dCGGradFunc;
functionPrecomputedExtrusion *uvInt2dFunc;
functionPrecomputedExtrusion *uvAv2dFunc;
functionPrecomputedExtrusion *u2dFunc;
functionPrecomputedExtrusion *uvAvGrad2dFunc;
functionPrecomputedExtrusion *rhoSurf2dFunc;
functionPrecomputedExtrusion *rhoSurfGrad2dFunc;
functionPrecomputedExtrusion *wMeshSurf2dPrecompFunc;
......@@ -177,6 +177,8 @@ public:
dgDofContainer *uDof,*vDof;
/** uv horizontal 3d velocity */
dgDofContainer *uvDof;
/** uv horizontal deviation velocity (u' = u - u2d) */
dgDofContainer *uvDevDof;
/** uv horizontal Gradient 3d Velocity */
dgDofContainer *uvHorDivDof;
/** uv in center of triangles for GOTM */
......@@ -187,8 +189,10 @@ public:
dgDofContainer *uvAvDof2d;
/** depth integrated uvDof */
dgDofContainer *uvIntDof2d;
/** for correcting uvDof to match uvAvDof2d */
/** for correcting uvDevDof */
dgDofContainer *uvCorrDof3d;
/** term source of 2d momentum */
dgDofContainer *uvCorrDof2d;
/** uvDof evaluated at the middle of bottom most element for bottom friction */
dgDofContainer *uvBotDof2d;
/** bottom friction velocity */
......@@ -197,9 +201,6 @@ public:
dgDofContainer *uvTauBDof2dCenter;
/** windStress in center of triangles for GOTM */
dgDofContainer *windStressDof2dCenter;
/** for depth integrated advection of uvDeviation */
dgDofContainer *uvDeviationProdDof3d;
dgDofContainer *uvDeviationProdIntDof2d;
/** w vertical velocity */
dgDofContainer *wDof3d;
/** w vertical velocity copy for jumps in w*/
......@@ -220,9 +221,9 @@ public:
/** r, baroclinic head */
dgDofContainer *rDof3d;
/** r integrated in vertical */
dgDofContainer *rIntDof2d;
//dgDofContainer *rIntDof2d;
/** r integrated in vertical */
dgDofContainer *rBotDof2d;
//dgDofContainer *rBotDof2d;
/** rho, density */
dgDofContainer *rhoDof3d;
dgDofContainer *rhoGradDof3d;
......@@ -279,7 +280,6 @@ public:
dgDofContainer *nbDof3d;
/** uvGrad, for smagorinsky diffusion */
dgDofContainer *uvGradDof;
dgDofContainer *uvGradDof2d;
/** density gradient with cg weak formulation */
dgDofContainer *rhoGradGMVelDof3d;
dgDofContainer * streamFuncDof3d;
......
modules/slim3d/slim3dFunctions.cpp
View file @ e259abb5
......@@ -202,26 +202,6 @@ void dragNormUb::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
bottomVelocityDeviation::bottomVelocityDeviation(const functor *bath, const functor *eta, const functor *uvBot2d, const functor *uvInt2d) {
_bath = bath;
_eta = eta;
_uvBot = uvBot2d;
_uvInt = uvInt2d;
}
void bottomVelocityDeviation::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &eta = cm.get(*_eta);
const fullMatrix<double> &uvBot = cm.get(*_uvBot);
const fullMatrix<double> &uvInt = cm.get(*_uvInt);
val.resize(cm.nEvaluationPoint(),2,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
double H = eta(i,0) + bath(i,0);
val(i,0) = uvBot(i,0) - uvInt(i,0)/H;
val(i,1) = uvBot(i,1) - uvInt(i,1)/H;
}
}
uvBotStar::uvBotStar(const functor *dragCoeff, const functor *uvBot) {
_dragCoeff = dragCoeff;
_uvBot = uvBot;
......@@ -237,109 +217,22 @@ void uvBotStar::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
uvCorrector::uvCorrector(const functor *uvInt, const functor *uvIntTarget, const functor *eta, const functor *bath, const functor *heightAboveVertBottom) {
_uvIntTarget = uvIntTarget;
uvCorrector::uvCorrector(const functor *uvInt, const functor *eta, const functor *bath, double dt3d) {
_uvInt = uvInt;
_eta = eta;
_bath = bath;
_heightAboveVertBottom = heightAboveVertBottom;
_dt3d = dt3d;
}
void uvCorrector::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &uvIntTarget = cm.get(*_uvIntTarget);
const fullMatrix<double> &uvInt = cm.get(*_uvInt);
const fullMatrix<double> &eta = cm.get(*_eta);
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &heightAboveVertBottom = cm.get(*_heightAboveVertBottom);
val.resize(cm.nEvaluationPoint(),2,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
//double H = eta(i,0) + bath(i,0);
//val(i,0) = uvIntTarget(i,0) - uvInt(i,0)/H; // uvTarget is depth average
//val(i,1) = uvIntTarget(i,1) - uvInt(i,1)/H;
if (heightAboveVertBottom(i,0) > 0){
double H = eta(i,0) + bath(i,0);
double HStar = eta(i,0) + heightAboveVertBottom(i,0);
val(i,0) = (H*uvIntTarget(i,0) - uvInt(i,0))/HStar; // uvTarget is depth average
val(i,1) = (H*uvIntTarget(i,1) - uvInt(i,1))/HStar;
}
else{
val(i,0) = 0.;
val(i,1) = 0.;
}
}
}
uvDeviation::uvDeviation(const functor *uv, const functor *uvInt, const functor *eta, const functor *bath) {
_uvInt = uvInt;
_uv = uv;
_eta = eta;
_bath = bath;
}
void uvDeviation::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &uvInt = cm.get(*_uvInt);
const fullMatrix<double> &uv = cm.get(*_uv);
const fullMatrix<double> &eta = cm.get(*_eta);
const fullMatrix<double> &bath = cm.get(*_bath);
val.resize(cm.nEvaluationPoint(),2,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
double d = eta(i,0) + bath(i,0);
val(i,0) = uv(i,0) - uvInt(i,0)/d;
val(i,1) = uv(i,1) - uvInt(i,1)/d;
}
}
uvDeviationProducts::uvDeviationProducts(const functor *uvDev) {
_uvDev = uvDev;
}
void uvDeviationProducts::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &uvDev = cm.get(*_uvDev);
val.resize(cm.nEvaluationPoint(),3,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
val(i,0) = uvDev(i,0) * uvDev(i,0);
val(i,1) = uvDev(i,0) * uvDev(i,1);
val(i,2) = uvDev(i,1) * uvDev(i,1);
}
}
uvDevAdvectionTerm::uvDevAdvectionTerm(const functor *uvDevProdIntGrad, const functor *eta, const functor *bath) {
_uvDevProdIntGrad = uvDevProdIntGrad;
_eta = eta;
_bath = bath;
}
void uvDevAdvectionTerm::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &uvDevProdIntGrad = cm.get(*_uvDevProdIntGrad);
const fullMatrix<double> &eta = cm.get(*_eta);
const fullMatrix<double> &bath = cm.get(*_bath);
val.resize(cm.nEvaluationPoint(),2,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
double H2 = pow_int(eta(i,0)+bath(i,0),2);
double duudx = uvDevProdIntGrad(i,0);
double duvdx = uvDevProdIntGrad(i,3);
double dvvdy = uvDevProdIntGrad(i,7);
double duvdy = uvDevProdIntGrad(i,4);
val(i,0) = - duudx/H2 - duvdy/H2;
val(i,1) = - duvdx/H2 - dvvdy/H2;
}
}
uvAv2uvIntFunction::uvAv2uvIntFunction(const functor *uvAv2d, const functor *eta, const functor *bath) {
_uvAv2d = uvAv2d;
_eta = eta;