From 5502a7c3d60c17425787953559b6d1627f1e23d7 Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Wed, 20 Dec 2017 16:15:39 +0100
Subject: [PATCH] [feature] port obstacletest box and cc
---
.../mpnc/implicit/CMakeLists.txt | 76 +++---
.../mpnc/implicit/obstacleproblem.hh | 177 +++++-------
.../mpnc/implicit/obstaclespatialparams.hh | 68 ++---
.../mpnc/implicit/test_boxmpnc.cc | 59 ----
.../mpnc/implicit/test_boxmpnc.input | 14 -
.../mpnc/implicit/test_ccmpnc.cc | 59 ----
.../{test_ccmpnc.input => test_mpnc.input} | 4 +-
.../mpnc/implicit/test_mpnc_obstacle_fv.cc | 252 ++++++++++++++++++
8 files changed, 377 insertions(+), 332 deletions(-)
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_boxmpnc.cc
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_boxmpnc.input
delete mode 100644 test/porousmediumflow/mpnc/implicit/test_ccmpnc.cc
rename test/porousmediumflow/mpnc/implicit/{test_ccmpnc.input => test_mpnc.input} (80%)
create mode 100644 test/porousmediumflow/mpnc/implicit/test_mpnc_obstacle_fv.cc
diff --git a/test/porousmediumflow/mpnc/implicit/CMakeLists.txt b/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
index 5ca0dbc762..880da4effd 100644
--- a/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/mpnc/implicit/CMakeLists.txt
@@ -1,47 +1,44 @@
add_input_file_links()
-add_dumux_test(test_boxmpnc test_boxmpnc test_boxmpnc.cc
- ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/obstaclebox-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/obstaclebox-00009.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnc")
+dune_add_test(NAME test_boxmpnc
+ SOURCES test_mpnc_obstacle_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=ObstacleBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/obstaclebox-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnc-00009.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnc test_mpnc.input -Problem.Name test_boxmpnc")
-add_dumux_test(test_ccmpnc test_ccmpnc test_ccmpnc.cc
- ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/obstaclecc-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/obstaclecc-00009.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccmpnc")
-
-# build target for the one-phase Forchheimer test problem
-add_dumux_test(test_forchheimer1p test_forchheimer1p test_forchheimer1p.cc
- ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/forchheimer1p-reference.vtp
- ${CMAKE_CURRENT_BINARY_DIR}/test_forchheimer1p-00001.vtp
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_forchheimer1p")
-
-# build target for the two-phase Forchheimer test problem
-add_dumux_test(test_forchheimer2p test_forchheimer2p test_forchheimer2p.cc
- ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/forchheimer2p-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_forchheimer2p-00008.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_forchheimer2p")
+dune_add_test(NAME test_ccmpnc
+ SOURCES test_mpnc_obstacle_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=ObstacleCCProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/obstaclecc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_ccmpnc-00009.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccmpnc test_mpnc.input -Problem.Name test_ccmpnc")
# build target for the full kinetic test problem
-add_dumux_test(test_boxmpnckinetic test_boxmpnckinetic test_boxmpnckinetic.cc
- ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/evaporationatmosphere-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/evaporationatmosphere-00005.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnckinetic"
- --zeroThreshold {"TnMinusTs":0.5,"TwMinusTn":0.5,"velocity_w_1":2e-12})
+dune_add_test(NAME test_boxmpnckinetic
+ SOURCES test_boxmpnckinetic.cc
+ COMPILE_DEFINITIONS TYPETAG=EvaporationAtmosphereBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/evaporationatmosphere-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnckinetic-00011.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpnckinetic test_boxmpnckinetic.input -Problem.Name test_boxmpnckinetic")
# build target for the energy kinetic test problem, two energy balance equations
-add_dumux_test(test_boxmpncthermalnonequil test_boxmpncthermalnonequil test_boxmpncthermalnonequil.cc
- "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpncthermalnonequil")
+dune_add_test(NAME test_boxmpncthermalnonequil
+ SOURCES test_boxmpncthermalnonequil.cc
+ COMPILE_DEFINITIONS TYPETAG=CombustionProblemOneComponentBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/combustion-reference.vtp
+ ${CMAKE_CURRENT_BINARY_DIR}/test_boxmpncthermalnonequil-00084.vtp
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxmpncthermalnonequil test_boxmpncthermalnonequil.input -Problem.Name test_boxmpncthermalnonequil")
+
+
#install sources
install(FILES
@@ -49,15 +46,10 @@ combustionproblem1c.hh
combustionspatialparams.hh
evaporationatmosphereproblem.hh
evaporationatmospherespatialparams.hh
-forchheimer1pproblem.hh
-forchheimer2pproblem.hh
-forchheimerspatialparams.hh
obstacleproblem.hh
obstaclespatialparams.hh
test_boxmpnc.cc
test_boxmpnckinetic.cc
test_boxmpncthermalnonequil.cc
test_ccmpnc.cc
-test_forchheimer1p.cc
-test_forchheimer2p.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/mpnc)
diff --git a/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh b/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
index 921f8925e9..561f42979c 100644
--- a/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
@@ -27,12 +27,15 @@
#include <dune/common/parametertreeparser.hh>
-#include <dumux/porousmediumflow/mpnc/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/mpnc/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
-#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/fluidstates/compositional.hh>
+#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include "obstaclespatialparams.hh"
@@ -46,7 +49,7 @@ namespace Properties
{
NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(MPNC, ObstacleSpatialParams));
NEW_TYPE_TAG(ObstacleBoxProblem, INHERITS_FROM(BoxModel, ObstacleProblem));
-NEW_TYPE_TAG(ObstacleCCProblem, INHERITS_FROM(CCModel, ObstacleProblem));
+NEW_TYPE_TAG(ObstacleCCProblem, INHERITS_FROM(CCTpfaModel, ObstacleProblem));
// Set the grid type
SET_TYPE_PROP(ObstacleProblem, Grid, Dune::YaspGrid<2>);
@@ -57,24 +60,9 @@ SET_TYPE_PROP(ObstacleProblem,
ObstacleProblem<TypeTag>);
// Set fluid configuration
-SET_PROP(ObstacleProblem, FluidSystem)
-{ private:
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
- using type = FluidSystems::H2ON2<Scalar, /*useComplexRelations=*/false>;
-};
-
-// Enable smooth upwinding?
-SET_BOOL_PROP(ObstacleProblem, ImplicitEnableSmoothUpwinding, true);
-
-// Enable molecular diffusion of the components?
-SET_BOOL_PROP(ObstacleProblem, EnableDiffusion, true);
-
-// Enable the re-use of the jacobian matrix whenever possible?
-SET_BOOL_PROP(ObstacleProblem, ImplicitEnableJacobianRecycling, true);
-
-// Reassemble the jacobian matrix only where it changed?
-SET_BOOL_PROP(ObstacleProblem, ImplicitEnablePartialReassemble, true);
+SET_TYPE_PROP(ObstacleProblem,
+ FluidSystem,
+ FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/false>);
// decide which type to use for floating values (double / quad)
SET_TYPE_PROP(ObstacleProblem, Scalar, double);
@@ -111,19 +99,28 @@ SET_TYPE_PROP(ObstacleProblem, Scalar, double);
*/
template <class TypeTag>
class ObstacleProblem
- : public ImplicitPorousMediaProblem<TypeTag>
+ : public PorousMediumFlowProblem<TypeTag>
{
- using ParentType = ImplicitPorousMediaProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
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 Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using ParameterCache = typename FluidSystem::ParameterCache;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ 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 MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using ParameterCache = typename FluidSystem::ParameterCache;
// world dimension
enum {dimWorld = GridView::dimensionworld};
@@ -137,13 +134,10 @@ class ObstacleProblem
enum {s0Idx = Indices::s0Idx};
enum {p0Idx = Indices::p0Idx};
- using Element = typename GridView::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using GlobalPosition = Dune::FieldVector<typename GridView::Grid::ctype, dimWorld>;
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
public:
/*!
@@ -152,8 +146,8 @@ public:
* \param timeManager The time manager
* \param gridView The grid view
*/
- ObstacleProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ ObstacleProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
temperature_ = 273.15 + 25; // -> 25°C
@@ -167,42 +161,7 @@ public:
int np = 1000;
FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
- }
-
- /*!
- * \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_ = getParam<std::string>("Problem.Name");
}
/*!
@@ -223,7 +182,7 @@ public:
*
* \param globalPos The global position
*/
- Scalar temperatureAtPos(const GlobalPosition &globalPos) const
+ Scalar temperature() const
{ return temperature_; }
/*!
@@ -240,37 +199,31 @@ public:
* \name Boundary conditions
*/
// \{
-
/*!
* \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary control volume.
+ * used for which equation on a given boundary segment
*
- * \param values The boundary types for the conservation equations
- * \param globalPos The position of the center of the finite volume
+ * \param globalPos The global position
*/
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
+ BoundaryTypes bcTypes;
if (onInlet_(globalPos) || onOutlet_(globalPos))
- values.setAllDirichlet();
+ bcTypes.setAllDirichlet();
else
- values.setAllNeumann();
+ bcTypes.setAllNeumann();
+ return bcTypes;
}
/*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
+ * \brief Evaluates 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 center of the finite volume which ought to be set.
- *
- * For this method, the \a values parameter stores primary variables.
+ * \param globalPos The global position
*/
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initial_(globalPos);
}
/*!
@@ -287,13 +240,13 @@ public:
*
* 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.; }
+ PrimaryVariables neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ return PrimaryVariables(0.0);
+ }
// \}
@@ -314,12 +267,13 @@ public:
*
* 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
+ //! \copydoc Dumux::ImplicitProblem::source()
+ PrimaryVariables source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
{
- values = Scalar(0.0);
+ return PrimaryVariables(0.0);
}
/*!
@@ -331,20 +285,18 @@ public:
* For this method, the \a values parameter stores primary
* variables.
*/
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
- Valgrind::CheckDefined(values);
+ return initial_(globalPos);
}
// \}
private:
// the internal method for the initial condition
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
+ PrimaryVariables values(0.0);
FluidState fs;
int refPhaseIdx;
@@ -389,7 +341,7 @@ private:
fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
// calulate the capillary pressure
- const MaterialLawParams &matParams =
+ const auto& matParams =
this->spatialParams().materialLawParamsAtPos(globalPos);
PhaseVector pc;
MaterialLaw::capillaryPressures(pc, matParams, fs);
@@ -422,6 +374,7 @@ private:
// first pressure
values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
+ return values;
}
bool onInlet_(const GlobalPosition &globalPos) const
diff --git a/test/porousmediumflow/mpnc/implicit/obstaclespatialparams.hh b/test/porousmediumflow/mpnc/implicit/obstaclespatialparams.hh
index e9565a2c89..16bb62b5b5 100644
--- a/test/porousmediumflow/mpnc/implicit/obstaclespatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/obstaclespatialparams.hh
@@ -29,7 +29,6 @@
#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>
@@ -73,26 +72,24 @@ template<class TypeTag>
class ObstacleSpatialParams : public FVSpatialParams<TypeTag>
{
using ParentType = FVSpatialParams<TypeTag>;
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- 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 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;
+ 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;
enum {dimWorld=GridView::dimensionworld};
- enum {wPhaseIdx = FluidSystem::wPhaseIdx};
-
- using Element = typename GridView::template Codim<0>::Entity;
- using DimWorldVector = Dune::FieldVector<CoordScalar, dimWorld>;
public:
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
- ObstacleSpatialParams(const GridView &gridView)
- : ParentType(gridView)
+
+ ObstacleSpatialParams(const Problem &problem)
+ : ParentType(problem)
{
// intrinsic permeabilities
coarseK_ = 1e-12;
@@ -118,31 +115,14 @@ public:
~ObstacleSpatialParams()
{}
- /*!
- * \brief Update the spatial parameters with the flow solution
- * after a timestep.
- *
- * \param globalSol The current solution vector
- */
- void update(const SolutionVector &globalSol)
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
{
- }
-
- /*!
- * \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
- {
- if (isFineMaterial_(fvGeometry.subContVol[scvIdx].global))
+ if (isFineMaterial_(scv.dofPosition()))
return fineK_;
- return coarseK_;
+ else
+ return coarseK_;
}
/*!
@@ -153,9 +133,9 @@ public:
* \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
+ Scalar porosity(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolutionVector &elemSol) const
{
return porosity_;
}
@@ -166,9 +146,9 @@ public:
* \param pos The global position of the sub-control volume.
* \return the material parameters object
*/
- const MaterialLawParams& materialLawParamsAtPos(const DimWorldVector &pos) const
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
- if (isFineMaterial_(pos))
+ if (isFineMaterial_(globalPos))
return fineMaterialParams_;
else
return coarseMaterialParams_;
@@ -179,7 +159,7 @@ private:
* \brief Returns whether a given global position is in the
* fine-permeability region or not.
*/
- static bool isFineMaterial_(const DimWorldVector &pos)
+ static bool isFineMaterial_(const GlobalPosition &pos)
{
return
10 - eps_ <= pos[0] && pos[0] <= 20 + eps_ &&
diff --git a/test/porousmediumflow/mpnc/implicit/test_boxmpnc.cc b/test/porousmediumflow/mpnc/implicit/test_boxmpnc.cc
deleted file mode 100644
index 64f79caa08..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_boxmpnc.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 MpNc box model.
- */
-#include <config.h>
-#include "obstacleproblem.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 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)
-{
- using ProblemTypeTag = TTAG(ObstacleBoxProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/mpnc/implicit/test_boxmpnc.input b/test/porousmediumflow/mpnc/implicit/test_boxmpnc.input
deleted file mode 100644
index 1c58b3803a..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_boxmpnc.input
+++ /dev/null
@@ -1,14 +0,0 @@
-[TimeManager]
-DtInitial = 250 # [s]
-TEnd = 1e4 # [s]
-
-[Grid]
-UpperRight = 60 40
-Cells = 24 16
-
-[LinearSolver]
-ResidualReduction = 1e-12
-
-[Problem]
-Name = obstaclebox
-
diff --git a/test/porousmediumflow/mpnc/implicit/test_ccmpnc.cc b/test/porousmediumflow/mpnc/implicit/test_ccmpnc.cc
deleted file mode 100644
index 6c0340f8d9..0000000000
--- a/test/porousmediumflow/mpnc/implicit/test_ccmpnc.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 MpNc CC model.
- */
-#include <config.h>
-#include "obstacleproblem.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 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)
-{
- using ProblemTypeTag = TTAG(ObstacleCCProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/mpnc/implicit/test_ccmpnc.input b/test/porousmediumflow/mpnc/implicit/test_mpnc.input
similarity index 80%
rename from test/porousmediumflow/mpnc/implicit/test_ccmpnc.input
rename to test/porousmediumflow/mpnc/implicit/test_mpnc.input
index dbb4c29d38..0ca290b167 100644
--- a/test/porousmediumflow/mpnc/implicit/test_ccmpnc.input
+++ b/test/porousmediumflow/mpnc/implicit/test_mpnc.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 250 # [s]
TEnd = 1e4 # [s]
@@ -10,5 +10,5 @@ Cells = 24 16
ResidualReduction = 1e-12
[Problem]
-Name = obstaclecc
+Name = obstacle
diff --git a/test/porousmediumflow/mpnc/implicit/test_mpnc_obstacle_fv.cc b/test/porousmediumflow/mpnc/implicit/test_mpnc_obstacle_fv.cc
new file mode 100644
index 0000000000..664daa0e1c
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/test_mpnc_obstacle_fv.cc
@@ -0,0 +1,252 @@
+// -*- 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 mpnc porousmedium box flow model
+ */
+#include <config.h>
+
+#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 "obstacleproblem.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/nonlinear/newtoncontroller.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.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 End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.LowerLeft Lower left corner coordinates\n"
+ "\t-Grid.UpperRight Upper right corner coordinates\n"
+ "\t-Grid.Cells Number of cells in respective coordinate directions\n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+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 SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.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 = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonController = NewtonController<TypeTag>;
+ using NewtonMethod = 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;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
--
GitLab