Commit 5a65f5e8 authored by Philippe Delandmeter's avatar Philippe Delandmeter Committed by Jonathan Lambrechts
Browse files

Lax-Friedrichs for the new equations

parent 087a703b
...@@ -54,6 +54,11 @@ class Slim3d_equations : ...@@ -54,6 +54,11 @@ class Slim3d_equations :
self._z0 = [0.005, 0.02] self._z0 = [0.005, 0.02]
self._linear2d = False self._linear2d = False
self._linear_solver_2D = True self._linear_solver_2D = True
self._LFF = 0.
self._additional_shear = None
self._areaFunc = None
self._nuh = None
self._add_shear = None
def set_implicit_vertical_diffusion(self, flag): def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver slimSolver = self._slimSolver
...@@ -80,6 +85,9 @@ class Slim3d_equations : ...@@ -80,6 +85,9 @@ class Slim3d_equations :
if slimSolver.getSolveT(): if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag) slimSolver.setFlagTLimiter(flag)
def set_lax_friedrichs(self, factor=1.):
self._LFF = factor
def set_conservative_ale(self, flag): def set_conservative_ale(self, flag):
self._slimSolver.setUseConservativeALE(flag) self._slimSolver.setUseConservativeALE(flag)
...@@ -91,6 +99,10 @@ class Slim3d_equations : ...@@ -91,6 +99,10 @@ class Slim3d_equations :
self._hor_visc_fact = factor self._hor_visc_fact = factor
self._hor_visc_max = maximum 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): def set_horizontal_diffusivity(self, mode, constant=1, factor=0.03, maximum=1e9):
if (mode != "okubo") and (mode != "constant"): if (mode != "okubo") and (mode != "constant"):
dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode) dgpy.Msg.Fatal("Unknown mode for horizontal_viscosity : "+mode)
...@@ -239,6 +251,7 @@ class Loop: ...@@ -239,6 +251,7 @@ class Loop:
self._export_eta_nc = False self._export_eta_nc = False
self._export_z_nc = False self._export_z_nc = False
self._export_rho = False self._export_rho = False
self._export_nuh = False
self._timeIntegrator = None self._timeIntegrator = None
self._export = [] self._export = []
self._ncExport = [] self._ncExport = []
...@@ -259,6 +272,7 @@ class Loop: ...@@ -259,6 +272,7 @@ class Loop:
def export_vertical_diffusivity(self): self._export_kappav = True def export_vertical_diffusivity(self): self._export_kappav = True
def export_vertical_viscosity(self): self._export_nuv = True def export_vertical_viscosity(self): self._export_nuv = True
def export_rho(self): self._export_rho = 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_uv_nc(self): self._export_uv_nc = True
def export_salinity_nc(self): self._export_S_nc = True def export_salinity_nc(self): self._export_S_nc = True
def export_elevation_nc(self): self._export_eta_nc = True def export_elevation_nc(self): self._export_eta_nc = True
...@@ -282,6 +296,7 @@ class Loop: ...@@ -282,6 +296,7 @@ class Loop:
if self._export_uv: self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", self._export_uv_vect)) 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_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_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_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_eta: self._export.append(dgpy.dgIdxExporter(d.etaDof2d, self._odir + "/eta"))
if self._export_S: if self._export_S:
...@@ -326,6 +341,8 @@ class Loop: ...@@ -326,6 +341,8 @@ class Loop:
if self._restart_ind > -1: if self._restart_ind > -1:
timeIntegrator.reStart(self._restart_dir, self._restart_ind) timeIntegrator.reStart(self._restart_dir, self._restart_ind)
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount()) if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
if self._export_nuh:
d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
for e in self._export : for e in self._export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time) e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
if self._ncWriter3d: if self._ncWriter3d:
...@@ -345,6 +362,9 @@ class Loop: ...@@ -345,6 +362,9 @@ class Loop:
def export_fields(self): def export_fields(self):
timeIntegrator = self._timeIntegrator timeIntegrator = self._timeIntegrator
if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount()) if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
d = self._slimSolver.getDofs()
if self._export_nuh:
d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
for e in self._export : for e in self._export :
e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time) e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
if self._ncWriter3d: if self._ncWriter3d:
...@@ -367,8 +387,6 @@ class Loop: ...@@ -367,8 +387,6 @@ class Loop:
def terminate(self): def terminate(self):
self._timeIntegrator.terminate() self._timeIntegrator.terminate()
def loop(self): def loop(self):
if not self._timeIntegrator: if not self._timeIntegrator:
self.setup() self.setup()
......
...@@ -91,11 +91,14 @@ def slim3d_setup(loop): ...@@ -91,11 +91,14 @@ def slim3d_setup(loop):
f.rhoFunc = dgpy.linearStateEquation(d.zDof.getFunction(), eq._linear_density[2], eq._linear_density[1]) f.rhoFunc = dgpy.linearStateEquation(d.zDof.getFunction(), eq._linear_density[2], eq._linear_density[1])
if eq._vertical_diffusivity == 'gotm': 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.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) dgpy.slim3dParameters.setRho0(eq._reference_density)
dgpy.slim3dParameters.setG(eq._gravity) dgpy.slim3dParameters.setG(eq._gravity)
eq._gravityFunc = dgpy.functionConstant(eq._gravity) eq._gravityFunc = dgpy.functionConstant(eq._gravity)
eq._areaFunc = None 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': if eq._horizontal_viscosity == 'constant':
f.nuh = dgpy.functionConstant(eq._hor_visc_const) f.nuh = dgpy.functionConstant(eq._hor_visc_const)
elif eq._horizontal_viscosity == 'smagorinsky': elif eq._horizontal_viscosity == 'smagorinsky':
...@@ -103,7 +106,9 @@ def slim3d_setup(loop): ...@@ -103,7 +106,9 @@ def slim3d_setup(loop):
slimSolver.copy2d3d.computeArea2d( area2d) slimSolver.copy2d3d.computeArea2d( area2d)
area2d.scatter() area2d.scatter()
eq._areaFunc = dgpy.functionPrecomputedExtrusion(slimSolver.extrusion(), 3, area2d.getFunction()) 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)
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': if eq._horizontal_diffusivity == 'constant':
f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const) f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const)
...@@ -262,7 +267,7 @@ def slim3d_setup(loop): ...@@ -262,7 +267,7 @@ def slim3d_setup(loop):
sw2DEq.addBoundaryCondition(openBnd.tag, sw2DEq.newOutsideValueBoundary(eq._sw2D_open[openBnd])) sw2DEq.addBoundaryCondition(openBnd.tag, sw2DEq.newOutsideValueBoundary(eq._sw2D_open[openBnd]))
horMomEq = e.horMomEq horMomEq = e.horMomEq
horMomEq.setLaxFriedrichsFactor(0.0) #horMomEq.setLaxFriedrichsFactor(eq._LFF)
if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc : if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc :
#horMomEq.setBoundary0Flux(topTags) # TODO CHECK WHICH !!! #horMomEq.setBoundary0Flux(topTags) # TODO CHECK WHICH !!!
horMomEq.setBoundarySymmetry(topTags) horMomEq.setBoundarySymmetry(topTags)
...@@ -290,6 +295,7 @@ def slim3d_setup(loop): ...@@ -290,6 +295,7 @@ def slim3d_setup(loop):
e.vertMomUEq.setBoundary0Flux(openBnd.tag) e.vertMomUEq.setBoundary0Flux(openBnd.tag)
wEq = e.wEq wEq = e.wEq
wEq.setLaxFriedrichsFactor(eq._LFF)
wEq.setBoundarySymmetry(topTags) wEq.setBoundarySymmetry(topTags)
wEq.setBoundary0Flux(bottomTags) wEq.setBoundary0Flux(bottomTags)
wEq.setBoundary0Flux(eq._boundary_coast) wEq.setBoundary0Flux(eq._boundary_coast)
...@@ -310,7 +316,7 @@ def slim3d_setup(loop): ...@@ -310,7 +316,7 @@ def slim3d_setup(loop):
#e.newRGradEq.addBoundaryCondition(['open', 'openSouth'], rGradBndOut) #e.newRGradEq.addBoundaryCondition(['open', 'openSouth'], rGradBndOut)
if slimSolver.getSolveS() : if slimSolver.getSolveS() :
e.SEq.setLaxFriedrichsFactor(0.0) e.SEq.setLaxFriedrichsFactor(eq._LFF)
e.SEq.setBoundary0Flux(topTags) e.SEq.setBoundary0Flux(topTags)
e.SEq.setBoundary0Flux(bottomTags) e.SEq.setBoundary0Flux(bottomTags)
e.SEq.setBoundary0Flux(eq._boundary_coast) e.SEq.setBoundary0Flux(eq._boundary_coast)
...@@ -322,7 +328,7 @@ def slim3d_setup(loop): ...@@ -322,7 +328,7 @@ def slim3d_setup(loop):
e.SEq.addBoundaryCondition(openBnd.tag, SBndOpen) e.SEq.addBoundaryCondition(openBnd.tag, SBndOpen)
if slimSolver.getSolveT() : if slimSolver.getSolveT() :
e.TEq.setLaxFriedrichsFactor(0.0) e.TEq.setLaxFriedrichsFactor(eq._LFF)
e.TEq.setBoundary0Flux(topTags) e.TEq.setBoundary0Flux(topTags)
e.TEq.setBoundary0Flux(bottomTags) e.TEq.setBoundary0Flux(bottomTags)
e.TEq.setBoundary0Flux(eq._boundary_coast) e.TEq.setBoundary0Flux(eq._boundary_coast)
...@@ -335,7 +341,7 @@ def slim3d_setup(loop): ...@@ -335,7 +341,7 @@ def slim3d_setup(loop):
if slimSolver.getAdvectTurbulentEnergy() : if slimSolver.getAdvectTurbulentEnergy() :
for turbEq in [e.tkeAdvEq, e.epsAdvEq] : for turbEq in [e.tkeAdvEq, e.epsAdvEq] :
turbEq.setLaxFriedrichsFactor(0.0) turbEq.setLaxFriedrichsFactor(eq._LFF)
turbEq.setBoundary0Flux(topTags) turbEq.setBoundary0Flux(topTags)
turbEq.setBoundary0Flux(bottomTags) turbEq.setBoundary0Flux(bottomTags)
turbEq.setBoundary0Flux(eq._boundary_coast) turbEq.setBoundary0Flux(eq._boundary_coast)
......
...@@ -36,6 +36,8 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do ...@@ -36,6 +36,8 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
const fullMatrix<double> &etaOut = cache.get(*_eta, 1); const fullMatrix<double> &etaOut = cache.get(*_eta, 1);
const fullMatrix<double> &UIn = cache.get(*_uv, 0); const fullMatrix<double> &UIn = cache.get(*_uv, 0);
const fullMatrix<double> &UOut = cache.get(*_uv, 1); 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); const fullMatrix<double> &bath = cache.get(*_bathymetryMinCG, 0);
double g = slim3dParameters::g; double g = slim3dParameters::g;
// determine which element is below // determine which element is below
...@@ -47,28 +49,24 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do ...@@ -47,28 +49,24 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
double unIn = nx*UIn(i,0) + ny*UIn(i,1); double unIn = nx*UIn(i,0) + ny*UIn(i,1);
double unOut = nx*UOut(i,0) + ny*UOut(i,1); double unOut = nx*UOut(i,0) + ny*UOut(i,1);
double un = (unIn+unOut)/2; double un = (unIn+unOut)/2;
double wn = 0; if (fabs(nz) < 1e-8) { // lateral face
if (fabs(nz) < 1e-8) { // vertical face if ( _laxFriedrichsFactor > 1e-8 )
//linear equations un = ( unIn+unOut + _laxFriedrichsFactor* (wOldIn(i,0) - wOldOut(i,0)) )/2;
double hMin = bath(i,0); else{
double etaM = 0.5*(etaIn(i,0)+etaOut(i,0)); double hMin = bath(i,0);
double sq_g_h = sqrt(g/(hMin+etaM)); double etaM = 0.5*(etaIn(i,0)+etaOut(i,0));
un = (unIn+unOut + (etaIn(i,0)-etaOut(i,0))*sq_g_h)/2; double sq_g_h = sqrt(g/(hMin+etaM));
//un = (unIn+unOut + 1*(_wOldIn(i,0) - _wOldOut(i,0)))/2; un = (unIn+unOut + (etaIn(i,0)-etaOut(i,0))*sq_g_h)/2;
//if (fabs(_wOldIn(i,0) - _wOldOut(i,0)) >1e-6) }
// printf("erreur!!!!+grand que 0 %g\n", _wOldIn(i,0) - _wOldOut(i,0)); } else { // top-bottom face
//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 (outIsTop == _propagateFromTop) { if (outIsTop == _propagateFromTop) {
un = unOut; un = unOut;
} else { } else {
un = unIn; 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 ...@@ -76,7 +74,7 @@ dgConservationLawSW3dContinuity::dgConservationLawSW3dContinuity( int propagateF
dgConservationLawFunction(1), _propagateFromTop(propagateFromTop){ dgConservationLawFunction(1), _propagateFromTop(propagateFromTop){
_uv = NULL; _uv = NULL;
_eta = NULL; _eta = NULL;
_bathymetry = NULL; _wOld = NULL;
} }
void dgConservationLawSW3dContinuity::setup() { void dgConservationLawSW3dContinuity::setup() {
_volumeTerm0[""] = NULL; _volumeTerm0[""] = NULL;
...@@ -86,6 +84,8 @@ void dgConservationLawSW3dContinuity::setup() { ...@@ -86,6 +84,8 @@ void dgConservationLawSW3dContinuity::setup() {
Msg::Fatal("dgConservationLawSW3dContinuity: uv function must be set"); Msg::Fatal("dgConservationLawSW3dContinuity: uv function must be set");
if (!_eta) if (!_eta)
Msg::Fatal("dgConservationLawSW3dContinuity: eta function must be set"); Msg::Fatal("dgConservationLawSW3dContinuity: eta function must be set");
if (!_wOld)
Msg::Fatal("dgConservationLawSW3dContinuity: wOld function must be set");
} }
dgConservationLawSW3dContinuity::~dgConservationLawSW3dContinuity() { dgConservationLawSW3dContinuity::~dgConservationLawSW3dContinuity() {
if (_volumeTerm0[""]) delete _volumeTerm0[""]; if (_volumeTerm0[""]) delete _volumeTerm0[""];
......
...@@ -12,8 +12,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction { ...@@ -12,8 +12,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
// u3d is horizontal 3d velocity // u3d is horizontal 3d velocity
// the 2d fields have been projected on 3d container for easier access // the 2d fields have been projected on 3d container for easier access
// they are thus constants in z // they are thus constants in z
const functor *_uv, *_eta, *_bathymetry, *_bathymetryMinCG;// *_w; const functor *_uv, *_eta, *_bathymetryMinCG, *_wOld;
bool _propagateFromTop; bool _propagateFromTop;
double _laxFriedrichsFactor;
void riemann(functorCache &cache, fullMatrix<double> &val) const; void riemann(functorCache &cache, fullMatrix<double> &val) const;
void advection(functorCache &cache, fullMatrix<double> &val) const; void advection(functorCache &cache, fullMatrix<double> &val) const;
void boundaryWall(functorCache &cache, fullMatrix<double> &val) const; void boundaryWall(functorCache &cache, fullMatrix<double> &val) const;
...@@ -25,8 +26,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction { ...@@ -25,8 +26,9 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
dgBoundaryCondition *newBoundaryWall(); dgBoundaryCondition *newBoundaryWall();
inline void setUV(const functor *f) { _uv = f; }; inline void setUV(const functor *f) { _uv = f; };
inline void setEta(const functor *f) { _eta = 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 setBathymetryMinCG(const functor *bathymetryMinCG) { _bathymetryMinCG = bathymetryMinCG;}
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
void finalize(const dgGroupOfElements &group, dgDofContainer &residual) const; void finalize(const dgGroupOfElements &group, dgDofContainer &residual) const;
//inline void setW(const functor *f) { _w = f; }; //inline void setW(const functor *f) { _w = f; };
}; };
......
...@@ -120,8 +120,6 @@ void dgConservationLawSW3dMomentum::source(functorCache &cm, fullMatrix<double> ...@@ -120,8 +120,6 @@ void dgConservationLawSW3dMomentum::source(functorCache &cm, fullMatrix<double>
} }
void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double> &val) const { void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &etaL = cm.get(*_eta, 0);
const fullMatrix<double> &etaR = cm.get(*_eta, 1);
const fullMatrix<double> &UDevL = cm.get(*function::getSolution(), 0); const fullMatrix<double> &UDevL = cm.get(*function::getSolution(), 0);
const fullMatrix<double> &UDevR = cm.get(*function::getSolution(), 1); const fullMatrix<double> &UDevR = cm.get(*function::getSolution(), 1);
const fullMatrix<double> &U2dL = cm.get(*_uv2d, 0); const fullMatrix<double> &U2dL = cm.get(*_uv2d, 0);
...@@ -130,111 +128,50 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double> ...@@ -130,111 +128,50 @@ void dgConservationLawSW3dMomentum::riemann(functorCache &cm, fullMatrix<double>
const fullMatrix<double> &wR = cm.get(*_w, 1); const fullMatrix<double> &wR = cm.get(*_w, 1);
const fullMatrix<double> &wMeshL = cm.get(*_wM, 0); const fullMatrix<double> &wMeshL = cm.get(*_wM, 0);
const fullMatrix<double> &wMeshR = cm.get(*_wM, 1); const fullMatrix<double> &wMeshR = cm.get(*_wM, 1);
const fullMatrix<double> &bath = cm.get(*_bathymetryMinCG);
const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0); const fullMatrix<double> &normals = cm.get(*function::getNormals(), 0);
const fullMatrix<double> *ipTerm =_ipTerm ? &cm.get(*_ipTerm) : NULL; const fullMatrix<double> *ipTerm =_ipTerm ? &cm.get(*_ipTerm) : NULL;
double g = slim3dParameters::g;
bool outIsTop = normals(0, 2) > 0; bool outIsTop = normals(0, 2) > 0;
val.resize(cm.nEvaluationPoint(), 2, false); val.resize(cm.nEvaluationPoint(), 2, false);
for(int i=0; i<cm.nEvaluationPoint(); i++) { for(int i=0; i<cm.nEvaluationPoint(); i++) {
double nx = normals(i,0); double nx = normals(i,0);
double ny = normals(i,1); double ny = normals(i,1);
double nz = normals(i,2); 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 u3dL = UDevL(i,0) + U2dL(i,0); double u3dL = UDevL(i,0) + U2dL(i,0);
double v3dL = UDevL(i,1) + U2dL(i,1); double v3dL = UDevL(i,1) + U2dL(i,1);
double u3dR = UDevR(i,0) + U2dR(i,0); double u3dR = UDevR(i,0) + U2dR(i,0);
double v3dR = UDevR(i,1) + U2dR(i,1); double v3dR = UDevR(i,1) + U2dR(i,1);
double unL = nx*u3dL + ny*v3dL, unR = nx*u3dR + ny*v3dR; 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 wL_ = wL(i,0), wR_ = wR(i,0); double wL_ = wL(i,0), wR_ = wR(i,0);
if (_wM) { if (_wM) {
wL_ += -wMeshL(i,0); wL_ += -wMeshL(i,0);
wR_ += -wMeshR(i,0); wR_ += -wMeshR(i,0);
} }
if (fabs(nz)<1e-8) { // lateral face double fluxu = 0., fluxv = 0.;
if ( _fullLaxFriedrichs ){ if (!_withoutAdvection){
eta = 0.5*(etaL(i,0)+etaR(i,0)); if (fabs(nz)<1e-8) { // lateral face
// Lax friedrichs fluxu = -((unL+unR)>0 ? unL*u3dL : unR * u3dR);
double fluxum = -(unL * u3dL + unR * u3dR)/2; fluxv = -((unL+unR)>0 ? unL*v3dL : unR * v3dR);
double fluxvm = -(unL * v3dL + unR * v3dR)/2; } else { // top-bottom face
double c = sqrt(g*(hMin+eta))/10; // down value (Aizinger & Dawson)
val(i,0) = fluxum - (UDevL(i,0) - UDevR(i,0))/2 * c; double un, u, v;
val(i,1) = fluxvm - (UDevL(i,1) - UDevR(i,1))/2 * c; if (outIsTop) {
if (ipTerm) { un = nx*u3dL + ny*v3dL + nz*wL_;
val(i,0) += (*ipTerm)(i,0); } else {
val(i,1) += (*ipTerm)(i,1); un = nx*u3dR + ny*v3dR + nz*wR_;
} }
continue; if ( un > 0 ) { // upwinding
} u = u3dL;
else { v = v3dL;
// should be symmetric for lateral bnd conditions } else {
double etaM = 0.5*(etaL(i,0)+etaR(i,0)); u = u3dR;
double sqgh = sqrt(g/(hMin+etaM)); v = v3dR;
}
// The commented lines below are not necessary, since etaR, u2dR and uR are correctly set up in the boundaryWall bnd condition fluxu = -u*un;
eta = etaM + ((un2dL-un2dR)/sqgh)/2; fluxv = -v*un;
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 {
// down value (Aizinger & Dawson)
if (outIsTop) {
un = nx*u3dL + ny*v3dL + nz*wL_;
eta = etaL(i,0);
} else {
un = nx*u3dR + ny*v3dR + nz*wR_;
eta = etaR(i,0);
}
if ( un > 0 ) { // upwinding
u = u3dL;
v = v3dL;
} else {
u = u3dR;
v = v3dR;
} }
} }
val(i,0) = 0; val(i,0) = fluxu;
val(i,1) = 0; val(i,1) = fluxv;
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) { if (ipTerm) {
val(i,0) += (*ipTerm)(i,0); val(i,0) += (*ipTerm)(i,0);
val(i,1) += (*ipTerm)(i,1); val(i,1) += (*ipTerm)(i,1);
...@@ -287,7 +224,6 @@ dgConservationLawSW3dMomentum::dgConservationLawSW3dMomentum() : dgConservationL ...@@ -287,7 +224,6 @@ dgConservationLawSW3dMomentum::dgConservationLawSW3dMomentum() : dgConservationL
_rGrad = _fzerov; _rGrad = _fzerov;
_maxDiffusivity = NULL; _maxDiffusivity = NULL;
_laxFriedrichsFactor = 0.0; _laxFriedrichsFactor = 0.0;
_fullLaxFriedrichs = false;
_withoutAdvection = false; _withoutAdvection = false;
} }
...@@ -903,37 +839,12 @@ void dgConservationLawSW2dEta::riemann(functorCache &cm, fullMatrix<double> &val ...@@ -903,37 +839,12 @@ void dgConservationLawSW2dEta::riemann(functorCache &cm, fullMatrix<double> &val
void dgConservationLawSW2dEta::boundaryWall(functorCache &cm, fullMatrix<double> &val) const { void dgConservationLawSW2dEta::boundaryWall(functorCache &cm, fullMatrix<double> &val) const {
val.resize(cm.nEvaluationPoint(), 1, true); //already set to zero in resize val.resize(cm.nEvaluationPoint(), 1, true); //already set to zero in resize
//for(int i=0; i<cm.nEvaluationPoint(); i++) {
// val(i,0) = 0.;
//}
} }
dgBoundaryCondition *dgConservationLawSW2dEta::newBoundaryWall(){ dgBoundaryCondition *dgConservationLawSW2dEta::newBoundaryWall(){
return newBoundaryConditionMember(*this, &dgConservationLawSW2dEta::boundaryWall); return newBoundaryConditionMember(*this, &dgConservationLawSW2dEta::boundaryWall);
} }
void dgConservationLawSW2dEta::boundaryOpenFree(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &normals = cm.get(*function::getNormals());
const fullMatrix<double> &eta = cm.get(*function::getSolution(),0);
const fullMatrix<double> &bath = cm.get(*_bathymetry,0);
fullMatrix<double> &externalUV = cm.set(*_UV,1);
externalUV.resize(cm.nEvaluationPoint(), 2, false);
cm.setDefaultToSide0(true);
double g = slim3dParameters::g;
for(int i=0; i<cm.nEvaluationPoint(); i++) {
double nx = normals