Commit 3eec101c authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

Revert "Interpolate solution to water content for grid adaptation"

This reverts commit 53f4b8a9.
Changes were falsely added to this branch.
parent 53f4b8a9
......@@ -20,7 +20,7 @@ public:
/// Do nothing.
/** \return Boolean if operators have to be updated
*/
virtual bool adapt_grid(Grid &grid, GV &gv, GFS &gfs, U &uold, U &unew, const RF time) { return false; }
virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { return false; }
};
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
......@@ -52,22 +52,9 @@ private:
/// Grid operator for Error LOP
using ESTGO = Dune::PDELab::GridOperator<GFS,AGFS,ESTLOP,MBE,RF,RF,RF,NoTrafo,NoTrafo>;
/// Type of the solution GFS
using GFSHelper = GridFunctionSpaceHelper<GV, RF,
Traits::fem_order, Traits::GridGeometryType>;
/// Solution DiscreteGridFunction
using UDGF = Dune::PDELab::DiscreteGridFunction<GFS, U>;
/// GridFunctionAdapter for water content (from matric head)
using WaterDGF = Dune::Dorie::WaterContentAdapter<Traits, Param, UDGF>;
/// GridFunctionAdapter for matric head (from water content)
using HeadDGF = Dune::Dorie::MatricHeadAdapter<Traits, UDGF, GV, Param>;
private:
const Dune::ParameterTree& inifile; //!< Parameter file parser
const Param& param; //!< Parameterization class
const Boundary& boundary; //!< Boundary Condition class
const int maxLevel; //!< Maximum local grid refinement level
const int minLevel; //!< Minimum local grid refinement level
const std::string strategy; //!< Refinement strategy
......@@ -76,50 +63,16 @@ private:
const float adaptivityThreshold; //!< Global error threshold below which no refinement is applied
const int verbose; //!< output verbosity of this object
/// Interpolate a (head) solution into water content space
U solution_to_water_content (const U& solution,
const GFS& gfs, const GV& gv)
{
U ret(gfs, 0.0);
UDGF udgf(gfs, solution);
WaterDGF waterdgf(gv, param, udgf);
Dune::PDELab::interpolate(waterdgf, gfs, ret);
return ret;
}
/// Interpolate a water content solution into matric head space
U water_content_to_solution(const U& water_content,
const GFS& gfs, const GV& gv)
{
U ret(gfs, 0.0);
UDGF waterdgf(gfs, water_content);
HeadDGF headdgf(gv, waterdgf, param);
Dune::PDELab::interpolate(headdgf, gfs, ret);
return ret;
}
void print_vector (const U& u, const std::string name="")
{
std::cout << name << std::endl;
const auto& vec = Dune::PDELab::Backend::native(u);
std::for_each(vec.begin(), vec.end(),
[](const auto val){ std::cout << val << std::endl; }
);
}
public:
/// Initialize members from config file parameters.
/** Disable grid closure for cubic grid.
/// Initialize members from config file parameters.
/** Disable grid closure for cubic grid.
* \param _inifile Parameter file parser
* \param grid Grid to adapt (reference is not saved)
*/
AdaptivityHandler(Dune::ParameterTree &_inifile, Grid &grid,
const Param& _param, const Boundary& _boundary)
: AdaptivityHandlerBase<Traits, GFS, Param, Boundary, U>(),
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>(),
inifile(_inifile),
param(_param),
boundary(_boundary),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
strategy(inifile.get<std::string>("adaptivity.markingStrategy")),
......@@ -151,9 +104,11 @@ public:
Grid& grid,
GV& gv,
GFS& gfs,
const Param& param,
const Boundary& boundary,
const RF time,
U& uold,
U& unew,
const RF time) override
U& unew) override
{
Dune::Timer timer2;
float t2;
......@@ -237,19 +192,7 @@ public:
t_mark = timer3.elapsed();
timer3.reset();
// interpolate solutions to ensure mass conservation
// print_vector(unew, "solution before adaptation");
auto wc_old = solution_to_water_content(uold, gfs, gv);
auto wc_new = solution_to_water_content(unew, gfs, gv);
// print_vector(wc_new, "water content before adaptation");
Dune::PDELab::adapt_grid(grid, gfs, wc_old, wc_new, 2*order);
// interpolate back
// print_vector(wc_new, "water content after adaptation");
uold = water_content_to_solution(wc_old, gfs, gv);
unew = water_content_to_solution(wc_new, gfs, gv);
// print_vector(unew, "solution after adaptation");
Dune::PDELab::adapt_grid(grid, gfs, uold, unew, 2*order);
t_adapt = timer3.elapsed();
timer3.reset();
......@@ -297,8 +240,6 @@ private:
Dune::ParameterTree& inifile;
Grid& grid;
const Param& param;
const Boundary& boundary;
public:
......@@ -307,11 +248,9 @@ public:
* \param _inifile Parameter file parser
* \param _grid Grid to adapt (reference is not saved)
*/
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid, const Param& _param, const Boundary& _boundary) :
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
grid(_grid),
param(_param),
boundary(_boundary)
grid(_grid)
{ }
/// Create a dummy AdaptivityHandler
......@@ -331,8 +270,7 @@ public:
std::is_same<Grid,Dune::UGGrid<3>>::value,
"Adaptivity is only implemented for UGGrid!");
std::unique_ptr<AHB> p;
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U>>
(inifile, grid, param, boundary);
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U>>(inifile,grid);
return p;
}
};
......
......@@ -4,6 +4,55 @@
namespace Dune{
namespace Dorie{
template<typename Traits, typename DGF, typename GV, typename Parameters>
class MatricHeadAdapter : public Dune::PDELab::GridFunctionBase
<Dune::PDELab::GridFunctionTraits<GV,typename Traits::RangeField,1,typename Traits::Scalar>,
MatricHeadAdapter<Traits,DGF,GV,Parameters> >
{
private:
using GFTraits = Dune::PDELab::GridFunctionTraits<GV,typename Traits::RangeField,1,typename Traits::Scalar>;
public:
/*
* \brief Constructor
* \param _gv GridView
* \param _dgf grid function to interpolate from
* \param _param Parameter object
*/
MatricHeadAdapter (const GV& _gv, const DGF& _dgf, const Parameters& _param) :
gv(_gv), dgf(_dgf), param(_param)
{}
/**
* \brief Evaluation of grid function at certain position
* \param e Element entity
* \param x Position in local entity coordinates
* \param y Return value
* \return Value of grid function
*/
void evaluate (const typename GFTraits::ElementType& e,
const typename GFTraits::DomainType& x,
typename GFTraits::RangeType& y) const
{
typename GFTraits::RangeType u;
dgf.evaluate(e,x,u);
// calculate matric head
const typename GFTraits::DomainType xGlobal = e.geometry().global(x);
const typename GFTraits::RangeType head = param.matricHead(u,xGlobal);
y = param.headRefToMiller(head,xGlobal);
}
const GV& getGridView () const
{
return gv;
}
private:
const GV& gv;
const Parameters& param;
const DGF& dgf;
};
template<typename Traits, typename SimTraits>
struct KnoFuInterfaceTraits
{
......
......@@ -48,7 +48,7 @@ Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid>
OutputWriterFactory output_fac(inifile,gv,*gfs,*param,*unew);
output = output_fac.create();
AdaptivityHandlerFactory adaptivity_fac(inifile,*grid, *param, *fboundary);
AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
adaptivity = adaptivity_fac.create();
if(verbose>1 && helper.rank()==0)
......@@ -171,8 +171,7 @@ void Simulation<Traits>::run ()
if(!compute_time_step()){
continue;
}
if (adaptivity->adapt_grid(*grid, gv, *gfs, *uold, *unew, t_start))
{ // reset operators if grid changes
if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary, t_start, *uold, *unew)){ // reset operators if grid changes
operator_setup();
}
output->write_vtk_output(controller->getTime());
......
......@@ -4,177 +4,133 @@
namespace Dune{
namespace Dorie{
template <typename Traits, typename DGF, typename GV, typename Parameters>
class MatricHeadAdapter : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<GV, typename Traits::RangeField, 1, typename Traits::Scalar>,
MatricHeadAdapter<Traits, DGF, GV, Parameters>>
{
private:
using GFTraits = Dune::PDELab::GridFunctionTraits<GV, typename Traits::RangeField, 1, typename Traits::Scalar>;
public:
/*
* \brief Constructor
* \param _gv GridView
* \param _dgf grid function to interpolate from
* \param _param Parameter object
*/
MatricHeadAdapter(const GV &_gv, const DGF &_dgf, const Parameters &_param) : gv(_gv), dgf(_dgf), param(_param)
{
}
/**
* \brief Evaluation of grid function at certain position
* \param e Element entity
* \param x Position in local entity coordinates
* \param y Return value
* \return Value of grid function
*/
void evaluate(const typename GFTraits::ElementType &e,
const typename GFTraits::DomainType &x,
typename GFTraits::RangeType &y) const
{
typename GFTraits::RangeType u;
dgf.evaluate(e, x, u);
// calculate matric head
const typename GFTraits::DomainType xGlobal = e.geometry().global(x);
const typename GFTraits::RangeType head = param.matricHead(u, xGlobal);
y = param.headRefToMiller(head, xGlobal);
}
const GV &getGridView() const
{
return gv;
}
private:
const GV &gv;
const Parameters &param;
const DGF &dgf;
};
/// VTK Plotting: Modified Dune::PDELab::DiscreteGridFunctionGradient class for calculating flux
template <typename T, typename X, typename P>
class GradientFluxAdapter
: public Dune::PDELab::GridFunctionInterface<
Dune::PDELab::GridFunctionTraits<
/// VTK Plotting: Modified Dune::PDELab::DiscreteGridFunctionGradient class for calculating flux
template<typename T, typename X, typename P>
class GradientFluxAdapter
: public Dune::PDELab::GridFunctionInterface<
Dune::PDELab::GridFunctionTraits<
typename T::Traits::GridViewType,
typename T::Traits::FiniteElementType::Traits::LocalBasisType ::Traits::RangeFieldType,
T::Traits::FiniteElementType::Traits::LocalBasisType::Traits ::dimDomain,
typename T::Traits::FiniteElementType::Traits::LocalBasisType
::Traits::RangeFieldType,
T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
::dimDomain,
FieldVector<
typename T::Traits::FiniteElementType::Traits ::LocalBasisType::Traits::RangeFieldType,
T::Traits::FiniteElementType::Traits::LocalBasisType::Traits ::dimDomain>>,
GradientFluxAdapter<T, X, P>>
{
typedef T GFS;
typedef typename GFS::Traits::FiniteElementType::Traits::
LocalBasisType::Traits LBTraits;
public:
typedef Dune::PDELab::GridFunctionTraits<
typename GFS::Traits::GridViewType,
typename LBTraits::RangeFieldType,
LBTraits::dimDomain,
FieldVector<
typename T::Traits::FiniteElementType::Traits
::LocalBasisType::Traits::RangeFieldType,
T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
::dimDomain> >,
GradientFluxAdapter<T,X,P> >
{
typedef T GFS;
typedef typename GFS::Traits::FiniteElementType::Traits::
LocalBasisType::Traits LBTraits;
public:
typedef Dune::PDELab::GridFunctionTraits<
typename GFS::Traits::GridViewType,
typename LBTraits::RangeFieldType,
LBTraits::dimDomain>>
Traits;
private:
typedef Dune::PDELab::GridFunctionInterface<
Traits,
GradientFluxAdapter<T, X, P>>
BaseT;
typedef Dune::PDELab::DiscreteGridFunction<T, X> DGF;
typedef typename GFS::Traits::FiniteElementType::Traits::
LocalBasisType::Traits::DomainFieldType DomainField;
enum
{
dim = GFS::Traits::GridViewType::dimension
};
typedef typename Dune::FieldVector<DomainField, dim> Domain;
public:
/**
LBTraits::dimDomain,
FieldVector<
typename LBTraits::RangeFieldType,
LBTraits::dimDomain> > Traits;
private:
typedef Dune::PDELab::GridFunctionInterface<
Traits,
GradientFluxAdapter<T,X,P> > BaseT;
typedef Dune::PDELab::DiscreteGridFunction<T,X> DGF;
typedef typename GFS::Traits::FiniteElementType::Traits::
LocalBasisType::Traits::DomainFieldType DomainField;
enum {dim = GFS::Traits::GridViewType::dimension};
typedef typename Dune::FieldVector<DomainField,dim> Domain;
public:
/**
* \param gfs The GridFunctionsSpace
* \param x_ The coefficients vector
* \param p_ The parameter class
* \param dgf_ The solution (head) GridFunction
*/
GradientFluxAdapter(const GFS &gfs, const X &x_, P &p_, const DGF &dgf_)
: pgfs(stackobject_to_shared_ptr(gfs)), lfs(gfs), lfs_cache(lfs), x_view(x_), xl(lfs.size()), p(p_), dgf(dgf_)
{
}
// Evaluate
void evaluate(const typename Traits::ElementType &e,
const typename Traits::DomainType &x,
typename Traits::RangeType &y) const
{
// get and bind local functions space
lfs.bind(e);
lfs_cache.update();
x_view.bind(lfs_cache);
// get local coefficients
xl.resize(lfs.size());
x_view.read(xl);
x_view.unbind();
// get Jacobian of geometry
const typename Traits::ElementType::Geometry::JacobianInverseTransposed
JgeoIT = e.geometry().jacobianInverseTransposed(x);
// get local Jacobians/gradients of the shape functions
std::vector<typename LBTraits::JacobianType> J(lfs.size());
lfs.finiteElement().localBasis().evaluateJacobian(x, J);
typename Traits::RangeType gradphi;
y = 0;
for (unsigned int i = 0; i < lfs.size(); ++i)
{
// compute global gradient of shape function i
gradphi = 0;
JgeoIT.umv(J[i][0], gradphi);
// sum up global gradients, weighting them with the appropriate coeff
y.axpy(xl[i], gradphi);
}
y -= p.gravity();
// get local conductivity
typename DGF::Traits::RangeType u;
dgf.evaluate(e, x, u);
const Domain xGlobal = e.geometry().global(x);
const typename DGF::Traits::RangeType head = p.headMillerToRef(u, xGlobal);
const typename DGF::Traits::RangeType sat = p.saturation(head, xGlobal);
const typename DGF::Traits::RangeType relCond = p.condFactor(sat, xGlobal);
const typename DGF::Traits::RangeType satCond = p.condRefToMiller(p.K(xGlobal), xGlobal);
const typename DGF::Traits::RangeType cond = relCond * satCond;
y *= -cond;
}
//! get a reference to the GridView
const typename Traits::GridViewType &getGridView() const
{
return pgfs->gridView();
}
private:
typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
typedef typename X::template ConstLocalView<LFSCache> XView;
shared_ptr<GFS const> pgfs;
mutable LFS lfs;
mutable LFSCache lfs_cache;
mutable XView x_view;
mutable std::vector<typename Traits::RangeFieldType> xl;
P &p;
const DGF &dgf;
GradientFluxAdapter (const GFS& gfs, const X& x_, P& p_, const DGF& dgf_)
: pgfs(stackobject_to_shared_ptr(gfs))
, lfs(gfs)
, lfs_cache(lfs)
, x_view(x_)
, xl(lfs.size())
, p(p_)
, dgf(dgf_)
{ }
// Evaluate
void evaluate (const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
// get and bind local functions space
lfs.bind(e);
lfs_cache.update();
x_view.bind(lfs_cache);
// get local coefficients
xl.resize(lfs.size());
x_view.read(xl);
x_view.unbind();
// get Jacobian of geometry
const typename Traits::ElementType::Geometry::JacobianInverseTransposed
JgeoIT = e.geometry().jacobianInverseTransposed(x);
// get local Jacobians/gradients of the shape functions
std::vector<typename LBTraits::JacobianType> J(lfs.size());
lfs.finiteElement().localBasis().evaluateJacobian(x,J);
typename Traits::RangeType gradphi;
y = 0;
for(unsigned int i = 0; i < lfs.size(); ++i)
{
// compute global gradient of shape function i
gradphi = 0;
JgeoIT.umv(J[i][0], gradphi);
// sum up global gradients, weighting them with the appropriate coeff
y.axpy(xl[i], gradphi);
}
y -= p.gravity();
// get local conductivity
typename DGF::Traits::RangeType u;
dgf.evaluate(e,x,u);
const Domain xGlobal = e.geometry().global(x);
const typename DGF::Traits::RangeType head = p.headMillerToRef(u, xGlobal);
const typename DGF::Traits::RangeType sat = p.saturation(head, xGlobal);
const typename DGF::Traits::RangeType relCond = p.condFactor(sat, xGlobal);
const typename DGF::Traits::RangeType satCond = p.condRefToMiller(p.K(xGlobal), xGlobal);
const typename DGF::Traits::RangeType cond = relCond * satCond;
y *= -cond;
}
//! get a reference to the GridView
const typename Traits::GridViewType& getGridView () const
{
return pgfs->gridView();
}
private:
typedef Dune::PDELab::LocalFunctionSpace<GFS> LFS;
typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
typedef typename X::template ConstLocalView<LFSCache> XView;
shared_ptr<GFS const> pgfs;
mutable LFS lfs;
mutable LFSCache lfs_cache;
mutable XView x_view;
mutable std::vector<typename Traits::RangeFieldType> xl;
P& p;
const DGF& dgf;
};
}
......@@ -280,7 +236,7 @@ namespace Dune{
* \param gf Grid function to probe
* \see RichardsEquationParameter
*/
WaterContentAdapter (const typename Traits::GridViewType& g_, const P& p_, const GF& gf_)
WaterContentAdapter (const typename Traits::GridViewType& g_, P& p_, const GF& gf_)
: g(g_), p(p_), gf(gf_)
{}
......@@ -316,7 +272,7 @@ namespace Dune{
private:
const typename Traits::GridViewType& g;
const P& p;
P& p;
const GF& gf;
};
}
......
......@@ -274,8 +274,8 @@ public:
continue;
}
if (this->adaptivity->adapt_grid(
*this->grid, this->gv, *this->gfs,
*this->uold, *this->unew, time))
*this->grid, this->gv, *this->gfs, *this->param, *this->fboundary,
time, *this->uold, *this->unew))
{
// reset operators if grid changes
this->operator_setup();
......@@ -343,10 +343,10 @@ int main (int argc, char** argv)
// run test
// inifile["adaptivity.useAdaptivity"] = "true";
// using TestSim = Dune::Dorie::TestSimulation<Cube<2,1>>;
// auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile, helper);
using TestSim = Dune::Dorie::TestSimulation<CubeAdaptive<2,1>>;
auto grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<2>>(inifile, helper);
using TestSim = Dune::Dorie::TestSimulation<Cube<2,1>>;
auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile, helper);
// using TestSim = Dune::Dorie::TestSimulation<CubeAdaptive<2,1>>;
// auto grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<2>>(inifile, helper);
TestSim sim(helper, grid, inifile);
sim.run_test();
......
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