Commit acf3c76a authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

modules/shallowWater : fix nudging term

parent cff75936
......@@ -28,11 +28,11 @@ class _merger_3d :
u1 = cm.get(self._u1)
u2 = cm.get(self._u2)
u3 = cm.get(self._u3)
nPt = sol.shape
nPt = u1.shape
val = np.zeros([nPt[0], 3])
val[:,0] = u1
val[:,1] = u2
val[:,2] = u3
val[:,0] = u1[:,0]
val[:,1] = u2[:,0]
val[:,2] = u3[:,0]
return val
class _merger_2d :
......@@ -43,10 +43,10 @@ class _merger_2d :
def call(self, cm) :
u1 = cm.get(self._u1)
u2 = cm.get(self._u2)
nPt = sol.shape
nPt = u1.shape
val = np.zeros([nPt[0], 2])
val[:,0] = u1
val[:,1] = u2
val[:,0] = u1[:,0]
val[:,1] = u2[:,0]
return val
class _assemble_solution_global :
......@@ -454,7 +454,8 @@ class ShallowWater2d:
else:
self._external_sol_local = _merger_3d(_external_sse, _external_ux, _external_uy).functor
self._ref += [_coeff_nudging, _external_sse, _external_ux, _external_uy, _external_uz, _external_vel_local_, _external_vel_local_x, _external_vel_local_y]
self._equation.setAbsLayer(_coeff_nudging, self._external_sol_local)
self._equation.setNudgingCoefficient(_coeff_nudging)
self._equation.setNudgingVelocity(self._external_sol_local)
def set_boundary_coast(self, physical_tag, slip=False):
"""
......
......@@ -217,14 +217,14 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
class dgConservationLawShallowWater2d::source: public function {
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _windStress, _nu, _xyz;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _absLayerCoef, _absLayerExtForc;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _nudgingCoeff, _nudgingVel;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer;
double _R, _g, _density;
public :
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF, const functor *windStressF,
const functor *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R, double g, double density,
const functor *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *absLayerCoefF, const functor *absLayerExtForcF) :function (3){
const functor *nudgingCoeffF, const functor *nudgingVelF) :function (3){
_linear = linear;
_R = R;
_g = g;
......@@ -240,9 +240,9 @@ class dgConservationLawShallowWater2d::source: public function {
setArgument(linearDissipation,linearDissipationF);
setArgument(quadraticDissipation,quadraticDissipationF);
_absLayer = false;
if (absLayerCoefF && absLayerExtForcF) {
setArgument(_absLayerCoef,absLayerCoefF);
setArgument(_absLayerExtForc,absLayerExtForcF);
if (nudgingVelF && nudgingCoeffF) {
setArgument(_nudgingCoeff,nudgingCoeffF);
setArgument(_nudgingVel,nudgingVelF);
_absLayer = true;
Msg::Info("AbsLayer is used.");
}
......@@ -319,20 +319,20 @@ class dgConservationLawShallowWater2d::source: public function {
val(i,1) *= J*J;
val(i,2) *= J*J;
if (_absLayer) {
//double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
//val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1));
//val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));
// Absorbing Layer of Lavelle & co (Pretty Good Sponge)
/*double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
val(i,0) -= (sigma[0]+sigma[1]) * (eta-_absLayerExtForc(i,0));
val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1));
val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));*/
// Absorbing Layer of Martinsen & co (Sponge Layer - FRS)
double sigma = _absLayerCoef(i,0);
double H = _bathymetry(i,0) + eta;
double sigma = _nudgingCoeff(i,0);
val(i,1) -= sigma * (u-_nudgingVel(i,0)/H);
val(i,2) -= sigma * (v-_nudgingVel(i,1)/H);
// Absorbing Layer of Lavelle & co (Pretty Good Sponge)
/*double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
val(i,0) -= (sigma[0]+sigma[1]) * (eta-_absLayerExtForc(i,0));
val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1));
val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));*/
// Absorbing Layer of Martinsen & co (Sponge Layer - FRS)
//double sigma = _absLayerCoef(i,0);
//val(i,0) -= sigma * (eta-_absLayerExtForc(i,0));
val(i,1) -= sigma * (u-_absLayerExtForc(i,1)/H);
val(i,2) -= sigma * (v-_absLayerExtForc(i,2)/H);
//val(i,1) -= sigma * (u-_absLayerExtForc(i,1));
//val(i,2) -= sigma * (v-_absLayerExtForc(i,2));
}
}
}
......@@ -599,6 +599,8 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_constantJacIMEX=false;
_linearFilter=false;
_laxFriedrichs = false;
_nudgingCoeff = NULL;
_nudgingVel = NULL;
}
void dgConservationLawShallowWater2d::setImexMode(imexMode mode) {
......@@ -632,11 +634,11 @@ void dgConservationLawShallowWater2d::setup() {
_diffusion = _nu ? new diffusivity(_nu) : NULL;
_ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL;
//_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear);
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel);
_gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz, _g);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false, _g);
_sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel);
_gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g);
_riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g);
......@@ -679,24 +681,6 @@ dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
}
}
void dgConservationLawShallowWater2d::setAbsLayer(const std::string tag, const functor *absCoef, const functor *extFunc) {
if (!absCoef || !extFunc)
Msg::Fatal("Dg/ShallowWater/AbsLayer: Bad 'absCoef' and/or 'extFunc'.");
if (_volumeTerm0[tag]) delete _volumeTerm0[tag];
_volumeTerm0[tag] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient,
_nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, absCoef, extFunc);
}
void dgConservationLawShallowWater2d::setAbsLayer(const functor *absCoef, const functor *extFunc) {
if (!absCoef || !extFunc)
Msg::Fatal("Dg/ShallowWater/AbsLayer: Bad 'absCoef' and/or 'extFunc'.");
if (_volumeTerm0[""]) delete _volumeTerm0[""];
_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient,
_nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, absCoef, extFunc);
}
/*==============================================================================
* Specific boundary conditions
*============================================================================*/
......
......@@ -41,7 +41,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_windStress, *_coriolisFactor, *_coordinatesF;
const functor *_movingBathFactor, *_movingBathFactorGradient;
const functor *_fzero, *_fzerov, *_nu, *_diffusiveFlux, *_diffusion, *_ipTerm, *_xyz;
const functor *_nudgingCoeff, *_nudgingVel;
const functor *_sourceTerm, *_sourceTermLin, *_gradPsiTerm, *_gradPsiTermLin,
*_riemannTerm, *_riemannTermLin, *_maxSpeed, *_maxSpeedFilt;
imexMode _imexMode;
......@@ -85,9 +85,10 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
inline void setBathymetryGradient(functor *bathymetryGradient) { _bathymetryGradient = bathymetryGradient;}
/**set the function to evaluate the diffusion term */
inline void setDiffusivity(const functor *nuFunction) { _nu = nuFunction;}
/**set an absorbing layer. Param: TAG, AbsCoef, ExtFields */
void setAbsLayer(const std::string, const functor *, const functor *);
void setAbsLayer(const functor *absCoeff, const functor *extFunc);
/**set the coefficient of the nudging layer. */
void setNudgingCoefficient(const functor *coeff){_nudgingCoeff = coeff;}
/**set the reference velocity for the nudging layer. */
void setNudgingVelocity(const functor *nudgingVel){_nudgingVel = nudgingVel;}
/**use the moving bath WD (Karna et al, 2011) and define the bathymetry and the moving bathymetry */
void setMovingBathWettingDrying(double alphaMovingBathWD);
/** a factor to smooth the if in Hv term in the roe riemann solver */
......
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