Commit f22f75c1 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts Committed by Valentin Vallaeys
Browse files

unified equation for uv

parent 6152d552
......@@ -296,30 +296,18 @@ def slim3d_setup(loop):
horMomEq = e.horMomEq
horMomEq.setLaxFriedrichsFactor(eq._LFF)
if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc :
horMomEq.setBoundary0Flux(topTags)
else :
if f.windStressFunc :
horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
else :
horMomEq.setBoundary0Flux(topTags)
for openBnd in eq._boundary_open:
horMomEq.setBoundaryNeumann(openBnd.tag, openBnd.horizontal_momentum_flux_f)
horMomBndWall = horMomEq.newBoundaryWall()
horMomBndWallVerticalBottom = horMomEq.newBoundaryWallVerticalBottom()
horMomEq.addBoundaryCondition(eq._boundary_coast, horMomBndWall)
horMomEq.addBoundaryCondition(bottomTags, horMomBndWall)
horMomEq.addBoundaryCondition('vertical_bottom', horMomBndWallVerticalBottom)
if slimSolver.getSolveUVImplicitVerticalDiffusion() :
if f.windStressFunc :
e.vertMomUEq.addBoundaryCondition(topTags, e.vertMomUEq.newBoundarySurface(f.windStressFunc))
else :
e.vertMomUEq.setBoundary0Flux(topTags)
e.vertMomUEq.addBoundaryCondition(bottomTags, e.vertMomUEq.newBoundaryBottom(f.bottomFriction2dPrecompFunc, f.zBotDist2dPrecompFunc))
e.vertMomUEq.setBoundary0Flux(eq._boundary_coast)
e.vertMomUEq.setBoundary0Flux('vertical_bottom')
for openBnd in eq._boundary_open:
e.vertMomUEq.setBoundary0Flux(openBnd.tag)
horMomEq.addBoundaryCondition('vertical_bottom', horMomBndWall)
horMomEq.addBoundaryCondition(bottomTags,horMomEq.newBoundaryBottom(f.z0BFunc, f.zBotDist2dPrecompFunc))
wEq = e.wEq
wEq.setBoundarySymmetry(topTags)
......
......@@ -2666,7 +2666,7 @@ bool dgSlopeLimiterUV::boundaryValues(const dgGroupOfFaces * faces, dgDofContain
}
}
dgSlopeLimiterUV::dgSlopeLimiterUV(dgConservationLaw* claw, dgGroupCollection* groups): dgLimiter(claw)
dgSlopeLimiterUV::dgSlopeLimiterUV(dgConservationLaw* claw, const dgGroupCollection* groups): dgLimiter(claw)
{
_groups = groups;
......
......@@ -83,12 +83,12 @@ public :
/** similar to original limiter. The difference is that for boundary values, the mean is computed on the interface */
class dgSlopeLimiterUV : public dgLimiter{
dgGroupCollection *_groups;
const dgGroupCollection *_groups;
dgDofContainer* _nodalVolume;
public :
virtual bool boundaryValues(const dgGroupOfFaces * faces, dgDofContainer * solution,
fullMatrix<double> & valueOnFaces, const functor * secondary=NULL);
dgSlopeLimiterUV (dgConservationLaw *claw, dgGroupCollection *groups);
dgSlopeLimiterUV (dgConservationLaw *claw, const dgGroupCollection *groups);
virtual ~dgSlopeLimiterUV();
virtual void apply (dgDofContainer *solution);
};
......
......@@ -36,8 +36,6 @@ 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> &bathL = cache.get(*_bathymetry, 0);
const fullMatrix<double> &bathR = cache.get(*_bathymetry, 1);
double g = slim3dParameters::g;
......@@ -70,7 +68,6 @@ dgConservationLawSW3dContinuity::dgConservationLawSW3dContinuity( int propagateF
_laxFriedrichsFactor = 0.0;
_uv = NULL;
_eta = NULL;
_wOld = NULL;
}
void dgConservationLawSW3dContinuity::setup() {
_volumeTerm0[""] = NULL;
......@@ -80,8 +77,6 @@ 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[""];
......
......@@ -12,7 +12,7 @@ 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, *_wOld;
const functor *_uv, *_eta, *_bathymetry;
bool _propagateFromTop;
double _laxFriedrichsFactor;
void riemann(functorCache &cache, fullMatrix<double> &val) const;
......@@ -28,7 +28,6 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
inline void setEta(const functor *f) { _eta = f; }
const functor *getEta()const {return _eta;}
const functor *getUV()const {return _uv;}
inline void setWOld(const functor *f) { _wOld = f; };
inline void setBathymetry(const functor *bathymetry) { _bathymetry = bathymetry;}
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
void finalize(const dgGroupOfElements &group, dgDofContainer &residual) const;
......
#ifndef _DG_CONSERVATION_LAW_SW_3D_MOMENTUM_H
#define _DG_CONSERVATION_LAW_SW_3D_MOMENTUM_H
#include "dgConservationLawFunction.h"
#include "dgConservationLawSW3dTracer.h"
#include "functor.h"
/**Horizontal momentum equation for the 3D shallow water. (u, v) */
class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
class dgConservationLawSW3dMomentum : public dgConservationLawSW3dIMEX {
void advection(functorCache &cm, fullMatrix<double> &val) const;
void source(functorCache &cm, fullMatrix<double> &val) const;
void riemann(functorCache &cm, fullMatrix<double> &val) const;
void diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const;
void maxConvectiveSpeed(functorCache &cm, fullMatrix<double> &val) const;
void boundarySurface(functorCache &cm, fullMatrix<double> &val, const functor *windF) const;
const functor *_uv2d, *_eta, *_rGrad, *_w, *_wM, *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry;
const functor *_uv2d, *_rGrad, *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry;
const functor *_fzero, *_fzerov;
const functor *_nuH, *_nuV, *_diffusiveFlux, *_ipTerm, *_maxDiffusivity;
const functor *_rhoDevInt, *_rhoDevMax;
double _laxFriedrichsFactor;
bool _withoutAdvection;
class boundaryWallArg;
void boundaryWall(functorCache &m, fullMatrix<double> &val, const boundaryWallArg arg) const;
class boundaryWallVerticalBottomArg;
void boundaryWallVerticalBottom(functorCache &m, fullMatrix<double> &val, const boundaryWallVerticalBottomArg arg) const;
public:
void setup();
void addIPTerm(functorCache &cm, fullMatrix<double> &val)const;
void diffusivityTensor(functorCache &cm, fullMatrix<double> &val) const;
/** set factor for scaling Lax-Friedrichs penalty term [default LFF = 0].
* 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;};
/**set eta function [eta] */
inline void setEta(const functor *eta){ _eta = eta;}
const functor * getEta()const{return _eta;}
/**set function for horizontal gradient of the baroclinic head [drdx drdy] */
inline void setRGrad(const functor *f){ _rGrad = f;}
inline void setRhoDevInt(const functor *f){ _rhoDevInt = f;}
......@@ -40,11 +35,6 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
const functor *getUV2d()const {return _uv2d;}
/**set depth averaged uvGrad function [dudx, dudy, dudz, dvdx, dvdy, dvdz] */
/**set vertical velocity function [w] */
inline void setW(const functor *w){ _w = w;}
inline const functor* getW(){ return _w;}
/**set vertical mesh velocity [wMesh] */
inline void setWMesh(const functor *wM){ _wM = wM;}
inline const functor* getWMesh(){ return _wM;}
/**Set vertical viscosity function (scalar) */
inline void setNuV(const functor *nuV){ _nuV = nuV;}
/**Set horizontal viscosity function (scalar) */
......@@ -64,57 +54,15 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
/**set the function to evaluate the bathymetry h (H = h+eta) */
inline void setBathymetry(const functor *bathymetry) { _bathymetry = bathymetry;}
/**set no to use the advection of momentum */
inline void setWithoutAdvection(bool b){ _withoutAdvection = b;}
void boundaryBottom(functorCache &cm, fullMatrix<double> &val, const functor *z0f, const functor *zBotD)const;
void boundaryWall(functorCache &cm, fullMatrix<double> &val)const;
dgConservationLawSW3dMomentum();
~dgConservationLawSW3dMomentum();
/** slip wall boundary */
dgBoundaryCondition *newBoundaryWall(const functor* extDiffusiveFlux=NULL, const functor* zBotDist=NULL);
/** same as boundary Wall, but for vertical bottom (uv2d is different) */
dgBoundaryCondition *newBoundaryWallVerticalBottom();
dgBoundaryCondition *newBoundaryWall();
dgBoundaryCondition *newBoundaryBottom(const functor* extDiffusiveFlux=NULL, const functor* zBotDist=NULL);
/**bnd condition for free surface */
dgBoundaryCondition *newBoundarySurface(const functor* wind);
};
/** Vertical advection and diffusion for 3d horizontal momentum equation for u or v.
* Since the system is the same for u and v, they are solved separately. */
class dgConservationLawSW3dMomWx : public dgConservationLawFunction {
void diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const;
void advectiveDiffusiveFlux(functorCache &cm, fullMatrix<double> &val) const;
void source(functorCache &cm, fullMatrix<double> &val) const;
void interfaceTerm(functorCache &cm, fullMatrix<double> &val) const;
void boundarySurface(functorCache &cm, fullMatrix<double> &val, const functor *wind) const;
const functor *_fzero, *_fzerov;
const functor *_nuV, *_diffusiveFlux, *_advectiveDiffusiveFlux, *_sourceTerm, *_interfaceTerm, *_ipTerm;
const functor *_w, *_wM, *_dwMdz;
bool _useConservativeALE;
double _laxFriedrichsFactor;
// determines whether u or v is being computed
int _comp;
public:
void diffusivityTensor(functorCache &cm, fullMatrix<double> &val) const;
void setup();
virtual bool isLinear() const {return true;}
/**set which component to solve 0 for u 1 for v [0] */
inline void setComponent(int c){ _comp = c;}
/**Set vertical viscosity function (scalar) */
inline void setNuV(const functor *nuV){ _nuV = nuV;}
dgConservationLawSW3dMomWx();
~dgConservationLawSW3dMomWx();
/**bnd condition for wall
* User may define a nonzero diffusive flux to create partial slip condition */
void boundaryBottom(functorCache &cm, fullMatrix<double> &val, const functor *extDiffFluxF, const functor *zBotDistF) const;
dgBoundaryCondition *newBoundaryBottom(const functor* extDiffusiveFlux = NULL,const functor* zBotDist = NULL);
/**bnd condition for free surface */
dgBoundaryCondition *newBoundarySurface(const functor* wind);
dgFaceTermAssembler *getFaceTerm(const dgGroupOfFaces &group) const;
/**set vertical velocity function [w] */
inline void setW(const functor *w){ _w = w;}
/**set vertical mesh velocity [wMesh] */
inline void setWMesh(const functor *wM){ _wM = wM;}
/**set derivative of vertical mesh velocity [DwMesh/Dz] */
inline void setWMeshDz(const functor *dwMdz){ _dwMdz = dwMdz;}
void setUseConservativeALE(bool b) { _useConservativeALE = b; };
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
};
#endif
......@@ -2,15 +2,19 @@
#include "functorMember.h"
#include "slim3dSolver.h"
void dgConservationLawSW3dTracer::setFullExplicit(){
_explicitC = NULL;
_explicitCGrad = NULL;
void dgConservationLawSW3dIMEX::setFullExplicit(){
_explicitSol = NULL;
_explicitSolGrad = NULL;
}
void dgConservationLawSW3dTracer::setImplicitVertical(const functor *explicitC, const functor *explicitCGrad)
void dgConservationLawSW3dIMEX::setImplicitVertical(const functor *explicitSol, const functor *explicitSolGrad)
{
_explicitC = explicitC;
_explicitCGrad = explicitCGrad;
_explicitSol = explicitSol;
_explicitSolGrad = explicitSolGrad;
}
dgConservationLawSW3dIMEX::dgConservationLawSW3dIMEX(int nbFields):dgConservationLawFunction(nbFields) {
setFullExplicit();
}
void dgConservationLawSW3dTracer::maxConvectiveSpeed(functorCache &cm, fullMatrix<double> &val) const {
......@@ -39,7 +43,7 @@ void dgConservationLawSW3dTracer::source(functorCache &cm, fullMatrix<double> &v
void dgConservationLawSW3dTracer::advection(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &c = cm.get(cm.solution());
const fullMatrix<double> &cE = cm.get(_explicitC ? *_explicitC : cm.solution());
const fullMatrix<double> &cE = cm.get(_explicitSol ? *_explicitSol : cm.solution());
const fullMatrix<double> &uv = cm.get(*_uv);
const fullMatrix<double> &w = cm.get(*_w);
const fullMatrix<double> *diffusiveFlux = _diffusiveFlux ? &cm.get(*_diffusiveFlux) : NULL;
......@@ -60,8 +64,8 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
const fullMatrix<double> &normals = cm.get(*functor::getNormals());
const fullMatrix<double> &cL = cm.get(cm.solution(), 0);
const fullMatrix<double> &cR = cm.get(cm.solution(), 1);
const fullMatrix<double> &cEL = cm.get(_explicitC ? *_explicitC : cm.solution(), 0);
const fullMatrix<double> &cER = cm.get(_explicitC ? *_explicitC : cm.solution(), 1);
const fullMatrix<double> &cEL = cm.get(_explicitSol ? *_explicitSol : cm.solution(), 0);
const fullMatrix<double> &cER = cm.get(_explicitSol ? *_explicitSol : cm.solution(), 1);
const fullMatrix<double> &bathL = cm.get(*_bathymetry, 0);
const fullMatrix<double> &bathR = cm.get(*_bathymetry, 1);
const fullMatrix<double> &etaL = cm.get(*_eta, 0);
......@@ -98,11 +102,20 @@ void dgConservationLawSW3dTracer::diffusiveFlux(functorCache &cm, fullMatrix<dou
const fullMatrix<double> &nuH = cm.get(*_kappaH);
const fullMatrix<double> &nuV = cm.get(*_kappaV);
int nQP = cm.nEvaluationPoint();
val.resize(nQP, 3, false);
for(int i=0; i< nQP; i++) {
val(i,0) = -solGradE(i,0)*nuH(i,0);
val(i,1) = -solGradE(i,1)*nuH(i,0);
val(i,2) = -solGrad(i,2)*nuV(i,0);
if (_kappaH) {
const fullMatrix<double> &nuH = cm.get(*_kappaH);
const fullMatrix<double> &solGradE = cm.get(_explicitSolGrad ? *_explicitSolGrad : cm.solutionGradient());
for(int i=0; i< nQP; i++) {
val(i,0) = -solGradE(i,0)*nuH(i,0);
val(i,1) = -solGradE(i,1)*nuH(i,0);
}
}
if (_kappaV) {
const fullMatrix<double> &solGrad = cm.get(cm.solutionGradient());
const fullMatrix<double> &nuV = cm.get(*_kappaV);
for(int i=0; i< nQP; i++) {
val(i,2) = -solGrad(i,2)*nuV(i,0);
}
}
}
......@@ -207,8 +220,8 @@ void dgConservationLawSW3dTracer::setup () {
void dgConservationLawSW3dTracer::addIPTerm(functorCache &cm, fullMatrix<double> &val)const {
const fullMatrix<double> &cL = cm.get(cm.solution(), 0);
const fullMatrix<double> &cR = cm.get(cm.solution(), 1);
const fullMatrix<double> &cEL = cm.get(_explicitC ? *_explicitC : cm.solution(), 0);
const fullMatrix<double> &cER = cm.get(_explicitC ? *_explicitC : cm.solution(), 1);
const fullMatrix<double> &cEL = cm.get(_explicitSol ? *_explicitSol : cm.solution(), 0);
const fullMatrix<double> &cER = cm.get(_explicitSol ? *_explicitSol : cm.solution(), 1);
const fullMatrix<double> &normals = cm.get(*functor::getNormals());
const fullMatrix<double> &diffusiveFluxL = cm.get(*_diffusiveFlux, 0);
const fullMatrix<double> &diffusiveFluxR = cm.get(*_diffusiveFlux, 1);
......@@ -246,7 +259,7 @@ dgConservationLawSW3dTracer::~dgConservationLawSW3dTracer() {
if (_diffusiveFlux) delete _diffusiveFlux;
}
dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawFunction(1)
dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawSW3dIMEX(1)
{
_sourceFunc = NULL;
_uv = NULL;
......@@ -259,5 +272,4 @@ dgConservationLawSW3dTracer::dgConservationLawSW3dTracer() : dgConservationLawFu
_ref=NULL;
_dt=0;
_T=1;
setFullExplicit();
}
......@@ -5,8 +5,23 @@
/**Advection and diffusion of a scalar, the advection and diffusion are provided by functions.
* \f[\frac{dc}{dt} + v. \nabla c - nu \Delta c = 0\f]
*/
class dgConservationLawSW3dTracer : public dgConservationLawFunction {
const functor *_uv, *_w, *_eta, *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_sourceFunc, *_explicitC, *_explicitCGrad;
class dgConservationLawSW3dIMEX : public dgConservationLawFunction {
protected:
const functor *_uv, *_w, *_eta, *_explicitSol, *_explicitSolGrad;
dgConservationLawSW3dIMEX(int nbFields);
public:
inline void setEta(const functor *eta) {_eta = eta;};
inline void setUV(const functor *f) {_uv = f;};
inline void setW(const functor *f) {_w = f;};
const functor *getUV()const {return _uv;}
const functor *getEta()const {return _eta;}
const functor *getW()const {return _w;}
void setFullExplicit();
void setImplicitVertical(const functor *explicitSol, const functor *explicitSolGrad);
};
class dgConservationLawSW3dTracer : public dgConservationLawSW3dIMEX {
const functor *_bathymetry, *_kappaH, *_kappaV, *_diffusiveFlux, *_sourceFunc;
int _dt, _T;
const functor *_ref;
void source(functorCache &cm, fullMatrix<double> &val) const ;
......@@ -22,19 +37,11 @@ class dgConservationLawSW3dTracer : public dgConservationLawFunction {
* Horizontal and vertical diffusivity are given by the functions kappaH and kappaV. */
/**A new advection-diffusion law */
dgConservationLawSW3dTracer();
/** set factor for scaling Lax-Friedrichs penalty term [default LFF = 0]. */
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 setKappaH(const functor *kap) {_kappaH = kap;};
inline void setKappaV(const functor *kap) {_kappaV = kap;};
inline void setEta(const functor *eta) {_eta = eta;};
inline void setBathymetry(const functor *bath) {_bathymetry = bath;};
void setFullExplicit();
void setImplicitVertical(const functor *explicitC, const functor *explicitGradC);
const functor *getUV()const {return _uv;}
const functor *getEta()const {return _eta;}
~dgConservationLawSW3dTracer();
/** Boundary condition that sets the in-flowing tracer concentration to a prescribed value */
void boundaryInFlux(functorCache &cm, fullMatrix<double> &val, const functor *outValue) const;
......
......@@ -53,9 +53,9 @@ static void setCp(double f){slim3dParameters::Cp = f;}
%include "slim3dTimeIntegrator.h"
%include "slim3dFunctions.h"
%include "slim3dUtils.h"
%include "dgConservationLawSW3dMomentum.h"
%include "dgConservationLawSW3dContinuity.h"
%include "dgConservationLawSW3dTracer.h"
%include "dgConservationLawSW3dMomentum.h"
%include "dgSW3dVerticalModel.h"
%include "dgSW3dTurbulenceGOTM.h"
%include "dgEddyTransportFlux.h"
......
......@@ -693,7 +693,6 @@ slim3dSolverFunctions::~slim3dSolverFunctions()
slim3dSolverEquations::~slim3dSolverEquations()
{
delete horMomEq;
delete vertMomUEq;
delete wEq;
delete newRGradEq;
#ifdef HAVE_GOTM
......@@ -722,7 +721,6 @@ slim3dSolverEquations::slim3dSolverEquations(slim3dSolver* solver) : _solver(sol
{
// create all equations with no arguments
horMomEq = new dgConservationLawSW3dMomentum();
vertMomUEq = new dgConservationLawSW3dMomWx();
wEq = new dgConservationLawSW3dContinuity(0);
newRGradEq = new dgConservationLawSW3dPressure(INTEGRATE_FROM_TOP);
#ifdef HAVE_GOTM
......@@ -829,15 +827,8 @@ void slim3dSolverEquations::initialize()
horMomEq->setBathymetry(funcs->bathFunc2d);
if (funcs->coriolisFunc)
horMomEq->setCoriolisFactor(funcs->coriolisFunc);
if (funcs->nuh)
horMomEq->setNuH(funcs->nuh);
if (!_solver->_solveUVDiff && funcs->nuvFunc) {
Msg::Warning("Vertical diffusion equation of UV disabled, assigning nuv to horizontal momentum equation");
horMomEq->setNuV(funcs->nuvFunc);
}
if ( _solver->_solveUVDiff ) {
vertMomUEq->setNuV(funcs->nuvFunc);
}
horMomEq->setNuH(funcs->nuh);
horMomEq->setNuV(funcs->nuvFunc);
wEq->setBathymetry(funcs->bathFunc2d);
......
......@@ -106,8 +106,6 @@ public:
// TODO rename
/** 3d horizontal momentum equation */
dgConservationLawSW3dMomentum *horMomEq;
/** Vertical diffusion for horizontal momentum equation */
dgConservationLawSW3dMomWx *vertMomUEq;
/** 2D shallow water equations */
dgConservationLawShallowWater2d *sw2DEq;
dgUVHorDivergence *uvHorDivEq;
......
......@@ -97,24 +97,30 @@ class slim3dMassMatrix {
_mult(dof,_matrices);
}
void addToSystem(dgDofManager &dof, double alpha = 1.) {
fullMatrix<double> m, alpham;
fullMatrix<double> m, m2;
int iG = 0;
for (auto &ms : _matrices) {
int nF = dof.nField(*dof.getGroups().getElementGroup(iG));
ms.getBlockProxy(0,m);
m2.resize(nF*m.size1(), nF*m.size2(), true);
for (size_t i = 0; i<ms.nBlock(); ++i) {
ms.getBlockProxy(i,m);
alpham.copy(m);
alpham.scale(alpha);
dof.assembleLHSMatrix(iG, i, iG, i, alpham);
for (int f = 0; f < nF; ++f) {
for (int a = 0; a <m.size1(); a++)
for (int b = 0; b < m.size2(); b++){
m2(f*m.size1()+a,f*m.size2()+b) = m(a,b)*alpha;
}
}
dof.assembleLHSMatrix(iG, i, iG, i, m2);
}
iG++;
}
}
};
class slim3dTracerSolver {
class slim3dIMEXSolver {
bool _implicitVertical, _applyLimiter;
dgConservationLawSW3dTracer &_law;
dgConservationLawSW3dIMEX &_law;
dgDofManager *_dof;
linearSystem<fullMatrix<double> > *_sys;
dgDofContainer &_c, _MCOld;
......@@ -140,16 +146,16 @@ class slim3dTracerSolver {
_dof->finalizeSparsity();
}
public :
slim3dTracerSolver(dgConservationLawSW3dTracer &law, dgDofContainer &c, bool implicitVertical, bool applyLimiter): _implicitVertical(implicitVertical), _applyLimiter(applyLimiter), _law(law), _c(c),_MCOld(*_c.getGroups(), 1){
slim3dIMEXSolver(dgConservationLawSW3dIMEX &law, dgDofContainer &c, bool implicitVertical, bool applyLimiter): _implicitVertical(implicitVertical), _applyLimiter(applyLimiter), _law(law), _c(c),_MCOld(*_c.getGroups(), law.getNbFields()){
if (_implicitVertical){
// _sys = new dgLinearSystemExtrusion(&law, _c.getGroups(), solver->columnInfo);
_sys = new linearSystemPETSc<fullMatrix<double> >();
//_sys->setParameter("petsc_solver_options", "-pc_type lu");
_dof = dgDofManager::newDGBlock (_c.getGroups(), 1, _sys);
_dof = dgDofManager::newDGBlock (_c.getGroups(), _law.getNbFields(), _sys);
_fillDofSparsity();
}
};
~slim3dTracerSolver() {
~slim3dIMEXSolver() {
if (_implicitVertical){
delete _sys;
delete _dof;
......@@ -164,10 +170,10 @@ class slim3dTracerSolver {
_law.setUV(uv);
_law.setW(w);
_law.setEta(eta);
dgDofContainer rhs(groups,1);
dgDofContainer rhs(groups,_law.getNbFields());
rhs.setAll(0.);
if(_implicitVertical) {
dgDofContainer SN(groups,1);
dgDofContainer SN(groups,_law.getNbFields());
SN.copy(_c);
_c.setAll(0.);
_law.setImplicitVertical(SN.getFunction(), SN.getFunctionGradient());
......@@ -197,14 +203,16 @@ class slim3dTracerSolver {
_c.axpy(rhs, dt);
massnext.multByInverse(_c);
}
if(_applyLimiter)
dgSlopeLimiterOriginal(&_law, &groups).apply(&_c);
_c.scatter();
if(_applyLimiter) {
if(_law.getNbFields() == 1)
dgSlopeLimiterOriginal(&_law, &groups).apply(&_c);
else
dgSlopeLimiterUV(&_law, &groups).apply(&_c);
}
}
};
class slim3dTimeIntegrator::tracerDeviationChecker {
dgDofContainer* _dof;
bool _initialized;
......@@ -238,11 +246,6 @@ slim3dTimeIntegrator::slim3dTimeIntegrator(slim3dSolver* solver) : _solver(solve
_exploded = false;
slim3dSolverEquations *e = solver->getEquations();
slim3dSolverDofs *d = solver->getDofs();
uSys = NULL;
momUNewton = NULL;
momVNewton = NULL;
momUNewtonSol = NULL;
momVNewtonSol = NULL;
_volConservationChecker = _SConservationChecker = _TConservationChecker = NULL;
_TDeviationChecker = _SDeviationChecker = NULL;
_sedDeviationChecker.resize(0);
......@@ -255,29 +258,19 @@ slim3dTimeIntegrator::slim3dTimeIntegrator(slim3dSolver* solver) : _solver(solve
rk = new dgDIRK(*e->sw2DEq, *dofManager2D, 1);
rk->getNewton().setAtol(1e-10);
sw2DSol = new dgDofContainer(*solver->groups2d, 3);
if ( _solver->getSolveUVImplicitVerticalDiffusion() || _solver->getSolveUVImplicitVerticalAdvection()) {
// uSys = new dgLinearSystemExtrusion(e->vertMomUEq, groups3d, solver->columnInfo);
uSys = new linearSystemPETSc<fullMatrix<double > >();
dofManagerU = dgDofManager::newDGBlock(solver->groups3d, 1, uSys);
momUNewton = new dgNewton(*dofManagerU, *e->vertMomUEq);
momUNewton->setMaxIt(4);
momUNewtonSol = new dgDofContainer(*solver->groups3d, 1);
momVNewton = new dgNewton(*dofManagerU, *e->vertMomUEq);
momVNewton->setMaxIt(4);
momVNewtonSol = new dgDofContainer(*solver->groups3d, 1);
}
_uvSolver = new slim3dIMEXSolver(*e->horMomEq, *d->uvDof, _solver->getSolveUVImplicitVerticalDiffusion(), _solver->getFlagUVLimiter());
if (_solver->getSolveS())
_tracerSolvers.push_back(new slim3dTracerSolver(*e->SEq, *d->SDof, _solver->getSolveSImplicitVerticalDiffusion(), _solver->getFlagSLimiter()));
_tracerSolvers.push_back(new slim3dIMEXSolver(*e->SEq, *d->SDof, _solver->getSolveSImplicitVerticalDiffusion(), _solver->getFlagSLimiter()));
if (_solver->getSolveT())
_tracerSolvers.push_back(new slim3dTracerSolver(*e->TEq, *d->TDof, _solver->getSolveTImplicitVerticalDiffusion(), _solver->getFlagTLimiter()));
_tracerSolvers.push_back(new slim3dIMEXSolver(*e->TEq, *d->TDof, _solver->getSolveTImplicitVerticalDiffusion(), _solver->getFlagTLimiter()));
#ifdef HAVE_GOTM
if (_solver->getSolveTurbulence() && _solver->getAdvectTurbulentEnergy()) {
_tracerSolvers.push_back(new slim3dTracerSolver(*e->tkeAdvEq, *d->tkeDof, false, true));
_tracerSolvers.push_back(new slim3dTracerSolver(*e->epsAdvEq, *d->epsDof, false, true));
_tracerSolvers.push_back(new slim3dIMEXSolver(*e->tkeAdvEq, *d->tkeDof, false, true));
_tracerSolvers.push_back(new slim3dIMEXSolver(*e->epsAdvEq, *d->epsDof, false, true));
}
#endif
//for (int iSedEq = 0; iSedEq < _solver->getNumSedEq(); ++iSedEq)
//_tracerSolvers.push_back(new slim3dTracerSolver(*e->sedEq[iSedEq], *d->sedDof[iSedEq],_solver->getSolveSedImplicitVerticalDiffusion(), false));
//_tracerSolvers.push_back(new slim3dIMEXSolver(*e->sedEq[iSedEq], *d->sedDof[iSedEq],_solver->getSolveSedImplicitVerticalDiffusion(), false));
if (_solver->getUseAdaptiveVerticalGrid()){
dgGroupCollection *groups3d = _solver->groups3d;
vJumpDof_K = new dgDofContainer(*groups3d, 1);
......@@ -288,14 +281,9 @@ slim3dTimeIntegrator::slim3dTimeIntegrator(slim3dSolver* solver) : _solver(solve
slim3dTimeIntegrator::~slim3dTimeIntegrator()
{
if ( uSys ) delete uSys;
if ( sw2DSys ) delete sw2DSys;
if ( rk ) delete rk;
if ( sw2DSol ) delete sw2DSol;
if( momUNewton ) delete momUNewton;
if( momVNewton ) delete momVNewton;
if( momUNewtonSol ) delete momUNewtonSol;
if( momVNewtonSol ) delete momVNewtonSol;
if (_volConservationChecker) delete _volConservationChecker;
if (_SConservationChecker) delete _SConservationChecker;
if (_TConservationChecker) delete _TConservationChecker;
......@@ -305,7 +293,7 @@ slim3dTimeIntegrator::~slim3dTimeIntegrator()
if ( vJumpDof_K ) delete vJumpDof_K;
for (auto ts : _tracerSolvers)
delete ts;