Commit 1d3b2e50 authored by Valentin Vallaeys's avatar Valentin Vallaeys
Browse files

interface : correct minor changes + add nudging in sw2d

parent 56ba37da
import slim
#slim.set_output_directory("output")
import dgpy
domain = slim.Domain("mesh.msh","bath_COMP_0.msh", g=9.81, order=1, solve_on_manifold=False)
#PREPRO domain.set_projection_UTM(55, a= 6370000, b= 6370000)
# PREPRO domain.set_projection("+proj=stere +ellps=WGS84 +lat_ts=69 +lat_0=90 +lon_0=85E +a=6370000 +b=6370000")
domain = slim.Domain("mesh.msh","bath_COMP_0.msh")
eq = slim.ShallowWater2d(domain, linear_equation = False, wetting_drying = None, export_every_sub_time_step = False)
eq.set_initial_condition("eta_init_COMP_0.msh")
eq = slim.ShallowWater2d(domain)
eq.set_initial_condition("eta_init_COMP_0.msh")
eq.set_viscosity("constant", constant_viscosity = 1e-6)
eq.set_dissipation("quadratic", quadratic_coefficient = 2.5e-3)
eq.set_coriolis("coriolis_COMP_0.msh")
......@@ -29,7 +24,7 @@ eq2.set_boundary_coast("Wall")
eq2.set_temporal_scheme("explicit")
loop=slim.Loop(initial_time="10 06 1991 22:00:00", final_time= "10 06 1991 23:00:00", maximum_time_step = 900, export_time = 3000)
loop=slim.Loop(initial_time="10 06 1991 22:00:00", final_time= "12 06 1991 23:00:00", maximum_time_step = 1000, export_time = 3600)
loop.add_equation(eq)
#loop.add_equation(eq2)
loop.set_output_directory("./output")
......
......@@ -145,8 +145,8 @@ class Domain:
density -- density of the liquid (default: 1025 kg/m^3)
solve_on_manifold -- boolean stating if you want to solve the shallow water equation on a 2d manifold instead of solving it on a plane (default: False)
"""
if order >10 ou order<1:
Msg.Fatal:("order of the finite elements functional basis must be larger than 1 and smaller than 11")
if order >10 or order<1:
Msg.Fatal("order of the finite elements functional basis must be larger than 1 and smaller than 11")
self._groups = dgGroupCollection(mesh, order)
self._bath = dgDofContainer(self._groups, 1)
_load(self._bath, bathy)
......@@ -195,7 +195,7 @@ class ShallowWater2d:
"""
if initial_time is None:
self._initial_time = 0.
elif isinstance(initial_time, str):
elif _is_string(initial_time):
_date = initial_time.split(' ')
_time = _date[3].split(':')
self._initial_time = calendar.timegm([int(_date[2]), int(_date[1]), int(_date[0]), int(_time[0]), int(_time[1]), float(_time[2])])
......@@ -203,7 +203,7 @@ class ShallowWater2d:
self._initial_time = float(initial_time)
if final_time is None:
self._final_time = 0.
elif isinstance(final_time, str):
elif _is_string(final_time):
_date = final_time.split(' ')
_time = _date[3].split(':')
self._final_time = calendar.timegm([int(_date[2]), int(_date[1]), int(_date[0]), int(_time[0]), int(_time[1]), float(_time[2])])
......@@ -373,7 +373,7 @@ class ShallowWater2d:
self._ref += [_wind_stress_x, _wind_stress_y, _wind_stress_z, _wind_stress_local_, _wind_stress_local_x, _wind_stress_local_y]
self._equation.setWindStress(self._wind_stress_local)
def set_wind_speed(self, density_air wind_speed_x, wind_speed_y, wind_speed_z= None):
def set_wind_speed(self, density_air, wind_speed_x, wind_speed_y, wind_speed_z= None):
"""
"""
......@@ -437,13 +437,33 @@ class ShallowWater2d:
self._ref += [_forcing_x, _forcing_y, _forcing_z, _forcing_local_, _forcing_local_x, _forcing_local_y]
self._equation.setSource(self._forcing_local)
def set_boundary_coast(self, physical_tag, slip):
def set_nudging(self, coefficient, external_sse, external_ux, external_uy, external_uz=None):
"""
Add a nudging term as a source term in the shallow water equation.
keyword argument:
coefficient -- netcdf or .msh file containing the nudging coefficient 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 s^-1]
external_sse -- netcdf or .msh file containing the external sea surface elevation 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]
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_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]
"""
self._coeff_nudging = _load_function(coefficient, self._domain._groups)
self._external_sse = _load_function(external_sse, self._domain._groups)
self._external_ux = _load_function(external_ux, self._domain._groups)
if self._domain._solve_on_manifold:
self._forcing_local_basis = _domain._space_transform.newFromGlobalFunctor(self._gravitational_potential)
self._equation.setSource(self._forcing_local_basis)
else:
self._equation.setSource(self._gravitational_potential)
def set_boundary_coast(self, physical_tag, slip=False):
"""
Boundary condition: normal speed set to zero and tangential speed set to zero or free (no slip or free slip condition)
keyword arguments:
physical_tag -- tag of the part of the boundary on which this boundary condition is applied.
slip -- tag for applying slip (True) or no slip (False) condition at the boundary
slip -- flag wheter applying slip (True) or no slip (False) condition at the boundary (Default: False)
"""
self._equation.addBoundaryCondition(physical_tag, self._equation.newBoundaryWall(slip))
......@@ -572,7 +592,7 @@ class ShallowWaterTracer2d:
"""
if initial_time is None:
self._initial_time = 0.
elif isinstance(initial_time, str):
elif _is_string(initial_time):
_date = initial_time.split(' ')
_time = _date[3].split(':')
self._initial_time = calendar.timegm([int(_date[2]), int(_date[1]), int(_date[0]), int(_time[0]), int(_time[1]), float(_time[2])])
......@@ -580,7 +600,7 @@ class ShallowWaterTracer2d:
self._initial_time = float(initial_time)
if final_time is None:
self._final_time = 0.
elif isinstance(final_time, str):
elif _is_string(final_time):
_date = final_time.split(' ')
_time = _date[3].split(':')
self._final_time = calendar.timegm([int(_date[2]), int(_date[1]), int(_date[0]), int(_time[0]), int(_time[1]), float(_time[2])])
......
......@@ -7,12 +7,8 @@
#include "dgMeshJacobian.h"
#include "slimMovingBathWettingDrying.h"
#include "functorMember.h"
#include "functor.h"
/*==============================================================================
* DG-terms
*============================================================================*/
......@@ -947,63 +943,63 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newForcedDischarge(dgDofCo
}
// BC : Coupling SW1D
class dgConservationLawShallowWater2d::boundaryCouplingSW1D:public dgBoundaryCondition {
class term:public function {
fullMatrix<double> normals, solExt, sol, bath, _xyz;
double _R, _g;
public:
term(const functor *solExtF, const functor *bathF,const double R, double g):function(3){
setArgument(solExt, solExtF);
setArgument(sol, function::getSolution());
setArgument(bath, bathF);
setArgument(normals, function::getNormals());
_R = R;
_g = g;
if (_R>0)
setArgument(_xyz, function::getCoordinates());
}
void call (dataCacheMap *map, fullMatrix<double> &val) {
double J = 1;
for (int i=0; i<sol.size1(); i++){
if (_R>0) {
double x = _xyz(i,0);
double y = _xyz(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
}
double nx = normals(i,0), ny = normals(i,1);
double uIn = nx * sol(i,1) + ny * sol(i,2);
double vIn = -ny * sol(i,1) + nx * sol(i,2);
double uExt = solExt(i,1);
double vExt = 0.0;
double etaIn = sol(i,0), etaExt = solExt(i,0);
double h = bath(i,0);
double HIn = h + etaIn;
double HExt = h + etaExt;
double uStar,vStar,HStar, AStar;
//printf("2D u = %.5e, %.5e H = %.5e %.5e eta %.5e %.5e \n",uIn, uExt, HIn, HExt, etaIn, etaExt);
roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, uStar, vStar, HStar, AStar, _g);
double etaStar = HStar - h;
double Feta = - HStar * uStar * J;
double Fun = - _g * etaStar * J - uStar * uStar * J;
double Fut = - uStar * vStar * J;
val(i,0) = Feta;
val(i,1) = Fun * nx - Fut * ny;
val(i,2) = Fun * ny + Fut * nx;
}
}
};
public:
term _term;
boundaryCouplingSW1D (dgConservationLawShallowWater2d *claw, const functor *solExtF) : _term(solExtF, claw->getBathymetry(), claw->getRadius(), claw->_g) {
_term0 = &_term;
}
};
dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryCouplingSW1D(const functor *solExtF){
return new boundaryCouplingSW1D(this, solExtF);
}
//
//class dgConservationLawShallowWater2d::boundaryCouplingSW1D:public dgBoundaryCondition {
// class term:public function {
// fullMatrix<double> normals, solExt, sol, bath, _xyz;
// double _R, _g;
// public:
// term(const functor *solExtF, const functor *bathF,const double R, double g):function(3){
// setArgument(solExt, solExtF);
// setArgument(sol, function::getSolution());
// setArgument(bath, bathF);
// setArgument(normals, function::getNormals());
// _R = R;
// _g = g;
// if (_R>0)
// setArgument(_xyz, function::getCoordinates());
// }
// void call (dataCacheMap *map, fullMatrix<double> &val) {
// double J = 1;
// for (int i=0; i<sol.size1(); i++){
// if (_R>0) {
// double x = _xyz(i,0);
// double y = _xyz(i,1);
// J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
// J = 1.0;
// }
// double nx = normals(i,0), ny = normals(i,1);
// double uIn = nx * sol(i,1) + ny * sol(i,2);
// double vIn = -ny * sol(i,1) + nx * sol(i,2);
// double uExt = solExt(i,1);
// double vExt = 0.0;
// double etaIn = sol(i,0), etaExt = solExt(i,0);
// double h = bath(i,0);
// double HIn = h + etaIn;
// double HExt = h + etaExt;
// double uStar,vStar,HStar, AStar;
// //printf("2D u = %.5e, %.5e H = %.5e %.5e eta %.5e %.5e \n",uIn, uExt, HIn, HExt, etaIn, etaExt);
// roeSolver(uIn, uExt, vIn, vExt, HIn, HExt, uStar, vStar, HStar, AStar, _g);
// double etaStar = HStar - h;
// double Feta = - HStar * uStar * J;
// double Fun = - _g * etaStar * J - uStar * uStar * J;
// double Fut = - uStar * vStar * J;
// val(i,0) = Feta;
// val(i,1) = Fun * nx - Fut * ny;
// val(i,2) = Fun * ny + Fut * nx;
// }
// }
// };
//public:
// term _term;
// boundaryCouplingSW1D (dgConservationLawShallowWater2d *claw, const functor *solExtF) : _term(solExtF, claw->getBathymetry(), claw->getRadius(), claw->_g) {
// _term0 = &_term;
// }
//};
//
//dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryCouplingSW1D(const functor *solExtF){
// return new boundaryCouplingSW1D(this, solExtF);
//}
/*==============================================================================
* IP Term
*============================================================================*/
......
......@@ -25,7 +25,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
class boundaryWall;
class boundaryForcedDischarge;
class boundaryCouplingSW1D;
//class boundaryCouplingSW1D;
bool _linear, _constantJac;
bool _useMovingBathWD;
......@@ -114,7 +114,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
/**slip wall boundary */
dgBoundaryCondition *newBoundaryWall(bool slip = false);
dgBoundaryCondition *newForcedDischarge(dgDofContainer *solution, const functor *discharge, std::string tag);
dgBoundaryCondition *newBoundaryCouplingSW1D(const functor *solExtF);
//dgBoundaryCondition *newBoundaryCouplingSW1D(const functor *solExtF);
};
/*==============================================================================
* IP Term
......
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