The TS-GitLab will have to shut down in the near future — please plan migrating your projects to GitLab.com or GitHub. Contact @yunus for more information.

Commit cbe9b125 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

Fixed timestep bug


Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 22b1cc1b
......@@ -34,6 +34,26 @@ TransportSimulation<Traits>::TransportSimulation(
controller = std::make_unique<CalculationController>(inifile,*sboundary,helper);
std::string ts_method = _inifile.get<std::string>("numerics.timestepMethod");
if (ts_method== "explicit_euler")
ts_param = std::make_unique<Dune::PDELab::ExplicitEulerParameter<RF>>();
else if (ts_method== "implicit_euler")
ts_param = std::make_unique<Dune::PDELab::ImplicitEulerParameter<RF>>();
else if (ts_method== "heun")
ts_param = std::make_unique<Dune::PDELab::HeunParameter<RF>>();
else if (ts_method== "shu3")
ts_param = std::make_unique<Dune::PDELab::Shu3Parameter<RF>>();
else if (ts_method== "rk4")
ts_param = std::make_unique<Dune::PDELab::RK4Parameter<RF>>();
else if (ts_method== "alex2")
ts_param = std::make_unique<Dune::PDELab::Alexander2Parameter<RF>>();
else if (ts_method== "fractional")
ts_param = std::make_unique<Dune::PDELab::FractionalStepParameter<RF>>();
else if (ts_method== "alex3")
ts_param = std::make_unique<Dune::PDELab::Alexander3Parameter<RF>>();
else
DUNE_THROW(Dune::NotImplemented,"Time method not supported!");
if (output_policy() != OutputPolicy::None)
{
// // --- Utility Class Setup --- //
......@@ -81,12 +101,12 @@ void TransportSimulation<Traits>::operator_setup ()
// pdesolver->setVerbosityLevel(verbose);
// --- Time Step Operators ---
Dune::PDELab::OneStepMethodResult osm_result;
if(osm){
osm_result = osm->result(); // cache old result
}
osm = std::make_unique<OSM>(alex2,*igo,*pdesolver);
osm->setResult(osm_result);
// Dune::PDELab::OneStepMethodResult osm_result;
// if(osm){
// osm_result = osm->result(); // cache old result
// }
osm = std::make_unique<OSM>(*ts_param,*igo,*pdesolver);
// osm->setResult(osm_result);
osm->setVerbosityLevel(verbose+1);
gfs->update();
......@@ -114,11 +134,12 @@ void TransportSimulation<Traits>::step ()
{
if (operator_setup_flag)
operator_setup();
bool succed = false;
bool succeed = false;
bool exception = false;
while (not succed)
while (not succeed)
{
const RF t = controller->getTime();
this->suggest_timestep(controller->getDT());
// Check whether time intervals of grid functions are valid
if (Dune::FloatCmp::gt(_waterFluxInterval.begin(),t+_dt))
......@@ -145,6 +166,9 @@ void TransportSimulation<Traits>::step ()
if (not parmetis_warning)
dwarn.pop();
u = unext;
if (not ts_param->implicit())
this->suggest_timestep(slop->suggestTimestep(_dt));
}
catch (Dune::ISTLError &e){
exception = true;
......@@ -154,7 +178,7 @@ void TransportSimulation<Traits>::step ()
std::cerr << "FATAL ERROR occured!" << std::endl;
DUNE_THROW(Dune::Exception,e);
}
succed = true;
succeed = true;
}
// controller reacts to outcome of solution
controller->validate(exception);
......
......@@ -194,7 +194,7 @@ protected:
using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
/// Time stepping scheme
using ExplicitTimeStepScheme = Dune::PDELab::ExplicitEulerParameter<RF>;
using TimeSteppingParameter = Dune::PDELab::TimeSteppingParameterInterface<RF>;
Dune::MPIHelper& helper;
std::shared_ptr<Grid> grid;
......@@ -218,7 +218,7 @@ protected:
std::unique_ptr<LS> ls;
std::unique_ptr<PDESOLVER> pdesolver;
std::unique_ptr<OSM> osm;
ExplicitTimeStepScheme alex2;
std::unique_ptr<TimeSteppingParameter> ts_param;
std::shared_ptr<U> u;
......
......@@ -54,7 +54,7 @@ class TransportFVSpatialOperator
GridFunctionWaterFlux,
GridFunctionSaturation
>
>
>
, public Dune::PDELab::NumericalJacobianSkeleton<
TransportFVSpatialOperator<
Boundary,
......@@ -80,10 +80,10 @@ public:
enum { doPatternSkeleton = true };
// Residual assembly flags
enum { doLambdaBoundary = true };
enum { doAlphaVolume = true }; //
enum { doAlphaBoundary = true };
enum { doAlphaSkeleton = true };
enum { doLambdaVolume = true }; //
private:
static_assert(std::is_same<typename GridFunctionWaterFlux::Traits::GridViewType,
typename GridFunctionSaturation::Traits::GridViewType>::value,
......@@ -129,29 +129,8 @@ public:
* @tparam R The type for r
*/
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{}
/*---------------------------------------------------------------------*//**
* @brief Jacobian of volume term depending on test and ansatz
* functions
*
* @param[in] eg The entity of the grid
* @param[in] lfsu The ansatz local function space
* @param[in] x The coefficients of the lfsu
* @param[in] lfsv The test local function space
* @param mat The view of the jacobian matrix w.r.t lfsu
*
* @tparam EG The type for eg
* @tparam LFSU The type for lfsu
* @tparam X The type for x
* @tparam LFSV The type for lfsv
* @tparam M The type for mat
*/
template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, M& mat) const
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{}
/*---------------------------------------------------------------------*//**
......@@ -227,7 +206,7 @@ public:
// Inside unknown value
RF u = x_i(lfsu_i,0);
#if 0
#if 1
// Upwinding
u = (vn>=0) ? x_i(lfsu_i,0) : x_o(lfsu_o,0);
#endif
......@@ -241,15 +220,15 @@ public:
auto distance = global_inside.two_norm();
// Finite difference of u between the two entities
RF dudn = (x_i(lfsu_i,0)-x_o(lfsu_o,0))/distance;
RF dudn = (x_o(lfsu_o,0)-x_i(lfsu_i,0))/distance;
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_inside;
_gfSaturation->evaluate(ig.inside(),entity_inside_local,sat_inside);
// Symmetric contribution to residual on inside element
r_i.accumulate(lfsu_i, 0, ( u*vn + D_eff_avg*dudn*sat_inside )*face_volume );
r_o.accumulate(lfsu_o, 0, -( u*vn + D_eff_avg*dudn*sat_inside )*face_volume );
r_i.accumulate(lfsu_i, 0, ( u*vn - D_eff_avg*dudn*sat_inside )*face_volume );
r_o.accumulate(lfsu_o, 0, -( u*vn - D_eff_avg*dudn*sat_inside )*face_volume );
if (_first_stage)
_dtmin = std::min(_dtmin,suggestTimestepIntersection(ig));
......@@ -257,19 +236,24 @@ public:
/*---------------------------------------------------------------------*//**
* @brief Residual contribution of boundary integral
* @brief Boundary integral depending on test and ansatz functions.
*
* @param[in] ig The intersection entity of the grid (inside + outside
* entities)
* @param[in] lfsu_i The inside ansatz local function space
* @param[in] x_i The coefficients of the lfsu_i
* @param[in] lfsv_i The inside test local function space
* @param r_i The view of the residual vector w.r.t lfsu_i
*
* @tparam IG The type for ig
* @tparam LFSV The type for lfsv_i and
* @tparam LFSU The type for lfsu_i
* @tparam X The type for x_i
* @tparam LFSV The type for lfsv_i
* @tparam R The type for r_i
*/
template<typename IG, typename LFSV, typename R>
void lambda_boundary (const IG& ig, const LFSV& lfsv_i, R& r_i) const
template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
const LFSV& lfsv_i, R& r_i) const
{
// Get geometries
auto geo = ig.geometry();
......@@ -280,6 +264,10 @@ public:
if (BoundaryCondition::isDirichlet(bc))
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// Inside entity center
auto entity_inside = ig.inside();
auto geo_inside = entity_inside.geometry();
......@@ -308,9 +296,6 @@ public:
auto normal_face = ig.unitOuterNormal(face_local);
auto vn = flux_inside*normal_face;
// Contribution to residual from Dirichlet boundary
r_i.accumulate(lfsv_i, 0, (g*vn)*face_volume);
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_inside;
_gfSaturation->evaluate(ig.inside(),entity_inside_local,sat_inside);
......@@ -318,7 +303,11 @@ public:
// Evaluate diffusion coefficient from inside
auto D_eff_inside = _D_eff; // TODO: from parametrization
r_i.accumulate(lfsv_i, 0, -(D_eff_inside*g*sat_inside/distance)*face_volume);
// Inside unknown value
RF u = x_i(lfsu_i,0);
// Contribution to residual from Dirichlet boundary
r_i.accumulate(lfsu_i, 0, ((g*vn)-(D_eff_inside*sat_inside*(g-u)/distance))*face_volume);
}
else if (BoundaryCondition::isNeumann(bc))
{
......@@ -338,78 +327,55 @@ public:
// Contribution to residual from Neumann boundary
auto j = _boundary.j(ig.intersection(),face_local,_time);
r_i.accumulate(lfsv_i, 0, (D_eff_inside*j*sat_inside)*face_volume );
r_i.accumulate(lfsv_i, 0, -(D_eff_inside*j*sat_inside)*face_volume );
}
if (_first_stage)
_dtmin = std::min(_dtmin,suggestTimestepIntersection(ig));
}
/*---------------------------------------------------------------------*//**
* @brief Residual contribution of volume integral
*
* @param[in] eg THe entity of the grid
* @param[in] lfsu The ansatz local function space
* @param r The view of the residual vector w.r.t lfsu
*
* @tparam EG The type for eg
* @tparam LFSV The type for lfsv and
* @tparam R The type for r
*/
template<typename EG, typename LFSV, typename R>
void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
{
// // get cell
// const auto& cell = eg.entity();
// // cell center
// auto geo = eg.geometry();
// auto ref_el = referenceElement(geo);
// auto local_inside = ref_el.position(0,0);
// evaluate source and sink term
// auto f = param.f(cell,local_inside);
// r.accumulate(lfsv,0,-10000.*geo.volume());
}
/*---------------------------------------------------------------------*//**
* @brief Boundary integral depending on test and ansatz functions.
* @brief Residual contribution of boundary integral
*
* @param[in] ig The intersection entity of the grid (inside + outside
* entities)
* @param[in] lfsu_i The inside ansatz local function space
* @param[in] x_i The coefficients of the lfsu_i
* @param[in] lfsv_i The inside test local function space
* @param r_i The view of the residual vector w.r.t lfsu_i
*
* @tparam IG The type for ig
* @tparam LFSU The type for lfsu_i
* @tparam X The type for x_i
* @tparam LFSV The type for lfsv_i
* @tparam LFSV The type for lfsv_i and
* @tparam R The type for r_i
*/
template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
const LFSV& lfsv_i, R& r_i) const
{
// Get geometries
auto geo = ig.geometry();
auto ref_el = referenceElement(geo);
auto face_local = ref_el.position(0,0);
auto bc = _boundary.bc(ig.intersection(),face_local,_time);
if (BoundaryCondition::isDirichlet(bc))
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// Inside entity center
auto entity_inside = ig.inside();
auto geo_inside = entity_inside.geometry();
auto global_inside = geo_inside.center();
// Face center in global coordinates
auto global_face = geo.center();
// Compute distance of these two points
global_inside -= global_face;
auto distance = global_inside.two_norm();
// Inside unknown value
RF u = x_i(lfsu_i,0);
// Face volume for integration
auto face_volume = geo.volume();
// Get geometry of intersection in local coord. of neighbor entity
auto geo_in_inside = ig.geometryInInside();
auto entity_inside_local = geo_in_inside.global(face_local);
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_inside;
_gfSaturation->evaluate(ig.inside(),entity_inside_local,sat_inside);
// Evaluate diffusion coefficient from inside
auto D_eff_inside = _D_eff; // TODO: from parametrization
r_i.accumulate(lfsv_i, 0, (D_eff_inside*u*sat_inside/distance)*face_volume);
}
if (_first_stage)
_dtmin = std::min(_dtmin,suggestTimestepIntersection(ig));
}
template<typename IG, typename LFSV, typename R>
void lambda_boundary (const IG& ig, const LFSV& lfsv_i, R& r_i) const
{}
/**
* @brief Set up CFL timestep calculation in first stage
......@@ -606,15 +572,24 @@ public:
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
RF u = x(lfsu,0);
// entity geometry
auto geo = eg.geometry();
// inside cell center
auto center = eg.geometry().center();
auto center = geo.center();
typename GridFunctionSaturation::Traits::RangeType saturation;
_gfSaturation->evaluate(eg.entity(),center,saturation);
// update residual
r.accumulate(lfsv ,0 , saturation * x(lfsu,0));
r.accumulate(lfsv ,0 , saturation*u*geo.volume());
}
/*---------------------------------------------------------------------*//**
......
......@@ -14,7 +14,7 @@ boundary.file = "{_asset_path}/bcs/solute_2d.dat"
parameters.arrayFile = "{_asset_path}/parfields/sand.h5"
parameters.molecularDifussion = 2E-9
time.end = 1E7
time.end = 1E8
time.maxTimestep = 1E7
time.startTimestep = 1E4
......@@ -23,11 +23,14 @@ adaptivity.useAdaptivity = false
grid.FEorder = 0
grid.gridType = rectangular
grid.initialLevel = 0
grid.cells = 1 160, 1 40, 1 20 | expand prec
grid.cells = 160 160, 40 40, 20 20 | expand prec
[NewtonParameters]
AbsoluteLimit = 1E-10
Reduction = 1E-10
[dg]
penaltyFactor = 10
\ No newline at end of file
penaltyFactor = 10
[numerics]
timestepMethod = rk4
\ No newline at end of file
......@@ -79,17 +79,17 @@ add_custom_target(test_param
COMMAND ctest --output-on-failure --tests-regex ^.+test-parameterization_.+$
)
# TEMPORAL trasport test
dorie_add_metaini_test(UNIT_TEST
SOURCE test-transport.cc
BASENAME test-transport
CREATED_TARGETS test-transport
METAINI test-transport.mini.in
SCRIPT
)
add_custom_target(test-transport
COMMAND ctest --output-on-failure --tests-regex ^.+test-parameterization_.+$
)
# # TEMPORAL trasport test
# dorie_add_metaini_test(UNIT_TEST
# SOURCE test-transport.cc
# BASENAME test-transport
# CREATED_TARGETS test-transport
# METAINI test-transport.mini.in
# SCRIPT
# )
# add_custom_target(test-transport
# COMMAND ctest --output-on-failure --tests-regex ^.+test-parameterization_.+$
# )
# dune excludes test targets from 'make all'; undo that here where applicable
set_property(TARGET dorie PROPERTY EXCLUDE_FROM_ALL 0)
......
......@@ -2,5 +2,6 @@ spatial_resolution_north 2 0.25 0.75
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 dirichlet 0 dirichlet 1 dirichlet 0 dirichlet 0
\ No newline at end of file
number_BC_change_times 2
0 dirichlet 0 dirichlet 1 dirichlet 0 dirichlet 0
1E7 dirichlet 0 dirichlet 0 dirichlet 0 dirichlet 0
\ No newline at end of file
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