Skip to content
GitLab
Commit 55d702d8 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

Lax-Friedrichs for the new equations

parent 5f41cb50
Pipeline #1308 passed with stage
in 26 minutes and 55 seconds
Hide whitespace changes
Inline Side-by-side
dgpy/scripts/slim3d.py
View file @ 55d702d8
......@@ -49,6 +49,9 @@ class Slim3d_equations :
self._boundary_open = []
self._gotm_option_file = None
self._z0 = [0.005, 0.02]
self._LFF = 0.
self._linear2d = False
self._additional_shear = None
def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver
......@@ -69,6 +72,12 @@ class Slim3d_equations :
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_linear_2d_equations(self, flag):
self._linear2d = flag
def set_lax_friedrichs(self, factor=1.):
self._LFF = factor
def set_conservative_ale(self, flag):
self._slimSolver.setUseConservativeALE(flag)
......@@ -80,6 +89,10 @@ class Slim3d_equations :
self._hor_visc_fact = factor
self._hor_visc_max = maximum
def set_additional_artificial_horizontal_shear(self, additional_shear):
self._additional_shear = additional_shear
def set_horizontal_diffusivity(self, mode, constant=1, factor=0.03, maximum=1e9):
if (mode != "okubo") and (mode != "constant"):
dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode)
......@@ -232,6 +245,7 @@ class Loop:
self._export_eta_nc = False
self._export_z_nc = False
self._export_rho = False
self._export_nuh = False
self._timeIntegrator = None
self._export = []
self._ncExport = []
......@@ -252,6 +266,7 @@ class Loop:
def export_vertical_diffusivity(self): self._export_kappav = True
def export_vertical_viscosity(self): self._export_nuv = True
def export_rho(self): self._export_rho = True
def export_nuh(self): self._export_nuh = True
def export_uv_nc(self): self._export_uv_nc = True
def export_salinity_nc(self): self._export_S_nc = True
def export_elevation_nc(self): self._export_eta_nc = True
......@@ -275,6 +290,7 @@ class Loop:
if self._export_uv: self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", self._export_uv_vect))
if self._export_w: self._export.append(dgpy.dgIdxExporter(d.wDof3d, self._odir + "/w"))
if self._export_rho: self._export.append(dgpy.dgIdxExporter(d.rhoDof3d, self._odir + "/rho"))
if self._export_nuh: self._export.append(dgpy.dgIdxExporter(d.nuhDof3d, self._odir + "/nuh"))
if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2d, self._odir + "/uvAv2d", self._export_uvAv2d_vect))
if self._export_eta: self._export.append(dgpy.dgIdxExporter(d.etaDof2d, self._odir + "/eta"))
if self._export_S:
......@@ -319,6 +335,8 @@ class Loop:
if self._restart_ind > -1:
timeIntegrator.reStart(self._restart_dir, self._restart_ind)
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
if self._export_nuh:
d.nuhDof3d.interpolate(f.nuh)
for e in self._export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
if self._ncWriter3d:
......@@ -338,6 +356,10 @@ class Loop:
def export_fields(self):
timeIntegrator = self._timeIntegrator
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
d = self._slimSolver.getDofs()
f = self._slimSolver.functions
if self._export_nuh:
d.nuhDof3d.interpolate(f.nuh)
for e in self._export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
if self._ncWriter3d:
......
dgpy/scripts/slim3d_private.py
View file @ 55d702d8
......@@ -91,9 +91,17 @@ def slim3d_setup(loop):
f.rhoFunc = dgpy.linearStateEquation(d.zDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
if eq._vertical_diffusivity == 'gotm':
dgpy.Msg.Fatal('GOTM does not work with linear density. S and T must be set, and the automatic Jackett-McDougall state equation will be computed')
#dgpy.Msg.Warning('GOTM does not work with linear density. What you are doing is weird..')
dgpy.slim3dParameters.setRho0(eq._reference_density)
eq._areaFunc = None
eq._nuh = None
def fun_add(cmap, val, f1, f2) :
val[:] = f1[:] + f2[:]
if eq._additional_shear and eq._horizontal_viscosity != 'smagorinsky':
dgpy.Msg.Fatal('Additional shear is currently only implemented with Smagorinsky formulation')
if eq._horizontal_viscosity == 'constant':
f.nuh = dgpy.functionConstant(eq._hor_visc_const)
elif eq._horizontal_viscosity == 'smagorinsky':
......@@ -101,7 +109,11 @@ def slim3d_setup(loop):
slimSolver.copy2d3d.computeArea2d( area2d)
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)
if eq._additional_shear:
eq._add_shear = slim_private._load_function(eq._additional_shear, slimSolver.groups3d)
else:
eq._add_shear = None
f.nuh = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof.getFunction(), eq._areaFunc, eq._add_shear)
if eq._horizontal_diffusivity == 'constant':
f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const)
......@@ -244,7 +256,7 @@ def slim3d_setup(loop):
horMomEq = e.horMomEq
horMomEq.setLaxFriedrichsFactor(0.0)
#horMomEq.setLaxFriedrichsFactor(eq._LFF)
if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc :
#horMomEq.setBoundary0Flux(topTags) # TODO CHECK WHICH !!!
horMomEq.setBoundarySymmetry(topTags)
......@@ -272,6 +284,7 @@ def slim3d_setup(loop):
e.vertMomUEq.setBoundary0Flux(openBnd.tag)
wEq = e.wEq
wEq.setLaxFriedrichsFactor(eq._LFF)
wEq.setBoundarySymmetry(topTags)
wEq.setBoundary0Flux(bottomTags)
wEq.setBoundary0Flux(eq._boundary_coast)
......@@ -292,7 +305,7 @@ def slim3d_setup(loop):
#e.newRGradEq.addBoundaryCondition(['open', 'openSouth'], rGradBndOut)
if slimSolver.getSolveS() :
e.SEq.setLaxFriedrichsFactor(0.0)
e.SEq.setLaxFriedrichsFactor(eq._LFF)
e.SEq.setBoundary0Flux(topTags)
e.SEq.setBoundary0Flux(bottomTags)
e.SEq.setBoundary0Flux(eq._boundary_coast)
......@@ -304,7 +317,7 @@ def slim3d_setup(loop):
e.SEq.addBoundaryCondition(openBnd.tag, SBndOpen)
if slimSolver.getSolveT() :
e.TEq.setLaxFriedrichsFactor(0.0)
e.TEq.setLaxFriedrichsFactor(eq._LFF)
e.TEq.setBoundary0Flux(topTags)
e.TEq.setBoundary0Flux(bottomTags)
e.TEq.setBoundary0Flux(eq._boundary_coast)
......@@ -317,7 +330,7 @@ def slim3d_setup(loop):
if slimSolver.getAdvectTurbulentEnergy() :
for turbEq in [e.tkeAdvEq, e.epsAdvEq] :
turbEq.setLaxFriedrichsFactor(0.0)
turbEq.setLaxFriedrichsFactor(eq._LFF)
turbEq.setBoundary0Flux(topTags)
turbEq.setBoundary0Flux(bottomTags)
turbEq.setBoundary0Flux(eq._boundary_coast)
......@@ -326,6 +339,9 @@ def slim3d_setup(loop):
turbEq.addBoundaryCondition(openBnd.tag, turbEq.newInFluxBoundary(f.tinyFunc))
eta2dEq = e.eta2dEq
eta2dEq.setIsLinear(eq._linear2d)
if eq._LFF > 1e-8:
eta2dEq.setLaxFriedrichs(True)
eta2dEq.setBoundary0Flux(eq._boundary_coast)
for openBnd in eq._boundary_open:
toReplaceEtaOpen = dgpy.VectorFunctorConst([uvAv2d, dgpy.function.getSolution()])
......@@ -334,6 +350,9 @@ def slim3d_setup(loop):
eta2dEq.addBoundaryCondition(openBnd.tag, etaBndOpen)
uv2dEq = e.uv2dEq
uv2dEq.setIsLinear(eq._linear2d)
if eq._LFF > 1e-8:
uv2dEq.setLaxFriedrichs(True)
uv2dBndWall = uv2dEq.newBoundaryWall()
uv2dEq.addBoundaryCondition(eq._boundary_coast,uv2dBndWall)
for openBnd in eq._boundary_open:
......
modules/slim3d/dgConservationLawSW3dContinuity.cpp
View file @ 55d702d8
......@@ -36,6 +36,8 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
const fullMatrix<double> &etaOut = cache.get(*_eta, 1);
const fullMatrix<double> &UIn = cache.get(*_uv, 0);
const fullMatrix<double> &UOut = cache.get(*_uv, 1);
const fullMatrix<double> &wOldIn = cache.get(*_wOld, 0);
const fullMatrix<double> &wOldOut = cache.get(*_wOld, 1);
const fullMatrix<double> &bath = cache.get(*_bathymetryMinCG, 0);
double g = slim3dParameters::g;
// determine which element is below
......@@ -47,28 +49,24 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
double unIn = nx*UIn(i,0) + ny*UIn(i,1);
double unOut = nx*UOut(i,0) + ny*UOut(i,1);
double un = (unIn+unOut)/2;
double wn = 0;
if (fabs(nz) < 1e-8) { // vertical face
//linear equations
double hMin = bath(i,0);
double etaM = 0.5*(etaIn(i,0)+etaOut(i,0));
double sq_g_h = sqrt(g/(hMin+etaM));
un = (unIn+unOut + (etaIn(i,0)-etaOut(i,0))*sq_g_h)/2;
//un = (unIn+unOut + 1*(_wOldIn(i,0) - _wOldOut(i,0)))/2;
//if (fabs(_wOldIn(i,0) - _wOldOut(i,0)) >1e-6)
// printf("erreur!!!!+grand que 0 %g\n", _wOldIn(i,0) - _wOldOut(i,0));
//printf("win, wout :%g %g\n", _wIn(i,0), _wOut(i,0));
val(i,0) = - (un + wn);
val(i,1) = -val(i,0);
} else { // horizontal face, top/bottom
if (fabs(nz) < 1e-8) { // lateral face
if ( _laxFriedrichsFactor > 1e-8 )
un = ( unIn+unOut + _laxFriedrichsFactor* (wOldIn(i,0) - wOldOut(i,0)) )/2;
else{
double hMin = bath(i,0);
double etaM = 0.5*(etaIn(i,0)+etaOut(i,0));
double sq_g_h = sqrt(g/(hMin+etaM));
un = (unIn+unOut + (etaIn(i,0)-etaOut(i,0))*sq_g_h)/2;
}
} else { // top-bottom face
if (outIsTop == _propagateFromTop) {
un = unOut;
} else {
un = unIn;
}
val(i,0) = -un;
val(i,1) = -val(i,0);
}
val(i,0) = -un;
val(i,1) = -val(i,0);
}
}
......@@ -76,7 +74,7 @@ dgConservationLawSW3dContinuity::dgConservationLawSW3dContinuity( int propagateF
dgConservationLawFunction(1), _propagateFromTop(propagateFromTop){
_uv = NULL;
_eta = NULL;
_bathymetry = NULL;
_wOld = NULL;
}
void dgConservationLawSW3dContinuity::setup() {
_volumeTerm0[""] = NULL;
......@@ -86,6 +84,8 @@ void dgConservationLawSW3dContinuity::setup() {
Msg::Fatal("dgConservationLawSW3dContinuity: uv function must be set");
if (!_eta)
Msg::Fatal("dgConservationLawSW3dContinuity: eta function must be set");
if (!_wOld)
Msg::Fatal("dgConservationLawSW3dContinuity: wOld function must be set");
}
dgConservationLawSW3dContinuity::~dgConservationLawSW3dContinuity() {
if (_volumeTerm0[""]) delete _volumeTerm0[""];
......
modules/slim3d/dgConservationLawSW3dContinuity.h
View file @ 55d702d8
......@@ -12,8 +12,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
// u3d is horizontal 3d velocity
// the 2d fields have been projected on 3d container for easier access
// they are thus constants in z
const functor *_uv, *_eta, *_bathymetry, *_bathymetryMinCG;// *_w;
const functor *_uv, *_eta, *_bathymetryMinCG, *_wOld;
bool _propagateFromTop;
double _laxFriedrichsFactor;
void riemann(functorCache &cache, fullMatrix<double> &val) const;
void advection(functorCache &cache, fullMatrix<double> &val) const;
void boundaryWall(functorCache &cache, fullMatrix<double> &val) const;
......@@ -25,8 +26,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
dgBoundaryCondition *newBoundaryWall();
inline void setUV(const functor *f) { _uv = f; };
inline void setEta(const functor *f) { _eta = f; };
inline void setBathymetry(const functor *f) { _bathymetry = f; };
inline void setWOld(const functor *f) { _wOld = f; };
inline void setBathymetryMinCG(const functor *bathymetryMinCG) { _bathymetryMinCG = bathymetryMinCG;}
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
void finalize(const dgGroupOfElements &group, dgDofContainer &residual) const;
//inline void setW(const functor *f) { _w = f; };
};
......
modules/slim3d/dgConservationLawSW3dMomentum.cpp
View file @ 55d702d8
......@@ -160,16 +160,13 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
double ny = normals(i,1);
double nz = normals(i,2);
double hMin = bath(i,0);
double un;//,vn; // in normal,tangential coordinates
double u,v; // in x,y,z coordinates
double un, u, v;
double u3dL = UDevL(i,0) + U2dL(i,0);
double v3dL = UDevL(i,1) + U2dL(i,1);
double u3dR = UDevR(i,0) + U2dR(i,0);
double v3dR = UDevR(i,1) + U2dR(i,1);
double unL = nx*u3dL + ny*v3dL, unR = nx*u3dR + ny*v3dR;
//double vnL = -ny*_UL(i,0) + nx*_UL(i,1), vnR = - ny*_UR(i,0) + nx*_UR(i,1);
double un2dL = nx*U2dL(i,0) + ny*U2dL(i,1), un2dR = nx*U2dR(i,0) + ny*U2dR(i,1);
//double vn2dL = -ny*_U2dL(i,0) + nx*_U2dL(i,1), vn2dR = - ny*_U2dR(i,0) + nx*_U2dR(i,1);
double eta;
double rhoSurf = 0.5*(rhoSurfL(i,0) + rhoSurfR(i,0));
double wL_ = wL(i,0), wR_ = wR(i,0);
......@@ -178,14 +175,18 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
wR_ += -wMeshR(i,0);
}
if (fabs(nz)<1e-8) { // lateral face
if ( _fullLaxFriedrichs ){
if ( _laxFriedrichsFactor > 1e-8 ){
eta = 0.5*(etaL(i,0)+etaR(i,0));
// Lax friedrichs
double fluxum = -(unL * u3dL + unR * u3dR)/2 - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *nx;
double fluxvm = -(unL * v3dL + unR * v3dR)/2 - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *ny;
double c = sqrt(g*(hMin+eta))/10;
val(i,0) = fluxum - (UDevL(i,0) - UDevR(i,0))/2 * c;
val(i,1) = fluxvm - (UDevL(i,1) - UDevR(i,1))/2 * c;
double fluxum = -((unL+unR)>0 ? unL*u3dL : unR * u3dR) - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *nx;
double fluxvm = -((unL+unR)>0 ? unL*v3dL : unR * v3dR) - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *ny;
val(i,0) = fluxum;
val(i,1) = fluxvm;
//double fluxum = -(unL * u3dL + unR * u3dR)/2 - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *nx;
//double fluxvm = -(unL * v3dL + unR * v3dR)/2 - g*( (rhoSurfL(i,0)/rho0)*etaL(i,0) + (rhoSurfR(i,0)/rho0)*etaR(i,0) )/2 *ny;
//double c = _laxFriedrichsFactor;// sqrt(g*(hMin+eta))/10;
//val(i,0) = fluxum - (UDevL(i,0) - UDevR(i,0))/2 * c;
//val(i,1) = fluxvm - (UDevL(i,1) - UDevR(i,1))/2 * c;
if (ipTerm) {
val(i,0) += (*ipTerm)(i,0);
val(i,1) += (*ipTerm)(i,1);
......@@ -196,34 +197,11 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
// should be symmetric for lateral bnd conditions
double etaM = 0.5*(etaL(i,0)+etaR(i,0));
double sqgh = sqrt(g/(hMin+etaM));
// The commented lines below are not necessary, since etaR, u2dR and uR are correctly set up in the boundaryWall bnd condition
eta = etaM + ((un2dL-un2dR)/sqgh)/2;
un = 0.5*(unL+unR) + 0.5*(etaL(i,0)-etaR(i,0))*sqgh;
//if (_xyz(i,2) > -hMin) {
// eta = etaM + ((un2dL-un2dR)/sqgh/(hMin+etaM))/2;
// un = 0.5*(unL+unR) + 0.5*(_etaL(i,0)-_etaR(i,0))*sqgh;
//}
//else {
// eta = etaM;
// un = 0;
//}
//vn = un > 0 ? vnL : vnR;
#if 0
un = 0;
vn = 0;
approxRoe(_etaL(i,0), _etaR(i,0), un2dL/hL, un2dR/hR, vn2dL/hL, vn2dR/hR, hL, hR, eta, un, vn);
#endif
}
if ( un > 0 ) { // upwinding
u = u3dL;
v = v3dL;
} else {
u = u3dR;
v = v3dR;
}
} else {
} else { // top-bottom face
// down value (Aizinger & Dawson)
if (outIsTop) {
un = nx*u3dL + ny*v3dL + nz*wL_;
......@@ -232,6 +210,10 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
un = nx*u3dR + ny*v3dR + nz*wR_;
eta = etaR(i,0);
}
}
val(i,0) = -g*(rhoSurf/rho0)*eta*nx;
val(i,1) = -g*(rhoSurf/rho0)*eta*ny;
if (!_withoutAdvection){
if ( un > 0 ) { // upwinding
u = u3dL;
v = v3dL;
......@@ -239,22 +221,9 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
u = u3dR;
v = v3dR;
}
}
val(i,0) = -g*(rhoSurf/rho0)*eta*nx;
val(i,1) = -g*(rhoSurf/rho0)*eta*ny;
if (!_withoutAdvection){
val(i,0) += -u*un;
val(i,1) += -v*un;
}
// Lax friedrichs
if (_laxFriedrichsFactor) {
//double uvwnL = nx*_UL(i,0) + ny*_UL(i,1) + nz*wL_;
//double uvwnR = nx*_UR(i,0) + ny*_UR(i,1) + nz*wR_;
double gamma = 1;//0.5 * ( fabs(uvwnL) + fabs(uvwnR) );
double uJ = (UDevL(i,0)-UDevR(i,0))/2, vJ = (UDevL(i,1)-UDevR(i,1))/2;
val(i,0) += -_laxFriedrichsFactor*gamma*uJ;
val(i,1) += -_laxFriedrichsFactor*gamma*vJ;
}
if (ipTerm) {
val(i,0) += (*ipTerm)(i,0);
val(i,1) += (*ipTerm)(i,1);
......@@ -308,7 +277,6 @@ dgConservationLawSW3dMomentum::dgConservationLawSW3dMomentum() : dgConservationL
_rhoSurf = _fzero;
_maxDiffusivity = NULL;
_laxFriedrichsFactor = 0.0;
_fullLaxFriedrichs = false;
_withoutAdvection = false;
}
......@@ -1076,7 +1044,7 @@ void dgConservationLawSW2dU::riemann(functorCache &cm, fullMatrix<double> &val)
const double c = sqrt(g * (hMin + etaM));
Fun = - g * (1+rhoSurf/rho0) * etaM -pa/rho0 + c * (uR - uL) / 2;
// TODO needed ? or 0 is ok ?
Fut = c * (vR - vL) / 2;
Fut = 0;// c * (vR - vL) / 2;
}
else{
linearRoe(etaL(i,0), etaR(i,0), uL, uR, h, eta, u);
......
modules/slim3d/dgConservationLawSW3dMomentum.h
View file @ 55d702d8
......@@ -16,7 +16,6 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
const functor *_fzero, *_fzerov;
const functor *_nuH, *_nuV, *_diffusiveFlux, *_ipTerm, *_maxDiffusivity;
double _laxFriedrichsFactor;
bool _fullLaxFriedrichs;
bool _useConservativeALE;
bool _withoutAdvection;
class boundaryWallArg;
......@@ -31,7 +30,6 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
* Larger value requires smaller explicit time step : CFL ~= 1/LFF *
* Value LFF=1 corresponds to Aizinger and Dawson CMAME 2007. */
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
inline void setFullLaxFriedrichs(bool b) {_fullLaxFriedrichs = b; };
void setUseConservativeALE(bool b) { _useConservativeALE = b; };
/**set eta function [eta] */
inline void setEta(const functor *eta){ _eta = eta;}
......
modules/slim3d/dgConservationLawSW3dTracer.cpp
View file @ 55d702d8
......@@ -136,6 +136,8 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
const fullMatrix<double> *uvR = _uv ? &cm.get(*_uv, 1) : NULL;
const fullMatrix<double> *wL = _w ? &cm.get(*_w, 0) : NULL;
const fullMatrix<double> *wR = _w ? &cm.get(*_w, 1) : NULL;
const fullMatrix<double> *wOldL = _wOld ? &cm.get(*_wOld, 0) : NULL;
const fullMatrix<double> *wOldR = _wOld ? &cm.get(*_wOld, 1) : NULL;
const fullMatrix<double> *wMeshL = _wM ? &cm.get(*_wM, 0) : NULL;
const fullMatrix<double> *wMeshR = _wM ? &cm.get(*_wM, 1) : NULL;
const fullMatrix<double> *uGML = _uGM ? &cm.get(*_uGM, 0) : NULL;
......@@ -158,12 +160,12 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
wR_ += (*wR)(i,0);
}
if(_uGM) {
uL += (*uGML)(i,0);
vL += (*uGML)(i,1);
wL_ += (*uGML)(i,2);
uR += (*uGMR)(i,0);
vR += (*uGMR)(i,1);
wR_ += (*uGMR)(i,2);
uL += (*uGML)(i,0);
vL += (*uGML)(i,1);
wL_ += (*uGML)(i,2);
uR += (*uGMR)(i,0);
vR += (*uGMR)(i,1);
wR_ += (*uGMR)(i,2);
}
if(_uv || _wM) {
double nx = normals(i,0);
......@@ -177,12 +179,14 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
if(_uv) {
double unL = nx*uL + ny*vL;
double unR = nx*uR + ny*vR;
un = (unL+unR)/2;
if(fabs(nz) < 1e-8) { // vertical face
if(_fullLaxFriedrichs){
double fluxM = -(cL(i,0) * unL + cR(i,0) * unR)/2.;
double cvel = 2.;
val(i,0) = fluxM - (cL(i,0) - cR(i,0))/2 * cvel ;
if(fabs(nz) < 1e-8) { // lateral face
if ( _laxFriedrichsFactor > 1e-8 ){
un = (unL+unR + _laxFriedrichsFactor*((*wOldL)(i,0) - (*wOldR)(i,0)))/2;
//double fluxM = - un * ( cL(i,0) + cR(i,0) ) /2;
//double cvel = _laxFriedrichsFactorTracer;
//val(i,0) = fluxM - (cL(i,0) - cR(i,0))/2 * cvel ;
double c = un>0 ? cL(i,0) : cR(i,0);
val(i,0) = - un * c;
if(ipTerm)
val(i,0) += (*ipTerm)(i,0);
continue;
......@@ -194,33 +198,22 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
double sq_g_h = sqrt(g/(hMin+etaM));
un = (unL+unR + (etaL(i,0)-etaR(i,0))*sq_g_h)/2;
}
} else { // horizontal face, top/bottom
} else { // top-bottom face
if (rightIsBelow)
un = unR + nz*wR_;
else
un = unL + nz*wL_;
}
}
else {
else
un = 0.;
}
double c = un>0 ? cL(i,0) : cR(i,0);
val(i,0) = -c*un;
if(_laxFriedrichsFactor) {
// Lax friedrichs
double uvwnL = nx*uL + ny*vL + nz*wL_;
double uvwnR = nx*uR + ny*vR + nz*wR_;
double gamma = 0.5*( fabs(uvwnL) + fabs(uvwnR) );
double cJ = (cL(i,0) - cR(i,0))/2;
val(i,0) += -_laxFriedrichsFactor*gamma*cJ;
}
}
if(ipTerm) {
if(ipTerm)
val(i,0) += (*ipTerm)(i,0);
}
if(ipTermIso){
if(ipTermIso)
val(i,0) += (*ipTermIso)(i,0);
}
}
}
......@@ -469,7 +462,6 @@ dgConservationLawSW3dTracer::dgConservationLawSW3dTracer(const functor *uv,const
_diffusiveFlux = NULL;
_ipTerm = NULL;
_useConservativeALE = false;
_fullLaxFriedrichs = false;
_uGM=NULL;
}
......@@ -498,7 +490,6 @@ dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawFu
_dt=0;
_T=1;
_useConservativeALE = false;
_fullLaxFriedrichs = false;
_uGM=NULL;
}
......
modules/slim3d/dgConservationLawSW3dTracer.h
View file @ 55d702d8
......@@ -7,10 +7,9 @@ class function;
* \f[\frac{dc}{dt} + v. \nabla c - nu \Delta c = 0\f]
*/
class dgConservationLawSW3dTracer : public dgConservationLawFunction {
const functor *_uv, *_w, *_wM, *_eta, *_bathymetryMinCG, *_dwMdz, *_kappaH, *_kappaV, *_diffusiveFlux, *_ipTerm, *_slope, *_Kedd, *_isopycnalDiffusion, *_interfaceIso, *_ipTermIso, *_sourceFunc, *_uGM;
const functor *_uv, *_w, *_wOld, *_wM, *_eta, *_bathymetryMinCG, *_dwMdz, *_kappaH, *_kappaV, *_diffusiveFlux, *_ipTerm, *_slope, *_Kedd, *_isopycnalDiffusion, *_interfaceIso, *_ipTermIso, *_sourceFunc, *_uGM;
const functor *_fzero, *_fzeroVec;
double _laxFriedrichsFactor;
bool _fullLaxFriedrichs;
int _dt,_T;
const functor *_ref;
bool _useConservativeALE;
......@@ -37,12 +36,12 @@ class dgConservationLawSW3dTracer : public dgConservationLawFunction {
dgConservationLawSW3dTracer();
/** set factor for scaling Lax-Friedrichs penalty term [default LFF = 0]. */
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;}
inline void setFullLaxFriedrichs(bool b) {_fullLaxFriedrichs = b; };
inline void setUseConservativeALE(bool b) {_useConservativeALE = b;};
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;};
......
modules/slim3d/slim3dEquations.cpp
View file @ 55d702d8
......@@ -12,7 +12,7 @@ slim3dSolverDofs::slim3dSolverDofs(slim3dSolver *solver): _solver(solver) {
uvIntDof2d = uvCorrDof2d = uvCorrDof3d = uvBotDof2d = NULL;
windStressDof2dCenter = NULL;
uvTauBDof2d = uvTauBDof2dCenter = NULL;
wDof3d = wSurfDof2d = NULL;
wDof3d = wOldDof3d = wSurfDof2d = NULL;
//wMeshDof3d = wMeshSurfDof2d = NULL;
//dwMeshDzDof3d = NULL;
etaDof2d = etaDof2dCG = etaDof3dCenter = NULL;
......@@ -29,6 +29,7 @@ slim3dSolverDofs::slim3dSolverDofs(slim3dSolver *solver): _solver(solver) {
bathDof2d = bathCGMinDof2d = NULL;
arbitraryMovingMeshDof3d = kappaZEqDof3d = NULL;
vJumpDof3d = errorDof3d = zStarDof3d = oneOverErrorDof3d = oneOverErrorIntDof2d = NULL;
nuhDof3d = NULL;
testDof3d = testDof2d = NULL;
_initialized = false;
}
......@@ -44,6 +45,7 @@ void slim3dSolverDofs::initialize()
uvDevDof = new dgDofContainer(*groups3d, 2);
uvHorDivDof = new dgDofContainer(*groups3d, 1);
wDof3d = new dgDofContainer(*groups3d, 1);
wOldDof3d = new dgDofContainer(*groups3d, 1);
//wDof3dCopy = new dgDofContainer(*groups3d, 1);
uvAvDof2d = new dgDofContainer(*groups2d, 2);
uvIntDof2d = new dgDofContainer(*groups2d, 2);
......@@ -149,6 +151,7 @@ void slim3dSolverDofs::initialize()
oneOverErrorDof3d = new dgDofContainer(*groups3d, 1);
oneOverErrorIntDof2d = new dgDofContainer(*groups2d, 1);
}
nuhDof3d = new dgDofContainer(*groups3d, 1);
testDof3d = new dgDofContainer(*groups3d, 1);
testDof2d = new dgDofContainer(*groups2d, 1);
_initialized = true;
......@@ -160,6 +163,7 @@ slim3dSolverDofs::~slim3dSolverDofs()
if ( uvDevDof ) delete uvDevDof;
if ( uvHorDivDof ) delete uvHorDivDof;
if ( wDof3d ) delete wDof3d;
if ( wOldDof3d ) delete wOldDof3d;
if ( uvAvDof2d ) delete uvAvDof2d;
if ( uvIntDof2d ) delete uvIntDof2d;
if ( uvCorrDof2d ) delete uvCorrDof2d;
......@@ -226,6 +230,7 @@ slim3dSolverDofs::~slim3dSolverDofs()
if ( zStarDof3d ) delete zStarDof3d;
if ( oneOverErrorDof3d ) delete oneOverErrorDof3d;
if ( oneOverErrorIntDof2d )delete oneOverErrorIntDof2d;
if ( nuhDof3d ) delete nuhDof3d;
if ( testDof3d) delete testDof3d;
if ( testDof2d) delete testDof2d;
......@@ -279,6 +284,7 @@ slim3dSolverFunctions::slim3dSolverFunctions(slim3dSolver *solver) : _solver(sol
wMeshSurf2dPrecompFunc = NULL;