Commit 63d78946 authored by Santiago Ospina's avatar Santiago Ospina

Moved cfl condition outside of the local operator

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 19288ec4
......@@ -141,18 +141,21 @@ void TransportSimulation<Traits>::step ()
bool exception = false;
while (not succeed)
{
const RF t = controller->getTime();
this->suggest_timestep(controller->getDT());
auto t = controller->getTime();
suggest_timestep(controller->getDT());
// for some reason dt!=this->dt
const auto _dt = this->dt;
// Check whether time intervals of grid functions are valid
if (Dune::FloatCmp::gt(water_flux_interval.begin,t+dt))
if (Dune::FloatCmp::gt(water_flux_interval.begin,t+this->dt))
DUNE_THROW(Dune::Exception," ");
if (Dune::FloatCmp::gt(saturation_interval.begin,t+dt))
if (Dune::FloatCmp::gt(saturation_interval.begin,t+this->dt))
DUNE_THROW(Dune::Exception," ");
if (Dune::FloatCmp::lt(water_flux_interval.end,t+dt))
DUNE_THROW(Dune::Exception,water_flux_interval.end<< " < "<<t+dt);
if (Dune::FloatCmp::lt(saturation_interval.end,t+dt))
if (Dune::FloatCmp::lt(water_flux_interval.end,t+this->dt))
DUNE_THROW(Dune::Exception,water_flux_interval.end<< " < "<<t+this->dt);
if (Dune::FloatCmp::lt(saturation_interval.end,t+this->dt))
DUNE_THROW(Dune::Exception," ");
const bool solver_warnings = verbose > 0 && helper.rank() == 0 ?
......@@ -164,13 +167,11 @@ void TransportSimulation<Traits>::step ()
if (not solver_warnings)
dwarn.push(false);
osm->apply(t, dt, *u, *unext);
osm->apply(t, this->dt, *u, *unext);
if (not solver_warnings)
dwarn.pop();
u = unext;
if (not ts_param->implicit())
this->suggest_timestep(slop->suggestTimestep(dt));
}
catch (Dune::ISTLError &e){
exception = true;
......
......@@ -24,6 +24,7 @@
#include "../solver/util_controller.hh"
#include "../solver/util_writer.hh"
#include "../solver/util_grid_creator.hh"
#include "../solver/util_cfl_condition.hh"
namespace Dune{
namespace Dorie{
......@@ -297,7 +298,11 @@ public:
void suggest_timestep(TimeField _dt) override
{
dt = std::min(_dt,controller->getDT());
auto cfl_dt = 0.3*Dune::Dorie::cfl_condition(*gf_water_flux); // FIXME
std::cout << "cfl_dt: " << cfl_dt << std::endl;
this->dt = std::min(_dt,controller->getDT());
this->dt = std::min(dt,cfl_dt);
std::cout << "dt: " << cfl_dt << std::endl;
}
// TODO: Define data exchange with richards
......
......@@ -238,9 +238,6 @@ public:
// Symmetric contribution to residual on inside element
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume );
r_o.accumulate(lfsv_o, 0, -normal_flux*face_volume );
if (_first_stage)
_dt_min = std::min(_dt_min,suggest_timestep_intersection(ig));
}
......@@ -347,8 +344,6 @@ public:
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume );
}
if (_first_stage)
_dt_min = std::min(_dt_min,suggest_timestep_intersection(ig));
}
/*---------------------------------------------------------------------*//**
......@@ -382,38 +377,6 @@ public:
void lambda_boundary (const IG& ig, const LFSV& lfsv_i, R& r_i) const
{}
/**
* @brief Set up CFL timestep calculation in first stage
*
* @param[in] time The time
* @param[in] r Stage
*/
template<class RF>
void preStage(RF time, int r)
{
if (r == 1)
{
_first_stage = true;
_dt_min = std::numeric_limits<RF>::max();
}
else
_first_stage = false;
}
/**
* @brief Suggest timestep based on CFL condition
*/
template<class RF>
RF suggestTimestep(RF dt) const
{
if (dt*_dt_min > 0.)
return _gf_flux->getGridView().comm().max(_dt_min);
else
return -_gf_flux->getGridView().comm().max(_dt_min);
}
/*---------------------------------------------------------------------*//**
* @brief Sets the time.
*
......@@ -427,95 +390,12 @@ public:
_time = t;
}
protected:
/*---------------------------------------------------------------------*//**
* @brief suggest a timestep for an intersection
*
* @param[in] ig The intersection entity of the grid (inside + outside
* entities)
*
* @tparam IG The type for ig
*
* @return Maximum admissible timestep for ig
*/
template<class IG>
auto suggest_timestep_intersection(IG ig) const
{
// Get entity inside entity from intersection
const auto& entity_i = ig.inside();
// Get geometries
auto geo = ig.geometry();
auto geo_i = entity_i.geometry();
// Compute minimum local mesh width h
const DF face_volume = geo.volume();
const DF entity_i_volume = geo_i.volume();
DF distance_eff = entity_i_volume/face_volume;
typename WaterFlux::value_type max_velocity = 0.;
int intorder = 1; //FIXME
WaterFlux water_flux;
Saturation sat;
// loop over quadrature points
for (const auto& point : quadratureRule(geo_i, intorder))
{
const auto& x = point.position();
_gf_flux->evaluate(entity_i, x, water_flux);
_gf_sat->evaluate(entity_i,x, sat);
// Water velocity
const typename WaterFlux::value_type entity_max_velocity = water_flux.two_norm()/sat;
max_velocity = std::max(max_velocity,entity_max_velocity);
}
if (ig.neighbor()) {
// Get entity outside entity from intersection
const auto& entity_o = ig.outside();
// Get geometry
auto geo_o = entity_o.geometry();
// Compute minimum local mesh width h
const DF entity_o_volume = geo_o.volume();
distance_eff = (2.0/(1.0/(entity_i_volume+ 1E-30)
+ 1.0/(entity_o_volume+1E-30)))/face_volume;
// loop over quadrature points
for (const auto& point : quadratureRule(geo_o, intorder))
{
const auto& x = point.position();
_gf_flux->evaluate(entity_o, x, water_flux);
_gf_sat->evaluate(entity_o, x, sat);
// Water velocity
const typename WaterFlux::value_type entity_max_velocity = water_flux.two_norm()/sat;
max_velocity = std::max(max_velocity,entity_max_velocity);
}
}
return distance_eff/max_velocity;
}
private:
Boundary& _boundary;
const std::shared_ptr<GridFunctionWaterFlux> _gf_flux;
const std::shared_ptr<GridFunctionSaturation> _gf_sat;
double _time;
double _diff_coeff;
mutable bool _first_stage;
mutable double _dt_min;
};
/**
......@@ -588,7 +468,7 @@ public:
_gfSaturation->evaluate(eg.entity(),p_center,saturation);
// update residual
r.accumulate(lfsv ,0 , saturation*u*geo.volume());
r.accumulate(lfsu ,0 , saturation*u*geo.volume());
}
/*---------------------------------------------------------------------*//**
......
#ifndef DUNE_DORIE_CFL_CONDITION_HH
#define DUNE_DORIE_CFL_CONDITION_HH
#include <dune/common/float_cmp.hh>
#include <limits>
namespace Dune{
namespace Dorie{
/**
* @brief Computes the CFL condition for a grid function representing a
* velocity field. In particular, the return value is the
* suggested timestep for the lowest CFL ratio for each cell:
*
* @f{eqnarray*}{
* CFL = \min_{T\in\mathcal{T}_h}{||\beta||^{-1}_{L^inf(T)h^T} \\
* \delta t \lt \rho CFL
* @f}
* where $\rho$ is the courant number and $\beta$ is the velocity
* field.
*
* @param[in] gf Grid function representing the veolcity field
* $\beta$.
* @param[in] intorder The integration order. This value determines the
* quadrature points to evaluate
* $||\beta||_{L^inf(T)h^T}$.
*
* @tparam GF The grid function type.
* @tparam TF The time field. Type to represent time values.
*
* @return { description_of_the_return_value }
*/
template<class GF, class TF = double>
TF cfl_condition(const GF& gf, unsigned int intorder = 1)
{
using Traits = typename GF::Traits;
using DomainField = typename Traits::DomainFieldType;
using RangeField = typename Traits::RangeFieldType;
const auto& gv = gf.getGridView();
TF suggested_dt = std::numeric_limits<TF>::max();
// loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
auto geo = cell.geometry();
// Check that cell is not a vertex
assert(geo.corners()>1);
DomainField h_T = 0.;
// Cell diameter: maximum distance between all vertices
for (int i = 1; i < geo.corners(); ++i)
for (int j = 0; j < i-1; ++j)
h_T = std::max(h_T,(geo.corner(i)-geo.corner(j)).two_norm());
std::cout << "h_T: " << h_T << std::endl;
// Check that h_T is not zero in debug mode, who knows.
assert(Dune::FloatCmp::gt(h_T,0.));
RangeField max_water_flux = 0.;
// Find the maximum velocity in the cell
for (const auto& point : quadratureRule(geo, intorder))
{
typename Traits::RangeType water_flux;
const auto& x = point.position();
gf.evaluate(cell,x,water_flux);
max_water_flux = std::max(max_water_flux,water_flux.two_norm());
}
std::cout << "max_water_flux: " << max_water_flux << std::endl;
// Calculate the cfl for this cell
if (Dune::FloatCmp::gt(max_water_flux,0.))
suggested_dt = std::min(suggested_dt,h_T/max_water_flux);
std::cout << "suggested_dt: " << suggested_dt << std::endl;
}
std::cout << "FINAL suggested_dt: " << suggested_dt << std::endl;
// Communicate the lowest cfl number to all procesors
suggested_dt = gv.comm().min(suggested_dt);
return suggested_dt;
}
}
}
#endif
\ 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