From 374053d5a0c3a988b49f89e5f73cddda8bfd9e47 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 21 Dec 2017 15:02:03 +0100
Subject: [PATCH] [mpnc] Remove forchheimer files temporiliy. Reintroduce
properly if needed
---
.../mpnc/implicit/forchheimer1pproblem.hh | 450 ------------------
.../mpnc/implicit/forchheimer2pproblem.hh | 444 -----------------
.../mpnc/implicit/forchheimerspatialparams.hh | 202 --------
.../mpnc/implicit/test_forchheimer1p.cc | 59 ---
.../mpnc/implicit/test_forchheimer1p.input | 13 -
.../mpnc/implicit/test_forchheimer2p.cc | 59 ---
.../mpnc/implicit/test_forchheimer2p.input | 13 -
7 files changed, 1240 deletions(-)
delete mode 100644 test/porousmediumflow/mpnc/implicit/forchheimer1pproblem.hh
delete mode 100644 test/porousmediumflow/mpnc/implicit/forchheimer2pproblem.hh
delete mode 100644 test/porousmediumflow/mpnc/implicit/forchheimerspatialparams.hh
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_forchheimer1p.cc
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_forchheimer1p.input
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_forchheimer2p.cc
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_forchheimer2p.input
diff --git a/test/porousmediumflow/mpnc/implicit/forchheimer1pproblem.hh b/test/porousmediumflow/mpnc/implicit/forchheimer1pproblem.hh
deleted file mode 100644
index e263c8ccc6..0000000000
--- a/test/porousmediumflow/mpnc/implicit/forchheimer1pproblem.hh
+++ /dev/null
@@ -1,450 +0,0 @@
-// -*- 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 Problem where liquid water is injected by means of a
- * dirchlet condition on the left of the domain.
- * Velocity according to Forchheimer.
- */
-#ifndef DUMUX_FORCHEIMER1P_HH
-#define DUMUX_FORCHEIMER1P_HH
-
-#include <dune/common/parametertreeparser.hh>
-
-#include <dumux/porousmediumflow/mpnc/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
-
-#include <dumux/material/fluidsystems/h2oair.hh>
-#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
-#include <dumux/material/fluidstates/compositional.hh>
-
-#include "forchheimerspatialparams.hh"
-
-namespace Dumux
-{
-
-template <class TypeTag>
-class Forchheimer1pProblem;
-
-namespace Properties
-{
-NEW_TYPE_TAG(Forchheimer1pProblem, INHERITS_FROM(BoxMPNC, ForchheimerSpatialParams));
-
-// Set the grid type
-SET_TYPE_PROP(Forchheimer1pProblem, Grid, Dune::YaspGrid<1>);
-
-// Set the problem property
-SET_TYPE_PROP(Forchheimer1pProblem,
- Problem,
- Forchheimer1pProblem<TypeTag>);
-
-
-// Set fluid configuration
-SET_PROP(Forchheimer1pProblem, FluidSystem)
-{ private:
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
- using type = FluidSystems::H2OAir<Scalar, SimpleH2O<Scalar>, /*useComplexRelations=*/false>;
-};
-
-// Enable molecular diffusion of the components?
-SET_BOOL_PROP(Forchheimer1pProblem, EnableDiffusion, false);
-
-// Enable gravity
-SET_BOOL_PROP(Forchheimer1pProblem, ProblemEnableGravity, false);
-
-// decide which type to use for floating values (double / quad)
-SET_TYPE_PROP(Forchheimer1pProblem, Scalar, double);
-
-// decide which to use for velocity calculation: Darcy / Forchheimer
-SET_TYPE_PROP(Forchheimer1pProblem, BaseFluxVariables, ImplicitForchheimerFluxVariables<TypeTag>);
-
-SET_BOOL_PROP(Forchheimer1pProblem, VtkAddVelocities, true);
-}
-
-
-/*!
- * \ingroup MPNCModel
- * \ingroup ImplicitTestProblems
- * \brief Problem where liquid water is injected by means of a
- * dirichlet condition on the left of the domain.
- * Velocity according to Forchheimer.
- *
- * The setup is for testing of the one phase Forchheimer relation.
- * The whole domain is filled with water at 1e5 Pa. On the left hand
- * side the pressure is raised to 42e5 Pa.
- *
- * The induced flow field is calculated by means of the Forchheimer relation.
- * This selection is chosen via the BaseFluxVariables property, which can also
- * be set to the Darcy relation.
- *
- * The setup is on purpose very simple allowing for the analytical calculation
- * of the correct velocity.
- * Just paste the following into e.g. matlab
- * (viscosity, forchheimer coefficient, density, permeability, pressure gradient)
- * -> mu = 1e-03 ; cE =0.55; rho = 1000; K = 1e-12; a = mu /( cE *rho* sqrt(K) );
- * -> gradP = -41e5; v = - a /2 + sqrt (a^2/4-sqrt(K)*gradP/(rho *cE))
- *
- * This problem uses the \ref MPNCModel.
- *
- * To run the simulation execute the following line in shell:
- * <tt>./test_forchheimer1p -parameterFile test_forchheimer1p.input</tt>
- */
-template <class TypeTag>
-class Forchheimer1pProblem
- : public ImplicitPorousMediaProblem<TypeTag>
-{
- using ParentType = ImplicitPorousMediaProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using ParameterCache = typename FluidSystem::ParameterCache;
- using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
- // Grid and world dimension
- enum {dim = GridView::dimension};
- enum {dimWorld = GridView::dimensionworld};
- enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
- enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
- enum {nPhaseIdx = FluidSystem::nPhaseIdx};
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
- enum {wCompIdx = FluidSystem::wCompIdx};
- enum {nCompIdx = FluidSystem::nCompIdx};
- enum {fug0Idx = Indices::fug0Idx};
- enum {s0Idx = Indices::s0Idx};
- enum {p0Idx = Indices::p0Idx};
-
- 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 GlobalPosition = Dune::FieldVector<typename GridView::Grid::ctype, dimWorld>;
- using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
-
-public:
- /*!
- * \brief The constructor
- *
- * \param timeManager The time manager
- * \param gridView The grid view
- */
- Forchheimer1pProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
- {
- pMax_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.pMax);
- pMin_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.pMin);
- outputName_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.outputName);
-
- temperature_ = 273.15 + 25; // -> 25°C
-
- // initialize the tables of the fluid system
- Scalar Tmin = temperature_ - 1.0;
- Scalar Tmax = temperature_ + 1.0;
- int nT = 3;
-
- Scalar pmin = 1.0e5 * 0.75;
- Scalar pmax = 2.0e5 * 1.25;
- int np = 1000;
-
- FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
- }
-
- /*!
- * \brief User defined output after the time integration
- *
- * Will be called diretly after the time integration.
- */
- void postTimeStep()
- {
- // Calculate storage terms of the individual phases
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- PrimaryVariables phaseStorage;
- 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";
- }
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief Returns the problem name
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const std::string name() const
- { return outputName_; }
-
- /*!
- * \brief Returns the temperature \f$ K \f$
- */
- Scalar temperatureAtPos(const GlobalPosition &globalPos) const
- { return temperature_; }
-
- /*!
- * \brief Write a restart file?
- */
- bool shouldWriteRestartFile() const
- {
- return ParentType::shouldWriteRestartFile();
- }
-
- // \}
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param values The boundary types for the conservation equations
- * \param globalPos The global position
- */
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
- {
- if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
- values.setAllDirichlet();
- else
- values.setAllNeumann();
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
- *
- * \param values Stores the Dirichlet values for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} ] \f$
- * \param globalPos The global position
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- initial_(values, globalPos);
- }
-
- /*!
- * \brief Evaluates the boundary conditions for a Neumann
- * boundary segment.
- *
- * \param values Stores the Neumann values for the conservation equations in
- * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
- * \param element The finite element
- * \param fvGeometry The finite volume geometry of the element
- * \param intersection The intersection between element and boundary
- * \param scvIdx The local index of the sub-control volume
- * \param boundaryFaceIdx The index of the boundary face
- *
- * Negative values mean influx.
- */
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- const unsigned int scvIdx,
- const unsigned int boundaryFaceIdx) const
- { values = 0.; }
-
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
- /*!
- * \brief Evaluate the source term for all balance equations within a given
- * sub-control-volume.
- *
- * \param values Stores the solution for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
- * \param element The finite element
- * \param fvGeometry The finite volume geometry of the element
- * \param scvIdx The local index of the sub-control volume
- *
- * Positive values mean that mass is created, negative ones mean that it vanishes.
- */
- void source(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- {
- values = Scalar(0.0);
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values Stores the solution for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} ] \f$
- * \param globalPos The global position
- */
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- initial_(values, globalPos);
- Valgrind::CheckDefined(values);
- }
-
- // \}
-
-private:
- // the internal method for the initial condition
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- FluidState fs;
-
- int refPhaseIdx;
- int otherPhaseIdx;
-
- // set the fluid temperatures
- fs.setTemperature(this->temperatureAtPos(globalPos));
-
-
- if (onLeftBoundary_(globalPos)) {
- // set pressure of the liquid phase
- fs.setPressure(wPhaseIdx, pMax_);
- }
- else {
- const Scalar interpolation = (globalPos[0] - this->bBoxMin()[0] ) / (this->bBoxMax()[0] - this->bBoxMin()[0]);
- const Scalar interpolatedPressure = pMax_ - interpolation * (pMax_ - pMin_) ;
-
- // set pressure of the gas phase
- fs.setPressure(wPhaseIdx, interpolatedPressure);
- }
-
- const Scalar belowSolubilityMassFraction= 1e-5;
- // set the liquid composition to equilibrium concentration
- fs.setMoleFraction(wPhaseIdx, nCompIdx, belowSolubilityMassFraction);
- fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.-belowSolubilityMassFraction);
- // only liquid on inlet
- refPhaseIdx = wPhaseIdx;
- otherPhaseIdx = nPhaseIdx;
-
-
- // set liquid saturation
- fs.setSaturation(wPhaseIdx, 1.);
-
- // set the other saturation
- fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
-
- // calulate the capillary pressure
- const MaterialLawParams &matParams =
- this->spatialParams().materialLawParamsAtPos(globalPos);
- PhaseVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fs);
- fs.setPressure(otherPhaseIdx,
- fs.pressure(refPhaseIdx)
- + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
-
- // make the fluid state consistent with local thermodynamic
- // equilibrium
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
-
- ParameterCache paramCache;
- ComputeFromReferencePhase::solve(fs,
- paramCache,
- refPhaseIdx,
- /*setViscosity=*/false,
- /*setEnthalpy=*/false);
-
- ///////////
- // assign the primary variables
- ///////////
-
- // all N component fugacities
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fs.fugacity(refPhaseIdx, compIdx);
-
- // first M - 1 saturations
- for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
-
- // first pressure
- values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
- }
-
-private:
- /*!
- * \brief Give back whether the testes position (input) is a specific region (left) in the domain
- */
- bool onLeftBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] < this->bBoxMin()[0] + eps_; }
-
- /*!
- * \brief Give back whether the testes position (input) is a specific region (right) in the domain
- */
- bool onRightBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] > this->bBoxMax()[0] - eps_; }
-
- /*!
- * \brief Give back whether the testes 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 testes 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_; }
-
- Scalar temperature_;
- static constexpr Scalar eps_ = 1e-6;;
- std::string outputName_;
- Scalar pMax_, pMin_ ;
-
-};
-} //end namespace
-
-#endif
diff --git a/test/porousmediumflow/mpnc/implicit/forchheimer2pproblem.hh b/test/porousmediumflow/mpnc/implicit/forchheimer2pproblem.hh
deleted file mode 100644
index c87cab3041..0000000000
--- a/test/porousmediumflow/mpnc/implicit/forchheimer2pproblem.hh
+++ /dev/null
@@ -1,444 +0,0 @@
-// -*- 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 Problem for testing the two-phase forchheimer relation.
- */
-#ifndef DUMUX_FORCHEIMER2P_HH
-#define DUMUX_FORCHEIMER2P_HH
-
-#include <dune/common/parametertreeparser.hh>
-
-#include <dumux/porousmediumflow/mpnc/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
-
-#include <dumux/material/fluidsystems/h2on2.hh>
-#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
-#include <dumux/material/fluidstates/compositional.hh>
-
-#include "forchheimerspatialparams.hh"
-
-namespace Dumux
-{
-
-template <class TypeTag>
-class Forchheimer2pProblem;
-
-namespace Properties
-{
-NEW_TYPE_TAG(Forchheimer2pProblem, INHERITS_FROM(BoxMPNC, ForchheimerSpatialParams));
-
-// Set the grid type
-SET_TYPE_PROP(Forchheimer2pProblem, Grid, Dune::YaspGrid<2>);
-
-// Set the problem property
-SET_TYPE_PROP(Forchheimer2pProblem,
- Problem,
- Forchheimer2pProblem<TypeTag>);
-
-
-// Set fluid configuration
-SET_PROP(Forchheimer2pProblem, FluidSystem)
-{ private:
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
- using type = FluidSystems::H2ON2<Scalar, /*useComplexRelations=*/false>;
-};
-
-
-// Enable molecular diffusion of the components?
-SET_BOOL_PROP(Forchheimer2pProblem, EnableDiffusion, false);
-
-// Enable gravity
-SET_BOOL_PROP(Forchheimer2pProblem, ProblemEnableGravity, true);
-
-// Enable the re-use of the jacobian matrix whenever possible?
-SET_BOOL_PROP(Forchheimer2pProblem, ImplicitEnableJacobianRecycling, true);
-
-// decide which type to use for floating values (double / quad)
-SET_TYPE_PROP(Forchheimer2pProblem, Scalar, double);
-
-// decide how to calculate velocity: Darcy / Forchheimer
-SET_TYPE_PROP(Forchheimer2pProblem, BaseFluxVariables, ImplicitForchheimerFluxVariables<TypeTag>);
-
-SET_BOOL_PROP(Forchheimer2pProblem, VtkAddVelocities, true);
-
-// set the linear solver
-SET_TYPE_PROP(Forchheimer2pProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag>);
-}
-
-
-/*!
- * \ingroup MPNCModel
- * \ingroup ImplicitTestProblems
- * \brief Problem where liquid water is injected by means of a
- * dirchlet condition on the left of the domain.
- * Velocity according to Forchheimer.
- *
- * The setup is for testing of the two phase Forchheimer relation.
- * The whole domain is filled with gas at 1e5 Pa. On the left hand
- * side the pressure is raised to 42e5 Pa and the saturation of the
- * water is set to one: Water phase invades the domain.
- *
- * The induced flow field is calculated by means of the Forchheimer relation.
- * This selection is chosen via the BaseFluxVariables property, which can also
- * be set to the Darcy relation.
- *
- * This problem uses the \ref MPNCModel.
- *
- * To run the simulation execute the following line in shell:
- * <tt>./test_forchheimer2p -parameterFile test_forchheimer2p.input</tt>
- */
-template <class TypeTag>
-class Forchheimer2pProblem
- : public ImplicitPorousMediaProblem<TypeTag>
-{
- using ParentType = ImplicitPorousMediaProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using ParameterCache = typename FluidSystem::ParameterCache;
- using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
- // Grid and world dimension
- enum {dim = GridView::dimension};
- enum {dimWorld = GridView::dimensionworld};
- enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
- enum {numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
- enum {nPhaseIdx = FluidSystem::nPhaseIdx};
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
- enum {wCompIdx = FluidSystem::wCompIdx};
- enum {nCompIdx = FluidSystem::nCompIdx};
- enum {fug0Idx = Indices::fug0Idx};
- enum {s0Idx = Indices::s0Idx};
- enum {p0Idx = Indices::p0Idx};
-
- 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 GlobalPosition = Dune::FieldVector<typename GridView::Grid::ctype, dimWorld>;
- using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
-
-public:
- /*!
- * \brief The constructor
- *
- * \param timeManager The time manager
- * \param gridView The grid view
- */
- Forchheimer2pProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
- {
- pMax_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.pMax);
- pMin_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.pMin);
- outputName_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.outputName);
-
- temperature_ = 273.15 + 25; // -> 25°C
-
- // initialize the tables of the fluid system
- Scalar Tmin = temperature_ - 1.0;
- Scalar Tmax = temperature_ + 1.0;
- int nT = 3;
-
- Scalar pmin = 1.0e5 * 0.75;
- Scalar pmax = 2.0e5 * 1.25;
- int np = 1000;
-
- FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
- }
-
- /*!
- * \brief User defined output after the time integration
- *
- * Will be called diretly after the time integration.
- */
- void postTimeStep()
- {
- // Calculate storage terms of the individual phases
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- PrimaryVariables phaseStorage;
- 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";
- }
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief Returns the problem name
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const std::string name() const
- { return outputName_; }
-
- /*!
- * \brief Returns the temperature \f$ K \f$
- */
- Scalar temperatureAtPos(const GlobalPosition &globalPos) const
- { return temperature_; }
-
- /*!
- * \brief Write a restart file?
- */
- bool shouldWriteRestartFile() const
- {
- return ParentType::shouldWriteRestartFile();
- }
-
- // \}
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param values The boundary types for the conservation equations
- * \param globalPos The global position
- */
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
- {
- if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
- values.setAllDirichlet();
- else
- values.setAllNeumann();
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
- *
- * \param values Stores the Dirichlet values for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} ] \f$
- * \param globalPos The global position
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- initial_(values, globalPos);
- }
-
- /*!
- * \brief Evaluates the boundary conditions for a Neumann
- * boundary segment.
- *
- * \param values Stores the Neumann values for the conservation equations in
- * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
- * \param element The finite element
- * \param fvGeometry The finite volume geometry of the element
- * \param intersection The intersection between element and boundary
- * \param scvIdx The local index of the sub-control volume
- * \param boundaryFaceIdx The index of the boundary face
- *
- * Negative values mean influx.
- */
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- const unsigned int scvIdx,
- const unsigned int boundaryFaceIdx) const
- { values = 0.; }
-
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
- /*!
- * \brief Evaluate the source term for all balance equations within a given
- * sub-control-volume.
- *
- * \param values Stores the solution for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
- * \param element The finite element
- * \param fvGeometry The finite volume geometry of the element
- * \param scvIdx The local index of the sub-control volume
- *
- * Positive values mean that mass is created, negative ones mean that it vanishes.
- */
- void source(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- {
- values = Scalar(0.0);
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values Stores the solution for the conservation equations in
- * \f$ [ \textnormal{unit of primary variable} ] \f$
- * \param globalPos The global position
- */
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- initial_(values, globalPos);
- Valgrind::CheckDefined(values);
- }
-
- // \}
-
-private:
- // the internal method for the initial condition
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- FluidState fs;
-
- int refPhaseIdx;
- int otherPhaseIdx;
-
- // set the fluid temperatures
- fs.setTemperature(this->temperatureAtPos(globalPos));
-
- // invasion over the left hand side domain border
- if (onLeftBoundary_(globalPos)) {
- // only liquid on inlet
- refPhaseIdx = wPhaseIdx;
- otherPhaseIdx = nPhaseIdx;
-
- // set liquid saturation
- fs.setSaturation(wPhaseIdx, 1.);
-
- // set pressure of the liquid phase
- fs.setPressure(wPhaseIdx, pMax_);
-
- // set the liquid composition to pure water
- fs.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
- fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);
- }
- else {
- // elsewhere, only gas
- refPhaseIdx = nPhaseIdx;
- otherPhaseIdx = wPhaseIdx;
-
- // set gas saturation
- fs.setSaturation(nPhaseIdx, 1.);
-
- // set pressure of the gas phase
- fs.setPressure(nPhaseIdx, pMin_);
-
- // set the gas composition to 99% nitrogen and 1% steam
- fs.setMoleFraction(nPhaseIdx, nCompIdx, 0.99);
- fs.setMoleFraction(nPhaseIdx, wCompIdx, 0.01);
- }
-
- // set the other saturation
- fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
-
- // calulate the capillary pressure
- const MaterialLawParams &matParams =
- this->spatialParams().materialLawParamsAtPos(globalPos);
- PhaseVector pc;
- MaterialLaw::capillaryPressures(pc, matParams, fs);
- fs.setPressure(otherPhaseIdx,
- fs.pressure(refPhaseIdx)
- + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
-
- // make the fluid state consistent with local thermodynamic
- // equilibrium
- using ComputeFromReferencePhase = ComputeFromReferencePhase<Scalar, FluidSystem>;
-
- ParameterCache paramCache;
- ComputeFromReferencePhase::solve(fs,
- paramCache,
- refPhaseIdx,
- /*setViscosity=*/false,
- /*setEnthalpy=*/false);
-
- ///////////
- // assign the primary variables
- ///////////
-
- // all N component fugacities
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- values[fug0Idx + compIdx] = fs.fugacity(refPhaseIdx, compIdx);
-
- // first M - 1 saturations
- for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
- values[s0Idx + phaseIdx] = fs.saturation(phaseIdx);
-
- // first pressure
- values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
- }
-
-private:
- /*!
- * \brief Give back whether the testes position (input) is a specific region (left) in the domain
- */
- bool onLeftBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] < this->bBoxMin()[0] + eps_; }
-
- /*!
- * \brief Give back whether the testes position (input) is a specific region (right) in the domain
- */
- bool onRightBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[0] > this->bBoxMax()[0] - eps_; }
-
- Scalar temperature_;
- static constexpr Scalar eps_ = 1e-6;
- std::string outputName_;
- Scalar pMax_, pMin_ ;
-
-};
-} //end namespace
-
-#endif
diff --git a/test/porousmediumflow/mpnc/implicit/forchheimerspatialparams.hh b/test/porousmediumflow/mpnc/implicit/forchheimerspatialparams.hh
deleted file mode 100644
index f8869ad466..0000000000
--- a/test/porousmediumflow/mpnc/implicit/forchheimerspatialparams.hh
+++ /dev/null
@@ -1,202 +0,0 @@
-// -*- 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 Definition of the spatial params properties for the Forchheimer problem
- *
- */
-
-#ifndef DUMUX_FORCHHEIMER_SPATIAL_PARAMS_HH
-#define DUMUX_FORCHHEIMER_SPATIAL_PARAMS_HH
-
-#include <dumux/material/spatialparams/fv.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-
-#include <dumux/porousmediumflow/mpnc/implicit/model.hh>
-#include <dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/mp/2padapter.hh>
-
-namespace Dumux
-{
-
-//forward declaration
-template<class TypeTag>
-class ForchheimerSpatialParams;
-
-namespace Properties
-{
-// The spatial parameters TypeTag
-NEW_TYPE_TAG(ForchheimerSpatialParams);
-
-// Set the spatial parameters
-SET_TYPE_PROP(ForchheimerSpatialParams, SpatialParams, ForchheimerSpatialParams<TypeTag>);
-
-// Set the material Law
-SET_PROP(ForchheimerSpatialParams, MaterialLaw)
-{
-private:
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
- // define the material law
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using EffMaterialLaw = RegularizedLinearMaterial<Scalar>;
- using TwoPMaterialLaw = EffToAbsLaw<EffMaterialLaw>;
-public:
- using type = TwoPAdapter<wPhaseIdx, TwoPMaterialLaw>;
-};
-}
-
-/**
- * \ingroup MPNCModel
- * \ingroup ImplicitTestProblems
- * \brief Definition of the spatial params properties for the Forchheimer problem
- *
- */
-template<class TypeTag>
-class ForchheimerSpatialParams : public FVSpatialParams<TypeTag>
-{
- using ParentType = FVSpatialParams<TypeTag>;
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using CoordScalar = typename Grid::ctype;
-
- enum {dimWorld=GridView::dimensionworld};
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
-
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
-
-public:
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename MaterialLaw::Params;
-
- ForchheimerSpatialParams(const GridView &gridView)
- : ParentType(gridView)
- {
- // intrinsic permeabilities
- K_ = 1e-12;
-
- // the porosity
- porosity_ = 0.3;
-
- // residual saturations
- materialParams_.setSwr(0.0);
- materialParams_.setSnr(0.0);
-
- // parameters for the linear law, i.e. minimum and maximum
- // pressures
- materialParams_.setEntryPc(0.0);
- materialParams_.setMaxPc(0.0);
- }
-
- ~ForchheimerSpatialParams()
- {}
-
- /*!
- * \brief Update the spatial parameters with the flow solution
- * after a timestep.
- *
- * \param globalSol The current solution vector
- */
- void update(const SolutionVector &globalSol)
- {
- }
-
- /*!
- * \brief Returns the intrinsic permeability tensor.
- *
- * \param element The current finite element
- * \param fvGeometry The current finite volume geometry of the element
- * \param scvIdx The index sub-control volume where the
- * intrinsic permeability is given.
- */
- Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- {
- return K_;
- }
-
- /*!
- * \brief Define the porosity \f$[-]\f$ of the soil
- *
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume where
- * the porosity needs to be defined
- */
- double porosity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- {
- return porosity_;
- }
-
- /*!
- * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
- *
- * \param pos The global position of the sub-control volume.
- * \return the material parameters object
- */
- const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition &pos) const
- {
- return materialParams_;
- }
-
-
- /*!
- * \brief Apply the Forchheimer coefficient for inertial forces
- * calculation.
- *
- * Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
- * Actually the Forchheimer coefficient is also a function of the dimensions of the
- * porous medium. Taking it as a constant is only a first approximation
- * (Nield, Bejan, Convection in porous media, 2006, p. 10)
- *
- * \param element The current finite element
- * \param fvGeometry The current finite volume geometry of the element
- * \param scvIdx The index sub-control volume face where the
- * intrinsic velocity ought to be calculated.
- *
- */
- Scalar forchCoeff(const Element &element,
- const FVElementGeometry &fvGeometry,
- const unsigned int scvIdx) const
- {
- // If there are better measures / estimates / values available than this default number:
- // here is the place to implement it
- return ParentType::forchCoeff(element,
- fvGeometry,
- scvIdx);
- }
-
- Scalar K_;
- Scalar porosity_;
- MaterialLawParams materialParams_;
-};
-
-}
-
-#endif
diff --git a/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.cc b/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.cc
deleted file mode 100644
index 551914055e..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.cc
+++ /dev/null
@@ -1,59 +0,0 @@
-// -*- 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 Forchheimer flux calculation with one phase.
- */
-#include <config.h>
-#include "forchheimer1pproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \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.
- */
-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 += "\n\nThe List of Mandatory arguments for this program is:\n"
- "\t-TimeManager.TEnd The end of the simulation. [s] \n"
- "\t-TimeManager.DtInitial The initial timestep size. [s] \n"
- "\t-Grid.File The 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)
-{
- using ProblemTypeTag = TTAG(Forchheimer1pProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.input b/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.input
deleted file mode 100644
index 0acf3f7e82..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_forchheimer1p.input
+++ /dev/null
@@ -1,13 +0,0 @@
-[TimeManager]
-DtInitial = 1e7 # [s]
-TEnd = 1e7 # [s]
-
-[Grid]
-UpperRight = 1
-Cells = 100
-
-[Problem]
-pMin = 1e5
-pMax = 42e5
-outputName = test_forchheimer1p
-
diff --git a/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.cc b/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.cc
deleted file mode 100644
index 51df73295c..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.cc
+++ /dev/null
@@ -1,59 +0,0 @@
-// -*- 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 Forchheimer flux calculation with two phases.
- */
-#include <config.h>
-#include "forchheimer2pproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \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.
- */
-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 += "\n\nThe List of Mandatory arguments for this program is:\n"
- "\t-TimeManager.TEnd The end of the simulation. [s] \n"
- "\t-TimeManager.DtInitial The initial timestep size. [s] \n"
- "\t-Grid.File The 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)
-{
- using ProblemTypeTag = TTAG(Forchheimer2pProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.input b/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.input
deleted file mode 100644
index 90028d6f98..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_forchheimer2p.input
+++ /dev/null
@@ -1,13 +0,0 @@
-[TimeManager]
-DtInitial = 330 # [s]
-TEnd = 3e3 # [s]
-
-[Grid]
-UpperRight = 60 40
-Cells = 24 16
-
-[Problem]
-pMin = 1e5
-pMax = 42e5
-outputName = test_forchheimer2p
-
--
GitLab