Allow adaptivity for simplices in coupled simulation.

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 9f452bf7
......@@ -113,6 +113,18 @@ struct RichardsSimulationTraits : public BaseTraits
/// Methods computing the time step
using OSM = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,U,U>;
/// Simulation state provider
struct State
{
using GridFunctionSpace = GFS;
using Coefficients = U;
// using TimeInterval = typename BaseTraits::TimeInterval;
std::shared_ptr<GridFunctionSpace> grid_function_space;
std::shared_ptr<Coefficients> coefficients;
// TimeInterval time_interval;
};
// -- Grid Functions of state variables -- //
/// Matric head
using GFMatricHead = Dune::PDELab::DiscreteGridFunction<GFS,U>;
......@@ -214,6 +226,9 @@ private:
/// Methods computing the time step
using OSM = typename Traits::OSM;
/// Simulation state provider
using State = typename Traits::State;
// -- Grid Functions of state valriables -- //
/// Matric head
using GFMatricHead = typename Traits::GFMatricHead;
......@@ -311,18 +326,18 @@ public:
/*------------------------------------------------------------------------*//**
* @brief Mark the grid in order to improve the current model.
*/
virtual void mark_grid() override;
void mark_grid() override;
/*------------------------------------------------------------------------*//**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
*/
virtual void adapt_grid() override;
void adapt_grid() override;
/*------------------------------------------------------------------------*//**
* @brief Setup operators to fit the new grid.
*/
virtual void post_adapt_grid() override;
void post_adapt_grid() override;
/*------------------------------------------------------------------------*//**
* @brief Method that provides the begin time of the model.
......@@ -345,6 +360,17 @@ public:
*/
double current_time() const final {return controller->getTime();}
/*------------------------------------------------------------------------*//**
* @brief Provides an object with all the information that describes
* the current state of the solved variable.
*
* @return Current state of the model.
*/
State current_state() {
return State{gfs,u/*,time_interval*/};
}
/*------------------------------------------------------------------------*//**
* @brief Suggest a time step to the model.
*
......
......@@ -169,7 +169,6 @@ public:
const auto& entity = eg.entity();
auto geo = entity.geometry();
const auto volume = geo.volume();
typename EG::Entity::Geometry::JacobianInverseTransposed jac;
// loop over quadrature points and integrate normal flux
......@@ -609,11 +608,6 @@ public:
for (unsigned int i = 0; i<lfsu_i.size(); i++)
u_i += x_i(lfsu_i,i) * phiu_i[i];
// evaluate velocity field assume H(div) velocity field => may choose any side
WaterFlux water_flux_i;
_gf_water_flux->evaluate(entity_i, position_i, water_flux_i);
auto water_flux_n = water_flux_i * normal_f;
if (BoundaryCondition::isDirichlet(bc))
{
......@@ -827,20 +821,20 @@ public:
using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits;
using LocalBasisTraitsV = typename LFSV::Traits::FiniteElementType::
Traits::LocalBasisType::Traits;
// using LocalBasisTraitsV = typename LFSV::Traits::FiniteElementType::
// Traits::LocalBasisType::Traits;
using RangeU = typename LocalBasisTraitsU::RangeType;
using RangeV = typename LocalBasisTraitsV::RangeType;
// using RangeV = typename LocalBasisTraitsV::RangeType;
using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
using RangeFieldV = typename LocalBasisTraitsV::RangeFieldType;
// using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
// using RangeFieldV = typename LocalBasisTraitsV::RangeFieldType;
using GradientU = Dune::FieldVector<RangeFieldU,dim>;
using GradientV = Dune::FieldVector<RangeFieldV,dim>;
// using GradientU = Dune::FieldVector<RangeFieldU,dim>;
// using GradientV = Dune::FieldVector<RangeFieldV,dim>;
using JacobianU = typename LocalBasisTraitsU::JacobianType;
using JacobianV = typename LocalBasisTraitsV::JacobianType;
// using JacobianU = typename LocalBasisTraitsU::JacobianType;
// using JacobianV = typename LocalBasisTraitsV::JacobianType;
// get polynomial degree
const int order = lfsu.finiteElement().localBasis().order();
......
......@@ -83,60 +83,35 @@ void RichardsTransportCouplingSimulation<Traits>::step()
write_data();
}
// template<typename Traits>
// void RichardsTransportCouplingSimulation<Traits>::adapt()
// {
// if constexpr (Traits::BaseTraits::adapt_grid)
// {
// if (!_grid)
// DUNE_THROW(Dune::Exception,"Dune grid is in invalid state!");
// // for some reason this is not being done in richards
// if(T::BaseTraits::GridGeometryType == Dune::GeometryType::BasicType::cube){
// auto closure(T::BaseTraits::Grid::ClosureType::NONE);
// _grid->setClosureType(closure);
// }
// _richards->mark_grid();
// using GFMatricHead = typename RichardsSimulationTraits::GFMatricHead;
// using GFSRichards = typename GFMatricHead::Traits::GridFunctionSpace;
// using XRichards = typename GFMatricHead::Traits::CoefficientVector;
// // if (_matricHeadStorage.size() > 2)
// // DUNE_THROW(Dune::Exception,
// // "DORiE do not handle more than two richards solutions when adapting.
// // This only happens when richards and transport simulations have similar timesteps.
// // Please wither increase timesteps in richards or decrease them in transport");
// auto it = _matricHeadStorage.begin();
// std::shared_ptr<GFMatricHead> gridFunctionPtr = it->second;
// std::shared_ptr<GFSRichards> gfs_richards = gridFunctionPtr->gridFunctionSpace();
// std::shared_ptr<XRichards> x_richards_1_ptr = gridFunctionPtr->coefficientVector(); // reference last solution
// std::shared_ptr<XRichards> x_richards_2_ptr;
// it++;
// if (it != _matricHeadStorage.end())
// gridFunctionPtr = it->second;
// // x_richards_2_ptr has to be adapted regartheless it contains real solution or not,
// // so here we create a fake one if it does not exist
// if (gridFunctionPtr)
// x_richards_2_ptr = gridFunctionPtr->coefficientVector();
// else
// x_richards_2_ptr = std::make_shared<XRichards>(*x_richards_1_ptr);
// int integration_order = 2;
// auto pack_richards = Dune::PDELab::transferSolutions(*gfs_richards,integration_order,*x_richards_1_ptr,*x_richards_2_ptr);
// auto solutePtr = _transport->getSolute();
// auto pack_transport = Dune::PDELab::transferSolutions(*(solutePtr->gridFunctionSpace()),integration_order,*(solutePtr->coefficientVector()));
// Dune::PDELab::adaptGrid(*_grid,pack_richards,pack_transport);
// }
// }
template<typename Traits>
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){
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
auto richards_state = _richards->current_state();
auto transport_state = _transport->current_state();
auto richards_int_order = RichardsSimulationTraits::fem_order + add_integration_order;
auto& richards_gfs = *(richards_state.grid_function_space);
auto& richards_coeff = *(richards_state.coefficients);
auto richards_pack = Dune::PDELab::transferSolutions(richards_gfs,richards_int_order,richards_coeff);
auto transport_int_order = TransportSimulationTraits::fem_order + add_integration_order;
auto& transport_gfs = *(transport_state.grid_function_space);
auto& transport_coeff = *(transport_state.coefficients);
auto transport_pack = Dune::PDELab::transferSolutions(transport_gfs,transport_int_order,transport_coeff);
Dune::PDELab::adaptGrid(*_grid,richards_pack,transport_pack);
}
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -80,7 +80,42 @@ public:
*/
void step() override;
// void adapt();
/**
* @brief Mark the grid in order to improve the current model.
*/
void mark_grid() override
{
if (adaptivity_policy() == AdaptivityPolicy::WaterFlux)
{
_richards->set_policy(adaptivity_policy());
_richards->mark_grid();
} else if (adaptivity_policy() != AdaptivityPolicy::None)
DUNE_THROW(Dune::NotImplemented,"Not known adaptivity policy!");
}
/**
* @brief Operations before adaptation of the grid
*/
void pre_adapt_grid() override
{
_richards->pre_adapt_grid();
_transport->pre_adapt_grid();
};
/**
* @brief Adapt the grid together it every dependency of the
* grid (e.g. solution vector and grid function spaces).
*/
void adapt_grid() override;
/**
* @brief Operations after adaptation of the grid
*/
void post_adapt_grid() override
{
_richards->post_adapt_grid();
_transport->post_adapt_grid();
};
/*------------------------------------------------------------------------*//**
* @brief Method that provides the begin time of the model.
......
......@@ -43,11 +43,11 @@ namespace Dorie{
* @tparam GFWaterFluxType Grid Function type of the water flux
* @tparam GFWaterContentType Grid Function type of the water content
*/
template<class BaseTraits, class GFWaterFluxType, class GFWaterContentType, unsigned int k = 0>
template<class BaseTraits, class GFWaterFluxType, class GFWaterContentType, unsigned int order = 0>
struct TransportSimulationTraits : public BaseTraits
{
static constexpr int dim = BaseTraits::dim;
static constexpr int order = k;
static constexpr int fem_order = order;
using RF = typename BaseTraits::RF;
using Grid = typename BaseTraits::Grid;
using GV = typename BaseTraits::GV;
......@@ -64,9 +64,9 @@ struct TransportSimulationTraits : public BaseTraits
/// GFS Helper
using GFSHelper =
std::conditional_t<order==0,
std::conditional_t<fem_order==0,
Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV,RF,BaseTraits::GridGeometryType>,
Dune::Dorie::GridFunctionSpaceHelper<GV,RF,order,BaseTraits::GridGeometryType>>;
Dune::Dorie::GridFunctionSpaceHelper<GV,RF,fem_order,BaseTraits::GridGeometryType>>;
/// Problem GFS
using GFS = typename GFSHelper::Type;
......@@ -82,7 +82,7 @@ struct TransportSimulationTraits : public BaseTraits
using SoluteInitial = Dune::Dorie::SoluteInitial<BaseTraits>;
/// Local spatial operator
using SLOP =
std::conditional_t<order==0,
std::conditional_t<fem_order==0,
Dune::Dorie::Operator::TransportFVSpatialOperator<SoluteBoundary,GFWaterFlux,GFWaterContent>,
Dune::Dorie::Operator::TransportDGSpatialOperator<SoluteBoundary,GFWaterFlux,GFWaterContent>>;
......@@ -111,6 +111,18 @@ struct TransportSimulationTraits : public BaseTraits
/// Methods computing the time step
using OSM = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,U,U>;
/// Simulation state provider
struct State
{
using GridFunctionSpace = GFS;
using Coefficients = U;
// using TimeInterval = typename BaseTraits::TimeInterval;
std::shared_ptr<GridFunctionSpace> grid_function_space;
std::shared_ptr<Coefficients> coefficients;
// TimeInterval time_interval;
};
// -- Grid Functions of state valriables -- //
/// Solute
using GFSolute = Dune::PDELab::DiscreteGridFunction<GFS,U>;
......@@ -147,7 +159,7 @@ public:
protected:
static constexpr int dim = Traits::dim;
static constexpr int order = Traits::fem_order; // TODO: FV when 0, DG if higher!
static constexpr int fem_order = Traits::fem_order;
using TimeField = typename Traits::TimeField;
using TimeInterval = typename Traits::TimeInterval;
using RF = typename Traits::RF;
......@@ -198,10 +210,13 @@ protected:
/// Methods computing the time step
using OSM = typename Traits::OSM;
/// Simulation state provider
using State = typename Traits::State;
// -- Grid Functions of state valriables -- //
/// Solute
/// Solute consetration
using GFSolute = typename Traits::GFSolute;
// Tota solute
using GFTotalSolute = typename Traits::GFTotalSolute;
// // -- Utility Class Definitions -- //
......@@ -345,17 +360,13 @@ public:
}
/*-----------------------------------------------------------------------*//**
* @brief Get a function for the solute for the current time of the
* model.
* @brief Provides an object with all the information that describes the
* current state of the solved variable.
*
* @return A pointer to a grid function of solute.
* @return Current state of the model.
*/
std::shared_ptr<GFSolute> get_solute()
{
if (!c_w)
DUNE_THROW(Dune::InvalidStateException,
"Pointer to c_w in invalid state!");
return c_w;
State current_state() {
return State{gfs,u/*,time_interval*/};
}
protected:
......
......@@ -11,7 +11,7 @@ except NameError:
pass
# paths set by cmake
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Release/dorie"
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
MPIEXEC = "/usr/bin/mpiexec"
MPIEXEC_NUMPROC_FLAG = "-n"
MPIEXEC_PREFLAG = ""
......
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