Commit 94b9f369 authored by Philippe Delandmeter's avatar Philippe Delandmeter

slim3d: Vertical adaptation. benchmarks + smoothing in sigma coords

parent 7c9d09ad
Pipeline #2418 failed with stage
in 1 minute and 32 seconds
......@@ -18,11 +18,11 @@ equations.set_implicit_vertical(True)
equations.set_boundary_open(['river'], u=(data_dir+'open_bnd.nc', 'u'), v=(data_dir+'open_bnd.nc', 'v'), salinity=(data_dir+'open_bnd.nc', 'Sriv'))
equations.set_boundary_open(['open'], salinity=(data_dir+'open_bnd.nc', 'Ssea'))
equations.set_bottom_friction(False)
#equations.set_lax_friedrichs_factor(0.05)
equations.set_lax_friedrichs_factor(0.05)
equations.set_limiter(True)
equations.set_initial_temperature('vertical_gradient', None, 20, 0)
equations.set_initial_salinity('vertical_gradient', None, 35, 0.)
equations.set_vertical_adaptation(tau=200, uv_factor=0, rho_factor=1, minimum_height=.1, maximum_height=15, resize_factor=3, background_error=1e-4)
equations.set_vertical_adaptation(tau=200, uv_factor=0, rho_factor=1, minimum_height=.1, maximum_height=15, resize_factor=3, background_error=1e-4, smoothing_factor=1)
#equations.set_vertical_adaptation(tau=5, minimum_height=.1, maximum_height=15, resize_factor=3, background_error=1e-3)
......
......@@ -30,7 +30,7 @@ equations.set_limiter(True)
#equations.set_initial_temperature('netcdf', (pre_data_dir+'initial_temperature.nc', 'T_init'))
equations.set_initial_temperature('netcdf', pre_data_dir_base+'T_init/T_init.idx')
#equations.set_temperature_boundary_relaxation('netcdf', pre_data_dir_base+'T_init/T_init.idx', tau=4*86400, belowBottom=100)
equations.set_temperature_boundary_relaxation('netcdf', (pre_data_dir+'surface_temperature.nc', 'temp'), tau=10*86400, belowBottom=100)
equations.set_temperature_boundary_relaxation('netcdf', (pre_data_dir+'surface_temperature.nc', 'temp'), tau=10*86400, depth=12)
equations.set_initial_salinity('vertical_gradient', None, 0, 0.)
equations.set_wind_stress('speed', (pre_data_dir+'wind2000-2005.nc','u'), (pre_data_dir+'wind2000-2005.nc','v') )
equations.set_coriolis((pre_data_dir+'coriolis.nc', 'coriolis'))
......
......@@ -176,8 +176,8 @@ class Slim3d_equations :
def set_initial_temperature(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0):
self._initial_temperature = (mode, temperature, surface_temperature, vertical_gradient)
def set_temperature_boundary_relaxation(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0, tau=0, belowBottom=0):
self._temperature_boundary_relaxation = (mode, temperature, surface_temperature, vertical_gradient, tau, belowBottom)
def set_temperature_boundary_relaxation(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0, tau=0, depth=0):
self._temperature_boundary_relaxation = (mode, temperature, surface_temperature, vertical_gradient, tau, depth)
# TODO check value of air_density
def set_wind_stress(self, mode, wind_x, wind_y, density_air=1.25):
......
......@@ -359,7 +359,7 @@ def slim3d_setup(loop):
e.TEq.setBoundary0Flux(eq._boundary_coast)
e.TEq.setBoundary0Flux('vertical_bottom')
if f.TRelaxFunc:
e.TEq.setRelaxation(f.TRelaxFunc, eq._temperature_boundary_relaxation[4])
e.TEq.setRelaxation(f.TRelaxFunc, eq._temperature_boundary_relaxation[4], eq._temperature_boundary_relaxation[5])
for openBnd in eq._boundary_open:
e.TEq.setBoundaryNeumann(openBnd.tag, openBnd.t_flux_f)
......
......@@ -34,7 +34,6 @@ void dgConservationLawSW3dTracer::source(functorCache &cm, fullMatrix<double> &v
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++) {
......@@ -43,8 +42,8 @@ void dgConservationLawSW3dTracer::source(functorCache &cm, fullMatrix<double> &v
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;
if (xyz(i,2) > -_relaxDepth)
val(i,0) += ((*relaxRef)(i,0)-cE(i,0)) * (_relaxDepth+xyz(i,2))/_relaxDepth / _relaxTau;
}
}
}
......@@ -279,4 +278,5 @@ dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawSW
_dt=0;
_T=1;
_relaxTau=0;
_relaxDepth = 0;
}
......@@ -26,7 +26,7 @@ class dgConservationLawSW3dTracer : public dgConservationLawSW3dIMEX {
const functor *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_sourceFunc;
int _dt, _T;
const functor *_ref;
double _relaxTau;
double _relaxTau, _relaxDepth;
const functor *_relaxRef;
void source(functorCache &cm, fullMatrix<double> &val) const ;
void advection(functorCache &cm, fullMatrix<double> &val) const ;
......@@ -42,7 +42,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 setRelaxation(const functor *f, double tau, double depth) {_relaxRef = f; _relaxTau = tau; _relaxDepth = depth;};
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;};
......
......@@ -29,6 +29,7 @@ slim3dSolverDofs::slim3dSolverDofs(slim3dSolver *solver): _solver(solver) {
//sedDof = sedBottomDof2d = sedGroundDof2d = sedBottomFluxDof2d = NULL;
sedDof.resize(0); sedBottomDof2d.resize(0); sedGroundDof2d.resize(0); sedBottomFluxDof2d.resize(0);
zDof = deltaZDof3d = deltaZDofSurf2d = zDofCenter = zBotDof2d = NULL;
zOrigDof = NULL;
zDofExport = NULL;
tmpzBotDof2d = NULL;
xyzOrigDof = nnDof = ssDof = nuvDof = nuvDofCenter = kappavDof = kappavDofCenter = NULL;
......@@ -62,6 +63,7 @@ void slim3dSolverDofs::initialize()
}
bathDof2d = new dgDofContainer(*groups2d, 1);
xyzOrigDof = new dgDofContainer(*groups3d, 3);
zOrigDof = new dgDofContainer(*groups3d, 1);
if ( _solver->_solveTurbulence ) {
nuvDof = new dgDofContainer(*groups3d, 1);
......@@ -242,6 +244,7 @@ slim3dSolverDofs::~slim3dSolverDofs()
//if ( wMeshDof3d ) delete wMeshDof3d;
if ( xyzOrigDof ) delete xyzOrigDof;
if ( zDof ) delete zDof;
if ( zOrigDof ) delete zOrigDof;
if ( zDofExport ) delete zDofExport;
if ( deltaZDof3d ) delete deltaZDof3d;
if ( deltaZDofSurf2d ) delete deltaZDofSurf2d;
......
......@@ -240,6 +240,7 @@ public:
std::vector<dgDofContainer *> sedDof, sedBottomDof2d, sedGroundDof2d, sedBottomFluxDof2d;
/** vertical coordinates of nodes for moving mesh */
dgDofContainer *zDof;
dgDofContainer *zOrigDof;
dgDofContainer *zDofExport;
dgDofContainer *deltaZDof3d;
dgDofContainer *deltaZDofSurf2d;
......
......@@ -411,6 +411,7 @@ void slim3dTimeIntegrator::_initializeFields()
zNewAdaptative newZFunc(f.eta2dCGFunc, f.eta2dCGFuncOld, s.verticalBottomRemover->heightAboveVertBottom->getFunction());
dgDofContainer zDof(*s.groups3d, 1);
zDof.interpolate(&newZFunc, _time);
d.zOrigDof->copy(zDof);
s.groups3d->updateCoordinates(&zDof);
s.newDepthIntegrator->update();
s.nodalVolume3d->update();
......@@ -656,17 +657,6 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
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);
//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);
functionPrecomputedExtrusion nodalVol2dFunc(s.extrusion(), 3);
nodalVol2dFunc.compute(*s.nodalVolume2d->get()->getFunction(), functorCache::NODE_GROUP_MODE);
......@@ -676,7 +666,7 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
dgCGMeanFilter dgCG3d(s.groups3d, &nodalVolume2d);
dgDofContainer meshVelocity(*s.groups3d,1);
verticalMeshAdapter->importDataFromDof(&vJump, &vJump, &zDof);
verticalMeshAdapter->importDataFromDof(&vJump, &zDof);
verticalMeshAdapter->computeVertMeshSize();
verticalMeshAdapter->exportZVel2Dof(&meshVelocity);
dgDofContainer newZDof(*s.groups3d,1);
......@@ -684,39 +674,19 @@ static void _computeAdaptiveMovingMeshVelocity(slim3dSolver &s, const dgDofConta
dgCG3d.apply(&newZDof, &newZDof);
dgDofContainer fixedNodesDof(*s.groups3d,1);
verticalMeshAdapter->exportFixedNodes2Dof(&fixedNodesDof);
//double dt = verticalMeshAdapter->getDt3d();
//zDof.axpy(meshVelocity, dt);
//zDof.exportMsh("newZ1");
slim3dSolverDofs &d = *s.getDofs();
newZDof.axpy(*d.zOrigDof, -1.);
dgDofContainer newZSmoothDof(*s.groups3d,1);
prismConverter.convertToPZeroHorizontal(&newZDof, &newZSmoothDof);
dgCG3d.apply(&newZSmoothDof, &newZSmoothDof);
newZSmoothDof.axpy(*d.zOrigDof, 1.);
newZDof.axpy(*d.zOrigDof, 1.);
double lambda = verticalMeshAdapter->getSmoothingFactor();
newZSmoothDof.scale(lambda);
newZSmoothDof.axpy(newZDof, 1-lambda);
verticalMeshAdapter->rebuildFixedNodes(&newZDof, &newZSmoothDof);
//slim3dForceZeroAtBnd forceZero(s.groups3d);
//forceZero.setup({"vertical_bottom"});
//forceZero.forceZeroAtBnd(&meshVelocity);
//dgDofContainer minDof(*s.groups3d,1);
//dgDofContainer maxDof(*s.groups3d,1);
//verticalMeshAdapter->limitVerticalMeshGradient(&zDof, &newZDof, &minDof, &maxDof);
//s.dgCG3d->dgCGMax(&minDof, &minDof);
////minDof.setAll(-1e10);
//s.dgCG3d->dgCGMin(&maxDof, &maxDof);
////maxDof.setAll(1e10);
//verticalMeshAdapter->limitVerticalMeshVelMinMax(&newZDof, &minDof, &maxDof);
//dgDofContainer minDof(*s.groups3d,1);
//dgDofContainer maxDof(*s.groups3d,1);
//verticalMeshAdapter->limitVerticalMeshGradient(&zDof, &meshVelocity, &minDof, &maxDof);
//s.dgCG3d->dgCGMax(&minDof, &minDof);
//s.dgCG3d->dgCGMin(&maxDof, &maxDof);
//verticalMeshAdapter->limitVerticalMeshVelMinMax(&meshVelocity, &minDof, &maxDof);
verticalMeshAdapter->limitVerticalMeshHeight(&zDof, &newZSmoothDof);
dgCG3d.apply(&newZSmoothDof, &newZSmoothDof);
verticalMeshAdapter->rebuildFixedNodes(&zDof, &newZSmoothDof, &fixedNodesDof);
......
......@@ -1352,29 +1352,26 @@ void slim3dVerticalAdaptiveMesh::computeJumps(const dgDofContainer *dof3d, dgDof
jumpDof3d->axpy(*_minDof,-1.);
}
void slim3dVerticalAdaptiveMesh::importDataFromDof(const dgDofContainer* dof3d, const dgDofContainer* dof3d2, const dgDofContainer* zdof3d)
void slim3dVerticalAdaptiveMesh::importDataFromDof(const dgDofContainer* dof3d, const dgDofContainer* zdof3d)
{
for (size_t i = 0; i < _2dCGNbNodes; ++i)
for (size_t j = 0; j < _vertJumps[i].size(); ++j)
_vertJumps[i][j] = 0;
dgFullMatrix<double> dataJump, dataJump2, dataZ;
dgFullMatrix<double> dataJump, dataZ;
std::vector<double> vertJumpLocal(4,0);
size_t iVertical, iLevel, iGroup2d, iElem2d, iNode2d;
size_t iGroup3d, iElem3d;
for (size_t iCol = 0; iCol < _columns->getNbColumns(); ++iCol) {
for (size_t iL = 0; iL < _columns->getNbElemsInColumn(iCol); ++iL) { //loop from the top
for (size_t iL = 0; iL < _columns->getNbElemsInColumn(iCol); ++iL) { //loop from top to bottom
_columns->mapColumnPosition2ElementID(iCol, iL, iGroup3d, iElem3d);
dof3d->getElementProxy(iGroup3d, iElem3d, dataJump);
dof3d2->getElementProxy(iGroup3d, iElem3d, dataJump2);
zdof3d->getElementProxy(iGroup3d, iElem3d, dataZ);
int nbNodes = _groups3d->getElementGroup(iGroup3d)->getNbNodes();
_columns->mapDofToVEdgeCG(iGroup3d, iElem3d, 0, iVertical, iLevel);
if (iL != iLevel)
Msg::Fatal("BUG HERE!! iL:%d, iLevel:%d\n", iL, iLevel);
_columns->mapVEdgeTo2dNode(iVertical,iGroup2d,iElem2d,iNode2d);
const std::vector<int> cgVertices = _HOMesh2d->vertices(iGroup2d, iElem2d);
if (iL == 0){ // TOP
if (iL == 0){ // top
for (int iNode = 0; iNode < nbNodes/2; iNode++){
_vertJumps[cgVertices[iNode]][iLevel] = -1;
_z[cgVertices[iNode]][iLevel] = dataZ(iNode, 0);
......@@ -1382,12 +1379,11 @@ void slim3dVerticalAdaptiveMesh::importDataFromDof(const dgDofContainer* dof3d,
}
else {
for (int iNode = 0; iNode < nbNodes/2; iNode++) {
_vertJumps[cgVertices[iNode]][iLevel] = std::max( (dataJump(iNode,0) * dataJump(iNode,0)),
(dataJump2(iNode,0) * dataJump2(iNode,0)) );
_vertJumps[cgVertices[iNode]][iLevel] = dataJump(iNode,0) * dataJump(iNode,0);
_z[cgVertices[iNode]][iLevel] = dataZ(iNode, 0);
}
}
if (iL == _columns->getNbElemsInColumn(iCol) -1){ // BOTTOM
if (iL == _columns->getNbElemsInColumn(iCol) -1){ // bottom
for (int iNode = 0; iNode < nbNodes/2; iNode++){
_vertJumps[cgVertices[iNode]][iLevel+1] = -1;
_z[cgVertices[iNode]][iLevel+1] = dataZ(iNode+nbNodes/2, 0);
......@@ -1395,12 +1391,6 @@ void slim3dVerticalAdaptiveMesh::importDataFromDof(const dgDofContainer* dof3d,
}
}
}
//for (size_t iCol = 0; iCol < _vertJumps.size(); ++iCol){
// for (size_t iLev = 0; iLev < _vertJumps[iCol].size(); ++iLev)
// printf("%g ", _vertJumps[iCol][iLev]);
// printf("\n");
//}
}
void _computeBackGroundError(double *z, double *backGroundError, size_t n)
......@@ -1457,80 +1447,6 @@ void slim3dVerticalAdaptiveMesh::computeVertMeshSizeOnOneStage(double *vJumps, d
std::vector<double> error(n-1,0);
std::vector<double> h(n-1,0);
#define byDiff
#ifndef byDiff
// Compute error and h
for (size_t i=0; i < n-1; ++i){
error[i] = vJumps[i] + vJumps[i+1];
h[i] = z[i] - z[i+1];
}
if (vJumps[0] < 0) error[0] = 2 * vJumps[1];
if (vJumps[n-1] < 0) error[n-2] = 2 * vJumps[n-2];
// Integrate 1/error
double oneOverErrInt = 0;
for (size_t i=0; i < n-1; ++i)
oneOverErrInt = oneOverErrInt + 1/error[i] * h[i];
double H = z[0] - z[n-1];
// Computing meshSize
newZ[n-1] = z[n-1];
zDiff[n-1] = 0;
for (int i=n-2; i >= 0; --i){
meshSizeP0[i] = h[i] * H/oneOverErrInt / error[i];
//meshSizeP0[i] = std::max(_hMin, std::min(_hMax, meshSizeP0[i]));
newZ[i] = newZ[i+1] + meshSizeP0[i];
zDiff[i] = newZ[i] - z[i];
}
//// Re-integrate 1/error and H where meshSizeP0 not saturated
//oneOverErrInt = 0;
//double hUnsat = 0;
//double hLost = 0;
//double tol = 1e-3;
//for (size_t i=0; i < n-1; ++i){
// bool saturated = (fabs(meshSizeP0[i]-_hMin) < tol) || (fabs(meshSizeP0[i]-_hMax) < tol);
// oneOverErrInt = oneOverErrInt + 1/error[i] * h[i] * !saturated;
// hUnsat = hUnsat + h[i] * !saturated;
// hLost = hLost + (h[i]-meshSizeP0[i]) * saturated;
//}
//// Re-computing meshSize and computing newZ
//newZ[n-1] = z[n-1];
//zDiff[n-1] = 0;
//for (int i=n-2; i >= 0; --i){
// bool saturated = (fabs(meshSizeP0[i]-_hMin) < tol) || (fabs(meshSizeP0[i]-_hMax) < tol);
// if (!saturated)
// meshSizeP0[i] = h[i] * hUnsat/oneOverErrInt / error[i] * (hUnsat+hLost)/hUnsat;
// newZ[i] = newZ[i+1] + meshSizeP0[i];
// zDiff[i] = newZ[i] - z[i];
//}
if ( fabs(newZ[0]-z[0]) > 1e-12 )
Msg::Error("BUG HERE!! newZ[0]: %.16e; z[0]: %.16e", newZ[0], z[0]);
newZ[0] = z[0];
zDiff[0] = 0;
// Simple limiting
for (size_t i=1; i < n-1; ++i){ // from top 2 bottom
double elemNewHeight = (z[i-1] + zDiff[i-1] / _tau * *_dt3d) - (z[i] + zDiff[i] / _tau * *_dt3d);
if (( elemNewHeight < _hMin ) & ( zDiff[i] > 0 ))
zDiff[i] = std::max(0., ((z[i-1] + zDiff[i-1] / _tau * *_dt3d) - _hMin - z[i]) * _tau / *_dt3d);
else if (( elemNewHeight > _hMax ) & ( zDiff[i] < 0 ))
zDiff[i] = std::min(0., ((z[i-1] + zDiff[i-1] / _tau * *_dt3d) - _hMax - z[i]) * _tau / *_dt3d);
}
for (size_t i=n-2; i > 0; --i){ // from bottom 2 top
double elemNewHeight = (z[i] + zDiff[i] / _tau * *_dt3d) - (z[i+1] + zDiff[i+1] / _tau * *_dt3d);
if (( elemNewHeight < _hMin ) & ( zDiff[i] < 0 ))
zDiff[i] = std::min(0., ((z[i+1] + zDiff[i+1] / _tau * *_dt3d) + _hMin - z[i]) * _tau / *_dt3d);
else if (( elemNewHeight > _hMax ) & ( zDiff[i] > 0 ))
zDiff[i] = std::min(0., ((z[i+1] + zDiff[i+1] / _tau * *_dt3d) + _hMax - z[i]) * _tau / *_dt3d);
}
#else
fixedNodes[0] = 1;
fixedNodes[n-1] = 1;
......@@ -1599,11 +1515,6 @@ void slim3dVerticalAdaptiveMesh::computeVertMeshSizeOnOneStage(double *vJumps, d
// zDiff[i] = std::min(0., ((z[i+1] + zDiff[i+1] / _tau * *_dt3d) + hMax - z[i]) * _tau / *_dt3d);
//}
//for (size_t i=0; i < n; ++i){
// newZ[i] = z[i] + zDiff[i] * (*_dt3d) / _tau;
//}
#endif
}
......
......@@ -197,7 +197,7 @@ public:
~slim3dVerticalAdaptiveMesh();
void initialize();
void computeJumps(const dgDofContainer *dof3d, dgDofContainer *jumpDof3d);
void importDataFromDof(const dgDofContainer* dof3d, const dgDofContainer* dof3d2, const dgDofContainer* zdof3d);
void importDataFromDof(const dgDofContainer* dof3d, const dgDofContainer* zdof3d);
void computeVertMeshSize();
void computeVertMeshSizeOnOneStage(double *vJumps, double *z, double *newZ, double *zDiff, double *backGroundError, double *meshSizeP0, double *fixedNodes, size_t length);
void exportCGData2Dof(const std::vector<std::vector<double> > vec, dgDofContainer* dof3d, double fact = 1);
......
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