diff --git a/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh b/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..905b6c534931e0a4777e6c043d3dab6bfa301c20
--- /dev/null
+++ b/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh
@@ -0,0 +1,367 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*
+ * \file
+ * \brief Test for the OnePTwoCModel in combination with the NI model for a convection problem:
+ * 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_1P2CNI_CONVECTION_TEST_PROBLEM_HH
+#define DUMUX_1P2CNI_CONVECTION_TEST_PROBLEM_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/porousmediumflow/1pnc/implicit/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/h2on2.hh>
+#include <dumux/material/components/h2o.hh>
+#include "1pnctestspatialparams.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class OnePTwoCNIConvectionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(OnePTwoCNIConvectionProblem, INHERITS_FROM(OnePNCNI));
+NEW_TYPE_TAG(OnePTwoCNIConvectionCCTpfaProblem, INHERITS_FROM(CCTpfaModel, OnePTwoCNIConvectionProblem));
+NEW_TYPE_TAG(OnePTwoCNIConvectionCCMpfaProblem, INHERITS_FROM(CCMpfaModel, OnePTwoCNIConvectionProblem));
+NEW_TYPE_TAG(OnePTwoCNIConvectionBoxProblem, INHERITS_FROM(BoxModel, OnePTwoCNIConvectionProblem));
+
+// Set the grid type
+#if HAVE_UG
+SET_TYPE_PROP(OnePTwoCNIConvectionProblem, Grid, Dune::UGGrid<2>);
+#else
+SET_TYPE_PROP(OnePTwoCNIConvectionProblem, Grid, Dune::YaspGrid<2>);
+#endif
+
+// Set the problem property
+SET_TYPE_PROP(OnePTwoCNIConvectionProblem, Problem, OnePTwoCNIConvectionProblem<TypeTag>);
+
+// Set fluid configuration
+SET_TYPE_PROP(OnePTwoCNIConvectionProblem,
+ FluidSystem,
+ FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), true>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(OnePTwoCNIConvectionProblem, SpatialParams, OnePNCTestSpatialParams<TypeTag>);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(OnePTwoCNIConvectionProblem, UseMoles, true);
+}
+
+
+/*!
+ * \ingroup OnePTwoCModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief Test for the OnePTwoCModel in combination with the NI model for a convection problem:
+ * The simulation domain is a tube where water with an elevated temperature is injected
+ * at a constant rate on the left hand side.
+ *
+ * Initially the domain is fully saturated with water at a constant temperature.
+ * On the left hand side water is injected at a constant rate and on the right hand side
+ * a Dirichlet boundary with constant pressure, saturation and temperature is applied.
+ *
+ * The results are compared to an analytical solution where a retarded front velocity is calculated as follows:
+ \f[
+ v_{Front}=\frac{q S_{water}}{\phi S_{total}}
+ \f]
+ *
+ * The result of the analytical solution is written into the vtu files.
+ *
+ * This problem uses the \ref OnePModel and \ref NIModel model.
+ *
+ * To run the simulation execute the following line in shell: <br>
+ * <tt>./test_box1p2cniconvection -ParameterFile ./test_box1p2cniconvection.input</tt> or <br>
+ * <tt>./test_cc1p2cniconvection -ParameterFile ./test_cc1p2cniconvection.input</tt>
+ */
+
+template <class TypeTag>
+class OnePTwoCNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using IapwsH2O = H2O<Scalar>;
+
+ // copy some indices for convenience
+ enum
+ {
+ // indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ massOrMoleFracIdx = Indices::firstMoleFracIdx,
+ temperatureIdx = Indices::temperatureIdx,
+
+ // indices of the equations
+ conti0EqIdx = Indices::conti0EqIdx,
+ transportEqIdx = Indices::firstTransportEqIdx,
+ energyEqIdx = Indices::energyEqIdx
+ };
+
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+ static const auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static const bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+ OnePTwoCNIConvectionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ //initialize fluid system
+ FluidSystem::init();
+
+ // stating in the console whether mole or mass fractions are used
+ if(useMoles)
+ std::cout<<"problem uses mole fractions"<<std::endl;
+ else
+ std::cout<<"problem uses mass fractions"<<std::endl;
+
+ temperatureExact_.resize(fvGridGeometry->numDofs(), 290.0);
+
+ darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
+
+ temperatureHigh_ = 291.;
+ temperatureLow_ = 290.;
+ pressureHigh_ = 2e5;
+ pressureLow_ = 1e5;
+ }
+
+ //! get the analytical temperature
+ const std::vector<Scalar>& getExactTemperature()
+ {
+ return temperatureExact_;
+ }
+
+ //! udpate the analytical temperature
+ void updateExactTemperature(const SolutionVector& curSol, Scalar time)
+ {
+ const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin());
+
+ ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry());
+ const auto someInitSol = initialAtPos(someElement.geometry().center());
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(someElement);
+ const auto someScv = *(scvs(fvGeometry).begin());
+
+ VolumeVariables volVars;
+ volVars.update(someElemSol, *this, someElement, someScv);
+
+ const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
+ const auto densityW = volVars.density(phaseIdx);
+ const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
+ const auto storageW = densityW*heatCapacityW*porosity;
+ const auto densityS = this->spatialParams().solidDensity(someElement, someScv, someElemSol);
+ const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol);
+ const auto storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
+ std::cout << "storage: " << storageTotal << '\n';
+
+ using std::max;
+ time = max(time, 1e-10);
+ const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
+ std::cout << "retarded velocity: " << retardedFrontVelocity << '\n';
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto dofIdxGlobal = scv.dofIndex();
+ auto dofPosition = scv.dofPosition();
+ temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
+ }
+ }
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain [K].
+ *
+ * This problem assumes a temperature of 20 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 20; } // in [K]
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param globalPos The position for which the bc type should be evaluated
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+
+ if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
+ {
+ values.setAllDirichlet();
+ }
+ else
+ {
+ values.setAllNeumann();
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * \param globalPos The position for which the bc type should be evaluated
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values = initial_(globalPos);
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param values The neumann values for the conservation equations in units of
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scvf The sub control volume face
+ *
+ * For this method, the \a values parameter stores the flux
+ * in normal direction of each phase. Negative values mean influx.
+ * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
+ */
+ ResidualVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ ResidualVector flux(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+
+ if(globalPos[0] < eps_)
+ {
+ flux[pressureIdx] = -darcyVelocity_*elemVolVars[scv].molarDensity(phaseIdx);
+ flux[massOrMoleFracIdx] = -darcyVelocity_*elemVolVars[scv].molarDensity(phaseIdx)*elemVolVars[scv].moleFraction(phaseIdx, massOrMoleFracIdx);
+ flux[temperatureIdx] = -darcyVelocity_*elemVolVars[scv].density(phaseIdx)
+ *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scv].pressure(phaseIdx));
+ }
+
+ return flux;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \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.
+ *
+ * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
+ */
+ PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+ { return PrimaryVariables(0.0); }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The position for which the initial condition should be evaluated
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ { return initial_(globalPos); }
+
+ // \}
+private:
+
+ // the internal method for the initial condition
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables priVars;
+ priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
+ priVars[massOrMoleFracIdx] = 1e-10; // initial condition for the N2 molefraction
+ priVars[temperatureIdx] = temperatureLow_;
+ return priVars;
+ }
+ static constexpr Scalar eps_ = 1e-6;
+ Scalar temperatureHigh_;
+ Scalar temperatureLow_;
+ Scalar pressureHigh_;
+ Scalar pressureLow_;
+ Scalar darcyVelocity_;
+
+ std::vector<Scalar> temperatureExact_;
+ };
+
+} //end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/1pnc/implicit/CMakeLists.txt b/test/porousmediumflow/1pnc/implicit/CMakeLists.txt
index 7c57dcd5235591ae3a251833c69ca437b973efda..dd29407b3220678ef64bb4b4512daa9ce4b8bfdb 100644
--- a/test/porousmediumflow/1pnc/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1pnc/implicit/CMakeLists.txt
@@ -32,6 +32,7 @@ dune_add_test(NAME test_1p2c_mpfa
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1p2c_mpfa test_1p2c_fv.input -Problem.Name 1pnctest_mpfa -Vtk.AddVelocity false")
# non-isothermal tests
+# conduction
dune_add_test(NAME test_1p2cni_conduction_box
SOURCES test_1p2cni_conduction_fv.cc
COMPILE_DEFINITIONS TYPETAG=OnePTwoCNIConductionBoxProblem
@@ -62,3 +63,35 @@ dune_add_test(NAME test_1p2cni_conduction_mpfa
--files ${CMAKE_SOURCE_DIR}/test/references/1p2cniccmpfaconduction-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/1p2cni_conductiontest_mpfa-00006.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1p2cni_conduction_mpfa test_1p2cni_conduction_fv.input -Problem.Name 1p2cni_conductiontest_mpfa -Vtk.AddVelocity false")
+
+# convection
+dune_add_test(NAME test_1p2cni_convection_box
+ SOURCES test_1p2cni_convection_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=OnePTwoCNIConvectionBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1p2cniboxconvection-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p2cni_convectiontest_box-00010.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p2cni_convection_box test_1p2cni_convection_fv.input -Problem.Name 1p2cni_convectiontest_box"
+ --zeroThreshold {"velocity_w \(m/s\)":1e-9})
+
+dune_add_test(NAME test_1p2cni_convection_tpfa
+ SOURCES test_1p2cni_convection_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=OnePTwoCNIConvectionCCTpfaProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1p2cniccconduction-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p2cni_convectiontest_tpfa-00006.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p2cni_convection_tpfa test_1p2cni_convection_fv.input -Problem.Name 1p2cni_convectiontest_tpfa"
+ --zeroThreshold {"velocity_w \(m/s\)_0":1e-9})
+
+#Note: The mpfa model does not support velocity output yet. We therefore deactivate it for the test and use a copy of the tpfa reference
+# solution where the velocities have been removed.
+dune_add_test(NAME test_1p2cni_convection_mpfa
+ SOURCES test_1p2cni_convection_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=OnePTwoCNIConvectionCCMpfaProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1p2cniccmpfaconduction-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p2cni_convectiontest_mpfa-00006.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p2cni_convection_mpfa test_1p2cni_convection_fv.input -Problem.Name 1p2cni_convectiontest_mpfa -Vtk.AddVelocity false")
diff --git a/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc b/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
index 8192151acf6ad08a4b3b6b7c944824d3fe7fba9e..10c8e66beb2d5a8c1c4043dd72c2c248fb9a47e3 100644
--- a/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
+++ b/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
@@ -120,7 +120,6 @@
vtkWriter.write(0.0);
// instantiate time loop
- // auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
@@ -137,10 +136,6 @@
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
NewtonMethod<NewtonController, Assembler, LinearSolver> nonLinearSolver(newtonController, assembler, linearSolver);
- // set some check points for the time loop
- // timeLoop->setPeriodicCheckPoint(1);
- // timeLoop->setPeriodicCheckPoint(tEnd/10.0);
-
// time loop
timeLoop->start(); do
{
@@ -174,8 +169,7 @@
timeLoop->advanceTimeStep();
// write vtk output
- // if (timeLoop->isCheckPoint())
- vtkWriter.write(timeLoop->time());
+ vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
diff --git a/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.cc b/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.cc
new file mode 100644
index 0000000000000000000000000000000000000000..71f424fcb6da78b7808fa1c1ce01b301a4ab2a89
--- /dev/null
+++ b/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.cc
@@ -0,0 +1,226 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the 1pnc model
+ */
+ #include <config.h>
+
+ #include "1p2cniconvectionproblem.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/discretization/methods.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/common/parameterparser.hh>
+
+ #include <dumux/nonlinear/newtoncontroller.hh>
+ #include <dumux/linear/seqsolverbackend.hh>
+ #include <dumux/nonlinear/newtonmethod.hh>
+
+ #include <dumux/assembly/fvassembler.hh>
+
+ #include <dumux/io/vtkoutputmodule.hh>
+
+ int main(int argc, char** argv) try
+ {
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TYPETAG);
+
+ ////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////
+
+ // 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);
+
+ // initialize parameter tree
+ Parameters::init(argc, argv);
+
+ //////////////////////////////////////////////////////////////////////
+ // 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);
+ auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+ auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+
+ // intialize the vtk output module
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.addField(problem->getExactTemperature(), "temperatureExact");
+ vtkWriter.write(0.0);
+ // output every vtkOutputInterval time step
+ const auto vtkOutputInterval = getParam<int>("Problem.OutputInterval");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, 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 = UMFPackBackend<TypeTag>;
+ using LinearSolver = ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::NewtonController<TypeTag>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod<NewtonController, Assembler, LinearSolver> 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.");
+ }
+
+ // update the exact time temperature
+ 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
+ if (timeLoop->timeStepIndex()==0 || timeLoop->timeStepIndex() % vtkOutputInterval == 0 || timeLoop->willBeFinished())
+ 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)
+ 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 (...)
+ {
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+ }
diff --git a/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.input b/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.input
new file mode 100644
index 0000000000000000000000000000000000000000..052d612b11e1929cf1d4f60bedbcc8b146a77dce
--- /dev/null
+++ b/test/porousmediumflow/1pnc/implicit/test_1p2cni_convection_fv.input
@@ -0,0 +1,18 @@
+[TimeLoop]
+DtInitial = 1 # [s]
+TEnd = 3e4 # [s]
+MaxTimeStepSize = 1e3
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 20 1
+Cells = 80 1
+
+[Problem]
+Name = 1p2cniconvection # name passed to the output routines
+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
diff --git a/test/references/1p2cniboxconvection-reference.vtu b/test/references/1p2cniboxconvection-reference.vtu
index 2af57a45df6bcbebebefb7460b0cecaf7a012cef..93d76411a549c1612f25e75cd835a78d02873aea 100644
--- a/test/references/1p2cniboxconvection-reference.vtu
+++ b/test/references/1p2cniboxconvection-reference.vtu
@@ -2,8 +2,8 @@
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="80" NumberOfPoints="162">
- <PointData Scalars="P" Vectors="velocity">
- <DataArray type="Float32" Name="P" NumberOfComponents="1" format="ascii">
+ <PointData Scalars="pressuere" Vectors="velocity_w (m/s)">
+ <DataArray type="Float32" Name="pressure" NumberOfComponents="1" format="ascii">
121579 121315 121579 121315 121051 121051 120787 120787 120522 120522 120258 120258
119994 119994 119730 119730 119465 119465 119201 119201 118936 118936 118672 118672
118407 118407 118141 118141 117876 117876 117609 117609 117343 117343 117075 117075
@@ -35,7 +35,7 @@
2170.43 2170.43 1899.12 1899.12 1627.82 1627.82 1356.52 1356.52 1085.21 1085.21 813.911 813.911
542.607 542.607 271.304 271.304 0 0
</DataArray>
- <DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
+ <DataArray type="Float32" Name="velocity_w (m/s)" NumberOfComponents="3" format="ascii">
0.0001 5.50706e-18 0 0.0001 -3.44192e-19 0 0.0001 5.50706e-18 0 0.0001 -3.44192e-19 0
0.0001 3.09773e-18 0 0.0001 3.09773e-18 0 0.0001 -1.06699e-17 0 0.0001 -1.06699e-17 0
0.0001 2.75349e-18 0 0.0001 2.75349e-18 0 0.0001 -1.37669e-17 0 0.0001 -1.37669e-17 0