Commit 15f75df2 authored by Santiago's avatar Santiago

Merge branch '101-couple-the-transportsimulation-with-richardssimulation' into...

Merge branch '101-couple-the-transportsimulation-with-richardssimulation' into 111-add-dg-local-operator-to-the-simulationtransport-class
parents 98380de2 8b220558
......@@ -46,6 +46,7 @@ dune_cmake_sphinx_doc(SPHINX_CONF ${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in
man-config-file.rst
man-parameter-file.rst
man-grid.rst
man-fluxes.rst
MODULE_ONLY)
if(TARGET sphinx_html)
......
......@@ -222,27 +222,39 @@ adding an empty line, make text **bold** or ``monospaced``.
<suggestion> 1 </suggestion>
<values> 1, 2, 3 </values>
</parameter>
</category>
<parameter name="fluxReconstruction">
<category name="fluxReconstruction">
<parameter name="enable">
<definition> Apply the flux reconstruction method to the solved matric
head and obtain conservative gradients. It always computes (internally)
the local lifting independent whether the ``lifting`` keyword is active
or not.
head and obtain conservative gradients. It always computes
(internally) the local lifting independent whether the ``lifting``
keyword is active or not.
</definition>
<suggestion> true </suggestion>
<values> true, false </values>
</parameter>
<parameter name="lifting">
<definition> Compute the local lifting for the discrete gradient. It serves
as a measure of the influence of symmetry term on the interior degrees
of freedom.
<parameter name="checkJumps">
<definition> Check that flux reconstruction engine is creating conforming
normal fluxes up to ``jumpTol``. ProTip: Setting warnings together
with a very low tolerance will let you track the changes on the
quality of the flux reconstruction.
</definition>
<suggestion> false </suggestion>
<values> true, false </values>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
<category name="NewtonParameters">
<parameter name="ForceIteration" hidden="true">
<definition> Force the Newton solver to calculate at least one iteration,
......
......@@ -266,4 +266,32 @@ adding an empty line, make text **bold** or ``monospaced``.
<suggestion> true </suggestion>
</parameter>
</category>
<category name="fluxReconstruction">
<parameter name="enable">
<definition> Apply the flux reconstruction method to the solved solute
concentration and obtain conservative fluxes.
</definition>
<suggestion> true </suggestion>
<values> true, false </values>
</parameter>
<parameter name="checkJumps">
<definition> Check that flux reconstruction engine is creating conforming
normal fluxes up to ``jumpTol``. ProTip: Setting warnings together
with a very low tolerance will let you track the changes on the
quality of the flux reconstruction.
</definition>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
</dorie>
......@@ -49,6 +49,7 @@ Manual
python-dorie-wrapper
man-config-file
man-grid
man-fluxes
man-bcfile
man-cookbook
man-parameter-file
......
......@@ -97,44 +97,6 @@ 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 (dg) 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
-------------------
The flux reconstruction is a projection of the fluxes of the dG method into a vector field function. Using correct elements, this procedure can ensure that fluxes in normal direction to the element are *equivalent* to those computed by the Discontinuous Galerkin method, and most importantly, it can also ensure the continuity of them. Hence, the resulting vector field is useful to compute other problems that rely on the fluxes of the water (i.e. solute transport).
The flux reconstruction technique used in DORiE always use Raviar Thomas finite elements of one degree less than the one set for the Richards model. It can be identified in the vtk file by the name ``flux_RT<k-1>``, where ``k`` is the finite element order set for the Richards model. Flux reconstruction is not available for non-conforming grids (i.e. Cube-adaptive grids).
+---------------------------+---+---+---+
| Richards FEorder | 1 | 2 | 3 |
+============+====+=========+===+===+===+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | ✓ |
| Non-adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | |
+------------+----+---------+---+---+---+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
| Adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
+------------+----+---------+---+---+---+
Examples
========
......
Fluxes
======
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.
.. DORiE allow finite volume computations under certain specific conditions. In such case, if generated, the raw flux output generated by DORiE has no meaning. The reason is that finite volumes are computed with finite elements of order 0 where gradients are 0.
Flux reconstruction
-------------------
The flux reconstruction is a projection of the fluxes used in the Discontinuous Galerkin method into a vector field function. Using correct elements, this procedure can ensure that fluxes in normal direction to the element are *equivalent* to those computed by the Discontinuous Galerkin method, and most importantly, it can also ensure the continuity of them. Hence, the resulting vector field is useful to compute other problems that rely on the fluxes of the water (i.e. solute transport).
The flux reconstruction technique always use Raviar Thomas finite elements of one degree less than the one set for the Richards model. It can be identified in the vtk file by the name ``flux_RT{k-1}``, where ``k`` is the finite element order set for the Richards model. Flux reconstruction is not available for non-conforming grids (i.e. Cube-adaptive grids).
+---------------------------+---+---+---+
| Richards FEorder | 1 | 2 | 3 |
+============+====+=========+===+===+===+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | ✓ |
| Non-adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | |
+------------+----+---------+---+---+---+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
| Adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
+------------+----+---------+---+---+---+
Usage
-----
To activate/deactivate flux reconstruction use the keyword ``richards.numeric.fluxReconstruction = true/false`` in the :doc:`config file<man-config-file>`.
......@@ -121,7 +121,8 @@ public:
const LocalAssembler & local_assembler_,
const GFSVVolume & gfsv_volume_,
const GFSVSkeleton & gfsv_skeleton_,
int int_order_add_ = 2
int int_order_add_ = 2,
int quadrature_factor_ = 2
)
: local_assembler(local_assembler_),
lop(local_assembler_.localOperator()),
......@@ -132,7 +133,8 @@ public:
lfsv_volume(*pgfsv_volume),
lfsv_skeleton(*pgfsv_skeleton),
lfsv_skeleton_n(pgfsv_skeleton),
int_order_add(int_order_add_)
int_order_add(int_order_add_),
quadrature_factor(quadrature_factor_)
{}
/*-----------------------------------------------------------------------*//**
......@@ -149,7 +151,8 @@ public:
const LocalAssembler & local_assembler_,
const std::shared_ptr<const GFSVVolume> & gfsv_volume_,
const std::shared_ptr<const GFSVSkeleton> & gfsv_skeleton_,
int int_order_add_ = 2
int int_order_add_ = 2,
int quadrature_factor_ = 2
)
: local_assembler(local_assembler_),
lop(local_assembler_.localOperator()),
......@@ -160,7 +163,8 @@ public:
pgfsv_skeleton(gfsv_skeleton_),
lfsv_skeleton(*pgfsv_skeleton),
lfsv_skeleton_n(pgfsv_skeleton),
int_order_add(int_order_add_)
int_order_add(int_order_add_),
quadrature_factor(quadrature_factor_)
{}
//! @name Query methods for the global grid assembler
......@@ -405,9 +409,6 @@ public:
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
using LFSVVolumeRange = typename BasisSwitchTest::Range;
// Calculate the mass matrix
int quadrature_factor = 2; // FIXME!
auto gt = eg.geometry();
const int intorder = int_order_add+quadrature_factor*order;
......@@ -469,13 +470,6 @@ public:
rl_view,rn_view);
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersection().indexInInside();
double weight = 1.;
if (ig.intersection().conforming()){
auto sub_entity_s_id = ig.intersection().indexInInside();
auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
weight = sub_entity_s.geometry().volume()/ig.geometry().volume();
}
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
......@@ -513,9 +507,6 @@ public:
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
using LFSVSkeletonRange = typename BasisSwitchTest::Range;
// Calculate the mass matrix
int quadrature_factor = 2; // FIXME!
const int intorder = int_order_add + quadrature_factor * order_s;
auto gtface = ig.geometry();
......@@ -606,8 +597,6 @@ public:
using LFSVSkeletonRange = typename BasisSwitchTest::Range;
// Calculate the mass matrix
int quadrature_factor = 2; // FIXME!
const int intorder = int_order_add + quadrature_factor * order_s;
auto gtface = ig.geometry();
......@@ -742,6 +731,7 @@ private:
bool non_conforming_element_warning;
int int_order_add;
int quadrature_factor;
};
} // namespace Dorie
......
......@@ -349,6 +349,13 @@ public:
*/
bool empty() const {return _list.empty();}
/*-----------------------------------------------------------------------*//**
* @brief Returns the number of elements in the queue.
*
* @return The number of elements in the underlying container.
*/
size_type size() const {return _list.size();}
/*-----------------------------------------------------------------------*//**
* @brief Evaluation of the grid function.
......
......@@ -47,7 +47,7 @@ public:
, _log(create_logger(log_name, helper, spdlog::level::from_str(log_level)))
{ }
~SimulationBase () = default;
virtual ~SimulationBase () = default;
/*-----------------------------------------------------------------------*//**
* @brief Sets the output policy.
*
......
......@@ -87,11 +87,6 @@ protected:
enum class LOPCase {DG,RTVolume,RTSkeleton}; //!< Local operator case
bool sym_term; //!< Symmetry term switch
bool cs_term; //!< Consistency term switch
bool pty_term; //!< Penalty term switch
bool diff_term; //!< Gradient term switch
public:
// pattern assembly flags
......@@ -139,11 +134,7 @@ public:
quadrature_factor(quadrature_factor_),
penalty_factor(config.get<RF>("numerics.penaltyFactor")),
time(0.0),
cache(20),
sym_term(true),
cs_term(true),
pty_term(true),
diff_term(true)
cache(20)
{
theta = 0.; // IIP
if (method == DGMethod::NIPG || method == DGMethod::OOB)
......@@ -169,9 +160,6 @@ public:
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
{
if (not diff_term)
return;
// Switches between local and global interface in test local space
using FETest = typename LFSV::Traits::FiniteElementType;
using BasisTest = typename FETest::Traits::LocalBasisType;
......@@ -219,7 +207,7 @@ public:
for (unsigned int i = 0; i<lfsu.size(); i++)
jac.mv(js[i][0],gradphiu[i]);
std::vector<Vector> gradphiv;
std::vector<Vector> gradphiv(lfsv.size());
// discrete gradient sign.
double dg_sign = 1.;
......@@ -402,7 +390,7 @@ public:
y_s /= detB_s;
BTransposed_s.mtv(y_s,gradphiv_s[i]);
}
auto BTransposed_n = ig.inside().geometry().jacobianTransposed(p_local_n);
auto BTransposed_n = ig.outside().geometry().jacobianTransposed(p_local_n);
auto detB_n = BTransposed_n.determinant();
for (unsigned int i=0; i<lfsv_n.size(); i++) {
Vector y_n(gradphiv_n[i]);
......@@ -474,12 +462,10 @@ public:
// compute numerical flux
RF numFlux_s(0.);
RF numFlux_n(0.);
if (pty_term)
numFlux_s = numFlux_n = penalty * jump;
if (cs_term) {
numFlux_s -= gradu_s * normal;
numFlux_n -= gradu_n * normal;
}
numFlux_s = numFlux_n = penalty * jump;
numFlux_s -= gradu_s * normal;
numFlux_n -= gradu_n * normal;
numFlux_s *= satCond_s;
numFlux_n *= satCond_n;
......@@ -543,8 +529,6 @@ public:
if constexpr (lopcase != LOPCase::RTSkeleton)
{
if (not sym_term)
continue;
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
......@@ -697,8 +681,6 @@ public:
// flux is given parallel to gravity vector
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );
if (not diff_term)
continue;
// update residual (flux term)
if constexpr (lopcase != LOPCase::RTVolume)
for (unsigned int i = 0; i<lfsv_s.size(); i++)
......@@ -734,12 +716,7 @@ public:
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux
RF numFlux(0.);
if (pty_term)
numFlux = penalty * jump;
if (cs_term)
numFlux -= gradu_s * normal;
numFlux *= satCond_s;
RF numFlux = satCond_s*(penalty * jump - gradu_s * normal);
// compute conductivity factor
RF relCond;
......@@ -770,8 +747,6 @@ public:
if constexpr (lopcase != LOPCase::RTSkeleton)
{
if (not sym_term)
continue;
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
......@@ -863,30 +838,6 @@ public:
{
time = t;
}
//! Switch symmetry term on/off in alpha methods.
void symmetry_term(bool t)
{
sym_term = t;
}
//! Switch diffusion term on/off in alpha methods.
void diffusion_term(bool t)
{
diff_term = t;
}
//! Switch consistency term on/off in alpha methods.
void consistency_term(bool t)
{
cs_term = t;
}
//! Switch penalty term on/off in alpha methods.
void penalty_term(bool t)
{
pty_term = t;
}
};
/**
......
......@@ -6,8 +6,8 @@ namespace Dorie{
template<typename Traits>
RichardsSimulation<Traits>::RichardsSimulation (
Dune::ParameterTree& _inifile,
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
const GridCreator& _grid_creator,
Dune::MPIHelper& _helper
):
SimulationBase(log_richards,
_inifile.get<std::string>("output.logLevel"),
......@@ -23,12 +23,11 @@ RichardsSimulation<Traits>::RichardsSimulation (
mbe_tlop(1),
time_before(0.0),
dt_before(0.0),
enable_fluxrc(_inifile.get<bool>("numerics.fluxReconstruction")),
enable_lifting(_inifile.get<bool>("numerics.lifting"))
enable_fluxrc(_inifile.get<bool>("fluxReconstruction.enable"))
{
// --- Grid Function Space ---
this->_log->trace("Setting up GridFunctionSpace");
gfs = std::make_unique<GFS>(GFSHelper::create(gv));
gfs = std::make_shared<GFS>(GFSHelper::create(gv));
gfs->name("matric_head");
cc = std::make_shared<CC>();
Dune::PDELab::constraints(*gfs,*cc,false);
......@@ -52,13 +51,13 @@ RichardsSimulation<Traits>::RichardsSimulation (
const auto upwinding = std::get<OP::DGUpwinding::Type>(settings);
const auto weights = std::get<OP::DGWeights::Type>(settings);
slop = std::make_unique<SLOP>(inifile, *fparam, *fboundary, *fsource,
slop = std::make_shared<SLOP>(inifile, *fparam, *fboundary, *fsource,
method, upwinding, weights);
#else
slop = std::make_unique<SLOP>(inifile, *fparam, *fboundary, *fsource);
slop = std::make_shared<SLOP>(inifile, *fparam, *fboundary, *fsource);
#endif // EXPERIMENTAL_DG_FEATURES
tlop = std::make_unique<TLOP>(inifile, *fparam);
tlop = std::make_shared<TLOP>(inifile, *fparam);
controller = std::make_unique<CalculationController>(
inifile, *fboundary, this->_log);
......@@ -66,23 +65,6 @@ RichardsSimulation<Traits>::RichardsSimulation (
u = std::make_shared<U>(*gfs,.0);
Dune::PDELab::interpolate(*finitial,*gfs,*u);
// --- Operator Setup ---
operator_setup();
update_adapters();
// --- Update Flux Reconstruction ---
if constexpr (enable_rt_engine)
if (enable_fluxrc)
try {
flux_rec_gf->update(*u);
} catch (Dune::Exception &e) {
this->_log->debug(" Flux Reconstruction failed");
this->_log->error(" Unexpected error in flux reconstruction: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in flux reconstruction");
}
// --- Utility Class Setup --- //
if (output_policy() != OutputPolicy::None)
{
......@@ -107,18 +89,9 @@ RichardsSimulation<Traits>::RichardsSimulation (
adaptivity = adaptivity_fac.create();
}
this->_log->info("Setup complete");
do_operator_setup = true;
if ((enable_fluxrc or enable_lifting) and not enable_rt_engine) {
if (enable_fluxrc)
this->_log->error(
" Flux reconstruction is not available for this configuration.");
if (enable_lifting)
this->_log->error(
" Local lifting is not available for this configuration.");
DUNE_THROW(Dune::NotImplemented,
"Flux reconstruction engine not implemented!");
}
this->_log->info("Setup complete");
}
template<typename Traits>
......@@ -158,17 +131,16 @@ void RichardsSimulation<Traits>::operator_setup()
gfs->update();
if constexpr (enable_rt_engine) {
if (enable_fluxrc)
flux_rec_gf = std::make_shared<GFWaterFluxReconstruction>(go0);
if (enable_lifting)
lifting_gf = std::make_shared<GFWaterFluxReconstruction>(go0);
}
do_operator_setup = false;
}
template<typename Traits>
void RichardsSimulation<Traits>::step()
{
// --- Operator Setup ---
if (do_operator_setup)
operator_setup();
// long time step log for level 'debug'
if (this->_log->should_log(spdlog::level::debug)) {
this->_log->info("Time Step {}:",
......@@ -178,6 +150,7 @@ void RichardsSimulation<Traits>::step()
bool step_succeed = false;
while(not step_succeed)
{
// Obtain time variables
const RF time = controller->getTime();
const RF dt = controller->getDT();
bool exception = false;
......@@ -194,6 +167,7 @@ void RichardsSimulation<Traits>::step()
_log->trace(" Allowed iterations for Newton solver: {}", iter);
pdesolver->setMaxIterations(iter);
// create a DOF vector for next step
std::shared_ptr<U> unext = std::make_shared<U>(*u);
dwarn.push(false);
......@@ -235,36 +209,6 @@ void RichardsSimulation<Traits>::step()
DUNE_THROW(Dune::Exception, "Critical error in Newton solver");
}
// --- Update Flux Reconstruction ---
if constexpr (enable_rt_engine) {
if (enable_fluxrc)
try {
// compute flux reconstruction
flux_rec_gf->update(*u);
} catch (Dune::Exception &e) {
this->_log->debug(" Flux Reconstruction failed");
this->_log->error(" Unexpected error in flux reconstruction: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in flux reconstruction");
}
if (enable_lifting)
try {
// compute lifting
slop->penalty_term(false);
slop->symmetry_term(false);
slop->diffusion_term(false);
lifting_gf->update(*u);
slop->penalty_term(true);
slop->symmetry_term(true);
slop->diffusion_term(true);
} catch (Dune::Exception &e) {
this->_log->debug(" Local lifting failed");
this->_log->error(" Unexpected error in local lifting: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in local lifting");
}
}
// saving last times for adaptivity
time_before = time;
dt_before = dt;
......@@ -272,6 +216,7 @@ void RichardsSimulation<Traits>::step()
// controller reacts to outcome of solution
step_succeed = controller->validate(exception);
}
if (this->output_policy() == OutputPolicy::EndOfStep)
write_data();
}
......@@ -301,17 +246,6 @@ template<typename Traits>
void RichardsSimulation<Traits>::post_adapt_grid()
{
operator_setup();
update_adapters();
}
template<typename Traits>
void RichardsSimulation<Traits>::update_adapters () const
{
matric_head_gf = std::make_shared<GFMatricHead>(gfs, u);
flux_gf = std::make_shared<GFWaterFlux>(gfs, u, fparam);
cond_gf = std::make_shared<GFConductivity>(gv, fparam);
water_cont_gf = std::make_shared<GFWaterContent>(matric_head_gf, gv, fparam);
sat_gf = std::make_shared<GFSaturation>(matric_head_gf, gv, fparam);
}
template<typename Traits>
......@@ -319,36 +253,40 @@ void RichardsSimulation<Traits>::write_data () const
{
if (vtkwriter)
{
update_adapters();