Commit 9936351b authored by Santiago Ospina's avatar Santiago Ospina

Merge branch '62-gradientfluxadapter-problem' into...

Merge branch '62-gradientfluxadapter-problem' into 101-couple-the-transportsimulation-with-richardssimulation
parents e58e7707 95d0a970
......@@ -97,6 +97,22 @@ Executing the Program
Analyzing Results
-----------------
Fluxes in DORiE
===============
Depending on the application for what you want to use DORiE, there might be the case that it is needed something more elaborated than the default fluxes that DORiE provides. Therefore, this section explains how to interpret the default fluxes and the flux reconstruction technique used to solve solute transport problems.
Understanding the water flux output
-----------------------------------
Firstly, we have to recall that DORiE solves a Discontinuous Galerking finite element problem with *matric head* as unknown. It means that the solution of the *matric head* (and therefore the *water flux*) is continuous only element-wise, or in other words, it is discontinuous on the intersections between elements. On the other hand, the dG method solves *numerical fluxes* on the intersections between elements composed together with a penalty term that increases with respect to the discontinuity of the *matric head*. This ensures that the local solution is conservative while keeps the discontinuity as low as possible.
From the description above one can infer that one has to distinguish between *water fluxes* at the interior of each element and at the intersections of all elements (we call these intersections skeleton of the grid). Unfortunately, there is no a standard form to write the skeleton fluxes on formats like VTK and that's the main reason why DORiE only provides the interior fluxes. However, assuming one can write both fluxes into some output format, they are still discontinuous (notice that direct use of discontinuous fluxes are useless for conservative computations since the transported quantities are very likely to get stagnated or over-transported in the nearby of intersections between elements). It means that it is needed some sort of post-processing that ensures that the *water mass* is still locally and globally conserved.
Flux reconstruction
-------------------
TODO
Examples
========
......
......@@ -5,6 +5,9 @@ if(dune-testtools_FOUND)
add_subdirectory("test")
endif()
add_executable("transport" transport.cc)
dune_target_link_libraries(transport dorie-impl ${DUNE_LIBS})
add_executable("dorie" dorie.cc)
dune_target_link_libraries(dorie dorie-impl ${DUNE_LIBS})
......
......@@ -152,12 +152,12 @@ int main(int argc, char** argv)
sim.run();
break;
}
case 3:{
Sim<Simplex<2,3>> sim(inifile, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 3:{ *
// Sim<Simplex<2,3>> sim(inifile, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -223,24 +223,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(inifile, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<Simplex<3,2>> sim(inifile, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<Simplex<3,3>> sim(inifile, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 1:{ *
// Sim<Simplex<3,1>> sim(inifile, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{ *
// Sim<Simplex<3,2>> sim(inifile, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 3:{ *
// Sim<Simplex<3,3>> sim(inifile, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -262,12 +262,12 @@ int main(int argc, char** argv)
sim.run();
break;
}
case 3:{
Sim<CubeAdaptive<3,3>> sim(inifile, grid_mapper, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 3:{ *
// Sim<CubeAdaptive<3,3>> sim(inifile, grid_mapper, helper);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -286,11 +286,11 @@ int main(int argc, char** argv)
sim.run();
break;
}
case 3:{
Sim<Cube<3,3>> sim(inifile, grid_mapper, helper);
sim.run();
break;
}
// case 3:{ *
// Sim<Cube<3,3>> sim(inifile, grid_mapper, helper);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1,0>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2,0>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
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 CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3,0>>;
// template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3,0>>; *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
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 CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1,0>>;
// template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1,0>>; *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
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 CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2,0>>;
// template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2,0>>; *
} // namespace Dorie
} // namespace Dune
......@@ -3,13 +3,13 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
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 CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3,0>>; *
// template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3,0>>; *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2>>;
template class CoupledSimulation<CoupledSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -3,12 +3,12 @@
#endif
#include <dune/dorie/impl/impl.hh>
#include <dune/dorie/interface/richards_simulation.cc>
#include <dune/dorie/interface/coupled_simulation.cc>
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>;
// template class CoupledsSimulation<CoupledsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3,0>>; *
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -7,7 +7,7 @@
namespace Dune{
namespace Dorie{
template<class BaseTraits, int RichardsOrder, int TransportOrder = RichardsOrder>
template<class BaseTraits, int RichardsOrder, int TransportOrder>
struct CoupledSimulationTraits : public BaseTraits
{
private:
......
......@@ -134,6 +134,8 @@ void RichardsSimulation<Traits>::operator_setup()
if (helper.rank()==0)
std::cout << " Total number of DOF: " << DOFnumber << std::endl;
}
waterfrgf = std::make_shared<GFWaterFluxReconstruction>(*go0); //FIXME!
}
template<typename Traits>
......@@ -215,6 +217,7 @@ 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);
}
template<typename Traits>
......@@ -223,19 +226,23 @@ void RichardsSimulation<Traits>::write_data () const
if (vtkwriter)
{
update_adapters();
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->addVertexData(udgf,"head");
vtkwriter->addVertexData(fluxdgf,"flux");
vtkwriter->addVertexData(condgf,"K_0");
vtkwriter->addVertexData(waterdgf,"theta_w");
vtkwriter->addVertexData(satdgf,"Theta");
vtkwriter->template addVertexData<GFMatricHead>(udgf,"head");
vtkwriter->template addVertexData<GFWaterFlux>(fluxdgf,"flux");
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)");
} else {
vtkwriter->addCellData(udgf,"head");
vtkwriter->addCellData(fluxdgf,"flux");
vtkwriter->addCellData(condgf,"K_0");
vtkwriter->addCellData(waterdgf,"theta_w");
vtkwriter->addCellData(satdgf,"Theta");
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)");
}
try{
vtkwriter->write(controller->getTime(),output_type);
}
......
......@@ -29,6 +29,7 @@
#include "../solver/adapters/saturation.hh"
#include "../solver/adapters/water_flux.hh"
#include "../solver/adapters/conductivity.hh"
#include "../solver/flux_reconstruction/rt_projection.hh"
namespace Dune{
......@@ -49,6 +50,7 @@ struct RichardsSimulationTraits : public BaseTraits
{
static constexpr int dim = BaseTraits::dim;
static constexpr int fem_order = order;
static constexpr int flux_order = order - 1; //FIXME!
using RF = typename BaseTraits::RF;
using Grid = typename BaseTraits::Grid;
using GV = typename BaseTraits::GV;
......@@ -118,7 +120,8 @@ struct RichardsSimulationTraits : public BaseTraits
using GFWaterContent = Dune::Dorie::WaterContentAdapter<BaseTraits,FlowParameters,GFMatricHead>;
/// Saturation
using GFSaturation = Dune::Dorie::SaturationAdapter<BaseTraits,FlowParameters,GFMatricHead>;
/// Water Flux reconstruction
using GFWaterFluxReconstruction = Dune::Dorie::RaviartThomasFluxReconstruction<GO0,BaseTraits::GridGeometryType,flux_order>;
// -- Utility Class Definitions -- //
/// Custom VTK output writer
......@@ -216,6 +219,8 @@ private:
using GFWaterContent = typename Traits::GFWaterContent;
/// Saturation
using GFSaturation = typename Traits::GFSaturation;
/// Water Flux reconstruction
using GFWaterFluxReconstruction = typename Traits::GFWaterFluxReconstruction;
// -- Utility Class Definitions -- //
......@@ -262,12 +267,12 @@ protected:
std::shared_ptr<U> u;
mutable std::shared_ptr<GFMatricHead> udgf;
mutable std::shared_ptr<GFWaterFlux> fluxdgf;
mutable std::shared_ptr<GFConductivity> condgf;
mutable std::shared_ptr<GFWaterContent> waterdgf;
mutable std::shared_ptr<GFSaturation> satdgf;
mutable std::shared_ptr<GFMatricHead> udgf;
mutable std::shared_ptr<GFWaterFlux> fluxdgf;
mutable std::shared_ptr<GFConductivity> condgf;
mutable std::shared_ptr<GFWaterContent> waterdgf;
mutable std::shared_ptr<GFSaturation> satdgf;
mutable std::shared_ptr<GFWaterFluxReconstruction> waterfrgf;
std::unique_ptr<Writer> vtkwriter;
std::unique_ptr<AdaptivityHandler> adaptivity;
......@@ -358,8 +363,8 @@ public:
std::shared_ptr<GFWaterFlux> get_waterflux()
{
update_adapters();
assert(fluxdgf);
return fluxdgf;
assert(waterfrgf);
return waterfrgf;
}
protected:
......
......@@ -211,21 +211,21 @@ void TransportSimulation<Traits>::write_data () const
update_adapters();
if (inifile.get<bool>("output.vertexData"))
{
vtkwriter->addVertexData(Cw,"solute");
vtkwriter->addVertexData(Ctotal,"solute (total)");
vtkwriter->template addVertexData<GFSolute>(Cw,"solute");
vtkwriter->template addVertexData<GFTotalSolute>(Ctotal,"solute (total)");
if (inifile.get<bool>("output.debugMode",false))
{
vtkwriter->addVertexData(_gfSaturation,"saturation");
vtkwriter->addVertexData(_gfWaterFlux,"flux");
vtkwriter->template addVertexData<GFSaturation>(_gfSaturation,"saturation");
vtkwriter->template addVertexData<GFWaterFlux>(_gfWaterFlux,"flux");
}
} else {
vtkwriter->addCellData(Cw,"solute");
vtkwriter->addCellData(Ctotal,"solute (total)");
vtkwriter->template addCellData<GFSolute>(Cw,"solute");
vtkwriter->template addCellData<GFTotalSolute>(Ctotal,"solute (total)");
if (inifile.get<bool>("output.debugMode",false))
{
vtkwriter->addCellData(_gfSaturation,"saturation");
vtkwriter->addCellData(_gfWaterFlux,"flux");
vtkwriter->template addCellData<GFSaturation>(_gfSaturation,"saturation");
vtkwriter->template addCellData<GFWaterFlux>(_gfWaterFlux,"flux");
}
}
......
......@@ -34,15 +34,15 @@ namespace Dune{
* @param gfs The GridFunctionsSpace
* @param x The coefficients vector
*/
DiscreteGridFunctionGradient(std::shared_ptr<GFS> gfs, std::shared_ptr<X> x)
DiscreteGridFunctionGradient(const std::shared_ptr<GFS> gfs, std::shared_ptr<X> x)
: Base(*gfs,*x)
, pgfs(gfs)
, px(x)
{}
private:
std::shared_ptr<GFS> pgfs;
std::shared_ptr<X> px;
const std::shared_ptr<GFS> pgfs;
const std::shared_ptr<X> px;
};
}
......
This diff is collapsed.
#ifndef DUNE_DORIE_VOLUME_RAVIART_THOMAS_COEFFICIENTS_HH
#define DUNE_DORIE_VOLUME_RAVIART_THOMAS_COEFFICIENTS_HH
#include <dune/localfunctions/common/localkey.hh>
#include <dune/geometry/type.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/power.hh>
#include <vector>
namespace Dune{
namespace Dorie{
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalCoefficients
{
static_assert(Dune::AlwaysFalse<DF>::value, "Volume Raviart Thomas coefficients is not implemented!");
};
template<class DF, class RF, int k, int dim>
class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::simplex>
{
static constexpr int volsize = (dim == 2) ? dim*k*(k+1)/2
: dim*k*(k+1)*(k+2)/6;
public:
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
for (unsigned int i = 0; i < size(); ++i)
li[i] = LocalKey(0,0,i);
}
constexpr std::size_t size() const
{
return volsize;
}
const LocalKey& localKey (std::size_t i) const
{
return li[i];
}
private:
std::vector<LocalKey> li;
};
template<class DF, class RF, int k, int dim>
class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::cube>
{
static constexpr int volsize = dim*Dune::StaticPower<k+1,dim-1>::power*k;
public:
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
for (unsigned int i = 0; i < size(); ++i)
li[i] = LocalKey(0,0,i);
}
constexpr std::size_t size() const
{
return volsize;
}
const LocalKey& localKey (std::size_t i) const
{
return li[i];
}
private:
std::vector<LocalKey> li;
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_VOLUME_RAVIART_THOMAS_COEFFICIENTS_HH
#ifndef DUNE_DORIE_VOLUME_RAVIART_THOMAS_INTERPOLATION_HH
#define DUNE_DORIE_VOLUME_RAVIART_THOMAS_INTERPOLATION_HH
#include <dune/common/exceptions.hh>
#include <vector>
namespace Dune{
namespace Dorie{
class VolumeRaviartThomasLocalInterpolation
{
public:
VolumeRaviartThomasLocalInterpolation()
{}