[Richards&Transport] Ensure constness of adapters & Add total solute GF

* This commit also fixes a mistake in creation of flux reconstruction, which could in certain cases share a cache set up with other times.
* Total solute adapter is essentially equal to product of grid functions defined in PDELab, however, this class get ownership of grid functions. This is necessary when sharing adapters outside of the main simulation.
parent b5b80e2c
......@@ -48,9 +48,9 @@ namespace Dune{
* @param p Parametrization class
* @see RichardsEquationParameter
*/
SaturationAdapter(const std::shared_ptr<GF> gf_h,
SaturationAdapter(std::shared_ptr<const GF> gf_h,
const typename Traits::GridViewType& gv,
const std::shared_ptr<P> p)
std::shared_ptr<const P> p)
: _gf_h(gf_h)
, _gv(gv)
, _p(p)
......@@ -88,9 +88,9 @@ namespace Dune{
}
private:
const std::shared_ptr<GF> _gf_h;
const std::shared_ptr<const GF> _gf_h;
const typename Traits::GridViewType _gv;
const std::shared_ptr<P> _p;
const std::shared_ptr<const P> _p;
};
}
......
......@@ -48,9 +48,9 @@ namespace Dune{
* @param p Parametrization class
* @see RichardsEquationParameter
*/
WaterContentAdapter(const std::shared_ptr<GF> gf_h,
WaterContentAdapter(std::shared_ptr<const GF> gf_h,
const typename Traits::GridViewType& gv,
const std::shared_ptr<P> p)
std::shared_ptr<const P> p)
: _gf_h(gf_h)
, _gv(gv)
, _p(p)
......@@ -89,9 +89,9 @@ namespace Dune{
}
private:
const std::shared_ptr<GF> _gf_h;
const std::shared_ptr<const GF> _gf_h;
const typename Traits::GridViewType _gv;
const std::shared_ptr<P> _p;
const std::shared_ptr<const P> _p;
};
}
}
......
......@@ -3,7 +3,7 @@
#include <cmath>
#include <vector>
#include <map>
#include <unordered_map>
#include <utility>
#include <functional>
#include <memory>
......@@ -71,12 +71,12 @@ private:
using ScaleAdapter = Dorie::ScalingAdapter<Traits>;
/// Storage for parameters and skaling factors
using ParameterStorage = std::vector<
std::tuple<std::shared_ptr<Parameterization>, Scaling>
using ParameterStorage = std::unordered_map<Index,
const std::tuple<const std::shared_ptr<const Parameterization>, const Scaling>
>;
/// Need this gruesome typedef because 'map::value_type' has const key
using Cache = typename ParameterStorage::value_type;
using Cache = std::tuple<std::shared_ptr<const Parameterization>, Scaling>;
/// Configuration file tree
const Dune::ParameterTree& _config;
/// Grid view of the coarsest grid configuration (level 0)
......@@ -97,7 +97,7 @@ private:
/// Check if an entity has been cached
void verify_cache () const
{
if (not std::get<std::shared_ptr<Parameterization>>(_cache)) {
if (not std::get<std::shared_ptr<const Parameterization>>(_cache)) {
_log->error("Parameterization cache is empty. Call 'bind' before "
"accessing the cache or the parameterization");
DUNE_THROW(Dune::InvalidStateException,
......@@ -153,7 +153,7 @@ public:
entity = entity.father();
}
const auto index = _mapper.index(entity);
_cache = _param[index];
_cache = _param.find(index)->second;
}
/// Return a scaled version of the Conductivity function
......@@ -165,7 +165,7 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<const Parameterization>>(_cache);
const auto cond_f = par->conductivity_f();
using Saturation = typename Parameterization::Saturation;
......@@ -185,7 +185,7 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<const Parameterization>>(_cache);
const auto sat_f = par->saturation_f();
using MatricHead = typename Parameterization::MatricHead;
......@@ -205,7 +205,7 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<const Parameterization>>(_cache);
// get water content function and apply the scaling
const auto& scale_por = std::get<Scaling>(_cache).scale_por;
......@@ -229,7 +229,6 @@ private:
// create the parameterization data
const auto parameterization_map = read_parameter_file(param_file);
_param.resize(element_index_map.size());
// build global scaling
using SAF = Dorie::ScalingAdapterFactory<Traits>;
......@@ -246,9 +245,9 @@ private:
const auto pos = elem.geometry().center();
// read values and insert
auto p = parameterization_map.at(element_index_map.at(index));
const auto p = parameterization_map.at(element_index_map.at(index));
const auto scale = scaling_adapter->evaluate(pos);
_param.at(index) = std::make_tuple(p, scale);
_param.emplace(index,std::make_tuple(p, scale));
}
// check that mapper can map to parameterization
......@@ -258,12 +257,12 @@ private:
/// Read the YAML parameter file and insert information into a map.
/** \param param_file The top YAML node to read from (file root)
*/
std::map<int, std::shared_ptr<Parameterization>> read_parameter_file (
std::map<int, std::shared_ptr<const Parameterization>> read_parameter_file (
const YAML::Node& param_file
)
{
// mapping: index on grid to parameterization object (will be returned)
std::map<int, std::shared_ptr<Parameterization>> ret;
std::map<int, std::shared_ptr<const Parameterization>> ret;
// set to check for duplicate indices
std::set<int> param_index_set;
......
......@@ -163,9 +163,9 @@ protected:
typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
const Parameter& param; //!< class providing equation parameters
const Boundary& boundary; //!< class providing boundary conditions
const SourceTerm& sourceTerm; //!< class providing source term information
std::shared_ptr<const Parameter> param; //!< class providing equation parameters
std::shared_ptr<const Boundary> boundary; //!< class providing boundary conditions
std::shared_ptr<const SourceTerm> sourceTerm; //!< class providing source term information
const RichardsDGMethod::Type method; //!< Interior penalty type
const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
const RichardsDGWeights::Type weights; //!< Gradient weighting
......@@ -205,9 +205,9 @@ public:
*/
RichardsDGSpatialOperator (
const Dune::ParameterTree& config,
const Parameter& param_,
const Boundary& boundary_,
const SourceTerm& sourceTerm_,
std::shared_ptr<const Parameter> param_,
std::shared_ptr<const Boundary> boundary_,
std::shared_ptr<const SourceTerm> sourceTerm_,
const RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
const RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
const RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
......@@ -273,9 +273,9 @@ public:
auto gt = eg.geometry();
// bind parameterization and retrieve functions
param.bind(eg.entity());
const auto saturation_f = param.saturation_f();
const auto conductivity_f = param.conductivity_f();
param->bind(eg.entity());
const auto saturation_f = param->saturation_f();
const auto conductivity_f = param->conductivity_f();
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
......@@ -330,7 +330,7 @@ public:
for (unsigned int i = 0; i<lfsu.size(); i++)
gradu.axpy(x(lfsu,i),gradphiu[i]);
gradu -= param.gravity();
gradu -= param->gravity();
// compute conductivity
const RF cond = conductivity_f(saturation_f(u));
......@@ -395,16 +395,16 @@ public:
auto gtface = ig.geometry();
// bind parameterization and retrieve functions
param.bind(ig.inside());
const auto saturation_f_s = param.saturation_f();
const auto conductivity_f_s = param.conductivity_f();
param->bind(ig.inside());
const auto saturation_f_s = param->saturation_f();
const auto conductivity_f_s = param->conductivity_f();
// cache saturated conductivity
const RF satCond_s = conductivity_f_s(1.0);
param.bind(ig.outside());
const auto saturation_f_n = param.saturation_f();
const auto conductivity_f_n = param.conductivity_f();
param->bind(ig.outside());
const auto saturation_f_n = param->saturation_f();
const auto conductivity_f_n = param->conductivity_f();
const RF satCond_n = conductivity_f_n(1.0);
......@@ -515,8 +515,8 @@ public:
gradu_n.axpy(x_n(lfsu_n,i),gradphiu_n[i]);
// update with gravity vector
gradu_s -= param.gravity();
gradu_n -= param.gravity();
gradu_s -= param->gravity();
gradu_n -= param->gravity();
// compute jump in solution
const RF jump = u_s - u_n;
......@@ -685,9 +685,9 @@ public:
auto gtface = ig.geometry();
// bind parameterization and retrieve functions
param.bind(ig.inside());
const auto saturation_f_s = param.saturation_f();
const auto conductivity_f_s = param.conductivity_f();
param->bind(ig.inside());
const auto saturation_f_s = param->saturation_f();
const auto conductivity_f_s = param->conductivity_f();
// cache saturated conductivity
const RF satCond_s = conductivity_f_s(1.0);
......@@ -747,14 +747,14 @@ public:
// query the type of boundary condition
// set the internal bcType flag
const auto bc = boundary.bc(ig.intersection(),it.position(),time);
const auto bc = boundary->bc(ig.intersection(),it.position(),time);
auto bcType = BCType::none;
RF normal_flux;
if (BoundaryCondition::isNeumann(bc))
{
// evaluate flux boundary condition
normal_flux = boundary.j(ig.intersection(),it.position(),time);
normal_flux = boundary->j(ig.intersection(),it.position(),time);
bcType = BCType::Neumann;
}
......@@ -772,7 +772,7 @@ public:
if (bcType == BCType::Neumann)
{
// flux is given parallel to gravity vector
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param->gravity() );
// update residual (flux term)
if constexpr (lopcase != LOPCase::RTVolume)
......@@ -790,13 +790,13 @@ public:
gradu_s.axpy(x_s(lfsu_s,i),gradphiu_s[i]);
// update with gravity vector
gradu_s -= param.gravity();
gradu_s -= param->gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = boundary.g(ig.intersection(),it.position(),time);
const RF g = boundary->g(ig.intersection(),it.position(),time);
const RF jump = u_s - g;
// conductivity factors
......@@ -857,9 +857,9 @@ public:
/* if (BoundaryCondition::isLimitedFlux(bc))
{
const RF satCond_s = param.K(*(ig.inside()),local_s);
normal_flux = boundary.j(ig.intersection(),it.position(),time);
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );
const RF satCond_s = param->K(*(ig.inside()),local_s);
normal_flux = boundary->j(ig.intersection(),it.position(),time);
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param->gravity() );
if (normal_flux < - satCond_s)
{
......@@ -910,7 +910,7 @@ public:
// scaled weight factor
const RF factor = it.weight() * gt.integrationElement(p_local);
const RF source = sourceTerm.q(eg.entity(),p_local,time) / gt.volume();
const RF source = sourceTerm->q(eg.entity(),p_local,time) / gt.volume();
// update residual
for (unsigned int i = 0; i<lfsv.size(); i++)
......@@ -966,7 +966,7 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
private:
const Parameter& param; //!< class providing equation parameters
std::shared_ptr<const Parameter> param; //!< class providing equation parameters
RF time;
const int intorderadd;
const int quadrature_factor;
......@@ -978,7 +978,7 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
*/
RichardsDGTemporalOperator (
const Dune::ParameterTree& config,
const Parameter& param_,
std::shared_ptr<const Parameter> param_,
int intorderadd_ = 0,
int quadrature_factor_ = 2
) :
......@@ -1002,9 +1002,9 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
auto gt = eg.geometry();
// bind parameterization and retrieve functions
param.bind(eg.entity());
const auto saturation_f = param.saturation_f();
const auto water_content_f = param.water_content_f();
param->bind(eg.entity());
const auto saturation_f = param->saturation_f();
const auto water_content_f = param->water_content_f();
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
......
......@@ -48,11 +48,11 @@ RichardsSimulation<Traits>::RichardsSimulation (
// --- Create the new parameter class
auto element_map = _grid_creator.element_index_map();
fparam = std::make_shared<FlowParameters>(inifile, grid, element_map);
fparam = std::make_shared<const FlowParameters>(inifile, grid, element_map);
// --- Operator Helper Classes ---
fboundary = std::make_unique<FlowBoundary>(inifile);
fsource = std::make_unique<FlowSource>(inifile);
fboundary = std::make_shared<const FlowBoundary>(inifile);
fsource = std::make_shared<const FlowSource>(inifile);
finitial = FlowInitialFactory::create(inifile, gv, this->_log);
// --- Local Operators ---
......@@ -64,13 +64,13 @@ RichardsSimulation<Traits>::RichardsSimulation (
const auto upwinding = std::get<OP::RichardsDGUpwinding::Type>(settings);
const auto weights = std::get<OP::RichardsDGWeights::Type>(settings);
slop = std::make_shared<SLOP>(inifile, *fparam, *fboundary, *fsource,
method, upwinding, weights);
slop = std::make_shared<SLOP>(inifile, fparam, fboundary, fsource,
method, upwinding, weights);
#else
slop = std::make_shared<SLOP>(inifile, *fparam, *fboundary, *fsource);
slop = std::make_shared<SLOP>(inifile, fparam, fboundary, fsource);
#endif // EXPERIMENTAL_DG_FEATURES
tlop = std::make_shared<TLOP>(inifile, *fparam);
tlop = std::make_shared<TLOP>(inifile, fparam);
controller = std::make_unique<CalculationController>(
inifile, *fboundary, this->_log);
......
......@@ -293,9 +293,9 @@ protected:
std::shared_ptr<LSGFS> lsgfs;
std::shared_ptr<LSCC> lscc;
std::shared_ptr<FlowParameters> fparam;
std::unique_ptr<FlowBoundary> fboundary;
std::unique_ptr<FlowSource> fsource;
std::shared_ptr<const FlowParameters> fparam;
std::shared_ptr<const FlowBoundary> fboundary;
std::shared_ptr<const FlowSource> fsource;
std::unique_ptr<FlowInitial> finitial;
std::unique_ptr<CalculationController> controller;
......@@ -424,9 +424,9 @@ public:
*
* @return Pointer to a matric head grid function
*/
std::shared_ptr<GFMatricHead> get_matric_head(ConstState state) const
std::shared_ptr<const GFMatricHead> get_matric_head(ConstState state) const
{
return std::make_shared<GFMatricHead>(state.grid_function_space,state.coefficients);
return std::make_shared<const GFMatricHead>(state.grid_function_space,state.coefficients);
}
/*------------------------------------------------------------------------*//**
......@@ -434,7 +434,7 @@ public:
*
* @return Pointer to a matric head grid function
*/
std::shared_ptr<GFMatricHead> get_matric_head() const
std::shared_ptr<const GFMatricHead> get_matric_head() const
{
return get_matric_head(current_state());
}
......@@ -444,9 +444,9 @@ public:
*
* @return Pointer to a conductivity grid function
*/
std::shared_ptr<GFConductivity> get_conductivity(ConstState state) const
std::shared_ptr<const GFConductivity> get_conductivity(ConstState state) const
{
return std::make_shared<GFConductivity>(gv, fparam);
return std::make_shared<const GFConductivity>(gv, fparam);
}
/*------------------------------------------------------------------------*//**
......@@ -454,7 +454,7 @@ public:
*
* @return Pointer to a conductivity grid function
*/
std::shared_ptr<GFConductivity> get_conductivity() const
std::shared_ptr<const GFConductivity> get_conductivity() const
{
return get_conductivity(current_state());
}
......@@ -464,9 +464,9 @@ public:
*
* @return Pointer to a water content grid function
*/
std::shared_ptr<GFWaterContent> get_water_content(ConstState state) const
std::shared_ptr<const GFWaterContent> get_water_content(ConstState state) const
{
return std::make_shared<GFWaterContent>(get_matric_head(state), gv, fparam);
return std::make_shared<const GFWaterContent>(get_matric_head(state), gv, fparam);
}
/*------------------------------------------------------------------------*//**
......@@ -474,7 +474,7 @@ public:
*
* @return Pointer to a water content grid function
*/
std::shared_ptr<GFWaterContent> get_water_content() const
std::shared_ptr<const GFWaterContent> get_water_content() const
{
return get_water_content(current_state());
}
......@@ -484,9 +484,9 @@ public:
*
* @return Pointer to a saturation grid function
*/
std::shared_ptr<GFSaturation> get_saturation(ConstState state) const
std::shared_ptr<const GFSaturation> get_saturation(ConstState state) const
{
return std::make_shared<GFSaturation>(get_matric_head(state), gv, fparam);
return std::make_shared<const GFSaturation>(get_matric_head(state), gv, fparam);
}
/*------------------------------------------------------------------------*//**
......@@ -494,7 +494,7 @@ public:
*
* @return Pointer to a saturation grid function
*/
std::shared_ptr<GFSaturation> get_saturation() const
std::shared_ptr<const GFSaturation> get_saturation() const
{
return get_saturation(current_state());
}
......@@ -504,9 +504,9 @@ public:
*
* @return Pointer to a water flux grid function
*/
std::shared_ptr<GFWaterFlux> get_water_flux(ConstState state) const
std::shared_ptr<const GFWaterFlux> get_water_flux(ConstState state) const
{
return std::make_shared<GFWaterFlux>(state.grid_function_space, state.coefficients, fparam);
return std::make_shared<const GFWaterFlux>(state.grid_function_space, state.coefficients, fparam);
}
/*------------------------------------------------------------------------*//**
......@@ -514,7 +514,7 @@ public:
*
* @return Pointer to a water flux grid function
*/
std::shared_ptr<GFWaterFlux> get_water_flux() const
std::shared_ptr<const GFWaterFlux> get_water_flux() const
{
return get_water_flux(current_state());
}
......@@ -527,16 +527,23 @@ public:
*
* @return Pointer to a (reconstructed) water flux grid function
*/
std::shared_ptr<GFFluxReconstruction> get_water_flux_reconstructed(ConstState state) const
std::shared_ptr<const GFFluxReconstruction> get_water_flux_reconstructed(ConstState state) const
{
std::shared_ptr<GFFluxReconstruction> gf_ptr;
auto& cache = cache_water_flux_gf_rt;
if constexpr (enable_rt_engine)
{
// use cache if possible
if (not cache_water_flux_gf_rt) {
// create a flux reconstruction grid function
cache_water_flux_gf_rt
= std::make_shared<GFFluxReconstruction>(this->_log,
gv,inifile.sub("fluxReconstruction"));
if (state.grid_function_space != gfs or
state.coefficients != u or
state.time != current_time())
{
// if state is different to current state, create flux from zero
gf_ptr = std::make_unique<GFFluxReconstruction>(
this->_log,gv,inifile.sub("fluxReconstruction"));
assert(fboundary);
assert(slop);
......@@ -544,18 +551,35 @@ public:
slop->setTime(state.time);
// update it with the state
cache_water_flux_gf_rt->update(*(state.coefficients),
*(state.grid_function_space),
*slop);
gf_ptr->update(*(state.coefficients),
*(state.grid_function_space),
*slop);
slop->setTime(current_time());
} else if (not cache) {
// if state is equal to current state, use cache.
cache = std::make_unique<GFFluxReconstruction>(
this->_log,gv,inifile.sub("fluxReconstruction"));
assert(fboundary);
assert(slop);
// update it with current state
cache->update(*u,*gfs,*slop);
gf_ptr = cache;
} else {
gf_ptr = cache;
}
} else {
DUNE_THROW(NotImplemented,
"Flux reconstruction is not implemented for the selected "
<< "configuration");
}
return cache_water_flux_gf_rt;
return gf_ptr;
}
/*------------------------------------------------------------------------*//**
......@@ -563,7 +587,7 @@ public:
*
* @return Pointer to a (reconstructed) water flux grid function
*/
std::shared_ptr<GFFluxReconstruction> get_water_flux_reconstructed() const
std::shared_ptr<const GFFluxReconstruction> get_water_flux_reconstructed() const
{
return get_water_flux_reconstructed(current_state());
}
......
#ifndef DUNE_DORIE_TOTAL_SOLUTE_ADAPTER_HH
#define DUNE_DORIE_TOTAL_SOLUTE_ADAPTER_HH
#include <dune/pdelab/common/function.hh>
namespace Dune{
namespace Dorie{
/*---------------------------------------------------------------------*//**
* @brief Adapter to obtain total solute
*
* @tparam GFSolute Solute concentration grid function
* @tparam GFWaterContent Water content grid function
*/
template<class GFSolute, class GFWaterContent>
class TotalSoluteAdapter
: public Dune::PDELab::GridFunctionBase<
typename GFSolute::Traits,
TotalSoluteAdapter<GFSolute,GFWaterContent>
>
{
static_assert(std::is_same_v<typename GFSolute::Traits,
typename GFWaterContent::Traits>,
"GFSolute and GFWaterContent must have the same traits!");
public:
using Traits = typename GFSolute::Traits;
public:
/*-------------------------------------------------------------------*//**
* @brief Constructor for the TotalSoluteAdapter
*
* @param[in] solute_gf The solute grid function
* @param[in] water_content_gf The water content grid function
*/
TotalSoluteAdapter(std::shared_ptr<const GFSolute> solute_gf,
std::shared_ptr<const GFWaterContent> water_content_gf)
: _solute_gf(solute_gf)
, _water_content_gf(water_content_gf)
{}
/*-------------------------------------------------------------------*//**
* @brief Evaluation of the total solute for a given entity e in an
* entity position x
*
* @param[in] e Entity of a grid
* @param[in] x Position in local coordinates with respect the entity
* @param y Total solute at position x
*/
void evaluate ( const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
typename Traits::RangeType solute;
typename Traits::RangeType water_content;
_solute_gf->evaluate(e, x, solute);
_water_content_gf->evaluate(e, x, water_content);
y = solute*water_content;
}
/*-------------------------------------------------------------------*//**
* @brief Function that returns a grid view valid for this grid
* function
*
* @return Grid view
*/
const typename Traits::GridViewType& getGridView() const
{
return _solute_gf->getGridView();
}
private:
const std::shared_ptr<const GFSolute> _solute_gf;
const std::shared_ptr<const GFWaterContent> _water_content_gf;
};
}
}
#endif
......@@ -109,7 +109,7 @@ public:
* @param[in] diff_coeff The diffusion coefficient
*/
TransportFVSpatialOperator(
Boundary& boundary,
std::shared_ptr<const Boundary> boundary,
std::shared_ptr<GridFunctionWaterFlux> gf_water_flux,
std::shared_ptr<GridFunctionWaterContent> gf_water_content,
double diff_coeff
......@@ -290,7 +290,7 @@ public:
const auto& ref_el_f = referenceElement(geo_f);
const auto& center_position_f = ref_el_f.position(0,0);
auto bc = _boundary.bc(entity_f,center_position_f,_time);
auto bc = _boundary->bc(entity_f,center_position_f,_time);
if (BoundaryCondition::isDirichlet(bc))
{
......@@ -305,7 +305,7 @@ public:
auto distance = center_position_i_g.two_norm();
// Evaluate Dirichlet condition
auto g = _boundary.g(entity_f,center_position_f,_time);
auto g = _boundary->g(entity_f,center_position_f,_time);
// Face volume for integration
auto volume_f = geo_f.volume();
......@@ -358,7 +358,7 @@ public:
auto diff_coeff_i = _diff_coeff;
// Contribution to residual from Neumann boundary
auto j = _boundary.j(entity_f,center_position_f,_time);
auto j = _boundary->j(entity_f,center_position_f,_time);
// Solute flux in normal direction w.r.t the intersection
auto soulte_flux_n = -diff_coeff_i*j*water_content_i;
......@@ -415,7 +415,7 @@ public:
}
private:
Boundary& _boundary;
const std::shared_ptr<const Boundary> _boundary;
const std::shared_ptr<GridFunctionWaterFlux> _gf_water_flux;
const std::shared_ptr<GridFunctionWaterContent> _gf_water_content;
double _time;
......
......@@ -45,7 +45,7 @@ TransportSimulation<Traits>::TransportSimulation(
Dune::PDELab::constraints(*gfs,*cc,false);
// --- Operator Helper Classes ---
sboundary = std::make_unique<SoluteBoundary>(inifile);
sboundary = std::make_shared<SoluteBoundary>(inifile);
sinitial = std::make_unique<SoluteInitial>(gv,0.);
controller = std::make_unique<CalculationController>(inifile,
......@@ -113,7 +113,7 @@ void TransportSimulation<Traits>::operator_setup()
// --- Local Operators ---
double diff_coeff = inifile.get<double>("parameters.molecularDiffusion");
slop = std::make_unique<SLOP>(
*sboundary,water_flux_gf,water_content_gf,diff_coeff);
sboundary,water_flux_gf,water_content_gf,diff_coeff);
tlop = std::make_unique<TLOP>(water_content_gf);
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
......@@ -217,7 +217,7 @@ void TransportSimulation<Traits>::step()
// invalidate cache for solute flux reconstruction
cache_solute_flux_gf_rt.reset();
time_steps++;
if (this->output_policy() == OutputPolicy::EndOfTransportStep)
......