diff --git a/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh b/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
index fadc4d0e5c76a553d12a4f710c2076b7db4bffb3..b3c8ff08d2c2fe61cc768fea384727d5711a73e6 100644
--- a/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
@@ -27,9 +27,11 @@
#include <math.h>
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/richards/implicit/model.hh>
-#include <dumux/porousmediumflow/richards/implicit/problem.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include "richardsnispatialparams.hh"
@@ -44,7 +46,7 @@ namespace Properties
{
NEW_TYPE_TAG(RichardsNIConductionProblem, INHERITS_FROM(RichardsNI, RichardsNISpatialParams));
NEW_TYPE_TAG(RichardsNIConductionBoxProblem, INHERITS_FROM(BoxModel, RichardsNIConductionProblem));
-NEW_TYPE_TAG(RichardsNIConductionCCProblem, INHERITS_FROM(CCModel, RichardsNIConductionProblem));
+NEW_TYPE_TAG(RichardsNIConductionCCProblem, INHERITS_FROM(CCTpfaModel, RichardsNIConductionProblem));
// Set the grid type
SET_TYPE_PROP(RichardsNIConductionProblem, Grid, Dune::YaspGrid<2>);
@@ -85,19 +87,23 @@ SET_TYPE_PROP(RichardsNIConductionProblem,
* <tt>./test_ccrichardsniconduction -ParameterFile ./test_ccrichardsniconduction.input</tt>
*/
template <class TypeTag>
-class RichardsNIConductionProblem : public RichardsProblem<TypeTag>
+class RichardsNIConductionProblem :public PorousMediumFlowProblem<TypeTag>
{
- typedef RichardsProblem<TypeTag> ParentType;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using IapwsH2O = H2O<Scalar>;
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
@@ -106,19 +112,14 @@ class RichardsNIConductionProblem : public RichardsProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dimWorld : 0 };
enum {
// indices of the primary variables
- pressureIdx = Indices::pwIdx,
+ pressureIdx = Indices::pressureIdx,
+ conti0EqIdx = Indices::conti0EqIdx,
+ wPhaseOnly = Indices::wPhaseOnly,
wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- temperatureIdx = Indices::temperatureIdx
- };
- enum {
- // index of the transport equation
- contiEqIdx = Indices::contiEqIdx,
+ temperatureIdx = Indices::temperatureIdx,
energyEqIdx = Indices::energyEqIdx
};
@@ -128,25 +129,15 @@ class RichardsNIConductionProblem : public RichardsProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
- RichardsNIConductionProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ RichardsNIConductionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
//initialize fluid system
FluidSystem::init();
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
- outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval);
+ name_ = getParam<std::string>("Problem.Name");
temperatureHigh_ = 300.;
- }
-
-
- bool shouldWriteOutput() const
- {
- return
- this->timeManager().timeStepIndex() == 0 ||
- this->timeManager().timeStepIndex() % outputInterval_ == 0 ||
- this->timeManager().episodeWillBeFinished() ||
- this->timeManager().willBeFinished();
+ temperatureExact_.resize(fvGridGeometry->numDofs());
}
/*!
@@ -154,57 +145,53 @@ public:
* from the solution of the current time step to the VTK
* writer.
*/
- void addOutputVtkFields()
+ //! get the analytical temperature
+ const std::vector<Scalar>& getExactTemperature()
{
- //Here we calculate the analytical solution
- typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
- unsigned numDofs = this->model().numDofs();
+ return temperatureExact_;
+ }
- //create required scalar fields
- ScalarField *temperatureExact = this->resultWriter().allocateManagedBuffer(numDofs);
+ //! udpate the analytical temperature
+ void updateExactTemperature(const SolutionVector& curSol, Scalar time)
+ {
+ const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin());
- FVElementGeometry fvGeometry;
- VolumeVariables volVars;
+ ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry());
+ const auto someInitSol = initialAtPos(someElement.geometry().center());
- const auto firstElement = *this->gridView().template begin<0>();
- fvGeometry.update(this->gridView(), firstElement);
- PrimaryVariables initialPriVars(0);
- GlobalPosition globalPos(0);
- initial_(initialPriVars, globalPos);
-
- //update the constant volume variables
- volVars.update(initialPriVars, *this, firstElement, fvGeometry, 0, false);
-
- Scalar porosity = this->spatialParams().porosity(firstElement, fvGeometry, 0);
- Scalar densityW = volVars.density(wPhaseIdx);
- Scalar heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
- Scalar densityS = this->spatialParams().solidDensity(firstElement, fvGeometry, 0);
- Scalar heatCapacityS = this->spatialParams().solidHeatCapacity(firstElement, fvGeometry, 0);
- Scalar storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
- Scalar effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
- firstElement, fvGeometry, 0);
- using std::max;
- Scalar time = max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(someElement);
+ const auto someScv = *(scvs(fvGeometry).begin());
- for (const auto& element : elements(this->gridView()))
+ VolumeVariables volVars;
+ volVars.update(someElemSol, *this, someElement, someScv);
+
+ const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
+ const auto densityW = volVars.density(wPhaseIdx);
+ const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
+ const auto densityS = this->spatialParams().solidDensity(someElement, someScv, someElemSol);
+ const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
+ const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
+ const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
+ someElement, fvGeometry, someScv);
+ using std::max;
+ time = max(time, 1e-10);
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
- fvGeometry.update(this->gridView(), element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int globalIdx = this->model().dofMapper().subIndex(element, scvIdx, dofCodim);
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = element.geometry().center();
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto globalIdx = scv.dofIndex();
+ const auto& globalPos = scv.dofPosition();
+ using std::erf;
+ using std::sqrt;
+ temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
+ *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
- using std::sqrt;
- using std::erf;
- (*temperatureExact)[globalIdx] = temperatureHigh_ + (initialPriVars[temperatureIdx] - temperatureHigh_)
- *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
}
}
- this->resultWriter().attachDofData(*temperatureExact, "temperatureExact", isBox);
}
/*!
* \name Problem parameters
@@ -231,13 +218,12 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
- * \param values The boundary types for the conservation equations
- * \param globalPos The position for which the bc type should be evaluated
+ * \param globalPos The position for which the boundary type is set
*/
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
- if(globalPos[0] < eps_ || globalPos[0] > this->bBoxMax()[0] - eps_)
+ BoundaryTypes values;
+ if(globalPos[0] < eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
{
values.setAllDirichlet();
}
@@ -245,6 +231,7 @@ public:
{
values.setAllNeumann();
}
+ return values;
}
/*!
@@ -256,14 +243,16 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
- void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ PrimaryVariables values(0.0);
+ values = initial_(globalPos);
if (globalPos[0] < eps_)
{
values[temperatureIdx] = temperatureHigh_;
}
+ return values;
}
/*!
@@ -274,14 +263,10 @@ public:
* in normal direction of each component. Negative values mean
* influx.
*/
- void neumann(PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- const int scvIdx,
- const int boundaryFaceIdx) const
+ PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
- priVars = 0;
+ PrimaryVariables values(0.0);
+ return values;
}
// \}
@@ -291,21 +276,6 @@ public:
*/
// \{
- /*!
- * \brief Evaluate the source term for all phases within a given
- * sub-control-volume.
- *
- * For this method, the \a priVars parameter stores the rate mass
- * of a component is generated or annihilate per volume
- * unit. Positive values mean that mass is created, negative ones
- * mean that it vanishes.
- */
- void sourceAtPos(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
- {
- priVars = Scalar(0.0);
- }
-
/*!
* \brief Returns the reference pressure [Pa] of the non-wetting
* fluid phase within a finite volume
@@ -318,9 +288,7 @@ public:
* \param scvIdx The sub control volume index inside the finite
* volume geometry
*/
- Scalar referencePressure(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar nonWettingReferencePressure() const
{ return 1e5; };
/*!
@@ -332,26 +300,29 @@ public:
* For this method, the \a values parameter stores primary
* variables.
*/
- void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initial_(globalPos);
}
// \}
private:
- // the internal method for the initial condition
- void initial_(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
+ PrimaryVariables priVars(0.0);
+ priVars.setState(wPhaseOnly);
priVars[pressureIdx] = 1e5; // initial condition for the pressure
+
priVars[temperatureIdx] = 290.;
+ return priVars;
}
Scalar temperatureHigh_;
static constexpr Scalar eps_ = 1e-6;
std::string name_;
int outputInterval_;
+ std::vector<Scalar> temperatureExact_;
};
} //end namespace
diff --git a/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh b/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
index 951edf57cac861ff0487131398b5f0670c145501..1472b7a585583d7c6b77c8a2dca260af91fa9e91 100644
--- a/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
@@ -23,14 +23,16 @@
* The simulation domain is a tube where water with an elevated temperature is injected
* at a constant rate on the left hand side.
*/
-#ifndef DUMUX_1PNI_CONVECTION_PROBLEM_HH
-#define DUMUX_1PNI_CONVECTION_PROBLEM_HH
+#ifndef DUMUX_RICHARDS_CONVECTION_PROBLEM_HH
+#define DUMUX_RICHARDS_CONVECTION_PROBLEM_HH
#include <math.h>
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/richards/implicit/model.hh>
-#include <dumux/porousmediumflow/richards/implicit/problem.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include "richardsnispatialparams.hh"
@@ -45,7 +47,7 @@ namespace Properties
{
NEW_TYPE_TAG(RichardsNIConvectionProblem, INHERITS_FROM(RichardsNI, RichardsNISpatialParams));
NEW_TYPE_TAG(RichardsNIConvectionBoxProblem, INHERITS_FROM(BoxModel, RichardsNIConvectionProblem));
-NEW_TYPE_TAG(RichardsNIConvectionCCProblem, INHERITS_FROM(CCModel, RichardsNIConvectionProblem));
+NEW_TYPE_TAG(RichardsNIConvectionCCProblem, INHERITS_FROM(CCTpfaModel, RichardsNIConvectionProblem));
// Set the grid type
SET_TYPE_PROP(RichardsNIConvectionProblem, Grid, Dune::YaspGrid<2>);
@@ -89,20 +91,24 @@ SET_TYPE_PROP(RichardsNIConvectionProblem,
* <tt>./test_ccrichardsniconvection -ParameterFile ./test_ccrichardsniconvection.input</tt>
*/
template <class TypeTag>
-class RichardsNIConvectionProblem : public RichardsProblem<TypeTag>
+class RichardsNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
{
- typedef RichardsProblem<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef H2O<Scalar> IapwsH2O;
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using IapwsH2O = H2O<Scalar>;
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
@@ -111,19 +117,14 @@ class RichardsNIConvectionProblem : public RichardsProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dimWorld : 0 };
enum {
// indices of the primary variables
- pressureIdx = Indices::pwIdx,
+ pressureIdx = Indices::pressureIdx,
+ conti0EqIdx = Indices::conti0EqIdx,
+ wPhaseOnly = Indices::wPhaseOnly,
wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- temperatureIdx = Indices::temperatureIdx
- };
- enum {
- // index of the transport equation
- contiEqIdx = Indices::contiEqIdx,
+ temperatureIdx = Indices::temperatureIdx,
energyEqIdx = Indices::energyEqIdx
};
@@ -133,86 +134,74 @@ class RichardsNIConvectionProblem : public RichardsProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
- RichardsNIConvectionProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ RichardsNIConvectionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
//initialize fluid system
FluidSystem::init();
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
- outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval);
- darcyVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, DarcyVelocity);
+ name_ = getParam<std::string>("Problem.Name");
+ darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
temperatureHigh_ = 291.;
temperatureLow_ = 290.;
pressureHigh_ = 2e5;
pressureLow_ = 1e5;
+ temperatureExact_.resize(fvGridGeometry->numDofs());
}
- bool shouldWriteOutput() const
- {
- return
- this->timeManager().timeStepIndex() == 0 ||
- this->timeManager().timeStepIndex() % outputInterval_ == 0 ||
- this->timeManager().episodeWillBeFinished() ||
- this->timeManager().willBeFinished();
- }
-
- /*!
+ /*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
* writer.
*/
- void addOutputVtkFields()
+ //! get the analytical temperature
+ const std::vector<Scalar>& getExactTemperature()
{
- //Here we calculate the analytical solution
- typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
- unsigned numDofs = this->model().numDofs();
+ return temperatureExact_;
+ }
- //create required scalar fields
- ScalarField *temperatureExact = this->resultWriter().allocateManagedBuffer(numDofs);
+ //! udpate the analytical temperature
+ void updateExactTemperature(const SolutionVector& curSol, Scalar time)
+ {
+ const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin());
- FVElementGeometry fvGeometry;
- VolumeVariables volVars;
+ ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry());
+ const auto someInitSol = initialAtPos(someElement.geometry().center());
- const auto firstElement = *this->gridView().template begin<0>();
- fvGeometry.update(this->gridView(), firstElement);
- PrimaryVariables initialPriVars(0);
- GlobalPosition globalPos(0);
- initial_(initialPriVars, globalPos);
-
- //update the constant volume variables
- volVars.update(initialPriVars, *this, firstElement, fvGeometry, 0, false);
-
- Scalar porosity = this->spatialParams().porosity(firstElement, fvGeometry, 0);
- Scalar densityW = volVars.density(wPhaseIdx);
- Scalar heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
- Scalar storageW = densityW*heatCapacityW*porosity;
- Scalar densityS = this->spatialParams().solidDensity(firstElement, fvGeometry, 0);
- Scalar heatCapacityS = this->spatialParams().solidHeatCapacity(firstElement, fvGeometry, 0);
- Scalar storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
- std::cout<<"storage: "<<storageTotal<<std::endl;
- using std::max;
- Scalar time = max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
- Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
- std::cout<<"retarded velocity: "<<retardedFrontVelocity<<std::endl;
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(someElement);
+ const auto someScv = *(scvs(fvGeometry).begin());
- for (const auto& element : elements(this->gridView()))
+ VolumeVariables volVars;
+ volVars.update(someElemSol, *this, someElement, someScv);
+
+ const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
+ const auto densityW = volVars.density(wPhaseIdx);
+ const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
+ const auto densityS = this->spatialParams().solidDensity(someElement, someScv, someElemSol);
+ const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
+ const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
+ const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
+ someElement, fvGeometry, someScv);
+ using std::max;
+ time = max(time, 1e-10);
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
- fvGeometry.update(this->gridView(), element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int globalIdx = this->model().dofMapper().subIndex(element, scvIdx, dofCodim);
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = element.geometry().center();
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto globalIdx = scv.dofIndex();
+ const auto& globalPos = scv.dofPosition();
+ using std::erf;
+ using std::sqrt;
+ temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
+ *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
- (*temperatureExact)[globalIdx] = globalPos[0] < retardedFrontVelocity*time ? temperatureHigh_ : temperatureLow_;
}
}
- this->resultWriter().attachDofData(*temperatureExact, "temperatureExact", isBox);
}
/*!
* \name Problem parameters
@@ -240,13 +229,12 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
- * \param values The boundary types for the conservation equations
- * \param globalPos The position for which the bc type should be evaluated
+ * \param globalPos The position for which the boundary type is set
*/
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
- if(globalPos[0] > this->bBoxMax()[0] - eps_)
+ BoundaryTypes values;
+ if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
{
values.setAllDirichlet();
}
@@ -254,6 +242,7 @@ public:
{
values.setAllNeumann();
}
+ return values;
}
/*!
@@ -265,50 +254,36 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
- void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initial_(globalPos);
}
/*!
- * \brief Evaluates the boundary conditions for a Neumann
- * boundary segment in dependency on the current solution.
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
*
- * \param values Stores the Neumann values for the conservation equations in
- * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param element The finite element
- * \param fvGeometry The finite volume geometry of the element
- * \param intersection The intersection between element and boundary
- * \param scvIdx The local index of the sub-control volume
- * \param boundaryFaceIdx The index of the boundary face
- * \param elemVolVars All volume variables for the element
- *
- * This method is used for cases, when the Neumann condition depends on the
- * solution and requires some quantities that are specific to the fully-implicit method.
- * The \a values store the mass flux of each phase normal to the boundary.
- * Negative values indicate an inflow.
+ * \param fvGeometry The finite-volume geometry in the box scheme
+ * \param elemVolVars The element volume variables
+ * \param scvf The subcontrolvolume face
+ * Negative values mean influx.
*/
- void solDependentNeumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- const int scvIdx,
- const int boundaryFaceIdx,
- const ElementVolumeVariables &elemVolVars) const
+ PrimaryVariables neumann(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
{
- values = 0;
- GlobalPosition globalPos(0);
- if (isBox)
- globalPos = fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
- else
- globalPos = fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
+ PrimaryVariables values(0.0);
+ const auto globalPos = scvf.ipGlobal();
if(globalPos[0] < eps_)
{
- values[pressureIdx] = -darcyVelocity_*elemVolVars[scvIdx].density(wPhaseIdx);
- values[temperatureIdx] = -darcyVelocity_*elemVolVars[scvIdx].density(wPhaseIdx)
- *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scvIdx].pressure(wPhaseIdx));
+ values[pressureIdx] = -darcyVelocity_*elemVolVars[scvf.insideScvIdx()].density(wPhaseIdx);
+ values[temperatureIdx] = -darcyVelocity_*elemVolVars[scvf.insideScvIdx()].density(wPhaseIdx)
+ *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scvf.insideScvIdx()].pressure(wPhaseIdx));
}
+ return values;
}
// \}
@@ -318,20 +293,6 @@ public:
*/
// \{
- /*!
- * \brief Evaluate the source term for all phases within a given
- * sub-control-volume.
- *
- * For this method, the \a priVars parameter stores the rate mass
- * of a component is generated or annihilate per volume
- * unit. Positive values mean that mass is created, negative ones
- * mean that it vanishes.
- */
- void sourceAtPos(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
- {
- priVars = Scalar(0.0);
- }
/*!
* \brief Returns the reference pressure [Pa] of the non-wetting
@@ -345,9 +306,7 @@ public:
* \param scvIdx The sub control volume index inside the finite
* volume geometry
*/
- Scalar referencePressure(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar nonWettingReferencePressure() const
{ return 1e5; };
/*!
@@ -359,20 +318,21 @@ public:
* For this method, the \a values parameter stores primary
* variables.
*/
- void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initial_(globalPos);
}
// \}
private:
- // the internal method for the initial condition
- void initial_(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
+ PrimaryVariables priVars(0.0);
+ priVars.setState(wPhaseOnly);
priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
priVars[temperatureIdx] = temperatureLow_;
+ return priVars;
}
Scalar temperatureHigh_;
@@ -382,7 +342,7 @@ private:
Scalar darcyVelocity_;
static constexpr Scalar eps_ = 1e-6;
std::string name_;
- int outputInterval_;
+ std::vector<Scalar> temperatureExact_;
};
} //end namespace
diff --git a/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh b/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
index ee40841c2f0f84bdbf1a2b2346deca39148813ad..3986b00d0dc43408243089e2461dc6b8c6dea83f 100644
--- a/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
+++ b/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
@@ -70,25 +70,33 @@ public:
template<class TypeTag>
class RichardsNISpatialParams : public ImplicitSpatialParams<TypeTag>
{
- typedef ImplicitSpatialParams<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GridView::template Codim<0>::Entity Element;
+ using ParentType = ImplicitSpatialParams<TypeTag>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using Element = typename GridView::template Codim<0>::Entity;
+ enum {
+ dimWorld=GridView::dimensionworld
+ };
+ using CoordScalar = typename Grid::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
//typedef LinearMaterial<Scalar> EffMaterialLaw;
public:
+ // export permeability type
+ using PermeabilityType = Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename MaterialLaw::Params MaterialLawParams;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
- RichardsNISpatialParams(const GridView &gridView)
- : ParentType(gridView)
+ RichardsNISpatialParams(const Problem& problem)
+ : ParentType(problem)
{
permeability_ = 1e-10;
porosity_ = 0.4;
@@ -113,16 +121,6 @@ public:
{}
- /*!
- * \brief Update the spatial parameters with the flow solution
- * after a timestep.
- *
- * \param globalSolution the global solution vector
- */
- void update(const SolutionVector &globalSolution)
- {
- }
-
/*!
* \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
*
@@ -130,9 +128,7 @@ public:
* \param fvGeometry The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
- const Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
return permeability_;
}
@@ -144,9 +140,7 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
*/
- double porosity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
return porosity_;
}
@@ -158,9 +152,7 @@ public:
* \param fvGeometry The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
- const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition &globalPos) const
{
return materialParams_;
}
@@ -170,13 +162,9 @@ public:
*
* This is only required for non-isothermal models.
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume
+ * \param globalPos The global position
*/
- Scalar solidHeatCapacity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{
return 790; // specific heat capacity of granite [J / (kg K)]
}
@@ -186,13 +174,9 @@ public:
*
* This is only required for non-isothermal models.
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume
+ * \param globalPos The global position
*/
- Scalar solidDensity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
{
return 2700; // density of granite [kg/m^3]
}
@@ -200,14 +184,11 @@ public:
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume where
- * the heat capacity needs to be defined
+ * This is only required for non-isothermal models.
+ *
+ * \param globalPos The global position
*/
- Scalar solidThermalConductivity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{
return lambdaSolid_;
}
diff --git a/test/porousmediumflow/richards/implicit/test_boxrichards.cc b/test/porousmediumflow/richards/implicit/test_boxrichards.cc
index 4d71dc885c67a23a56134cffe7b90ce9cb86f212..5805bb62baea7c9561600731e9a71703bd0f702c 100644
--- a/test/porousmediumflow/richards/implicit/test_boxrichards.cc
+++ b/test/porousmediumflow/richards/implicit/test_boxrichards.cc
@@ -19,7 +19,7 @@
/*!
* \file
*
- * \brief Test for the Richards CC model.
+ * \brief Test for the Richards box model.
*/
#include <config.h>
diff --git a/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.cc b/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.cc
index f14cc608144274cf04ea0e5e0b8ebe0695be4e39..014685dda5ce55739d91cd5f0668be873c847c29 100644
--- a/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.cc
+++ b/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.cc
@@ -19,12 +19,35 @@
/*!
* \file
*
- * \brief test for the RichardsNI box model
+ * \brief Test for the Richards box model.
*/
#include <config.h>
+
#include "richardsniconductionproblem.hh"
-#include <dumux/common/start.hh>
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -45,13 +68,183 @@ void usage(const char *progName, const std::string &errorMsg)
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
+
std::cout << errorMessageOut
<< "\n";
}
}
-int main(int argc, char** argv)
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(RichardsNIConductionBoxProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = Dumux::ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::RichardsNewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+ // compute the new analytical temperature field for the output
+ problem->updateExactTemperature(x, timeLoop->time()+timeLoop->timeStepSize());
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
+
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
{
- typedef TTAG(RichardsNIConductionBoxProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
}
diff --git a/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.input b/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.input
index afed3ddd3cde6462f0baad1c93be5f6c65387179..65a74ff7aab30a415329a78b1515707db06bfa61 100644
--- a/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.input
+++ b/test/porousmediumflow/richards/implicit/test_boxrichardsniconduction.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1 # [s]
TEnd = 1e5 # [s]
MaxTimeStepSize = 1e10 # [s]
@@ -12,5 +12,5 @@ Name = test_boxrichardsniconduction # name passed to the output routines
OutputInterval = 5 # every 5th timestep an output file is written
EnableGravity= 0 # disable gravity
-[Vtk]
-AddVelocity = 1 # enable velocity output
+[Newton]
+EnableChop = false # chop for better convergence
diff --git a/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.cc b/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.cc
index 10becacae6b3d1e0b84eeb8f5354ee59dd312833..c6b60e4770210543de8d2218ae600a5d3f97c4db 100644
--- a/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.cc
+++ b/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.cc
@@ -19,12 +19,35 @@
/*!
* \file
*
- * \brief test for the RichardsNI box model
+ * \brief Test for the Richards box model.
*/
#include <config.h>
+
#include "richardsniconvectionproblem.hh"
-#include <dumux/common/start.hh>
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -45,13 +68,183 @@ void usage(const char *progName, const std::string &errorMsg)
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
+
std::cout << errorMessageOut
<< "\n";
}
}
-int main(int argc, char** argv)
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(RichardsNIConvectionBoxProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = Dumux::ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::RichardsNewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+ // compute the new analytical temperature field for the output
+ problem->updateExactTemperature(x, timeLoop->time()+timeLoop->timeStepSize());
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
+
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
{
- typedef TTAG(RichardsNIConvectionBoxProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
}
diff --git a/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.input b/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.input
index 98fd0d6c70570f9963ab0ee8c7b28f75c8d8c914..be59a892ccb35ed43722f57c706d2f07064079c8 100644
--- a/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.input
+++ b/test/porousmediumflow/richards/implicit/test_boxrichardsniconvection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1 # [s]
TEnd = 3e4 # [s]
MaxTimeStepSize = 1e3 # [s]
@@ -13,5 +13,5 @@ OutputInterval = 5 # every 5th timestep an output file is written
DarcyVelocity = 1e-4 # [m/s] inflow at the left boundary
EnableGravity = 0 # disable gravity
-[Vtk]
-AddVelocity = 1 # enable velocity output
+[Newton]
+EnableChop = false # chop for better convergence
diff --git a/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.cc b/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.cc
index 212f9d00847d999b175aaa5ddf1ebbed6f7bc733..07812769c37c38df271c776bd727c7a9a218124e 100644
--- a/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.cc
+++ b/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.cc
@@ -19,12 +19,35 @@
/*!
* \file
*
- * \brief test for the RichardsNI cc model
+ * \brief Test for the Richards CC model.
*/
#include <config.h>
+
#include "richardsniconductionproblem.hh"
-#include <dumux/common/start.hh>
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -45,13 +68,183 @@ void usage(const char *progName, const std::string &errorMsg)
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
+
std::cout << errorMessageOut
<< "\n";
}
}
-int main(int argc, char** argv)
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(RichardsNIConductionCCProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = Dumux::ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::RichardsNewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+ // compute the new analytical temperature field for the output
+ problem->updateExactTemperature(x, timeLoop->time()+timeLoop->timeStepSize());
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
+
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
{
- typedef TTAG(RichardsNIConductionCCProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
}
diff --git a/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.input b/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.input
index 9f1d1ae0be5cd13fd00c17a9484ef04e71cd921f..fa4803bb64ad929e831eb4079febe74dea1dc96f 100644
--- a/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.input
+++ b/test/porousmediumflow/richards/implicit/test_ccrichardsniconduction.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1 # [s]
TEnd = 1e5 # [s]
MaxTimeStepSize = 1e10 # [s]
@@ -12,5 +12,5 @@ Name = test_ccrichardsniconduction # name passed to the output routines
OutputInterval = 5 # every 5th timestep an output file is written
EnableGravity= 0 # disable gravity
-[Vtk]
-AddVelocity = 1 # enable velocity output
+[Newton]
+EnableChop = false # chop for better convergence
diff --git a/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.cc b/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.cc
index f9406ad117b5014363abaed929ac637a795ec45e..47b703e730d8b33830e57c115736f2b9608364ec 100644
--- a/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.cc
+++ b/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.cc
@@ -19,12 +19,35 @@
/*!
* \file
*
- * \brief test for the RichardsNI cc model
+ * \brief Test for the Richards CC model.
*/
#include <config.h>
+
#include "richardsniconvectionproblem.hh"
-#include <dumux/common/start.hh>
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -45,13 +68,183 @@ void usage(const char *progName, const std::string &errorMsg)
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
+
std::cout << errorMessageOut
<< "\n";
}
}
-int main(int argc, char** argv)
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(RichardsNIConvectionCCProblem);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = Dumux::ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::RichardsNewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+ // compute the new analytical temperature field for the output
+ problem->updateExactTemperature(x, timeLoop->time()+timeLoop->timeStepSize());
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
+
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
{
- typedef TTAG(RichardsNIConvectionCCProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
}
diff --git a/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.input b/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.input
index a8d961715fd7da74e3a207179802805a25312375..8d1e9d7bc29d4fe9e0a2f777c89c6fec0abab165 100644
--- a/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.input
+++ b/test/porousmediumflow/richards/implicit/test_ccrichardsniconvection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1 # [s]
TEnd = 3e4 # [s]
MaxTimeStepSize = 1e3 # [s]
@@ -13,5 +13,5 @@ OutputInterval = 5 # every 5th timestep an output file is written
DarcyVelocity = 1e-4 # [m/s] inflow at the left boundary
EnableGravity = 0 # disable gravity
-[Vtk]
-AddVelocity = 1 # enable velocity output
+[Newton]
+EnableChop = false # chop for better convergence