Commit 68676ca3 authored by David Vincent's avatar David Vincent
Browse files

gravity is a functor

parent 9685342e
......@@ -73,4 +73,4 @@ def tideDynColat(lonlat, t):
dObldColat = obl*p_const/R*np.cos(2.*colat)*np.cos(lon)*np.cos(omega*t) # 1/R *dphi_ob/dcolat
return -(dObldColat + dEccdColat)
return -(dObldColat + dEccdColat)*-1# *-1 because lat and colat are anti parallel
......@@ -4,7 +4,7 @@ import shutil
import functions
import param
do = ["forcing", "coriolis"]
do = ["gravity", "bathy", "forcing", "coriolis"]
check_result = True
......@@ -79,5 +79,13 @@ if "forcing" in do:
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fx', param.data_dir+'Fx',50*param.dt)
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fy', param.data_dir+'Fy',50*param.dt)
if "gravity" in do:
lonlat = lonlat_global.coordinates
nPt = len(lonlat)
g = np.ones(nPt) * param.g
slimPre.write_file(param.data_dir+'gravity.nc', region=region_global, time=None, data=[('gravity', g)])
if check_result:
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'gravity.nc', 'gravity', param.data_dir+'gravity')
print( 'Preprocessing done ' )
slimPre.exit(0)
import slim
import param
domain = slim.Domain(param.mesh,(param.data_dir+"bathymetry.nc","bathymetry"), g=param.g, density=param.rho)
domain = slim.Domain(param.mesh,(param.data_dir+"bathymetry.nc","bathymetry"), g= (param.data_dir+"gravity.nc","gravity"), density=param.rho)
eq = slim.ShallowWater2d(domain, initial_time=param.initial_time, final_time=param.final_time, export_every_sub_time_step = False,wetting_drying = 0.5)
eq = slim.ShallowWater2d(domain, initial_time=param.initial_time, final_time=param.final_time, export_every_sub_time_step = False, wetting_drying = 0.5)
eq.set_viscosity("smagorinsky")
eq.set_dissipation("manning")
eq.set_coriolis((param.data_dir+"coriolis.nc","coriolis"))
......@@ -12,8 +12,16 @@ eq.set_boundary_coast("Wall", slip=False)
eq.set_temporal_scheme("implicit")
eq.compute_mass("mass.dat")
eq2 = slim.ShallowWaterTracer2d(domain,name="density", offline = False, initial_time=param.initial_time, final_time=param.final_time)
eq2.set_diffusivity("okubo")
eq2.set_hydro_solution(equation=eq)
eq2.set_boundary_coast("Wall")
eq2.set_temporal_scheme("implicit")
eq2.compute_mass("tracer.dat")
loop=slim.Loop(maximum_time_step = param.dt, export_time = param.dt,path="./output")
loop.add_equation(eq)
#loop.add_equation(eq2)
#loop.export_on_structured_grid(["sw2d"], 4.7e5, 5e5, -3.23e6, -3.2e6, 0.2e5, 0.2e5)
loop.run()
......
......@@ -4,7 +4,7 @@ import shutil
import functions
import param
do = ["bathy", "forcing", "coriolis"]
do = ["gravity", "bathy", "forcing", "coriolis"]
check_result = True
......@@ -75,6 +75,13 @@ if "forcing" in do:
if check_result:
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fx', param.data_dir+'Fx',50*param.dt)
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fy', param.data_dir+'Fy',50*param.dt)
if "gravity" in do:
lonlat = lonlat_global.coordinates
nPt = len(lonlat)
g = np.ones(nPt) * param.g
slimPre.write_file(param.data_dir+'gravity.nc', region=region_global, time=None, data=[('gravity', g)])
if check_result:
slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'gravity.nc', 'gravity', param.data_dir+'gravity')
print( 'Preprocessing done ' )
slimPre.exit(0)
import slim
import param
domain = slim.Domain(param.mesh,(param.data_dir+"bathymetry.nc","bathymetry"), g=param.g, density=param.rho, solve_on_sphere = True)
domain = slim.Domain(param.mesh,(param.data_dir+"bathymetry.nc","bathymetry"), g= (param.data_dir+"gravity.nc","gravity"), density=param.rho, solve_on_sphere = True)
eq = slim.ShallowWater2d(domain, initial_time=param.initial_time, final_time=param.final_time, export_every_sub_time_step = False)
eq.set_viscosity("smagorinsky")
eq.set_dissipation("manning")
eq.set_coriolis((param.data_dir+"coriolis.nc","coriolis"))
eq.set_forcing(forcing_x=(param.data_dir+"forcing.nc","Fx"), forcing_y=(param.data_dir+"forcing.nc","Fy"))
eq.set_forcing(forcing_x=(param.data_dir+"forcing.nc","Fx"), forcing_y=(param.data_dir+"forcing.nc","Fy"), forcing_z=(param.data_dir+"forcing.nc","Fz"))
eq.set_boundary_coast("Wall", slip=False)
eq.set_temporal_scheme("implicit")
eq.compute_mass("mass.dat")
......
......@@ -49,7 +49,7 @@ import dgpy
class Domain:
"""Create the numerical domain"""
def __init__(self, mesh, bathy, g=9.81, order=1, density = 1025, solve_on_sphere = False):
def __init__(self, mesh, bathy, g=None, order=1, density = 1025, solve_on_sphere = False):
"""keyword arguments:
* mesh
......@@ -57,7 +57,7 @@ class Domain:
* bathy
bathymetric file [in meters, positive] (.msh, .idx or .nc format)
* g
mean gravitational acceleration (default: 9.81 m/s^2)
map of the mean gravitational acceleration (.msh, .idx or .nc format) (default: 9.81 m/s^2)
* order
order of the finite elements functional basis (default: 1, maximum = 10)
* density
......@@ -80,12 +80,16 @@ class Domain:
slim_private._load(self._bath, bathy)
self._bath_PC = dgpy.functionPrecomputed(self._groups, 3)
self._bath_PC.compute(self._bath.getFunction())
self._bath_PC.compute(self._bath.getFunction(),dgpy.dataCacheMap.NODE_GROUP_MODE)
self._bath_gradient = self._bath.getFunctionGradient()
self._bath_gradient_PC = dgpy.functionPrecomputed(self._groups, 3)
self._bath_gradient_PC.compute(self._bath.getFunctionGradient())
self._bath_function = self._bath_PC
self._bath_gradient_function = self._bath_gradient_PC
self._g = g
if g is None:
self._g = dgpy.functionConstant(9.81)
else:
self._g = slim_private._load_function(g,self._groups)
self._density = density
class ShallowWater2d:
......@@ -836,7 +840,8 @@ class Loop:
if equation._export_every_sub_time_step:
self._export_dofs_full.append([dgpy.dgIdxExporter(equation._solution, self._output_directory +'/export_for_offline_'+equation._name), equation._initial_time, equation._final_time])
self._number_of_equations += 1
self._ref += [equation]
self._ref += [equation]
def set_time_step_offline(self, dt, periodic = False, index_start=1, n_index_per_period=1):
"""Set the time step for the tracer equation if the hydrodynamics has already been computed
......
......@@ -102,18 +102,17 @@ class dgConservationLawShallowWater2d::maxDiffusiveSpeed: public function {
};
class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
fullMatrix<double> bath, sol, _xyz;
fullMatrix<double> bath, sol, _xyz, _g;
bool _linear, _filterLinear;
double _R;
double _g;
public:
maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear, double g):function(1){
maxConvectiveSpeed (const functor *bathymetry, bool linear, double R, bool filterLinear, const functor *g):function(1){
setArgument(bath,bathymetry);
setArgument(sol,function::getSolution());
_linear = linear;
_filterLinear = filterLinear;
_R = R;
_g = g;
setArgument(_g, g);
if(_R>0) {
setArgument(_xyz, function::getCoordinates());
}
......@@ -124,17 +123,17 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
val.setAll(0);
else{
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(bath(i,0)*_g);
val(i,0) = sqrt(bath(i,0)*_g(i,0));
}
}else
if (!_filterLinear){
for (int i=0; i<val.size1(); i++) {
val(i,0) = sqrt((bath(i,0)+sol(i,0))*_g)+hypot(sol(i,1),sol(i,2));
val(i,0) = sqrt((bath(i,0)+sol(i,0))*_g(i,0))+hypot(sol(i,1),sol(i,2));
}
}
else{
for (int i=0; i<val.size1(); i++)
val(i,0) = sqrt(sol(i,0)*_g)+hypot(sol(i,1),sol(i,2)); }
val(i,0) = sqrt(sol(i,0)*_g(i,0))+hypot(sol(i,1),sol(i,2)); }
/*if (_R>0)
......@@ -148,11 +147,11 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public function {
};
class dgConservationLawShallowWater2d::gradPsiTerm: public function {
fullMatrix<double> bath, sol, diff, _xyz, movingBathFactor;
fullMatrix<double> bath, sol, diff, _xyz, movingBathFactor, _g;
bool _linear, _isDiffusive;
double _R, _g;
double _R;
public:
gradPsiTerm (const functor *bathymetryF, const functor *diffusiveFluxF, bool linear, double R, const functor *movingBathFactorF, const functor *xyz, double g):function(9){
gradPsiTerm (const functor *bathymetryF, const functor *diffusiveFluxF, bool linear, double R, const functor *movingBathFactorF, const functor *xyz, const functor *g):function(9){
setArgument(bath,bathymetryF);
setArgument(sol,function::getSolution());
_isDiffusive = diffusiveFluxF;
......@@ -161,7 +160,7 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
}
_linear = linear;
_R = R;
_g = g;
setArgument(_g, g);
if(_R>0) {
setArgument(_xyz, xyz ? xyz : function::getCoordinates());
}
......@@ -170,7 +169,6 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
void call(dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = val.size1();
double wd_g = _g;
double J = 1.0;
for(size_t i=0; i< nQP; i++) {
if (_R > 0){
......@@ -188,12 +186,12 @@ class dgConservationLawShallowWater2d::gradPsiTerm: public function {
H+=eta;
// flux_x
val(i,0) = H*u * J/A;
val(i,1) = wd_g*eta*J;
val(i,1) = _g(i,0)*eta*J;
val(i,2) = 0;
// flux_y
val(i,3) = H*v * J/A;
val(i,4) = 0;
val(i,5) = wd_g*eta*J;
val(i,5) = _g(i,0)*eta*J;
if (!_linear) {
val(i,1) += u*u*J;
val(i,2) += u*v*J;
......@@ -217,17 +215,17 @@ 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, _nudgingCoeff, _nudgingVel;
fullMatrix<double> _bathymetry, _bathymetryGradient, _movingBathFactor, _movingBathFactorGradient, _nudgingCoeff, _nudgingVel, _g;
bool _linear, _isDiffusive, _useMovingBathWD, _absLayer, _nudgingVelIsTransport;
double _R, _g, _density;
double _R, _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 *bathymetryF, const functor *bathymetryGradientF, const functor *nuF, bool linear, double R, const functor *gravityF, double density,
const functor *movingBathFactorF, const functor *movingBathFactorGradientF, const functor *xyz,
const functor *nudgingCoeffF, const functor *nudgingVelF, bool nudgingVelIsTransport) :function (3){
_linear = linear;
_R = R;
_g = g;
setArgument(_g, gravityF);
_density = density;
_isDiffusive = nuF;
_useMovingBathWD = movingBathFactorF;
......@@ -281,8 +279,8 @@ class dgConservationLawShallowWater2d::source: public function {
val (i,1) = coriolisFactor(i,0)*v - linearDissipation(i,0) * u + _source(i,0) + _windStress(i,0)/((eta+_bathymetry(i,0))*_density);
val (i,2) = -coriolisFactor(i,0)*u - linearDissipation(i,0) * v + _source(i,1) + _windStress(i,1)/((eta+_bathymetry(i,0))*_density);
if (_R>0) {
val (i,1) -= _g*(eta)*_xyz(i,0)/2.0/_R/_R;
val (i,2) -= _g*(eta)*_xyz(i,1)/2.0/_R/_R;
val (i,1) -= _g(i,0)*(eta)*_xyz(i,0)/2.0/_R/_R;
val (i,2) -= _g(i,0)*(eta)*_xyz(i,1)/2.0/_R/_R;
}
if(_isDiffusive){
double H = _bathymetry(i,0);
......@@ -407,6 +405,7 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
const fullMatrix<double> &solR = cache.get(cache.solution(), 1);
const fullMatrix<double> &bathL = cache.get(*_bathymetry,0);
const fullMatrix<double> &bathR = cache.get(*_bathymetry,1);
const fullMatrix<double> &g = cache.get(*_g,0);
const fullMatrix<double> &normalsL = cache.get(*function::getNormals(),0);
const fullMatrix<double> &normalsR = cache.get(*function::getNormals(),1);
const fullMatrix<double> &movingBathFactorL = cache.get(*_movingBathFactor,0);
......@@ -471,8 +470,8 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
uR = uR * hR/hMin;
uL = uL * hL/hMin;
const double unM = (uR + uL) /2;
const double c = sqrt(_g * hMin);
Fun = J * (- _g * etaM ) + c * (uR - uL) / 2;
const double c = sqrt(g(i,0) * hMin);
Fun = J * (- g(i,0) * etaM ) + c * (uR - uL) / 2;
Feta = J * (- hMin * unM) + c * (etaR - etaL) / 2;
Fut = 0;
}
......@@ -487,8 +486,8 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
vL = vL * hL/hMin;
const double unM = (uR + uL) /2;
const double utM = (vR + vL) /2;
const double c = sqrt(_g * (hMin + etaM) + fabs(unM));
Fun = J * (- _g * etaM - unM * unM) + c * (uR - uL) / 2;
const double c = sqrt(g(i,0) * (hMin + etaM) + fabs(unM));
Fun = J * (- g(i,0) * etaM - unM * unM) + c * (uR - uL) / 2;
Feta = J * (- (hMin + etaM) * unM) + c * (etaR - etaL) / 2;
Fut = J * (-unM * utM) + c * (vR - vL) / 2;
}
......@@ -499,10 +498,10 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
double etaR = solR(i,0);
double etaL = solL(i,0);
const double hMin = std::min(hL, hR);
double sq_g_h = sqrt(_g/hMin);
double sq_g_h = sqrt(g(i,0)/hMin);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -_g*etaRiemann*J;
Fun = -g(i,0)*etaRiemann*J;
Fut = 0;
Feta = -(hMin)*unRiemann*J;
}else{
......@@ -512,9 +511,9 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
Msg::Error("Discontinuous Bathymetry is not implemented for non-linear equation with Roe Solver. Please use the linear equations and/or the LaxFriedrichs Flux.(left %g/%g right)", hL, hR);
if(HR<0 || HL<0) printf(" HR = %e HL =%e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g, movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
double eta = H-hMin;
Fun = -_g*eta*J - u*u*J;
Fun = -g(i,0)*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
}
......@@ -589,7 +588,7 @@ dgConservationLawShallowWater2d::dgConservationLawShallowWater2d() : dgConservat
_constantJac = false;
_coordinatesF = NULL;
_R = -1.0;
_g = 9.80616;
_g = new functionConstant(9.80616);
_density = 1025;
_nu = NULL;
_useMovingBathWD = false;
......@@ -670,6 +669,7 @@ void dgConservationLawShallowWater2d::setMovingBathWettingDrying(double alphaMov
}
dgConservationLawShallowWater2d::~dgConservationLawShallowWater2d() {
delete _g;
delete _fzero;
delete _fzerov;
if (_volumeTerm0[""]) delete _volumeTerm0[""];
......@@ -868,11 +868,11 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(bool slip)
class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBoundaryCondition {
class term : public function {
fullMatrix<double> normals, sol, _xyz, bath, discharge, mbFact;
fullMatrix<double> normals, sol, _xyz, bath, discharge, mbFact, _g;
bool _linear;
double _R, S0, _g;
double _R, S0;
public:
term(bool linear, double R, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, double g):function(3){
term(bool linear, double R, const functor *bathF, dgDofContainer *solution, const functor *dischargeF, std::string tag, const functor *movingBathF, const functor *g):function(3){
setArgument(sol,function::getSolution());
setArgument(_xyz,function::getCoordinates());
setArgument(normals,function::getNormals());
......@@ -881,7 +881,7 @@ class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBounda
setArgument(mbFact,movingBathF);
_linear = linear;
_R = R;
_g = g;
setArgument(_g, g);
dgFunctionIntegratorInterface integrator(solution->getGroups(), bathF, solution);
fullMatrix<double> Smat;
integrator.compute(tag, Smat);
......@@ -905,19 +905,19 @@ class dgConservationLawShallowWater2d::boundaryForcedDischarge : public dgBounda
double etaL = sol(i,0);
double etaR = etaL;
if (_linear) {
double sq_g_h = sqrt(_g/h);
double sq_g_h = sqrt(_g(i,0)/h);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -_g*etaRiemann*J;
Fun = -_g(i,0)*etaRiemann*J;
Fut = 0;
Feta = -(h)*unRiemann*J;
} else {
double HR = etaR + h, HL = etaL + h;
if(HR<0 || HL<0) printf("forced Discharge HR = %e HL =%e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g, mbFact(i,0), mbFact(i,0), -1);
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g(i,0), mbFact(i,0), mbFact(i,0), -1);
double eta = H-h;
Fun = -_g*eta*J - u*u*J;
Fun = -_g(i,0)*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -(h+eta)*u*J/Amb;
}
......
......@@ -34,12 +34,11 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
bool _laxFriedrichs;
bool _nudgingVelIsTransport;
double _R;
double _g;
double _density;
const functor *_bathymetry, *_bathymetryGradient;
const functor *_originalBathymetry, *_originalBathymetryGradient;
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_windStress, *_coriolisFactor, *_coordinatesF;
const functor *_linearDissipation, *_quadraticDissipation, *_source, *_windStress, *_coriolisFactor, *_coordinatesF, *_g;
const functor *_movingBathFactor, *_movingBathFactorGradient;
const functor *_fzero, *_fzerov, *_nu, *_diffusiveFlux, *_diffusion, *_ipTerm, *_xyz;
const functor *_nudgingCoeff, *_nudgingVel;
......@@ -65,13 +64,14 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
virtual bool isConstantJac() const {return (_constantJac || _constantJacIMEX);}
/**if this flag is true, a spherical version of the equations are solved*/
inline void setIsSpherical(double R) {_R = R;}
void setGravity(double g) {_g = g;}
void setDensity(double density) {_density = density;}
inline double getRadius() {return _R;}
void setCoordinatesFunction(functor *xyz);
/**set the function to evaluate the coriolis factor \f{eqnarray*} \frac{du}{dt} &= &-f v \\ \frac{dv}{dt} &=& f u\f} */
inline void setCoriolisFactor(const functor *coriolisFactor) { _coriolisFactor = coriolisFactor;}
/*set the gravity as a functor */
inline void setGravity(const functor *gravityF) { _g = gravityF;}
/**set the function to evaluate the linear dissipation \f{eqnarray*}\frac{du}{dt} &=& -\gamma u\\ \frac {dv}{dt} &=& -\gamma v \f}*/
inline void setLinearDissipation(const functor *linearDissipation) { _linearDissipation = linearDissipation;}
/**set the function to evaluate the quadratic dissipation \f{eqnarray*}\frac{du}{dt} &=& -c_d u||u||\\ \frac{dv}{dt}&= &-c_d v \f}*/
......
......@@ -72,6 +72,7 @@ void dgConservationLawShallowWaterTracer2d::riemannRoeCB(functorCache &m, fullMa
const fullMatrix<double> &hydroSolL = m.get(*_hydroSolution, 0);
const fullMatrix<double> &hydroSolR = m.get(*_hydroSolution, 1);
const fullMatrix<double> &bath = m.get(*_bathymetry, 0);
const fullMatrix<double> &g = m.get(*_g, 0);
const fullMatrix<double> &normalsL = m.get(*function::getNormals(), 0);
const fullMatrix<double> &normalsR = m.get(*function::getNormals(),1);
const fullMatrix<double> &movingBathFactorL = m.get(*_movingBathFactor, 0);
......@@ -108,7 +109,7 @@ void dgConservationLawShallowWaterTracer2d::riemannRoeCB(functorCache &m, fullMa
double HL = hydroSolL(i,0) + bath(i,0);
if(HR<0 || HL<0) printf("in sw2d tracer Riemann: HR = %e HL = %e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, _g, movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactorRiemann);
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, g(i,0), movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactorRiemann);
double Feta = -H*u/J/Amb;
val(i,0) = Feta * (uL > 0 ? cL : cR) ;
val(i,1) = -val(i,0);
......@@ -140,6 +141,7 @@ void dgConservationLawShallowWaterTracer2d::forcedBoundaryCB(functorCache &m, fu
const fullMatrix<double> &tracerIn = m.get(*function::getSolution(), 0);
const fullMatrix<double> &tracerExt = m.get(*tracerExtF);
const fullMatrix<double> &bath = m.get(*_bathymetry);
const fullMatrix<double> &g = m.get(*_g);
const fullMatrix<double> &movingBathFactor = m.get(*_movingBathFactor);
val.resize(m.nEvaluationPoint(), 1);
for (int i=0; i<m.nEvaluationPoint(); i++){
......@@ -154,7 +156,7 @@ void dgConservationLawShallowWaterTracer2d::forcedBoundaryCB(functorCache &m, fu
double HIn= etaIn + bath(i,0) , HExt = etaExt + bath(i,0);
if(HIn<0 || HExt<0) printf(" HIn = %e HExt =%e\n", HIn, HExt);
double u,v,H,Amb;
roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, u, v, H, Amb, _g, movingBathFactor(i,0), movingBathFactor(i,0), _upwindFactorRiemann);
roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, u, v, H, Amb, g(i,0), movingBathFactor(i,0), movingBathFactor(i,0), _upwindFactorRiemann);
val(i,0) = -H*u * ( uIn > 0 ? cIn : cExt );
}
}
......@@ -195,7 +197,7 @@ dgConservationLawShallowWaterTracer2d::dgConservationLawShallowWaterTracer2d() :
_ipTerm = NULL;
_diffusiveFlux = NULL;
_R = -1.0;
_g = 9.80616;
_g = new functionConstant(9.80616);
_alphaMovingBathWD = 0.0;
_useMovingBathWD = false;
_bathymetry = NULL;
......@@ -210,6 +212,7 @@ dgConservationLawShallowWaterTracer2d::~dgConservationLawShallowWaterTracer2d ()
delete fzero;
delete fzerov;
delete _movingBathFactor;
delete _g;
if (_volumeTerm0[""]) delete _volumeTerm0[""];
if (_volumeTerm1[""]) delete _volumeTerm1[""];
if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
......
......@@ -12,7 +12,8 @@ class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
const functor *_ipTerm;
const functor *_movingBathFactor;
const functor *fzero, *fzerov;
double _R, _g;
const functor *_g;
double _R;
functor *_bathymetry, *_originalBathymetry;
double _upwindFactorRiemann;
bool _useMovingBathWD;
......@@ -33,9 +34,10 @@ class dgConservationLawShallowWaterTracer2d : public dgConservationLawFunction {
bool isLinear() const {return _isLinear;}
void setIsLinear(bool isLinear) {_isLinear = isLinear;}
void setIsConstantJac(bool isConstantJac) {_isConstantJac = isConstantJac;}
void setGravity(double g) {_g = g;}
double gravity()const {return _g;}
/*set the gravity as a functor */
inline void setGravity(const functor *gravityF) { _g = gravityF;}
/**set the function to evaluate the source in the tracer equation */
inline void setSource(const functor *source1){_source = source1;}
/**set the hydrodynamic solution function (section, velocity) computed by the corresponding hydrodynamic equations */
......
......@@ -22,7 +22,7 @@ mesh.setSpaceTransform(st, groups)
alpha=0.
R = 6371220.
rho = 1000.
g = 9.80616
g = functionConstant(9.80616)
omega = 7.292e-5
u0 = 2.*np.pi*R/(12.*24.*3600.)
......@@ -85,12 +85,13 @@ vel_local = st.newFromGlobalFunctor(vel_CI) #computing the velocity in the local
def initialCondition(cm):
coord = cm.get(lonColat)
velocity =cm.get(vel_local)
gravity=cm.get(g)
r = coord[:,0]
lat = math.pi/2.-coord[:,1]*math.pi/180.
lon = coord[:,2]*math.pi/180.
nPt = coord.shape
val = np.zeros([nPt[0],3])
val[:,0]=(2.94e4-(R*omega*u0+u0*u0*0.5)*(-np.cos(lon)*np.cos(lat)*np.sin(alpha)+np.sin(lat)*np.cos(alpha))**2)/g
val[:,0]=(2.94e4-(R*omega*u0+u0*u0*0.5)*(-np.cos(lon)*np.cos(lat)*np.sin(alpha)+np.sin(lat)*np.cos(alpha))**2)/gravity[:,0]
val[:,1]=velocity[:,0]
val[:,2]=velocity[:,1]
return val
......
......@@ -1107,7 +1107,7 @@ void quadraticBottomFriction::operator()(functorCache &cm, fullMatrix<double> &v
}
}
manningBottomFriction::manningBottomFriction(double coeff, double gravity, const functor *bathy, const functor *solution) {
manningBottomFriction::manningBottomFriction(double coeff, const functor *gravity, const functor *bathy, const functor *solution) {
_coeff = coeff;
_gravity = gravity;
_bathy = bathy;
......@@ -1117,9 +1117,10 @@ manningBottomFriction::manningBottomFriction(double coeff, double gravity, const
void manningBottomFriction::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &bathy = cm.get(*_bathy);
const fullMatrix<double> &solution = cm.get(*_solution);
const fullMatrix<double> &gravity = cm.get(*_gravity);
val.resize(cm.nEvaluationPoint(),1,false);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
val(i,0) = _gravity * _coeff * _coeff/ pow((bathy(i,0)+solution(i,0)), 4./3.);
val(i,0) = gravity(i,0) * _coeff * _coeff/ pow((bathy(i,0)+solution(i,0)), 4./3.);
}
}
......
......@@ -294,11 +294,11 @@ class quadraticBottomFriction: public functor {
};
class manningBottomFriction: public functor {
const functor *_bathy, *_solution;
double _coeff, _gravity;
const functor *_bathy, *_solution, *_gravity;
double _coeff;
public:
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
manningBottomFriction(double coeff, double gravity, const functor *bathy, const functor *solution);
manningBottomFriction(double coeff, const functor *gravity, const functor *bathy, const functor *solution);
};
class mergeVector: public functor {
......
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