From 63425fe02a5cbd51d9b7aa0cdd10a06d4ab7179d Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Wed, 20 Dec 2017 16:17:36 +0100
Subject: [PATCH] [feature] port evaporation nonequilibrium test
---
.../implicit/evaporationatmosphereproblem.hh | 440 ++++-------------
.../evaporationatmospherespatialparams.hh | 443 ++++++++----------
.../mpnc/implicit/test_boxmpnckinetic.cc | 233 ++++++++-
.../mpnc/implicit/test_boxmpnckinetic.input | 23 +-
4 files changed, 522 insertions(+), 617 deletions(-)
diff --git a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
index e2d1721c97..de1fc1aced 100644
--- a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
@@ -37,9 +37,6 @@
#ifndef DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
#define DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
-// set this to one for looking at a different concept of mass transfer (see mpnclocalresidualmasskinetic)
-#define FUNKYMASSTRANSFER 0
-
// this sets that the relation using pc_max is used.
// i.e. - only parameters for awn, ans are given,
// - the fit for ans involves the maximum value for pc, where Sw, awn are zero.
@@ -48,19 +45,14 @@
#include <dune/common/parametertreeparser.hh>
-#include <dumux/porousmediumflow/mpnc/implicit/modelkinetic.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
-#include <dumux/porousmediumflow/mpnc/implicit/velomodelnewtoncontroller.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/porousmediumflow/mpnc/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2kinetic.hh>
#include <dumux/io/gnuplotinterface.hh>
#include "plotoverline2d.hh"
-#include <dumux/material/fluidstates/nonequilibrium.hh>
-//#include <dumux/material/fluidstates/nonequilibriumenergy.hh>
-//#include <dumux/material/fluidstates/nonequilibriummass.hh>
-//#include <dumux/material/fluidstates/compositional.hh>
-
#include "evaporationatmospherespatialparams.hh"
@@ -73,7 +65,8 @@ class EvaporationAtmosphereProblem;
namespace Properties
{
NEW_TYPE_TAG(EvaporationAtmosphereProblem,
- INHERITS_FROM(BoxMPNCKinetic, EvaporationAtmosphereSpatialParams));
+ INHERITS_FROM(MPNCNonequil, EvaporationAtmosphereSpatialParams));
+NEW_TYPE_TAG(EvaporationAtmosphereBoxProblem, INHERITS_FROM(BoxModel, EvaporationAtmosphereProblem));
// Set the grid type
SET_TYPE_PROP(EvaporationAtmosphereProblem, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
@@ -81,20 +74,12 @@ SET_TYPE_PROP(EvaporationAtmosphereProblem, Grid, Dune::YaspGrid<2, Dune::Tensor
// Set the problem property
SET_TYPE_PROP(EvaporationAtmosphereProblem,
Problem,
- EvaporationAtmosphereProblem<TTAG(EvaporationAtmosphereProblem)>);
+ EvaporationAtmosphereProblem<TypeTag>);
// Set fluid configuration
-SET_PROP(EvaporationAtmosphereProblem, FluidSystem)
-{
-private: using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public: using type = FluidSystems::H2ON2Kinetic<Scalar, /*useComplexRelations=*/false>;
-};
-
-// Set the newton controller
SET_TYPE_PROP(EvaporationAtmosphereProblem,
- NewtonController,
- VeloModelNewtonController<TypeTag>);
-
+ FluidSystem,
+ FluidSystems::H2ON2Kinetic<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/false>);
//! Set the default pressure formulation: either pw first or pn first
SET_INT_PROP(EvaporationAtmosphereProblem,
@@ -103,55 +88,6 @@ SET_INT_PROP(EvaporationAtmosphereProblem,
// Set the type used for scalar values
SET_TYPE_PROP(EvaporationAtmosphereProblem, Scalar, double);
-
-
-// Specify whether diffusion is enabled
-SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableDiffusion, true);
-
-// do not use a chopped newton method in the beginning
-SET_BOOL_PROP(EvaporationAtmosphereProblem, NewtonEnableChop, false);
-
-// Enable the re-use of the jacobian matrix whenever possible?
-SET_BOOL_PROP(EvaporationAtmosphereProblem, ImplicitEnableJacobianRecycling, true);
-
-//#################
-// Which Code to compile
-// Specify whether there is any energy equation
-SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableEnergy, true );
-// Specify whether the kinetic energy module is used
-SET_INT_PROP(EvaporationAtmosphereProblem, NumEnergyEquations, 3);
-// Specify whether the kinetic mass module is use
-SET_BOOL_PROP(EvaporationAtmosphereProblem, EnableKinetic, true);
-//#################
-
-/*!
- * \brief The fluid state which is used by the volume variables to
- * store the thermodynamic state. This has to be chosen
- * appropriately for the model ((non-)isothermal, equilibrium, ...).
- */
-SET_PROP(EvaporationAtmosphereProblem, FluidState){
- private: using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- private: using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-// public: using type = NonEquilibriumEnergyFluidState<Scalar, FluidSystem>;
-// public: using type = NonEquilibriumMassFluidState<Scalar, FluidSystem>;
- public: using type = NonEquilibriumFluidState<Scalar, FluidSystem>;
-// public: using type = CompositionalFluidState<Scalar, FluidSystem>;
-};
-
-SET_BOOL_PROP(EvaporationAtmosphereProblem, UseMaxwellDiffusion, false);
-
-// Output variables
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddVelocities, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddReynolds, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddPrandtl, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddNusselt, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddBoundaryTypes, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddInterfacialArea, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddTemperatures, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddxEquil, true);
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VtkAddDeltaP, true);
-
-SET_BOOL_PROP(EvaporationAtmosphereProblem, VelocityAveragingInModel, true);
}
/*!
@@ -161,24 +97,35 @@ SET_BOOL_PROP(EvaporationAtmosphereProblem, VelocityAveragingInModel, true);
* a porous medium sub-domain into a gas filled "quasi-freeflow" sub-domain.
*/
template <class TypeTag>
-class EvaporationAtmosphereProblem
- : public ImplicitPorousMediaProblem<TypeTag>
+class EvaporationAtmosphereProblem: public PorousMediumFlowProblem<TypeTag>
{
- using ParentType = ImplicitPorousMediaProblem<TypeTag>;
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ 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 FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using ParameterCache = typename FluidSystem::ParameterCache;
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
/*!
* \brief The fluid state which is used by the volume variables to
* store the thermodynamic state. This should be chosen
* appropriately for the model ((non-)isothermal, equilibrium, ...).
* This can be done in the problem.
*/
- using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
- using ParameterCache = typename FluidSystem::ParameterCache;
enum { dim = GridView::dimension}; // Grid and world dimension
enum { dimWorld = GridView::dimensionworld};
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
@@ -187,17 +134,17 @@ class EvaporationAtmosphereProblem
enum { p0Idx = Indices::p0Idx};
enum { conti00EqIdx = Indices::conti0EqIdx };
enum { energyEq0Idx = Indices::energyEqIdx};
- enum { numEnergyEqs = Indices::numPrimaryEnergyVars};
enum { wPhaseIdx = FluidSystem::wPhaseIdx};
enum { nPhaseIdx = FluidSystem::nPhaseIdx};
enum { sPhaseIdx = FluidSystem::sPhaseIdx};
enum { wCompIdx = FluidSystem::H2OIdx};
enum { nCompIdx = FluidSystem::N2Idx};
- enum { enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)};
- enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
- enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
- enum { velocityOutput = GET_PROP_VALUE(TypeTag, VtkAddVelocities)};
- enum { velocityAveragingInModel = GET_PROP_VALUE(TypeTag, VelocityAveragingInModel)};
+ enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid)};
+ enum { numEnergyEqSolid = GET_PROP_VALUE(TypeTag, NumEnergyEqSolid)};
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+
+ static constexpr bool enableChemicalNonEquilibrium = GET_PROP_VALUE(TypeTag, EnableChemicalNonEquilibrium);
+ static constexpr bool enableThermalNonEquilibrium = GET_PROP_VALUE(TypeTag, EnableChemicalNonEquilibrium);
// formulations
enum {
@@ -206,21 +153,7 @@ class EvaporationAtmosphereProblem
leastWettingFirst = MpNcPressureFormulation::leastWettingFirst
};
- using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
- using Intersection = typename GridView::Intersection;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
-
- using GlobalPosition = Dune::FieldVector<typename GridView::Grid::ctype, dimWorld>;
-
- using ScalarBuffer = std::vector<Dune::FieldVector<Scalar, 1>>;
- using PhaseBuffer = std::array<ScalarBuffer, numPhases>;
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using VelocityField = Dune::BlockVector<DimVector>;
- using PhaseVelocityField = std::array<VelocityField, numPhases>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
/*!
@@ -229,134 +162,33 @@ public:
* \param timeManager The time manager
* \param gridView The grid view
*/
- EvaporationAtmosphereProblem(TimeManager &timeManager,
- const GridView &gridView)
- : ParentType(timeManager, gridView), gnuplot_(false)
- {
- this->timeManager().startNextEpisode(24.* 3600.);
- }
-
- void episodeEnd()
+ EvaporationAtmosphereProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
- // each day is one episode, result file at the end of each day: comparison
- this->timeManager().startNextEpisode(24.* 3600.);
- }
+ percentOfEquil_ = getParam<Scalar>("BoundaryConditions.percentOfEquil");
+ nTemperature_ = getParam<Scalar>("FluidSystem.nTemperature");
+ nPressure_ = getParam<Scalar>("FluidSystem.nPressure");
+ outputName_ = getParam<std::string>("Problem.Name");
+ TInitial_ = getParam<Scalar>("InitialConditions.TInitial");
+ SwPMInitial_ = getParam<Scalar>("InitialConditions.SwPMInitial");
+ SwFFInitial_ = getParam<Scalar>("InitialConditions.SwFFInitial");
+ pnInitial_ = getParam<Scalar>("InitialConditions.pnInitial");
+ pnInjection_ = getParam<Scalar>("InitialConditions.pnInjection");
+ TInject_ = getParam<Scalar>("BoundaryConditions.TInject");
+ massFluxInjectedPhase_ = getParam<Scalar>("BoundaryConditions.massFluxInjectedPhase");
- void init()
- {
- percentOfEquil_ = GET_RUNTIME_PARAM(TypeTag, Scalar,BoundaryConditions.percentOfEquil);
- nTemperature_ = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.nTemperature);
- nPressure_ = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.nPressure);
- outputName_ = GET_RUNTIME_PARAM(TypeTag, std::string, Constants.outputName);
- nRestart_ = GET_RUNTIME_PARAM(TypeTag, int, Constants.nRestart);
- TInitial_ = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.TInitial);
- SwPMInitial_ = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.SwPMInitial);
- SwFFInitial_ = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.SwFFInitial);
- pnInitial_ = GET_RUNTIME_PARAM(TypeTag, Scalar,InitialConditions.pnInitial);
- pnInjection_ = GET_RUNTIME_PARAM(TypeTag, Scalar,InitialConditions.pnInjection);
- TInject_ = GET_RUNTIME_PARAM(TypeTag, Scalar,BoundaryConditions.TInject);
- massFluxInjectedPhase_ = GET_RUNTIME_PARAM(TypeTag, Scalar,BoundaryConditions.massFluxInjectedPhase);
-
- this->spatialParams().setInputInitialize();
-
- initFluidSystem_();
-
- ParentType::init();
- this->model().initVelocityStuff();
- }
-
- void initFluidSystem_()
- {
// initialize the tables of the fluid system
FluidSystem::init(TInitial_ - 15.0, 453.15, nTemperature_, // T_min, T_max, n_T
0.75*pnInitial_, 2.25*pnInitial_, nPressure_); // p_min, p_max, n_p
}
- /*!
- * \brief User defined output after the time integration
- *
- * Will be called diretly after the time integration.
- */
- void postTimeStep()
- {
- // write two plot Over Lines to text files
- // each output gets it's writer object in order to allow for different header files
- PlotOverLine2D<TypeTag> writerInterface;
- PlotOverLine2D<TypeTag> writerDepth;
-
- // writer points: on the pm--ff Interface to file
- GlobalPosition pointInterfaceOne ;
- pointInterfaceOne[0] = this->bBoxMin()[0] ;
- pointInterfaceOne[1] = (this->bBoxMax()[1]) /2 ;
-
- GlobalPosition pointInterfaceTwo ;
- pointInterfaceTwo[0] = this->bBoxMax()[0] ;
- pointInterfaceTwo[1] = (this->bBoxMax()[1]) /2 ;
-
- writerInterface.write(*this,
- pointInterfaceOne,
- pointInterfaceTwo,
- "plotOverLineInterface");
-
- // writer points: cut through the domain
- GlobalPosition pointDepthOne ;
- pointDepthOne[0] = (this->bBoxMax()[0])/2 ;
- pointDepthOne[1] = this->bBoxMin()[1] ;
-
- GlobalPosition pointDepthTwo ;
- pointDepthTwo[0] = (this->bBoxMax()[0])/2 ;
- pointDepthTwo[1] = (this->bBoxMax()[1])/2 ;
-
- writerDepth.write(*this,
- pointDepthOne,
- pointDepthTwo,
- false, /* do not append data*/
- "plotOverLineIntoDepth");
-
- // Calculate storage terms of the individual phases
- for (int phaseIdx = 0; phaseIdx < numEnergyEqs; ++phaseIdx) {
- PrimaryVariables phaseStorage;
- /* How much conserved quantity is in the system (in the unit of the balance equation)
- * Summed up storage term for the whole system.
- */
- this->model().globalPhaseStorage(phaseStorage, phaseIdx);
-
- if (this->gridView().comm().rank() == 0) {
- std::cout
- <<"Storage in "
- << FluidSystem::phaseName(phaseIdx)
- << "Phase: ["
- << phaseStorage
- << "]"
- << "\n";
- }
- }
-
- // Calculate total storage terms
- PrimaryVariables storage;
- this->model().globalStorage(storage);
-
- // Write mass balance information for rank 0
- if (this->gridView().comm().rank() == 0) {
- std::cout
- <<"Storage total: [" << storage << "]"
- << "\n";
- }
+ void setGridVariables(std::shared_ptr<GridVariables> gridVariables)
+ { gridVariables_ = gridVariables; }
- // use gnuplot for plotting the line data
- gnuplot_.resetPlot();
- gnuplot_.setXlabel("xN2w [-]");
- gnuplot_.setYlabel("y [m]");
- std::ostringstream stream;
- stream << this->timeManager().time();
- gnuplot_.setOption("set label 'at time " + stream.str() + "' at screen 0.15,0.90 left");
- gnuplot_.setDatafileSeparator(' ');
- std::string fileName = outputName_ + "plotOverLineIntoDepth.dat";
- gnuplot_.addFileToPlot(fileName, "u 11:4 w l");
- gnuplot_.plot("plot");
- }
+ const GridVariables& gridVariables() const
+ { return *gridVariables_; }
- /*!
+ /*!
* \name Problem parameters
*/
// \{
@@ -369,31 +201,6 @@ public:
const std::string name() const
{ return outputName_ ; }
- /*!
- * \brief Returns the temperature \f$ K \f$
- */
- Scalar boxTemperature(const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- { return TInitial_; }
-
- /*!
- * \brief Write a restart file?
- * yes of course:
- * - every nRestart_'th steps
- */
- bool shouldWriteRestartFile() const
- {
- return this->timeManager().timeStepIndex() > 0 &&
- (this->timeManager().timeStepIndex() % nRestart_ == 0);
- }
-
- /*!
- * \brief Write a output file?
- */
- bool shouldWriteOutput() const
- { return true; }
-
// \}
/*!
@@ -408,32 +215,18 @@ public:
* \param bTypes The boundarytypes types for the conservation equations
* \param globalPos The global position
*/
- void boundaryTypesAtPos(BoundaryTypes &bTypes,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
+ BoundaryTypes bcTypes;
// Default: Neumann
- bTypes.setAllNeumann();
-
- // To the right: let out what wants out
- if(onRightBoundary_(globalPos) && this->spatialParams().inFF_(globalPos) )
- {
- bTypes.setAllOutflow();
- }
+ bcTypes.setAllNeumann();
// Put a dirichlet somewhere: we need this for convergence
- if(onUpperBoundary_(globalPos) )
+ if(onRightBoundary_(globalPos) && globalPos[1] > this->spatialParams().heightPM() + eps_)
{
- bTypes.setAllDirichlet();
- }
-
- // In the porous part the *temperature* is fixed on the boundary.
- // Mass however, is not allowed to pass (default neumann=0)
- if (this->spatialParams().inPM_(globalPos)
- && (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos) || onLowerBoundary_(globalPos)))
- {
- for (int energyEqIdx=0; energyEqIdx< numEnergyEqs; ++energyEqIdx)
- bTypes.setDirichlet(energyEq0Idx+energyEqIdx);
+ bcTypes.setAllDirichlet();
}
+ return bcTypes;
}
/*!
@@ -445,10 +238,9 @@ public:
* \param globalPos The global position
*
*/
- void dirichletAtPos(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(priVars, globalPos);
+ return initial_(globalPos);
}
/*!
@@ -468,21 +260,13 @@ public:
* The \a values store the mass flux of each phase normal to the boundary.
* Negative values indicate an inflow.
*/
- void neumann(PrimaryVariables & priVars,
- const Element & element,
- const FVElementGeometry & fvGeometry,
- const Intersection & intersection,
- const unsigned int scvIdx,
- const unsigned int boundaryFaceIdx) const
+ PrimaryVariables neumann(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
{
- // this is the coordinate of the vertex
- // distinction via vertex works better in this case, because this is also how the
- // permeabilities & co are defined. This way there is only injection in the free flow and
- // not also in the last porous medium node.
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global;
-
- priVars = 0.0;
-
+ PrimaryVariables priVars(0.0);
+ const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
const Scalar massFluxInjectedPhase = massFluxInjectedPhase_ ;
ParameterCache dummyCache;
@@ -500,7 +284,7 @@ public:
FluidSystem::calculateEquilibriumMoleFractions(fluidState, dummyCache) ;
// Now let's make the air phase less than fully saturated with water
- fluidState.setMoleFraction(nPhaseIdx, wCompIdx, fluidState.moleFraction(nPhaseIdx, wCompIdx) * percentOfEquil_ ) ;
+ fluidState.setMoleFraction(nPhaseIdx, wCompIdx, fluidState.moleFraction(nPhaseIdx, wCompIdx)*percentOfEquil_ ) ;
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1.-fluidState.moleFraction(nPhaseIdx, wCompIdx) ) ;
// compute density of injection phase
@@ -509,19 +293,12 @@ public:
nPhaseIdx);
fluidState.setDensity(nPhaseIdx, density);
- if(numEnergyEquations){
- for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
+ for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
+ {
const Scalar h = FluidSystem::enthalpy(fluidState,
dummyCache,
phaseIdx);
fluidState.setEnthalpy(phaseIdx, h);
- }
- }
- else{
- const Scalar h = FluidSystem::enthalpy(fluidState,
- dummyCache,
- nPhaseIdx);
- fluidState.setEnthalpy(nPhaseIdx, h);
}
const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(nPhaseIdx);
@@ -534,19 +311,12 @@ public:
priVars[conti00EqIdx + nPhaseIdx * numComponents + nCompIdx]
= -molarFlux * fluidState.moleFraction(nPhaseIdx, nCompIdx);
// energy equations are specified mass specifically
- if(numEnergyEquations){
- priVars[energyEq0Idx + nPhaseIdx] = - massFluxInjectedPhase
- * fluidState.enthalpy(nPhaseIdx) ;
- }
- else if(enableEnergy)
- priVars[energyEq0Idx] = - massFluxInjectedPhase
- * fluidState.enthalpy(nPhaseIdx) ;
+ priVars[energyEq0Idx + nPhaseIdx] = - massFluxInjectedPhase
+ * fluidState.enthalpy(nPhaseIdx) ;
}
+ return priVars;
}
- // \}
-
-
/*!
* \name Volume terms
*/
@@ -559,10 +329,9 @@ public:
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
- void initialAtPos(PrimaryVariables &priVars,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- initial_(priVars, globalPos);
+ return initial_(globalPos);
}
/*!
@@ -577,22 +346,21 @@ public:
*
* Positive values mean that mass is created, negative ones mean that it vanishes.
*/
- void source(PrimaryVariables & priVars,
- const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ PrimaryVariables source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
{
- priVars = Scalar(0.0);
- }
-
- // \}
+ PrimaryVariables priVars(0.0);
+ return priVars;
+ }
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);
const Scalar T = TInitial_;
Scalar S[numPhases];
@@ -615,17 +383,13 @@ private:
FluidState equilibriumFluidState;
//set saturation to inital values, this needs to be done in order for the fluidState to tell me pc
- for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx) {
+ for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx)
+ {
equilibriumFluidState.setSaturation(phaseIdx, S[phaseIdx]);
-
- if(numEnergyEquations){
- equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
- }
- else
- equilibriumFluidState.setTemperature(TInject_ );
+ equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
}
- const MaterialLawParams &materialParams =
+ const auto &materialParams =
this->spatialParams().materialLawParamsAtPos(globalPos);
Scalar capPress[numPhases];
@@ -658,9 +422,7 @@ private:
}
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << pressureFormulation << " is invalid.");
- // temperature
- if (enableEnergy || numEnergyEquations)
- for (int energyEqIdx=0; energyEqIdx< numEnergyEqs; ++energyEqIdx)
+ for (int energyEqIdx=0; energyEqIdx< numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx)
priVars[energyEq0Idx + energyEqIdx] = T;
for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
@@ -673,12 +435,12 @@ private:
FluidState dryFluidState(equilibriumFluidState);
// Now let's make the air phase less than fully saturated with vapor
dryFluidState.setMoleFraction(nPhaseIdx, wCompIdx, dryFluidState.moleFraction(nPhaseIdx, wCompIdx) * percentOfEquil_ ) ;
- dryFluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1.-dryFluidState.moleFraction(nPhaseIdx, wCompIdx) ) ;
+ dryFluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1.0-dryFluidState.moleFraction(nPhaseIdx, wCompIdx) ) ;
/* Difference between kinetic and MPNC:
* number of component related primVar and how they are calculated (mole fraction, fugacities, resp.)
*/
- if(enableKinetic){
+ if(enableChemicalNonEquilibrium){
for(int phaseIdx=0; phaseIdx < numPhases; ++ phaseIdx)
{
for(int compIdx=0; compIdx <numComponents; ++compIdx){
@@ -695,7 +457,8 @@ private:
}
}
}
- else{ //enableKinetic
+ else
+ {
// in the case I am using the "standard" mpnc model, the variables to be set are the "fugacities"
const Scalar fugH2O = FluidSystem::H2O::vaporPressure(T) ;
const Scalar fugN2 = p[nPhaseIdx] - fugH2O ;
@@ -716,31 +479,26 @@ private:
for (int i = 0; i < numComponents; ++i)
priVars[conti00EqIdx + i] = xl[i]*beta[i]; // this should be really fug0Idx but the compiler only knows one or the other
}
+ return priVars;
}
/*!
* \brief Give back whether the tested position (input) is a specific region (left) in the domain
*/
bool onLeftBoundary_(const GlobalPosition & globalPos) const
- { return globalPos[0] < this->bBoxMin()[0] + eps_; }
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
/*!
* \brief Give back whether the tested position (input) is a specific region (right) in the domain
*/
bool onRightBoundary_(const GlobalPosition & globalPos) const
- { return globalPos[0] > this->bBoxMax()[0] - eps_; }
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
/*!
* \brief Give back whether the tested position (input) is a specific region (down, (gravityDir)) in the domain
*/
bool onLowerBoundary_(const GlobalPosition & globalPos) const
- { return globalPos[dimWorld-1] < this->bBoxMin()[dimWorld-1] + eps_; }
-
- /*!
- * \brief Give back whether the tested position (input) is a specific region (up, (gravityDir)) in the domain
- */
- bool onUpperBoundary_(const GlobalPosition & globalPos) const
- { return globalPos[dimWorld-1] > this->bBoxMax()[dimWorld-1] - eps_; }
+ { return globalPos[dimWorld-1] < this->fvGridGeometry().bBoxMin()[dimWorld-1] + eps_; }
private:
static constexpr Scalar eps_ = 1e-6;
@@ -748,7 +506,6 @@ private:
int nTemperature_;
int nPressure_;
std::string outputName_;
- int nRestart_ ;
Scalar heatIntoSolid_;
Scalar TInitial_ ;
Scalar SwPMInitial_ ;
@@ -758,15 +515,12 @@ private:
Scalar pnInjection_;
Dune::ParameterTree inputParameters_;
Scalar x_[numPhases][numComponents] ;
- GnuplotInterface<Scalar> gnuplot_;
Scalar TInject_;
Scalar massFluxInjectedPhase_ ;
- PhaseVelocityField volumeDarcyVelocity_;
- PhaseBuffer volumeDarcyMagVelocity_ ;
- ScalarBuffer boxSurface_;
+ std::shared_ptr<GridVariables> gridVariables_;
public:
diff --git a/test/porousmediumflow/mpnc/implicit/evaporationatmospherespatialparams.hh b/test/porousmediumflow/mpnc/implicit/evaporationatmospherespatialparams.hh
index 1e7a5784b9..a74ca409f9 100644
--- a/test/porousmediumflow/mpnc/implicit/evaporationatmospherespatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/evaporationatmospherespatialparams.hh
@@ -43,8 +43,6 @@
#include <dune/common/parametertreeparser.hh>
-#include <dumux/porousmediumflow/mpnc/implicit/modelkinetic.hh>
-
namespace Dumux
{
@@ -64,7 +62,6 @@ SET_TYPE_PROP(EvaporationAtmosphereSpatialParams, SpatialParams, EvaporationAtmo
SET_PROP(EvaporationAtmosphereSpatialParams, MaterialLaw)
{
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-
enum {wPhaseIdx = FluidSystem::wPhaseIdx};
private:
@@ -92,16 +89,15 @@ private:
#endif
#if RegVGofT
- using EffectiveLaw = RegularizedVanGenuchtenOfTemperature<Scalar>;
+ using EffectiveLaw = RegularizedVanGenuchtenOfTemperature<Scalar>;
#endif
#if VG
- using EffectiveLaw = VanGenuchten<Scalar>;
+ using EffectiveLaw = VanGenuchten<Scalar> ;
#endif
using TwoPMaterialLaw = EffToAbsLaw<EffectiveLaw>;
public:
-// using type = TwoPOfTAdapter<wPhaseIdx, TwoPMaterialLaw>; // adapter for incorporating temperature effects on pc-S
using type = TwoPAdapter<wPhaseIdx, TwoPMaterialLaw>;
};
@@ -111,9 +107,9 @@ private:
SET_PROP(EvaporationAtmosphereSpatialParams, AwnSurface)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
using EffectiveIALaw = AwnSurfacePcMaxFct<Scalar>;
- // using EffectiveIALaw = AwnSurfacePolynomial2ndOrder<Scalar>;
public:
using type = EffToAbsLawIA<EffectiveIALaw, MaterialLawParams>;
};
@@ -123,7 +119,8 @@ public:
SET_PROP(EvaporationAtmosphereSpatialParams, AwsSurface)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
using EffectiveIALaw = AwnSurfacePolynomial2ndOrder<Scalar>;
public:
using type = EffToAbsLawIA<EffectiveIALaw, MaterialLawParams>;
@@ -133,7 +130,8 @@ public:
SET_PROP(EvaporationAtmosphereSpatialParams, AnsSurface)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
using EffectiveIALaw = AwnSurfaceExpSwPcTo3<Scalar>;
public:
using type = EffToAbsLawIA<EffectiveIALaw, MaterialLawParams>;
@@ -148,8 +146,15 @@ template<class TypeTag>
class EvaporationAtmosphereSpatialParams : public FVSpatialParams<TypeTag>
{
using ParentType = FVSpatialParams<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
enum {dim=GridView::dimension };
@@ -157,168 +162,128 @@ class EvaporationAtmosphereSpatialParams : public FVSpatialParams<TypeTag>
enum {wPhaseIdx = FluidSystem::wPhaseIdx};
enum {nPhaseIdx = FluidSystem::nPhaseIdx};
enum {sPhaseIdx = FluidSystem::sPhaseIdx};
- enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
+ enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid)};
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
- enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)};
-
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using LocalPosition = Dune::FieldVector<Scalar, dim>;
+ enum { enableThermalNonEquilibrium = GET_PROP_VALUE(TypeTag,EnableThermalNonEquilibrium)};
- using DimVector = Dune::FieldVector<Scalar, dim>;
+ using DimVector = Dune::FieldVector<Scalar,dim>;
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
public:
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename MaterialLaw::Params;
using AwnSurface = typename GET_PROP_TYPE(TypeTag, AwnSurface);
- using AwnSurfaceParams = typename AwnSurface::Params;
-
+ using AwnSurfaceParams = typename AwnSurface::Params;
using AwsSurface = typename GET_PROP_TYPE(TypeTag, AwsSurface);
- using AwsSurfaceParams = typename AwsSurface::Params;
-
+ using AwsSurfaceParams = typename AwsSurface::Params;
using AnsSurface = typename GET_PROP_TYPE(TypeTag, AnsSurface);
- using AnsSurfaceParams = typename AnsSurface::Params;
-
+ using AnsSurfaceParams = typename AnsSurface::Params;
+ using PermeabilityType = Scalar;
- EvaporationAtmosphereSpatialParams(const GridView &gridView)
- : ParentType(gridView)
+ EvaporationAtmosphereSpatialParams(const Problem &problem)
+ : ParentType(problem)
{
- }
-
- ~EvaporationAtmosphereSpatialParams()
- {}
-
- // load parameters from input file and initialize parameter values
- void setInputInitialize()
- {
- heightPM_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.InterfacePosY);
- heightDomain_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, Grid.UpperRight)[1];
- // BEWARE! First the input values have to be set, than the material parameters can be set
-
- // this is the parameter value from file part
- porosityPM_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PorousMedium.porosity);
- intrinsicPermeabilityPM_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PorousMedium.permeability);
-
- porosityFF_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.FreeFlow.porosity);
- intrinsicPermeabilityFF_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.FreeFlow.permeability);
-
- solidDensity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.density);
- solidThermalConductivity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.thermalConductivity);
- solidHeatCapacity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.heatCapacity);
-
- aWettingNonWettingA1_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aWettingNonWettingA1);
- aWettingNonWettingA2_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aWettingNonWettingA2);
- aWettingNonWettingA3_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aWettingNonWettingA3);
-
- aNonWettingSolidA1_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aNonWettingSolidA1);
- aNonWettingSolidA2_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aNonWettingSolidA2);
- aNonWettingSolidA3_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.aNonWettingSolidA3);
-
- BCPd_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.BCPd);
- BClambda_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.BClambda);
- Swr_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.Swr);
- Snr_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.Snr);
-
- characteristicLengthFF_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.FreeFlow.meanPoreSize);
-
- characteristicLengthPM_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PorousMedium.meanPoreSize);
- factorEnergyTransfer_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PorousMedium.factorEnergyTransfer);
- factorMassTransfer_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PorousMedium.factorMassTransfer);
-
- // residual saturations
- materialParamsFF_.setSwr(0.0);
- materialParamsFF_.setSnr(0.00);
-
- materialParamsPM_.setSwr(Swr_);
- materialParamsPM_.setSnr(Snr_);
-
- // pc / kr parameters
- materialParamsPM_.setLambda(BClambda_);
- materialParamsPM_.setPe(BCPd_);
-
- // for making pc == 0 in the FF
- materialParamsFF_.setLambda(42.);
- materialParamsFF_.setPe(0.);
-
- {//scope it
- // capillary pressure parameters
- FluidState fluidState ;
- Scalar S[numPhases] ;
- const MaterialLawParams &materialParams = materialParamsPM_;
- Scalar capPress[numPhases];
- //set saturation to inital values, this needs to be done in order for the fluidState to tell me pc
- for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx) {
- // set saturation to zero for getting pcmax
- S[phaseIdx] = 0. ;
- Scalar TInitial = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.TInitial);
-
- fluidState.setSaturation(phaseIdx, S[phaseIdx]);
- if(enableEnergy){
- if (numEnergyEquations>1)
- fluidState.setTemperature(phaseIdx,TInitial);
- else
- fluidState.setTemperature(TInitial);
- }
- }
-
- //obtain pc according to saturation
- MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
- using std::abs;
- pcMax_ = abs(capPress[0]);
-
- // set pressures from capillary pressures
- aWettingNonWettingSurfaceParams_.setPcMax(pcMax_);
+ heightPM_ = getParam<std::vector<Scalar>>("Grid.Positions1")[1];
+ heightDomain_ = getParam<std::vector<Scalar>>("Grid.Positions1")[2];
+
+ porosityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.porosity");
+ intrinsicPermeabilityPM_ = getParam<Scalar>("SpatialParams.PorousMedium.permeability");
+
+ porosityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.porosity");
+ intrinsicPermeabilityFF_ = getParam<Scalar>("SpatialParams.FreeFlow.permeability");
+
+ solidDensity_ = getParam<Scalar>("SpatialParams.soil.density");
+ solidThermalConductivity_ = getParam<Scalar>("SpatialParams.soil.thermalConductivity");
+ solidHeatCapacity_ = getParam<Scalar>("SpatialParams.soil.heatCapacity");
+
+ aWettingNonWettingA1_ = getParam<Scalar>("SpatialParams.soil.aWettingNonWettingA1");
+ aWettingNonWettingA2_ = getParam<Scalar>("SpatialParams.soil.aWettingNonWettingA2");
+ aWettingNonWettingA3_ = getParam<Scalar>("SpatialParams.soil.aWettingNonWettingA3");
+
+ aNonWettingSolidA1_ = getParam<Scalar>("SpatialParams.soil.aNonWettingSolidA1");
+ aNonWettingSolidA2_ = getParam<Scalar>("SpatialParams.soil.aNonWettingSolidA2");
+ aNonWettingSolidA3_ = getParam<Scalar>("SpatialParams.soil.aNonWettingSolidA3");
+
+ BCPd_ = getParam<Scalar>("SpatialParams.soil.BCPd");
+ BClambda_ = getParam<Scalar>("SpatialParams.soil.BClambda");
+ Swr_ = getParam<Scalar>("SpatialParams.soil.Swr");
+ Snr_ = getParam<Scalar>("SpatialParams.soil.Snr");
+
+ characteristicLengthFF_ = getParam<Scalar>("SpatialParams.FreeFlow.meanPoreSize");
+ characteristicLengthPM_ = getParam<Scalar>("SpatialParams.PorousMedium.meanPoreSize");
+
+ factorEnergyTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorEnergyTransfer");
+ factorMassTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorMassTransfer");
+
+ // residual saturations
+ materialParamsFF_.setSwr(0.0);
+ materialParamsFF_.setSnr(0.00);
+
+ materialParamsPM_.setSwr(Swr_);
+ materialParamsPM_.setSnr(Snr_);
+
+ // pc / kr parameters
+ materialParamsPM_.setLambda(BClambda_);
+ materialParamsPM_.setPe(BCPd_);
+
+ // for making pc == 0 in the FF
+ materialParamsFF_.setLambda(42.);
+ materialParamsFF_.setPe(0.);
+
+ {//scope it
+ // capillary pressure parameters
+ FluidState fluidState ;
+ Scalar S[numPhases] ;
+ const auto &materialParams = materialParamsPM_;
+ Scalar capPress[numPhases];
+ //set saturation to inital values, this needs to be done in order for the fluidState to tell me pc
+ for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx) {
+ // set saturation to zero for getting pcmax
+ S[phaseIdx] = 0. ;
+ Scalar TInitial = getParam<Scalar>("InitialConditions.TInitial");
+ fluidState.setSaturation(phaseIdx, S[phaseIdx]);
+ fluidState.setTemperature(phaseIdx,TInitial);
}
- // wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
- aWettingNonWettingSurfaceParams_.setA1(aWettingNonWettingA1_);
- aWettingNonWettingSurfaceParams_.setA2(aWettingNonWettingA2_);
- aWettingNonWettingSurfaceParams_.setA3(aWettingNonWettingA3_);
-
- // non-wetting-solid
- aNonWettingSolidSurfaceParams_.setA1(aNonWettingSolidA1_);
- aNonWettingSolidSurfaceParams_.setA2(aNonWettingSolidA2_);
- aNonWettingSolidSurfaceParams_.setA3(aNonWettingSolidA3_);
+ //obtain pc according to saturation
+ MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
+ using std::abs;
+ pcMax_ = abs(capPress[0]);
- // dummys for free flow: no interface where there is only one phase
- aWettingNonWettingSurfaceParamsFreeFlow_.setA1(0.);
- aWettingNonWettingSurfaceParamsFreeFlow_.setA2(0.);
- aWettingNonWettingSurfaceParamsFreeFlow_.setA3(0.);
- aWettingNonWettingSurfaceParamsFreeFlow_.setPcMax(42.); // not needed because it is anyways zero; silencing valgrind
+ // set pressures from capillary pressures
+ aWettingNonWettingSurfaceParams_.setPcMax(pcMax_);
+ }
- // dummys for free flow: no interface where there is only one phase
- aNonWettingSolidSurfaceParamsFreeFlow_.setA1(0.);
- aNonWettingSolidSurfaceParamsFreeFlow_.setA2(0.);
- aNonWettingSolidSurfaceParamsFreeFlow_.setA3(0.);
+ // wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
+ aWettingNonWettingSurfaceParams_.setA1(aWettingNonWettingA1_);
+ aWettingNonWettingSurfaceParams_.setA2(aWettingNonWettingA2_);
+ aWettingNonWettingSurfaceParams_.setA3(aWettingNonWettingA3_);
+
+ // non-wetting-solid
+ aNonWettingSolidSurfaceParams_.setA1(aNonWettingSolidA1_);
+ aNonWettingSolidSurfaceParams_.setA2(aNonWettingSolidA2_);
+ aNonWettingSolidSurfaceParams_.setA3(aNonWettingSolidA3_);
+
+ // dummys for free flow: no interface where there is only one phase
+ aWettingNonWettingSurfaceParamsFreeFlow_.setA1(0.);
+ aWettingNonWettingSurfaceParamsFreeFlow_.setA2(0.);
+ aWettingNonWettingSurfaceParamsFreeFlow_.setA3(0.);
+ aWettingNonWettingSurfaceParamsFreeFlow_.setPcMax(42.); // not needed because it is anyways zero; silencing valgrind
+
+ // dummys for free flow: no interface where there is only one phase
+ aNonWettingSolidSurfaceParamsFreeFlow_.setA1(0.);
+ aNonWettingSolidSurfaceParamsFreeFlow_.setA2(0.);
+ aNonWettingSolidSurfaceParamsFreeFlow_.setA3(0.);
+ }
- // dummy for aws not necessary: aws = ans - as, ans=0, as=ans(0,pcmax)=0 -> aws=0
- }
+ ~EvaporationAtmosphereSpatialParams()
+ {}
- /*! \brief Update the spatial parameters with the flow solution
- * after a timestep.
- */
- void update(const SolutionVector & globalSolutionFn)
- { }
- /*! \brief Returns the intrinsic permeability
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub control volume.
- * \param element The current finite element
- * \param fvGeometry The current finite volume geometry of the element
- * \param scvIdx The index sub-control volume */
- const Scalar intrinsicPermeability(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
{
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
+ const auto & globalPos = scv.dofPosition();
if (inFF_(globalPos) )
return intrinsicPermeabilityFF_ ;
else if (inPM_(globalPos))
@@ -334,13 +299,11 @@ public:
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume */
- const Scalar porosity(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ Scalar porosity(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
{
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
-// const LocalPosition & localPos = fvGeometry.subContVol[scvIdx].localCenter ;
-// const GlobalPosition & globalPos = element.geometry().global(localPos) ;
+ const auto& globalPos = scv.dofPosition();
if (inFF_(globalPos) )
return porosityFF_ ;
@@ -350,28 +313,21 @@ public:
DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
- /*!
- * \brief Return a reference to the material parameters of the material law.
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub control volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume */
- const MaterialLawParams & materialLawParams(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
{
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
- const MaterialLawParams & materialParams = materialLawParamsAtPos(globalPos) ;
- return materialParams ;
+ const auto& globalPos = scv.dofPosition();
+ if (inFF_(globalPos) )
+ return materialParamsFF_ ;
+ else if (inPM_(globalPos))
+ return materialParamsPM_ ;
+ else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
/*!
* \brief Return a reference to the material parameters of the material law.
- *
- * \param globalPos The position in global coordinates.
- */
+ * \param globalPos The position in global coordinates. */
const MaterialLawParams & materialLawParamsAtPos(const GlobalPosition & globalPos) const
{
if (inFF_(globalPos) )
@@ -389,11 +345,11 @@ public:
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume */
- const AwnSurfaceParams & aWettingNonWettingSurfaceParams(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ const AwnSurfaceParams & aWettingNonWettingSurfaceParams(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
{
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
+ const auto& globalPos = scv.dofPosition();
if (inFF_(globalPos) )
return aWettingNonWettingSurfaceParamsFreeFlow_ ;
else if (inPM_(globalPos))
@@ -410,10 +366,10 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume */
const AnsSurfaceParams & aNonWettingSolidSurfaceParams(const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
{
- const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
+ const auto& globalPos = scv.dofPosition();
if (inFF_(globalPos) )
return aNonWettingSolidSurfaceParamsFreeFlow_ ;
else if (inPM_(globalPos))
@@ -433,11 +389,12 @@ public:
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume */
- const Scalar pcMax(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ const Scalar pcMax(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
{ return aWettingNonWettingSurfaceParams_.pcMax() ; }
+
/*!\brief Return the characteristic length for the mass transfer.
*
* The position is determined based on the coordinate of
@@ -446,10 +403,14 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub control volume */
const Scalar characteristicLength(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
- { const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
- return characteristicLengthAtPos(globalPos); }
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
+
+ {
+ const auto& globalPos = scv.dofPosition();
+ return characteristicLengthAtPos(globalPos);
+ }
+
/*!\brief Return the characteristic length for the mass transfer.
* \param globalPos The position in global coordinates.*/
const Scalar characteristicLengthAtPos(const GlobalPosition & globalPos) const
@@ -468,14 +429,17 @@ public:
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub control volume */
- const Scalar factorEnergyTransfer(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
- { const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
- return factorEnergyTransferAtPos(globalPos); }
+ const Scalar factorEnergyTransfer(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
+ {
+ const auto& globalPos = scv.dofPosition();
+ return factorEnergyTransferAtPos(globalPos);
+ }
+
/*!\brief Return the pre factor the the energy transfer
* \param globalPos The position in global coordinates.*/
- const Scalar factorEnergyTransferAtPos(const GlobalPosition & globalPos) const
+ const Scalar factorEnergyTransferAtPos(const GlobalPosition & globalPos) const
{
if (inFF_(globalPos) )
return factorEnergyTransfer_ ;
@@ -484,21 +448,25 @@ public:
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
- /*!\brief Return the pre factor for the mass transfer
+ /*!\brief Return the pre factor the the mass transfer
*
* The position is determined based on the coordinate of
* the vertex belonging to the considered sub controle volume.
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub control volume */
- const Scalar factorMassTransfer(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
- { const GlobalPosition & globalPos = fvGeometry.subContVol[scvIdx].global ;
- return factorMassTransferAtPos(globalPos); }
+ const Scalar factorMassTransfer(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
+ {
+ const auto& globalPos = scv.dofPosition();
+ return factorMassTransferAtPos(globalPos);
+
+ }
+
/*!\brief Return the pre factor the the mass transfer
* \param globalPos The position in global coordinates.*/
- const Scalar factorMassTransferAtPos(const GlobalPosition & globalPos) const
+ const Scalar factorMassTransferAtPos(const GlobalPosition & globalPos) const
{
if (inFF_(globalPos) )
return factorMassTransfer_ ;
@@ -508,32 +476,34 @@ public:
}
- /*! \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
- * \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 */
- const Scalar solidHeatCapacity(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
+ /*!
+ * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param globalPos The global position
+ */
+ Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{ return solidHeatCapacity_ ;} // specific heat capacity of solid [J / (kg K)]
- /*!\brief Returns the density \f$[kg/m^3]\f$ of the rock matrix.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume */
- const Scalar solidDensity(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx) const
- { return solidDensity_ ;} // density of solid [kg/m^3]
+ /*!
+ * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param globalPos The global position
+ */
+ Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
+ {return solidDensity_ ;} // density of solid [kg/m^3]
- /*!\brief Returns the thermal conductivity \f$[W/(m K)]\f$ of the rock matrix.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume */
- const Scalar solidThermalConductivity(const Element & element,
- const FVElementGeometry & fvGeometry,
- const unsigned int scvIdx)const
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param globalPos The global position
+ */
+ Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{ return solidThermalConductivity_ ;} // conductivity of solid [W / (m K ) ]
/*!\brief Give back whether the tested position (input) is a specific region (porous medium part) in the domain
@@ -547,7 +517,7 @@ public:
* -> be careful with neumannAtPos
*/
bool inPM_(const GlobalPosition & globalPos) const
- { return ( (globalPos[dimWorld-1] > 0. - eps_) and (globalPos[dimWorld-1] < (heightPM_ + eps_) ) ); }
+ { return ( (globalPos[dimWorld-1] > 0. - eps_) and (globalPos[dimWorld-1] < (heightPM_ + eps_) ) ); }
/*!
* \brief Give back whether the tested position (input) is a specific region (above PM / "free flow") in the domain
@@ -561,21 +531,7 @@ public:
* -> be careful with neumannAtPos
*/
bool inFF_(const GlobalPosition & globalPos) const
- { return ( (globalPos[dimWorld-1] < heightDomain_ + eps_) and (globalPos[dimWorld-1] > (heightPM_ + eps_) ) ); }
-
- /*!
- * \brief Give back whether the tested position (input) is a specific region (above PM / "free flow") in the domain
- *
- * This setting ensures, that the boundary between the two domains has porous medium properties.
- * This is desirable, because I want to observe the boundary of the porous domain.
- * However, the position has to be the coordinate of the vertex and not the integration point
- * of the boundary flux. If the position is the ip of the neumann flux this leads to a situation
- * where the vertex belongs to porous medium and there is nonetheless injection over the boundary.
- * That does not work.
- * -> be careful with neumannAtPos
- */
- bool inInjection_(const GlobalPosition & globalPos) const
- { return ( (globalPos[dimWorld-1] < heightDomain_ - 0.25*heightDomain_ + eps_) and (globalPos[dimWorld-1] > (heightPM_ + eps_) ) ); }
+ { return ( (globalPos[dimWorld-1] < heightDomain_ + eps_) and (globalPos[dimWorld-1] > (heightPM_ + eps_) ) ); }
/*! \brief access function for the depth / height of the porous medium */
const Scalar heightPM() const
@@ -627,6 +583,7 @@ private:
Scalar BClambda_ ;
Scalar Swr_ ;
Scalar Snr_ ;
+ std::vector<Scalar> gridVector_;
};
}
diff --git a/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.cc b/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.cc
index 11a6ecd5e0..9478439173 100644
--- a/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.cc
+++ b/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.cc
@@ -16,43 +16,238 @@
* 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 kinetic modules of the mpnc box model.
+ * \brief Test for the three-phase box model
*/
#include <config.h>
-
-#include <dumux/common/start.hh>
#include "evaporationatmosphereproblem.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/properties.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/porousmediumflow/nonequilibrium/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
/*!
- * \brief Print a usage string for simulations.
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
*
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
*/
-void printUsage(const char *progName, const std::string &errorMsg)
+void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\nAn uncomplete list of mandatory options for this program is:\n"
- "[Grid]\n"
- "InterfacePosY Vertical position of the interface [m]\n"
- "\n";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\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)
+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);
+
+ // 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 GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(leafGridView.size(GridView::dimension));
+ 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);
+ problem->setGridVariables(gridVariables);
+
+ // 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.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::AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonController = Dumux::NonEquilibriumNewtonController<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.");
+ }
+
+ // 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 (...)
{
- using ProblemTypeTag = TTAG(EvaporationAtmosphereProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, printUsage);
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
}
diff --git a/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.input b/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.input
index edde273ddf..36b9448b59 100644
--- a/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.input
+++ b/test/porousmediumflow/mpnc/implicit/test_boxmpnckinetic.input
@@ -1,24 +1,18 @@
-[TimeManager]
-DtInitial = 1.5 # [s]
+[TimeLoop]
+DtInitial = 0.05 # [s]
TEnd = 10 # [s]
-MaxTimeStepSize = 2e20 # maximum allowed timestep size
[Grid]
Cells0 = 14
Cells1 = 15 15
Grading0 = 1.0
-Grading1 = 1.2 0.833333
+Grading1 = -0.833333 0.833333
Positions0 = 0.0 1.0
Positions1 = 0.0 0.25 0.5
-Cells = 14 30
-UpperRight = 1.0 0.5
-InterfacePosY = 0.25
-RefineTop = false
-
[InitialConditions]
-SwFFInitial = 1e-5 # - 0.001#
-SwPMInitial = 0.999 # -
+SwFFInitial = 1e-4 # - 0.001#
+SwPMInitial = 0.8 # -
pnInitial = 1e5 # 6.8e6 # 1e5 # 101475 # Pa
TInitial = 293
pnInjection = 100003
@@ -68,10 +62,15 @@ aNonWettingSolidA3 = 1.063e-09 #
heatIntoSolid = 0 #
[Constants]
-outputName = evaporationatmosphere
nRestart = 100 # after so many timesteps a restart file should be written
+[Problem]
+Name = evaporationatmosphere
+
[FluidSystem]
nTemperature = 100 # for the tabulation of water: so many interpolation points
nPressure = 100 # for the tabulation of water: so many interpolation points
hammer = 1e4
+
+[Vtk]
+AddVelocity = 1 # enable velocity output
--
GitLab