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

slim3d: adaptivity as function of drho/dz

compute thermocline
compute Thermocline depth
parent 92608567
......@@ -76,7 +76,7 @@ class Slim3d_equations :
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_vertical_adaptation(self, tau, uv_factor, rho_factor, minimum_height=1, maximum_height=1e10, background_error=0, resize_factor=1e10, smoothing_factor=1) :
def set_vertical_adaptation(self, tau, uv_factor, rho_factor, minimum_height=1, maximum_height=1e10, background_error=0, resize_factor=1e10, smoothing_factor=1, vertical_gradient_factor=0) :
self._slimSolver.setUseAdaptiveVerticalGrid(True)
self._verticalAdaptation['tau'] = tau
self._verticalAdaptation['uv_factor'] = uv_factor
......@@ -86,6 +86,7 @@ class Slim3d_equations :
self._verticalAdaptation['background_error'] = background_error
self._verticalAdaptation['resize_factor'] = resize_factor
self._verticalAdaptation['smoothing_factor'] = smoothing_factor
self._verticalAdaptation['vertical_gradient_factor'] = vertical_gradient_factor
def set_lax_friedrichs_factor(self, factor=1.):
self._LFF = factor
......
......@@ -269,6 +269,7 @@ def slim3d_setup(loop):
adaptVertMesh.setSmoothingFactor(eq._verticalAdaptation['smoothing_factor'])
adaptVertMesh.setUVFactor(eq._verticalAdaptation['uv_factor'])
adaptVertMesh.setRhoFactor(eq._verticalAdaptation['rho_factor'])
adaptVertMesh.setVertGradFactor(eq._verticalAdaptation['vertical_gradient_factor'])
bottomTags = eq._domain._bottomTags
topTags = eq._domain._topTags
......
......@@ -115,7 +115,7 @@ class dgExtrusion {
maxLev = std::max(maxLev,getNbDGPointsInVertical(i));
size_t nLevMax;
#ifdef HAVE_MPI
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_UNSIGNED_LONG, MPI_MAX, MPI_COMM_WORLD);
#else
nLevMax = maxLev;
#endif
......@@ -128,7 +128,7 @@ class dgExtrusion {
maxLev = std::max(maxLev,getNbCGPointsInVertical(i));
size_t nLevMax;
#ifdef HAVE_MPI
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_UNSIGNED_LONG, MPI_MAX, MPI_COMM_WORLD);
#else
nLevMax = maxLev;
#endif
......
......@@ -278,7 +278,6 @@ void dgConservationLawSW3dMomentum::boundaryBottom(functorCache &m, fullMatrix<d
val(i,0) = -Cd * uBot;
val(i,1) = -Cd * vBot;
}
val.setAll(0.);
}
void dgConservationLawSW3dMomentum::boundaryWall(functorCache &m, fullMatrix<double> &val) const
......
......@@ -622,12 +622,28 @@ class uvNormFunc : public functor {
}
};
class vertGradFunc : public functor {
const functor *_fGrad;
public:
vertGradFunc(const functor *fGrad){
_fGrad = fGrad;
}
void operator()(functorCache &cm, fullMatrix<double> &val) const{
const fullMatrix<double> &fGrad = cm.get(*_fGrad);
val.resize(cm.nEvaluationPoint(), 1, false);
for (int i = 0; i < cm.nEvaluationPoint();++i) {
val(i,0) = fabs(fGrad(i,2));
}
}
};
static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofContainer &rhoDof, const dgDofContainer &uvDof, dgDofContainer &zDof)
{
slim3dVerticalAdaptiveMesh* verticalMeshAdapter = s.getVerticalMeshAdapter();;
dgDofContainer vJump(*s.groups3d,1);
dgDofContainer processedField(*s.groups3d,1);
dgDofContainer processedFieldVertGrad(*s.groups3d,1);
dgDofContainer uvNorm(*s.groups3d,1);
uvNormFunc uvNormF(uvDof.getFunction());
uvNorm.interpolate(&uvNormF);
......@@ -636,6 +652,10 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
processedField.axpy(uvNorm, verticalMeshAdapter->getUVFactor());
verticalMeshAdapter->computeJumps(&processedField, &vJump);
vertGradFunc vertGradF(processedField.getFunctionGradient());
processedFieldVertGrad.interpolate(&vertGradF);
vJump.axpy(processedFieldVertGrad, verticalMeshAdapter->getVertGradFactor());
//double DT = _dt;
//d->testDof3d->copy(*d->rhoDof3d);
////s->copyField3d->copy(d->uvDof,d->testDof3d,0,0);
......@@ -648,6 +668,12 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
//dgDofContainer *vJump2 = new dgDofContainer(*s->groups3d,1);
//s->getVerticalMeshAdapter()->computeJumps(d->testDof3d, vJump2);
functionPrecomputedExtrusion nodalVol2dFunc(s.extrusion(), 3);
nodalVol2dFunc.compute(*s.nodalVolume2d->get()->getFunction(), functorCache::NODE_GROUP_MODE);
dgDofContainer nodalVolume2d(*s.groups3d, 1);
nodalVolume2d.interpolate(&nodalVol2dFunc);
dgPrismDofConverter prismConverter(&nodalVolume2d);
dgCGMeanFilter dgCG3d(s.groups3d, &nodalVolume2d);
dgDofContainer meshVelocity(*s.groups3d,1);
verticalMeshAdapter->importDataFromDof(&vJump, &vJump, &zDof);
......@@ -655,6 +681,7 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
verticalMeshAdapter->exportZVel2Dof(&meshVelocity);
dgDofContainer newZDof(*s.groups3d,1);
verticalMeshAdapter->exportNewZ2Dof(&newZDof);
dgCG3d.apply(&newZDof, &newZDof);
dgDofContainer fixedNodesDof(*s.groups3d,1);
verticalMeshAdapter->exportFixedNodes2Dof(&fixedNodesDof);
//double dt = verticalMeshAdapter->getDt3d();
......@@ -662,13 +689,7 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
//zDof.exportMsh("newZ1");
dgDofContainer newZSmoothDof(*s.groups3d,1);
functionPrecomputedExtrusion nodalVol2dFunc(s.extrusion(), 3);
nodalVol2dFunc.compute(*s.nodalVolume2d->get()->getFunction(), functorCache::NODE_GROUP_MODE);
dgDofContainer nodalVolume2d(*s.groups3d, 1);
nodalVolume2d.interpolate(&nodalVol2dFunc);
dgPrismDofConverter prismConverter(&nodalVolume2d);
prismConverter.convertToPZeroHorizontal(&newZDof, &newZSmoothDof);
dgCGMeanFilter dgCG3d(s.groups3d, &nodalVolume2d);
dgCG3d.apply(&newZSmoothDof, &newZSmoothDof);
double lambda = verticalMeshAdapter->getSmoothingFactor();
newZSmoothDof.scale(lambda);
......@@ -697,6 +718,7 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
//verticalMeshAdapter->limitVerticalMeshVelMinMax(&meshVelocity, &minDof, &maxDof);
verticalMeshAdapter->limitVerticalMeshHeight(&zDof, &newZSmoothDof);
dgCG3d.apply(&newZSmoothDof, &newZSmoothDof);
verticalMeshAdapter->rebuildFixedNodes(&zDof, &newZSmoothDof, &fixedNodesDof);
zDof.copy(newZSmoothDof);
......@@ -914,6 +936,7 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
ts->solve(uvStar.getFunction(), wStarFull.getFunction(), &etaFunc, _time, _dt, massnext);
timerTracerEq.stop();
timerRGradEq.start();
dgDofContainer rGradDof3d(*s.groups3d, 2);
functionPrecomputedExtrusion rFunc(s.extrusion(), 3);
......
......@@ -1269,7 +1269,7 @@ slim3dVerticalAdaptiveMesh::slim3dVerticalAdaptiveMesh(dgExtrusion* columns)
_topError = 0;
_topLayer = 0;
_smoothingFactor = 1;
_uvFactor = _rhoFactor = 0;
_uvFactor = _rhoFactor = _vertGradFactor = 0;
_minDof = new dgDofContainer(*_groups3d, 1);
_maxDof = new dgDofContainer(*_groups3d, 1);
_initialized = false;
......@@ -1865,3 +1865,36 @@ void slim3dVerticalAdaptiveMesh::rebuildFixedNodes(const dgDofContainer* zDof, d
}
}
}
void computeThermocline(dgExtrusion* columns, const dgDofContainer* TDof, const dgDofContainer* zDof, dgDofContainer* thermoclineDof, double TLim)
{
dgFullMatrix<double> TData, zData, thermoData;
thermoclineDof->setAll(999);
for (size_t iCol = 0; iCol < columns->getNbColumns(); ++iCol) {
for (size_t iL = 0; iL < columns->getNbElemsInColumn(iCol); ++iL) {
size_t iGroup, iElem;
columns->mapColumnPosition2ElementID(iCol, iL, iGroup, iElem); // we loop from the top
int nbNodes = columns->getGroups3d().getElementGroup(iGroup)->getNbNodes();
TDof->getElementProxy(iGroup,iElem,TData);
zDof->getElementProxy(iGroup,iElem,zData);
for (int iNode = 0; iNode < nbNodes/2; iNode++) {
size_t iGroup2, iElem2, iNode2;
columns->map3d2d(iGroup, iElem, iNode, iGroup2, iElem2, iNode2);
thermoclineDof->getElementProxy(iGroup2,iElem2,thermoData);
if ((TData(iNode,0) < TLim) && (iL == 0)){
thermoData(iNode2,0) = -999;
continue;
}
if ((TData(iNode+nbNodes/2,0) > TLim) || (thermoData(iNode2,0) < 998))
continue;
if (TData(iNode,0) < TLim)
thermoData(iNode2,0) = -zData(iNode,0);
else{
double xsi = (TLim-TData(iNode,0)) / (TData(iNode+nbNodes/2,0)-TData(iNode,0));
double z = zData(iNode,0) + xsi * (zData(iNode+nbNodes/2,0)-zData(iNode,0));
thermoData(iNode2,0) = -z;
}
}
}
}
}
......@@ -189,7 +189,7 @@ private:
size_t _topLayer;
int _nb3dIter; // to advance the rhoDof
double *_dt3d;
double _smoothingFactor, _uvFactor, _rhoFactor;
double _smoothingFactor, _uvFactor, _rhoFactor, _vertGradFactor;
dgDofContainer *_maxDof, *_minDof;
bool _initialized;
public:
......@@ -244,6 +244,10 @@ public:
inline double getUVFactor() { return _uvFactor;};
inline void setRhoFactor(double f) { _rhoFactor = f;};
inline double getRhoFactor() { return _rhoFactor;};
inline void setVertGradFactor(double f) { _vertGradFactor = f;};
inline double getVertGradFactor() { return _vertGradFactor;};
};
void computeThermocline(dgExtrusion* columns, const dgDofContainer* TDof, const dgDofContainer* zDof, dgDofContainer* thermoclineDof, double TLim);
#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