Enable simulations with no flux reconstruction.

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent d84f7025
......@@ -215,6 +215,13 @@ adding an empty line, make text **bold** or ``monospaced``.
<suggestion> 1 </suggestion>
<values> 1, 2, 3 </values>
</parameter>
<parameter name="fluxReconstruction">
<definition> Apply the flux reconstruction method to the solved matric
head to obtain conservative gradients for each element. </definition>
<suggestion> true </suggestion>
<values> true, false </values>
</parameter>
</category>
<category name="NewtonParameters">
......
......@@ -12,6 +12,25 @@
namespace Dune{
namespace Dorie{
template<class GO, Dune::GeometryType::BasicType gt, unsigned int order>
struct RaviartThomasFluxReconstructionHelper
{
private:
using GV = typename GO::Traits::TrialGridFunctionSpace::Traits::GridViewType;
static constexpr int dim = GV::dimension;
static_assert(dim==2||dim==3);
static_assert(gt==GeometryType::simplex||gt==GeometryType::cube);
// hardcoded bool helpers
static constexpr bool value_2d_simplex = (order <= 1);
static constexpr bool value_3d_simplex = false;
static constexpr bool value_2d_cube = (order <= 2);
static constexpr bool value_3d_cube = (order <= 1);
static constexpr bool value_simplex = (dim==2) ? value_2d_simplex : value_3d_simplex;
static constexpr bool value_cube = (dim==2) ? value_2d_cube : value_3d_cube;
public:
static constexpr bool enable = (gt==GeometryType::simplex) ? value_simplex : value_cube;
};
/*-------------------------------------------------------------------------*//**
* @brief Provide RT flux vector field from DG discretization. This code
* was adapted from Marian's code
......@@ -20,7 +39,7 @@ namespace Dorie{
* @tparam gt Geometry type of the grid.
* @tparam order Order of the RT reconstruction.
*/
template<class GO, Dune::GeometryType::BasicType gt, int order>
template<class GO, Dune::GeometryType::BasicType gt, unsigned int order>
class RaviartThomasFluxReconstruction
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<
......
......@@ -13,8 +13,8 @@ namespace Dorie{
template<class GV, class DF, class RF, std::size_t k>
class RaviartThomasSimplexLocalFiniteElementMap
: public Dune::PDELab::SimpleLocalFiniteElementMap<
Dune::Dorie::RaviartThomasSimplexLocalFiniteElement<GV::dimension,DF,RF>,
GV::dimension>
Dune::Dorie::RaviartThomasSimplexLocalFiniteElement<GV::dimension,DF,RF>,
GV::dimension>
{
static constexpr unsigned int dim = GV::dimension;
static constexpr unsigned int size_face = Dune::PB::PkSize<k,dim-1>::value;
......@@ -34,7 +34,8 @@ class RaviartThomasSimplexLocalFiniteElementMap
static_assert(dim == 2 || dim == 3,
"Raviart Thomas for simplices are only implemented for dimension 2 and 3!");
using FE = Dune::Dorie::RaviartThomasSimplexLocalFiniteElement<GV::dimension,DF,RF>;
using FE = Dune::Dorie::
RaviartThomasSimplexLocalFiniteElement<GV::dimension,DF,RF>;
using Base = Dune::PDELab::SimpleLocalFiniteElementMap<FE,GV::dimension>;
public:
RaviartThomasSimplexLocalFiniteElementMap(const GV& gv) : Base(_fe) {}
......@@ -80,14 +81,16 @@ RaviartThomasSimplexLocalFiniteElementMap<GV,DF,RF,k>::_fe
= Dune::Dorie::RaviartThomasSimplexLocalFiniteElement<GV::dimension,DF,RF>
(Dune::GeometryTypes::simplex(GV::dimension),k);
template<typename GV, int dim, typename DF, typename RF, std::size_t k, GeometryType::BasicType gt>
template<typename GV, int dim, typename DF, typename RF,
std::size_t k, GeometryType::BasicType gt>
struct RaviartThomasLocalFiniteElementMapBaseSelector
{
using type = Dune::PDELab::RaviartThomasLocalFiniteElementMap<GV,DF,RF,k,gt>;
};
template<typename GV, typename DF, typename RF>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,3,DF,RF,0,GeometryType::simplex>
struct RaviartThomasLocalFiniteElementMapBaseSelector<
GV,3,DF,RF,0,GeometryType::simplex>
{
using type = RaviartThomasSimplexLocalFiniteElementMap<GV,DF,RF,0>;
};
......@@ -107,8 +110,7 @@ struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,3,DF,RF,0,GeometryType:
* methods
*
* Available elements:
*
* @todo list of Available elements
*
* @warning This map should not be used to solve finite element problems. For
* that purpose use the
......@@ -136,9 +138,13 @@ public:
template<class GV, class DF, class RF, int k>
class RaviartThomasLocalFiniteElementMap<GV,DF,RF,k,Dune::GeometryType::simplex>
: public Dune::Dorie::RaviartThomasLocalFiniteElementMapBaseSelector<GV,GV::dimension,DF,RF,k,Dune::GeometryType::simplex>::type
: public Dune::Dorie::RaviartThomasLocalFiniteElementMapBaseSelector<
GV,GV::dimension,DF,RF,k,Dune::GeometryType::simplex>::type
{
using Base = typename Dune::Dorie::RaviartThomasLocalFiniteElementMapBaseSelector<GV,GV::dimension,DF,RF,k,Dune::GeometryType::simplex>::type;
using Base = typename Dune::Dorie::
RaviartThomasLocalFiniteElementMapBaseSelector<
GV,GV::dimension,DF,RF,k,Dune::GeometryType::simplex
>::type;
public:
RaviartThomasLocalFiniteElementMap(const GV& gv) : Base(gv) {}
};
......
......@@ -37,8 +37,7 @@ public:
/**
* @brief Local key for a given index.
*
* @param[in] i The index that identifies a shape function (see
* Dune::Dorie::VolumeRaviartThomasLocalBasis).
* @param[in] i The index that identifies a shape function.
*
* @return The local key associated to the index i. (see Dune:.LocalKey in
* dune-localfunctions)
......@@ -54,22 +53,6 @@ private:
std::vector<Dune::LocalKey> li;
};
/**
* \brief Raviart-Thomas local finite elements of arbitrary order
* for simplices of arbitrary dimension.
*
* These generic local finite elements are only available for
* simplex geometry types. The extension to cube and prism
* elements could be added.
*
* \ingroup RaviartThomas
*
* \tparam dimDomain dimension of reference elements
* \tparam D domain for basis functions
* \tparam R range for basis functions
* \tparam SF storage field for basis matrix
* \tparam CF compute field for basis matrix
*/
template<unsigned int dim, class D, class R,
class SF=R, class CF=SF>
class RaviartThomasSimplexLocalFiniteElement
......
......@@ -157,12 +157,12 @@ int main(int argc, char** argv)
sim.run();
break;
}
// case 3:{ *
// Sim<Simplex<2,3>> sim(richards_config, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 3:{ // No flux reconstruction available
Sim<Simplex<2,3>> sim(richards_config, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -228,24 +228,24 @@ int main(int argc, char** argv)
Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
auto grid_mapper = grid_creator.get_mapper();
switch(FEorder){
// case 1:{ *
// Sim<Simplex<3,1>> sim(richards_config, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{ *
// Sim<Simplex<3,2>> sim(richards_config, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 3:{ *
// Sim<Simplex<3,3>> sim(richards_config, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 1:{ // No flux reconstruction available
Sim<Simplex<3,1>> sim(richards_config, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{ // No flux reconstruction available
Sim<Simplex<3,2>> sim(richards_config, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{ // No flux reconstruction available
Sim<Simplex<3,3>> sim(richards_config, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -267,12 +267,12 @@ int main(int argc, char** argv)
sim.run();
break;
}
// case 3:{ *
// Sim<CubeAdaptive<3,3>> sim(richards_config, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 3:{ // No flux reconstruction available
Sim<CubeAdaptive<3,3>> sim(richards_config, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -291,11 +291,11 @@ int main(int argc, char** argv)
sim.run();
break;
}
// case 3:{ *
// Sim<Cube<3,3>> sim(richards_config, grid_mapper, helper);
// sim.run();
// break;
// }
case 3:{ // No flux reconstruction available
Sim<Cube<3,3>> sim(richards_config, grid_mapper, helper);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......
......@@ -12,4 +12,6 @@ using Geo = Dune::GeometryType::BasicType;
}
}
// A star "*" in the impl file is refers that it does not have flux reconstruction available.
#endif // DUNE_DORIE_RICHARDS_IMPL_HH
\ No newline at end of file
......@@ -9,7 +9,7 @@ namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>; *
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>; // *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/richards/impl/impl.hh>
#include <dune/dorie/model/richards/richards.cc>
......@@ -8,7 +9,7 @@ namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1>>; *
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1>>; // *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -9,7 +9,7 @@ namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>; *
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>; // *
} // namespace Dorie
} // namespace Dune
......@@ -8,8 +8,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3>>; *
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3>>; *
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3>>; //*
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3>>; // *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>; *
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>; // *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -133,7 +133,19 @@ void RichardsSimulation<Traits>::operator_setup()
std::cout << " Total number of DOF: " << DOFnumber << std::endl;
}
waterfrgf = std::make_shared<GFWaterFluxReconstruction>(*go0); //FIXME!
if constexpr (enable_fluxrc)
{
if (inifile.get<bool>("numerics.fluxReconstruction"))
{
waterfrgf = std::make_shared<GFWaterFluxReconstruction>(*go0);
}
} else {
if (inifile.get<bool>("numerics.fluxReconstruction"))
{
DUNE_THROW(Dune::NotImplemented,
"Flux reconstruction is not implemented for the selected configuration");
}
}
}
template<typename Traits>
......@@ -205,6 +217,7 @@ template<typename Traits>
void RichardsSimulation<Traits>::post_adapt_grid()
{
operator_setup();
update_adapters();
}
template<typename Traits>
......@@ -215,7 +228,9 @@ void RichardsSimulation<Traits>::update_adapters () const
condgf = std::make_shared<GFConductivity>(gv, fparam);
waterdgf = std::make_shared<GFWaterContent>(udgf, gv, fparam);
satdgf = std::make_shared<GFSaturation>(udgf, gv, fparam);
waterfrgf->update(*u);
if constexpr (enable_fluxrc)
if (inifile.get<bool>("numerics.fluxReconstruction"))
waterfrgf->update(*u);
}
template<typename Traits>
......@@ -223,6 +238,7 @@ void RichardsSimulation<Traits>::write_data () const
{
if (vtkwriter)
{
bool fluxrc = inifile.get<bool>("numerics.fluxReconstruction");
update_adapters();
if (inifile.get<bool>("output.vertexData")) {
......@@ -231,14 +247,18 @@ void RichardsSimulation<Traits>::write_data () const
vtkwriter->template addVertexData<GFConductivity>(condgf,"K_0");
vtkwriter->template addVertexData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addVertexData<GFSaturation>(satdgf,"Theta");
vtkwriter->template addVertexData<GFWaterFluxReconstruction>(waterfrgf,"flux (reconstruction)");
if constexpr (enable_fluxrc)
if (fluxrc)
vtkwriter->template addVertexData<GFWaterFluxReconstruction>(waterfrgf,"flux RT" + std::to_string(flux_order));
} else {
vtkwriter->template addCellData<GFMatricHead>(udgf,"head");
vtkwriter->template addCellData<GFWaterFlux>(fluxdgf,"flux");
vtkwriter->template addCellData<GFConductivity>(condgf,"K_0");
vtkwriter->template addCellData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addCellData<GFSaturation>(satdgf,"Theta");
vtkwriter->template addCellData<GFWaterFluxReconstruction>(waterfrgf,"flux (reconstruction)");
if constexpr (enable_fluxrc)
if (fluxrc)
vtkwriter->template addCellData<GFWaterFluxReconstruction>(waterfrgf,"flux RT" + std::to_string(flux_order));
}
try{
......
......@@ -43,21 +43,21 @@ namespace Dorie{
*
* @tparam BaseTraits Traits defining domain and range field properties of
* the simulation.
* @tparam order Order of the polynomial degree used for the basis
* @tparam k Order of the polynomial degree used for the basis
* functions in the DG method
*/
template<class BaseTraits, int order>
template<class BaseTraits, int k>
struct RichardsSimulationTraits : public BaseTraits
{
static constexpr int dim = BaseTraits::dim;
static constexpr int fem_order = order;
static constexpr int flux_order = order - 1; //FIXME!
static constexpr int order = k;
static constexpr int flux_order = order - 1;
using RF = typename BaseTraits::RF;
using Grid = typename BaseTraits::Grid;
using GV = typename BaseTraits::GV;
/// GFS Helper
using GFSHelper = Dune::Dorie::GridFunctionSpaceHelper<GV,RF,fem_order,BaseTraits::GridGeometryType>;
using GFSHelper = Dune::Dorie::GridFunctionSpaceHelper<GV,RF,order,BaseTraits::GridGeometryType>;
/// Problem GFS
using GFS = typename GFSHelper::Type;
/// GFS Constraints Container
......@@ -122,7 +122,8 @@ struct RichardsSimulationTraits : public BaseTraits
/// Saturation
using GFSaturation = Dune::Dorie::SaturationAdapter<BaseTraits,FlowParameters,GFMatricHead>;
/// Water Flux reconstruction
using GFWaterFluxReconstruction = Dune::Dorie::RaviartThomasFluxReconstruction<GO0,BaseTraits::GridGeometryType,flux_order>;
static constexpr bool enable_fluxrc = Dune::Dorie::RaviartThomasFluxReconstructionHelper<GO0,BaseTraits::GridGeometryType,flux_order>::enable;
using GFWaterFluxReconstruction = std::conditional_t<enable_fluxrc,Dune::Dorie::RaviartThomasFluxReconstruction<GO0,BaseTraits::GridGeometryType,flux_order>,void>;
// -- Utility Class Definitions -- //
/// Custom VTK output writer
......@@ -150,7 +151,9 @@ public:
private:
static constexpr int dim = Traits::dim;
static constexpr int order = Traits::fem_order;
static constexpr int order = Traits::order;
static constexpr int flux_order = Traits::flux_order;
static constexpr bool enable_fluxrc = Traits::enable_fluxrc;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
......
......@@ -9,7 +9,7 @@ except NameError:
pass
# paths set by cmake
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Release/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