Authentication method changed. UCLouvain users, please use the "UCLouvain SSO" button to connect on the website and use ssh + ssh key (https://git.immc.ucl.ac.be/-/profile/keys) or https + personnal access token (https://git.immc.ucl.ac.be/-/profile/personal_access_tokens) in your git clients.

Commit 8b4b8d58 authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

nudging : modified term

parent 8937fd61
...@@ -465,17 +465,18 @@ class ShallowWater2d: ...@@ -465,17 +465,18 @@ class ShallowWater2d:
self._ref += [_forcing_x, _forcing_y] self._ref += [_forcing_x, _forcing_y]
self._equation.setSource(self._forcing_local) self._equation.setSource(self._forcing_local)
def set_nudging(self, coefficient, external_sse, external_ux, external_uy, external_uz=None, ramp_period=0): def set_nudging(self, coefficient, external_sse, external_ux, external_uy, external_uz=None, ramp_period=0, transport_flux=False):
""" """
Add a nudging term as a source term in the shallow water equation. Add a nudging term as a source term in the shallow water equation.
keyword arguments: keyword arguments:
coefficient -- netcdf or .msh file containing the nudging coefficient for the whole domain [in s^-1] coefficient -- netcdf or .msh file containing the nudging coefficient for the whole domain [in s^-1]
external_sse -- netcdf or .msh file containing the external sea surface elevation for the whole domain [in m] external_sse -- netcdf or .msh file containing the external sea surface elevation for the whole domain [in m]
external_ux -- netcdf or .msh file containing the external velocity along the x-axis expressed in the local basis (except for equation on a manifold for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in m/s] external_ux -- netcdf or .msh file containing the external velocity along the x-axis expressed in the local basis (except for equation on a manifold for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in m/s]
external_uy -- netcdf or .msh file containing the external velocity along the y-axis expressed in the local basis (except for equation on a manifold for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in m/s] external_uy -- netcdf or .msh file containing the external velocity along the y-axis expressed in the local basis (except for equation on a manifold for which it has to be expressed in the global basis (x,y,z)) for the whole domain [in m/s]
external_uz -- netcdf or .msh file containing the external velocity along the z-axis expressed in the global basis (x,y,z) for the whole domain [in m/s] external_uz -- netcdf or .msh file containing the external velocity along the z-axis expressed in the global basis (x,y,z) for the whole domain [in m/s]
ramp_period -- period for the linear ramp before applying the nudging [in s] (Default: 0s, no ramp) ramp_period -- period for the linear ramp before applying the nudging [in s] (Default: 0s, no ramp)
transport_flux -- set whether ux and uy (and uz) are velocities (False) or transport (True) (Default: False)
""" """
_coeff_nudging = _load_function(coefficient, self._domain._groups) _coeff_nudging = _load_function(coefficient, self._domain._groups)
_external_sse = _load_function(external_sse, self._domain._groups) _external_sse = _load_function(external_sse, self._domain._groups)
...@@ -491,16 +492,16 @@ class ShallowWater2d: ...@@ -491,16 +492,16 @@ class ShallowWater2d:
_external_vel_local_x = functorNumpy(lambda cm : cm.get(_external_vel_local_)[:,0]) _external_vel_local_x = functorNumpy(lambda cm : cm.get(_external_vel_local_)[:,0])
_external_vel_local_y = functorNumpy(lambda cm : cm.get(_external_vel_local_)[:,1]) _external_vel_local_y = functorNumpy(lambda cm : cm.get(_external_vel_local_)[:,1])
_external_sol_local = _merger_3d(_external_sse, _external_vel_local_x, _external_vel_local_y ).functor _external_sol_local = _merger_3d(_external_sse, _external_vel_local_x, _external_vel_local_y ).functor
self._ref += [_coeff_nudging, _external_sse, _external_ux, _external_uy, _external_uz, _external_vel_local_, _external_vel_local_x, _external_vel_local_y, _external_sol_local, ramp_period] self._ref += [_external_uz, _external_vel_local_x, _external_vel_local_y]
else: else:
_external_sol_local = _merger_3d(_external_sse, _external_ux, _external_uy).functor _external_sol_local = _merger_3d(_external_sse, _external_ux, _external_uy).functor
self._ref += [_coeff_nudging, _external_sse, _external_ux, _external_uy, _external_vel_local_, _external_sol_local, ramp_period] self._ref += [_coeff_nudging, _external_sse, _external_ux, _external_uy, _external_vel_local_, _external_sol_local, ramp_period, transport_flux]
self._equation.setNudgingCoefficient(_coeff_nudging) self._equation.setNudgingCoefficient(_coeff_nudging)
if (ramp_period > 0.): if (ramp_period > 0.):
self._external_sol_local = ramp3d(_external_sol_local, function.getTime(), self._initial_time, ramp_period) self._external_sol_local = ramp3d(_external_sol_local, function.getTime(), self._initial_time, ramp_period)
else: else:
self._external_sol_local = _external_sol_local self._external_sol_local = _external_sol_local
self._equation.setNudgingVelocity(self._external_sol_local) self._equation.setNudgingVelocity(self._external_sol_local, transport_flux)
def set_boundary_coast(self, physical_tag, slip=False): def set_boundary_coast(self, physical_tag, slip=False):
""" """
......
...@@ -218,20 +218,19 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function { ...@@ -218,20 +218,19 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
class dgConservationLawShallowWater2d::source: public function { class dgConservationLawShallowWater2d::source: public function {
fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _windStress, _nu, _xyz; fullMatrix<double> sol, solGradient, coriolisFactor, linearDissipation, quadraticDissipation, _source, _windStress, _nu, _xyz;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _nudgingCoeff, _nudgingVel; fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _nudgingCoeff, _nudgingVel;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer; bool _linear, _isDiffusive, _useMovingBathWD, _absLayer, _nudgingVelIsTransport;
double _R, _g, _density; double _R, _g, _density;
public : public :
source(const functor *linearDissipationF, const functor *quadraticDissipationF, const functor *coriolisFactorF, const functor *sourceF, const functor *windStressF, 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 *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 *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *nudgingCoeffF, const functor *nudgingVelF) :function (3){ const functor *nudgingCoeffF, const functor *nudgingVelF, bool nudgingVelIsTransport) :function (3){
_linear = linear; _linear = linear;
_R = R; _R = R;
_g = g; _g = g;
_density = density; _density = density;
_isDiffusive = nuF; _isDiffusive = nuF;
_useMovingBathWD = movingBathFactorF; _useMovingBathWD = movingBathFactorF;
if(nuF && bathymetryF && !bathymetryGradientF){ if(nuF && bathymetryF && !bathymetryGradientF){
Msg::Fatal("The gradient of the bathymetry is missing since we want to use diffusion with a bathymetry"); Msg::Fatal("The gradient of the bathymetry is missing since we want to use diffusion with a bathymetry");
} }
...@@ -243,8 +242,8 @@ class dgConservationLawShallowWater2d::source: public function { ...@@ -243,8 +242,8 @@ class dgConservationLawShallowWater2d::source: public function {
if (nudgingVelF && nudgingCoeffF) { if (nudgingVelF && nudgingCoeffF) {
setArgument(_nudgingCoeff,nudgingCoeffF); setArgument(_nudgingCoeff,nudgingCoeffF);
setArgument(_nudgingVel,nudgingVelF); setArgument(_nudgingVel,nudgingVelF);
_nudgingVelIsTransport = nudgingVelIsTransport;
_absLayer = true; _absLayer = true;
Msg::Info("AbsLayer is used.");
} }
setArgument(_source,sourceF); setArgument(_source,sourceF);
setArgument(_windStress,windStressF); setArgument(_windStress,windStressF);
...@@ -319,20 +318,26 @@ class dgConservationLawShallowWater2d::source: public function { ...@@ -319,20 +318,26 @@ class dgConservationLawShallowWater2d::source: public function {
val(i,1) *= J*J; val(i,1) *= J*J;
val(i,2) *= J*J; val(i,2) *= J*J;
if (_absLayer) { if (_absLayer) {
double uref = _nudgingVel(i,1);
double vref = _nudgingVel(i,2);
if (_nudgingVelIsTransport){
uref = uref/_bathymetry(i,0);
vref = vref/_bathymetry(i,0);
}
//double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)}; //double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
//val(i,1) -= sigma[0] * (u-_nudgingVel(i,1)/H); double sigma = _nudgingCoeff(i,0);
//val(i,2) -= sigma[1] * (v-_nudgingVel(i,2)/H); val(i,1) -= sigma * (u-uref);
val(i,2) -= sigma * (v-vref);
// Absorbing Layer of Lavelle & co (Pretty Good Sponge) // Absorbing Layer of Lavelle & co (Pretty Good Sponge)
/*double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)}; /*double sigma[2] = {_absLayerCoef(i,0), _absLayerCoef(i,1)};
val(i,0) -= (sigma[0]+sigma[1]) * (eta-_absLayerExtForc(i,0)); val(i,0) -= (sigma[0]+sigma[1]) * (eta-_absLayerExtForc(i,0));
val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1)); val(i,1) -= sigma[0] * (u-_absLayerExtForc(i,1));
val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));*/ val(i,2) -= sigma[1] * (v-_absLayerExtForc(i,2));*/
// Absorbing Layer of Martinsen & co (Sponge Layer - FRS) // Absorbing Layer of Martinsen & co (Sponge Layer - FRS)
double H = _bathymetry(i,0) + eta; //double sigma = _nudgingCoeff(i,0);
double sigma = _nudgingCoeff(i,0); //val(i,0) -= sigma * (eta-_nudgingVel(i,0));
val(i,0) -= sigma * (eta-_nudgingVel(i,0)); //val(i,1) -= sigma * (u-uref);
val(i,1) -= sigma * (u-_nudgingVel(i,1)/H); //val(i,2) -= sigma * (v-vref);
val(i,2) -= sigma * (v-_nudgingVel(i,2)/H);
} }
} }
} }
...@@ -601,6 +606,7 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat ...@@ -601,6 +606,7 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_laxFriedrichs = false; _laxFriedrichs = false;
_nudgingCoeff = NULL; _nudgingCoeff = NULL;
_nudgingVel = NULL; _nudgingVel = NULL;
_nudgingVelIsTransport = false;
} }
void dgConservationLawShallowWater2d::setImexMode(imexMode mode) { void dgConservationLawShallowWater2d::setImexMode(imexMode mode) {
...@@ -634,11 +640,11 @@ void dgConservationLawShallowWater2d::setup() { ...@@ -634,11 +640,11 @@ void dgConservationLawShallowWater2d::setup() {
_diffusion = _nu ? new diffusivity(_nu) : NULL; _diffusion = _nu ? new diffusivity(_nu) : NULL;
_ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL; _ipTerm = _nu ? dgNewIpTermIsotropicOnSphere(3, _diffusiveFlux, _diffusion) : NULL;
//_volumeTerm0[""] = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _linear); //_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, _nudgingCoeff, _nudgingVel); _sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel, _nudgingVelIsTransport);
_gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz, _g); _gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz, _g);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann); _riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false, _g); _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, _nudgingCoeff, _nudgingVel); _sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _windStress, _bathymetry, _bathymetryGradient, _nu, true, _R, _g, _density, _movingBathFactor, _movingBathFactorGradient, _xyz, _nudgingCoeff, _nudgingVel, _nudgingVelIsTransport);
_gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g); _gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz, _g);
_riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann); _riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g); _maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true, _g);
......
...@@ -32,6 +32,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction { ...@@ -32,6 +32,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
double _alphaMovingBathWD; double _alphaMovingBathWD;
double _upwindFactorRiemann; double _upwindFactorRiemann;
bool _laxFriedrichs; bool _laxFriedrichs;
bool _nudgingVelIsTransport;
double _R; double _R;
double _g; double _g;
double _density; double _density;
...@@ -88,7 +89,10 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction { ...@@ -88,7 +89,10 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
/**set the coefficient of the nudging layer. */ /**set the coefficient of the nudging layer. */
void setNudgingCoefficient(const functor *coeff){_nudgingCoeff = coeff;} void setNudgingCoefficient(const functor *coeff){_nudgingCoeff = coeff;}
/**set the reference velocity for the nudging layer. */ /**set the reference velocity for the nudging layer. */
void setNudgingVelocity(const functor *nudgingVel){_nudgingVel = nudgingVel;} void setNudgingVelocity(const functor *nudgingVel, bool transport=false){
_nudgingVel = nudgingVel;
_nudgingVelIsTransport = transport;
}
/**use the moving bath WD (Karna et al, 2011) and define the bathymetry and the moving bathymetry */ /**use the moving bath WD (Karna et al, 2011) and define the bathymetry and the moving bathymetry */
void setMovingBathWettingDrying(double alphaMovingBathWD); void setMovingBathWettingDrying(double alphaMovingBathWD);
/** a factor to smooth the if in Hv term in the roe riemann solver */ /** 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