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

slim3d: vertical adaptive mesh

parent f746cc20
......@@ -59,6 +59,7 @@ class Slim3d_equations :
self._areaFunc = None
self._nuh = None
self._add_shear = None
self._verticalAdaptation = {}
def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver
......@@ -85,6 +86,14 @@ class Slim3d_equations :
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_vertical_adaptation(self, tau, minimum_height=1, maximum_height=1e10, minimum_error=0, resize_factor=1e10) :
self._slimSolver.setUseAdaptiveVerticalGrid(True)
self._verticalAdaptation['tau'] = tau
self._verticalAdaptation['minimum_height'] = minimum_height
self._verticalAdaptation['maximum_height'] = maximum_height
self._verticalAdaptation['minimum_error'] = minimum_error
self._verticalAdaptation['resize_factor'] = resize_factor
def set_lax_friedrichs(self, factor=1.):
self._LFF = factor
......
......@@ -156,7 +156,8 @@ def slim3d_setup(loop):
if eq._initial_salinity:
if eq._initial_salinity[0] == 'netcdf':
f.SInitFunc = slim_private._load_function(eq._initial_salinity[1], slimSolver.groups3d)
eq._SInitFunc = slim_private._load_function(eq._initial_salinity[1], slimSolver.groups3d)
f.SInitFunc = eq._SInitFunc
elif eq._initial_salinity[0] == 'vertical_gradient':
f.SInitFunc = dgpy.functorNumpy(lambda cm : eq._initial_salinity[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._initial_salinity[3])
else:
......@@ -166,7 +167,8 @@ def slim3d_setup(loop):
if eq._initial_temperature:
if eq._initial_temperature[0] == 'netcdf':
f.TInitFunc = slim_private._load_function(eq._initial_temperature[1], slimSolver.groups3d)
eq._TInitFunc = slim_private._load_function(eq._initial_temperature[1], slimSolver.groups3d)
f.TInitFunc = eq._TInitFunc
elif eq._initial_temperature[0] == 'vertical_gradient':
f.TInitFunc = dgpy.functorNumpy(lambda cm : eq._initial_temperature[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._initial_temperature[3])
else:
......@@ -250,6 +252,15 @@ def slim3d_setup(loop):
e = slimSolver.createEquations()
d = slimSolver.getDofs()
if eq._verticalAdaptation:
adaptVertMesh = slimSolver.getVerticalMeshAdapter()
adaptVertMesh.setMinHeight(eq._verticalAdaptation['minimum_height'])
adaptVertMesh.setMaxHeight(eq._verticalAdaptation['maximum_height'])
adaptVertMesh.setTau(eq._verticalAdaptation['tau'])
adaptVertMesh.setMinError(eq._verticalAdaptation['minimum_error'])
adaptVertMesh.setResizeFact(eq._verticalAdaptation['resize_factor'])
bottomTags = eq._domain._bottomTags
topTags = eq._domain._topTags
......
......@@ -47,7 +47,7 @@ slim3dSolver::slim3dSolver(const std::string meshFile3d, std::vector<std::string
nodalVolume2d = NULL;
prismConverter = NULL;
movingMesh = NULL;
zDiff1dEq = NULL;
//zDiff1dEq = NULL;
uvLimiter = SLimiter = TLimiter = tkeLimiter = epsLimiter = sedLimiter= NULL;
jumpDiffTracerS = jumpDiffTracerT = jumpDiffTracerVertS = jumpDiffTracerVertT = NULL;
dgCG3d = dgCG2d = NULL;
......@@ -120,7 +120,7 @@ slim3dSolver::~slim3dSolver()
if ( copyField3d ) delete copyField3d;
if ( movingMesh ) delete movingMesh;
if ( _verticalAdaptMesh ) delete _verticalAdaptMesh;
if ( zDiff1dEq ) delete zDiff1dEq;
// if ( zDiff1dEq ) delete zDiff1dEq;
if ( uvLimiter ) delete uvLimiter;
if ( SLimiter ) delete SLimiter;
if ( TLimiter ) delete TLimiter;
......@@ -161,8 +161,8 @@ slim3dSolverEquations* slim3dSolver::createEquations()
_verticalAdaptMesh = new slim3dVerticalAdaptiveMesh(columnInfo);
_verticalAdaptMesh->setDt3d(functions->getDt3d());
}
if (_useAdaptiveVerticalGrid || _useAdaptiveVerticalGridWithDiffEq )
zDiff1dEq = new slim3dZDiff1dEq(columnInfo, functions->getDt3d());
// if (_useAdaptiveVerticalGrid || _useAdaptiveVerticalGridWithDiffEq )
// zDiff1dEq = new slim3dZDiff1dEq(columnInfo, functions->getDt3d());
if (_useUVLimiter){
//uvLimiter = new dgSlopeLimiterVertex (_equations->horMomEq,groups3d, nodalVolume3d->get(), _useOptimalLimiter);
//uvLimiter = new dgSlopeLimiter4 (_equations->horMomEq,groups3d, _dofs->uvDof);
......
......@@ -70,7 +70,7 @@ public:
dgPrismDofConverter *prismConverter;
dgFieldCopy *copyField2d, *copyField3d;
dgSW3dMovingMesh *movingMesh;
slim3dZDiff1dEq *zDiff1dEq;
//slim3dZDiff1dEq *zDiff1dEq;
dgSW3dVerticalBottomRemover *verticalBottomRemover;
dgLimiter *uvLimiter, *SLimiter, *TLimiter, *tkeLimiter, *epsLimiter, *sedLimiter, *sedLimiterGround2d;
dgJumpDiffusionTracer *jumpDiffTracerS, *jumpDiffTracerT, *jumpDiffTracerVertS, *jumpDiffTracerVertT;
......
......@@ -922,6 +922,46 @@ class uvMeanAboveVertBottom : public functor {
}
};
static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofContainer &rhoDof, const dgDofContainer &zDof, dgDofContainer &meshVelocity)
{
if ( s.getUseAdaptiveVerticalGrid() )
{
dgDofContainer vJump(*s.groups3d,1);
s.getVerticalMeshAdapter()->computeJumps(&rhoDof, &vJump);
//double DT = _dt;
//d->testDof3d->copy(*d->rhoDof3d);
////s->copyField3d->copy(d->uvDof,d->testDof3d,0,0);
//for (int i = 0; i < s->getVerticalMeshAdapter()->getNb3dIter(); ++i){
// e->vJumpEq->computeAllTerms(_time, *d->testDof3d, *vJumpDof_K, false);
// vJumpDof_K->multiplyByInvMassMatrix();
// d->testDof3d->axpy(*vJumpDof_K, DT * s->getVerticalMeshAdapter()->getDtPredictorFact());
// //s->SLimiter->apply(d->testDof3d);
//}
//dgDofContainer *vJump2 = new dgDofContainer(*s->groups3d,1);
//s->getVerticalMeshAdapter()->computeJumps(d->testDof3d, vJump2);
s.getVerticalMeshAdapter()->importDataFromDof(&vJump, &vJump, &zDof);
s.getVerticalMeshAdapter()->computeVertMeshSize();
s.getVerticalMeshAdapter()->exportZVel2Dof(&meshVelocity);
s.prismConverter->convertToPZeroHorizontal(&meshVelocity, &meshVelocity);
s.dgCG3d->apply(&meshVelocity, &meshVelocity);
dgDofContainer minDof(*s.groups3d,1);
dgDofContainer maxDof(*s.groups3d,1);
s.getVerticalMeshAdapter()->limitVerticalMeshGradient(&zDof, &meshVelocity, &minDof, &maxDof);
s.dgCG3d->dgCGMax(&minDof, &minDof);
s.dgCG3d->dgCGMin(&maxDof, &maxDof);
s.getVerticalMeshAdapter()->limitVerticalMeshVelMinMax(&meshVelocity, &minDof, &maxDof);
s.getVerticalMeshAdapter()->limitVerticalMeshHeight(&zDof, &meshVelocity);
}
else{
meshVelocity.setAll(0.);
}
}
void slim3dTimeIntegratorPC::advanceOneTimeStep()
{
slim3dSolver &s = *_solver;
......@@ -1042,7 +1082,7 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
dgDofContainer etaStarCG(*s.groups2d,1);
s.dgCG2d->apply(&etaStar, &etaStarCG);
{
f.eta2dCGFunc->compute(*etaStarCG.getFunction(),functorCache::NODE_GROUP_MODE);
//f.eta2dCGFunc->compute(*etaStarCG.getFunction(),functorCache::NODE_GROUP_MODE);
f.eta2dCGFuncOld->compute(*eta.getFunction(),functorCache::NODE_GROUP_MODE);
zNewAdaptative newZFunc(d.xyzOrigDof->getFunction(), &etaFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction(), f.zeroFunc, _dt, _time);
dgDofContainer zDof(*s.groups3d, 1);
......@@ -1090,7 +1130,12 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
f.eta2dCGFunc->compute(*etaNext.getFunction(),functorCache::NODE_GROUP_MODE);
f.eta2dCGFuncOld->compute(*etaStar.getFunction(),functorCache::NODE_GROUP_MODE);
dgDofContainer zDof(*s.groups3d, 1);
zNewAdaptative newZFunc(d.xyzOrigDof->getFunction(), f.eta2dCGFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction(), f.zeroFunc, _dt, _time);
dgDofContainer meshVelocity(*s.groups3d, 1);
s.copyField3d->copy(&s.groups3d->coordinatesDof(),&zDof,2,0);
dgDofContainer rhoForMesh(*s.groups3d, 1);
rhoForMesh.interpolate(f.rhoFunc);
_computeAdaptiveMovingMeshVelocity(s, rhoForMesh, zDof, meshVelocity);
zNewAdaptative newZFunc(d.xyzOrigDof->getFunction(), f.eta2dCGFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction(), meshVelocity.getFunction(), _dt, _time);
zDof.interpolate(&newZFunc, _time);
slim3dMassMatrix massnext(*s.groups3d, zDof);
slim3dMassMatrix imassnext(massnext);
......@@ -1187,7 +1232,7 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
e.horMomEq->setUV2d(&uvMeanStarFunc);
//e.horMomEq->setRhoDevMax(&rhopMaxFunc);
e.horMomEq->setRhoDevInt(&rFunc);
e.horMomEq->setLaxFriedrichsFactor(1.);
e.horMomEq->setLaxFriedrichsFactor(.05);
d.tmpuvGradDof->interpolate(uvStar.getFunctionGradient());
e.horMomEq->computeAllTerms(_time,uvStar,*uvDof_K,false);
d.uvDof->copy(uvStar);
......@@ -1229,7 +1274,10 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
f.eta2dCGFunc->compute(*etaNext.getFunction(),functorCache::NODE_GROUP_MODE);
f.eta2dCGFuncOld->compute(*etaStar.getFunction(),functorCache::NODE_GROUP_MODE);
dgDofContainer zDof(*s.groups3d, 1);
zNewAdaptative newZFunc(d.xyzOrigDof->getFunction(), f.eta2dCGFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction(), f.zeroFunc, _dt, _time);
dgDofContainer meshVelocity(*s.groups3d, 1);
s.copyField3d->copy(&s.groups3d->coordinatesDof(),&zDof,2,0);
_computeAdaptiveMovingMeshVelocity(s, rhoForMesh, zDof, meshVelocity);
zNewAdaptative newZFunc(d.xyzOrigDof->getFunction(), f.eta2dCGFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction(), meshVelocity.getFunction(), _dt, _time);
zDof.interpolate(&newZFunc, _time);
s.groups3d->updateCoordinates(&zDof);
s.newDepthIntegrator->update();
......
This diff is collapsed.
......@@ -129,33 +129,33 @@ public:
/** Class to derive vertical coordinates */
class slim3dZDiff1dEq {
private:
dgGroupCollection* _groups;
dgExtrusion* _columns;
double _minHeight;
double _maxHeight;
double *_dt3d;
public:
slim3dZDiff1dEq(dgExtrusion* columns, double *dt3d);
~slim3dZDiff1dEq();
void diffuseZ(const dgDofContainer* zDof, const dgDofContainer* kappaDof, dgDofContainer* zUpdateDof);
/** limits the vertical mesh arbitrary vel if the mesh is already too resolved */
void limitVerticalMeshHeight(const dgDofContainer* zDof, dgDofContainer* velDof);
/** limit the velocities such that mesh does not become bad conditionned.
* velMinMax has two fields. velDof must stay between velMinDof and velMaxDof.
* outputs are velMinDof and velMaxDof*/
void limitVerticalMeshGradient(const dgDofContainer* zDof, const dgDofContainer *velDof, dgDofContainer *velMinDof, dgDofContainer *velMaxDof);
/** limit velDof between min and max allowed */
void limitVerticalMeshVelMinMax(dgDofContainer* velDof, const dgDofContainer* velMinDof, const dgDofContainer* velMaxDof);
void computeVerticalJumps(const dgDofContainer* dof3d, dgDofContainer* vJumpDof3d);
void computeError(const dgDofContainer* vJumpDof3d, const dgDofContainer* vJumpShiftDof3d, const dgDofContainer* zDof3d, dgDofContainer* errorDof3d, dgDofContainer* oneOverErrorDof3d);
void computeZStar(const dgDofContainer* errorDof3d, const dgDofContainer* oneOverErrorIntDof2d, const dgDofContainer* zDof3d, const dgDofContainer* bathDof2d, const dgDofContainer* etaDof2dCG, dgDofContainer *newZDof3d);
inline void setMinHeight(double h) { _minHeight = h; };
inline double getMinHeight() { return _minHeight; };
inline void setMaxHeight(double h) { _maxHeight = h; };
inline double getMaxHeight() { return _maxHeight; };
};
//class slim3dZDiff1dEq {
//private:
// dgGroupCollection* _groups;
// dgExtrusion* _columns;
// double _minHeight;
// double _maxHeight;
// double *_dt3d;
//public:
// slim3dZDiff1dEq(dgExtrusion* columns, double *dt3d);
// ~slim3dZDiff1dEq();
// void diffuseZ(const dgDofContainer* zDof, const dgDofContainer* kappaDof, dgDofContainer* zUpdateDof);
// /** limits the vertical mesh arbitrary vel if the mesh is already too resolved */
// void limitVerticalMeshHeight(const dgDofContainer* zDof, dgDofContainer* velDof);
// /** limit the velocities such that mesh does not become bad conditionned.
// * velMinMax has two fields. velDof must stay between velMinDof and velMaxDof.
// * outputs are velMinDof and velMaxDof*/
// void limitVerticalMeshGradient(const dgDofContainer* zDof, const dgDofContainer *velDof, dgDofContainer *velMinDof, dgDofContainer *velMaxDof);
// /** limit velDof between min and max allowed */
// void limitVerticalMeshVelMinMax(dgDofContainer* velDof, const dgDofContainer* velMinDof, const dgDofContainer* velMaxDof);
// void computeVerticalJumps(const dgDofContainer* dof3d, dgDofContainer* vJumpDof3d);
// void computeError(const dgDofContainer* vJumpDof3d, const dgDofContainer* vJumpShiftDof3d, const dgDofContainer* zDof3d, dgDofContainer* errorDof3d, dgDofContainer* oneOverErrorDof3d);
// void computeZStar(const dgDofContainer* errorDof3d, const dgDofContainer* oneOverErrorIntDof2d, const dgDofContainer* zDof3d, const dgDofContainer* bathDof2d, const dgDofContainer* etaDof2dCG, dgDofContainer *newZDof3d);
// inline void setMinHeight(double h) { _minHeight = h; };
// inline double getMinHeight() { return _minHeight; };
// inline void setMaxHeight(double h) { _maxHeight = h; };
// inline double getMaxHeight() { return _maxHeight; };
//};
class slim3dForceZeroAtBnd {
private:
......
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