[Coupling] Allow adapt the grid depending on the policy

parent 0b5884e6
......@@ -29,7 +29,7 @@ enum class OutputPolicy {None,EndOfStep,All};
* @brief Enum class for output policy. It defines which variable
* is the target in order to mark the grid
*/
enum class AdaptivityPolicy {None,WaterFlux};
enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
/// Return the estimation of entries per matrix row for the spatial GridOperator
/** This supposedly decreases matrix assembly time.
......
#include <dune/dorie/model/transport/richards_coupling.hh>
// #include <dune/pdelab/adaptivity/adaptivity.hh>
#include <dune/pdelab/adaptivity/adaptivity.hh>
namespace Dune{
namespace Dorie{
......@@ -42,12 +42,10 @@ RichardsTransportCouplingSimulation<Traits>::RichardsTransportCouplingSimulation
else
DUNE_THROW(NotImplemented,"not known output!");
// Create richards and transport simulations
// create richards simulations
_richards = std::make_unique<RichardsSimulation>(_inifile_richards,grid_creator,helper);
_transport = std::make_unique<TransportSimulation>(_inifile_transport,grid_creator,helper);
_richards->set_policy(richards_op);
_transport->set_policy(transport_op);
auto time_begin = _richards->begin_time();
......@@ -61,8 +59,9 @@ RichardsTransportCouplingSimulation<Traits>::RichardsTransportCouplingSimulation
auto gf_water_flux = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux,time_begin);
_transport->set_water_flux(igf_water_flux);
_transport->set_water_content(igf_water_content);
// create transport simulations
_transport = std::make_unique<TransportSimulation>(_inifile_transport,grid_creator,helper,igf_water_flux,igf_water_content);
_transport->set_policy(transport_op);
}
template<typename Traits>
......@@ -70,32 +69,32 @@ void RichardsTransportCouplingSimulation<Traits>::step()
{
TimeField time_begin = _richards->current_time();
richards_state_before_step = _richards->current_state();
// store richards sate in case of intermediate adaptivity
if (adaptivity_policy() == AdaptivityPolicy::SoluteFlux)
richards_state_before_step = _richards->current_state();
// Create instationary grid function containers
// create instationary grid function containers
auto igf_water_content = std::make_shared<GFWaterContentContainer>();
auto igf_water_flux = std::make_shared<GFWaterFluxContainer>();
// Set initial state of the water content to its container
auto gf_water_content_begin = _richards->get_water_content(richards_state_before_step);
// set initial state of the water content to container
auto gf_water_content_begin = _richards->get_water_content();
igf_water_content->push(gf_water_content_begin,time_begin);
// Set initial state of the water flux to its container
// set initial state of the water flux to container
auto gf_water_flux_begin = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_begin,time_begin);
// Always do an step of richards
// always do an step of richards
_richards->step();
TimeField time_end = _richards->current_time();
richards_state_after_step = _richards->current_state();
// Set final state of the water content to its container
auto gf_water_content_end = _richards->get_water_content(richards_state_after_step);
// set final state of the water content to container
auto gf_water_content_end = _richards->get_water_content();
igf_water_content->push(gf_water_content_end,time_end);
// Set final state of the water flux to its container
// set final state of the water flux to container
auto gf_water_flux_end = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_end,time_end);
......@@ -108,20 +107,21 @@ void RichardsTransportCouplingSimulation<Traits>::step()
return;
}
// Set instatonary grid functions to the transport solver
// set instatonary grid functions to the transport solver
_transport->set_water_flux(igf_water_flux);
_transport->set_water_content(igf_water_content);
// solve transport until synchronization with richards
while (Dune::FloatCmp::lt(_time_transport,_time_richards))
{
_time_transport = _transport->current_time();
_time_richards = _richards->current_time();
_transport->suggest_timestep(_time_richards - _time_transport);
_transport->step();
if (adaptivity_policy() == AdaptivityPolicy::SoluteFlux)
// intermediate adapt if policy targets solute
if (adaptivity_policy() == AdaptivityPolicy::SoluteFlux and
Dune::FloatCmp::lt(_transport->current_time(),_transport->end_time()))
{
// if ( Dune::FloatCmp::lt( _transport->current_time(), _time_richards) )
// {
mark_grid();
adapt_grid();
post_adapt_grid();
......@@ -129,7 +129,7 @@ void RichardsTransportCouplingSimulation<Traits>::step()
}
_current_time = _time_richards;
// Drop the old grid functions that the transport simulation holds.
// drop the old grid functions that the transport simulation have.
igf_water_content->pop();
igf_water_flux->pop();
if (_transport->output_policy() == OutputPolicy::EndOfStep)
......@@ -141,37 +141,47 @@ void RichardsTransportCouplingSimulation<Traits>::adapt_grid()
{
if (!_grid)
DUNE_THROW(Dune::InvalidStateException,"Dune grid is in invalid state!");
// for some reason this is not being done in richards
if(Traits::BaseTraits::GridGeometryType == Dune::GeometryType::BasicType::cube){
if(Traits::BaseTraits::GridGeometryType
== Dune::GeometryType::BasicType::cube){
DUNE_THROW(Dune::NotImplemented,"Coupling between richards and transport "
<< "is not implemented for grids with cubes due to the Raviart Thomas "
<< "flux reconstruction!");
}
int add_integration_order = 2; //FIXME
int add_int_order = 2;
// get current state of both models
auto richards_state = _richards->current_state();
auto transport_state = _transport->current_state();
auto richards_int_order = RichardsSimulationTraits::order + add_integration_order;
// get gfs and DOF vector of richards
auto richards_int_order = RichardsSimulationTraits::order + add_int_order;
auto& richards_gfs = *(richards_state.grid_function_space);
auto& richards_coeff = *(richards_state.coefficients);
auto transport_int_order = TransportSimulationTraits::order + add_integration_order;
// get gfs and DOF vector of transport
auto transport_int_order = TransportSimulationTraits::order + add_int_order;
auto& transport_gfs = *(transport_state.grid_function_space);
auto& transport_coeff = *(transport_state.coefficients);
// pack transport state
auto transport_pack = Dune::PDELab::transferSolutions(transport_gfs,transport_int_order,transport_coeff);
if (adaptivity_policy() == AdaptivityPolicy::WaterFlux)
{
{ // adaptivity on water flux assumes that simulations are synchronized
// pack richards state
auto richards_pack = Dune::PDELab::transferSolutions(richards_gfs,richards_int_order,richards_coeff);
// adapt grid and packs
Dune::PDELab::adaptGrid(*_grid,richards_pack,transport_pack);
}
else if (adaptivity_policy() == AdaptivityPolicy::SoluteFlux)
{
{ // adaptivity on solute flux assumes that simulations are not synchronized
// get DOF vector richards at begging of step
auto& richards_coeff_before = *(richards_state_before_step.coefficients);
// pack richards current and before state
auto richards_pack = Dune::PDELab::transferSolutions(richards_gfs,richards_int_order,richards_coeff,richards_coeff_before);
// adapt grid and packs
Dune::PDELab::adaptGrid(*_grid,richards_pack,transport_pack);
}
}
......@@ -179,31 +189,34 @@ void RichardsTransportCouplingSimulation<Traits>::adapt_grid()
template<typename Traits>
void RichardsTransportCouplingSimulation<Traits>::post_adapt_grid()
{
// if water state was adapted, containers have to be reset
if (adaptivity_policy() == AdaptivityPolicy::SoluteFlux)
{
// Create instationary grid function containers
// create instationary grid function containers
auto igf_water_content = std::make_shared<GFWaterContentContainer>();
auto igf_water_flux = std::make_shared<GFWaterFluxContainer>();
// Set initial state of the water content to its container
// set initial state of the water content to container
auto gf_water_content_begin = _richards->get_water_content(richards_state_before_step);
igf_water_content->push(gf_water_content_begin,richards_state_before_step.time);
// Set initial state of the water flux to its container
auto gf_water_flux_begin = _richards->get_water_flux_reconstructed(richards_state_before_step);
// set initial state of the water flux to container
auto gf_water_flux_begin = _richards->get_water_flux_reconstructed(richards_state_before_step);
igf_water_flux->push(gf_water_flux_begin,richards_state_before_step.time);
// Set final state of the water content to its container
// set final state of the water content to container
auto gf_water_content_end = _richards->get_water_content(richards_state_after_step);
igf_water_content->push(gf_water_content_end,richards_state_after_step.time);
// Set final state of the water flux to its container
auto gf_water_flux_end = _richards->get_water_flux_reconstructed(richards_state_after_step);
// set final state of the water flux to container
auto gf_water_flux_end = _richards->get_water_flux_reconstructed(richards_state_after_step);
igf_water_flux->push(gf_water_flux_end,richards_state_after_step.time);
_transport->set_water_flux(igf_water_flux);
_transport->set_water_content(igf_water_content);
}
// post adapt in richards and transport (e.g. operator setup)
_richards->post_adapt_grid();
_transport->post_adapt_grid();
}
......
......@@ -326,8 +326,10 @@ public:
Dune::MPIHelper& _helper,
std::shared_ptr<GFWaterFlux> _water_flux_gf,
std::shared_ptr<GFWaterContent> _water_content_gf
) : TransportSimulation(_inifile,_grid_creator,helper)
) : TransportSimulation(_inifile,_grid_creator,_helper)
{
set_water_flux(_water_flux_gf);
set_water_content(_water_content_gf);
operator_setup();
}
......
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