From 06beaace09f27ae939e048dade5990b8c7c917bc Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 27 Jul 2017 17:53:36 +0200
Subject: [PATCH] [fringe][tut] Improve example for fringe / infiltration
richards / 2p
---
tutorial/capillaryfringe/2pfringeproblem.hh | 54 ++++---
tutorial/capillaryfringe/CMakeLists.txt | 2 +-
tutorial/capillaryfringe/fringe_2p.input | 14 --
.../capillaryfringe/fringe_richards.input | 14 --
.../capillaryfringe/fringespatialparams.hh | 139 +++++++++++++++---
.../capillaryfringe/richardsfringeproblem.hh | 57 ++++---
tutorial/capillaryfringe/soil.input | 42 ++++++
7 files changed, 232 insertions(+), 90 deletions(-)
delete mode 100644 tutorial/capillaryfringe/fringe_2p.input
delete mode 100644 tutorial/capillaryfringe/fringe_richards.input
create mode 100644 tutorial/capillaryfringe/soil.input
diff --git a/tutorial/capillaryfringe/2pfringeproblem.hh b/tutorial/capillaryfringe/2pfringeproblem.hh
index d6eb06e162..f02d97dbb9 100644
--- a/tutorial/capillaryfringe/2pfringeproblem.hh
+++ b/tutorial/capillaryfringe/2pfringeproblem.hh
@@ -142,7 +142,7 @@ public:
TwoPFringeProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name) + "_2p";
}
/*!
@@ -195,6 +195,9 @@ public:
bcTypes.setAllDirichlet();
else
bcTypes.setAllNeumann();
+ static const bool bottomDirichlet = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Problem, BottomDirichlet);
+ if (bottomDirichlet && onLowerBoundary_(globalPos))
+ bcTypes.setAllDirichlet();
return bcTypes;
}
@@ -208,20 +211,23 @@ public:
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- return initialAtPos(globalPos);
- }
+ PrimaryVariables values(0.0);
+ if (onUpperBoundary_(globalPos))
+ {
+ static const Scalar sw = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TopSaturation);
+ const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
+ values[pwIdx] = nonWettingReferencePressure() - pc;
+ values[snIdx] = 1.0-sw;
+ }
+ else if(onLowerBoundary_(globalPos))
+ {
+ static const Scalar sw = 1.0;
+ const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
+ values[pwIdx] = nonWettingReferencePressure() - pc;
+ values[snIdx] = 0.0;
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
- *
- * \param globalPos The position for which the Dirichlet value is set
- *
- * For this method, the \a values parameter stores primary variables.
- */
- PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
- {
- return PrimaryVariables(0.0);
+ }
+ return values;
}
/*!
@@ -239,23 +245,30 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
+ static const Scalar levelFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, WaterLevelFactor);
+ static const Scalar swBottom = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, SwBottomInitial);
+ static const Scalar swTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, SwTopInitial);
+
PrimaryVariables values(0.0);
- if (globalPos[1] > 0.5*(this->bBoxMax()[1] - this->bBoxMin()[1]))
+ if (globalPos[1] > levelFactor*(this->bBoxMax()[1] - this->bBoxMin()[1]))
{
- const Scalar sw = 0.05;
+ const Scalar sw = swTop;
const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
values[pwIdx] = nonWettingReferencePressure() - pc;
values[snIdx] = 1.0-sw;
}
else
{
- const Scalar sw = 1.0;
+ const Scalar sw = swBottom;
const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
values[pwIdx] = nonWettingReferencePressure() - pc;
values[snIdx] = 1.0-sw;
}
return values;
- };
+ }
+
+ bool shouldWriteRestartFile() const
+ { return false; }
// \}
@@ -266,6 +279,11 @@ private:
return globalPos[1] > this->bBoxMax()[1] - eps_;
}
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bBoxMin()[1] + eps_;
+ }
+
static constexpr Scalar eps_ = 1.5e-7;
GlobalPosition lensLowerLeft_;
diff --git a/tutorial/capillaryfringe/CMakeLists.txt b/tutorial/capillaryfringe/CMakeLists.txt
index deb451555e..aeefd0843b 100644
--- a/tutorial/capillaryfringe/CMakeLists.txt
+++ b/tutorial/capillaryfringe/CMakeLists.txt
@@ -1,3 +1,3 @@
-add_input_file_links()
+dune_symlink_to_source_files(FILES soil.input)
dune_add_test(SOURCES fringe_richards.cc)
dune_add_test(SOURCES fringe_2p.cc)
diff --git a/tutorial/capillaryfringe/fringe_2p.input b/tutorial/capillaryfringe/fringe_2p.input
deleted file mode 100644
index 6e19e1f3a7..0000000000
--- a/tutorial/capillaryfringe/fringe_2p.input
+++ /dev/null
@@ -1,14 +0,0 @@
-[TimeManager]
-DtInitial = 3600 # [s]
-TEnd = 36000000 # [s]
-
-[Grid]
-UpperRight = 2 2
-Cells = 1 100
-
-[Problem]
-Name = 2p_fringe
-EnableGravity = true # enable gravity
-
-[Newton]
-EnableChop = true # chop for better convergence
diff --git a/tutorial/capillaryfringe/fringe_richards.input b/tutorial/capillaryfringe/fringe_richards.input
deleted file mode 100644
index 0e6e4c5228..0000000000
--- a/tutorial/capillaryfringe/fringe_richards.input
+++ /dev/null
@@ -1,14 +0,0 @@
-[TimeManager]
-DtInitial = 3600 # [s]
-TEnd = 36000000 # [s]
-
-[Grid]
-UpperRight = 2 2
-Cells = 1 100
-
-[Problem]
-Name = richards_fringe
-EnableGravity = true # enable gravity
-
-[Newton]
-EnableChop = true # chop for better convergence
diff --git a/tutorial/capillaryfringe/fringespatialparams.hh b/tutorial/capillaryfringe/fringespatialparams.hh
index f4a15dc476..187a257a72 100644
--- a/tutorial/capillaryfringe/fringespatialparams.hh
+++ b/tutorial/capillaryfringe/fringespatialparams.hh
@@ -19,24 +19,46 @@
/*!
* \file
*
- * \brief spatial parameters for the RichardsLensProblem
+ * \brief spatial parameters for the TwoPTwoCTestProblem
*/
#ifndef DUMUX_FRINGE_SPATIAL_PARAMETERS_HH
#define DUMUX_FRINGE_SPATIAL_PARAMETERS_HH
#include <dumux/material/spatialparams/implicit.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-#include <dumux/porousmediumflow/richards/implicit/model.hh>
-
namespace Dumux
{
+template<class TypeTag>
+class FringeSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(FringeSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(FringeSpatialParams, SpatialParams, FringeSpatialParams<TypeTag>);
+
+// Set the material law
+SET_PROP(FringeSpatialParams, MaterialLaw)
+{
+private:
+ // define the material law which is parameterized by effective
+ // saturations
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+ // define the material law parameterized by absolute saturations
+ using type = EffToAbsLaw<RegularizedVanGenuchten<Scalar>>;
+};
+} // end namespace Properties
/*!
- * \ingroup RichardsModel
+ * \ingroup TwoPTwoCModel
* \ingroup ImplicitTestProblems
- * \brief The spatial parameters for the RichardsLensProblem
+ * \brief The spatial parameters for the TwoPTwoCTestProblem
*/
template<class TypeTag>
class FringeSpatialParams : public ImplicitSpatialParams<TypeTag>
@@ -68,18 +90,29 @@ public:
FringeSpatialParams(const Problem& problem, const GridView& gridView)
: ParentType(problem, gridView)
{
- // residual saturations
- outerMaterialParams_.setSwr(0.05);
- outerMaterialParams_.setSnr(0.0);
-
- // parameters for the Van Genuchten law
- // alpha and n
- outerMaterialParams_.setVgAlpha(0.0037);
- outerMaterialParams_.setVgn(4.7);
- outerMaterialParams_.setPcLowSw(0.001);
-
- outerK_ = 5e-12;
+ lensLowerLeft_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, LensSpatialParams.LowerLeft);
+ lensUpperRight_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, LensSpatialParams.UpperRight);
+
+ // residual saturations and parameters for the Van Genuchten law
+ materialParams_.setSwr(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Swr));
+ materialParams_.setSnr(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Snr));
+ materialParams_.setVgAlpha(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.VanGenAlpha));
+ materialParams_.setVgn(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.VanGenN));
+ materialParams_.setPcLowSw(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.RegThresholdSw));
+
+ // same for lens: residual saturations and parameters for the Van Genuchten law
+ materialParamsLens_.setSwr(GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.Swr));
+ materialParamsLens_.setSnr(GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.Snr));
+ materialParamsLens_.setVgAlpha(GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.VanGenAlpha));
+ materialParamsLens_.setVgn(GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.VanGenN));
+ materialParamsLens_.setPcLowSw(GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.RegThresholdSw));
+
+ permeability_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Permeability);
+ lensPermeability_ = GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.Permeability);
+ porosity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Porosity);
+ lensPorosity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, LensSpatialParams.Porosity);
+ hasLens_ = GET_RUNTIME_PARAM(TypeTag, bool, LensSpatialParams.HasLens);
}
/*!
@@ -89,7 +122,11 @@ public:
*/
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
- return outerK_;
+ if (isInLens_(globalPos))
+ return lensPermeability_;
+ else
+ return permeability_;
+
}
/*!
@@ -98,27 +135,81 @@ public:
* \param globalPos The global position where we evaluate
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
- { return 0.4; }
+ {
+ if (isInLens_(globalPos))
+ return lensPorosity_;
+ else
+ return porosity_;
+ }
/*!
* \brief Returns the parameters for the material law at a given location
*
- * This method is not actually required by the Richards model, but provided
- * for the convenience of the RichardsLensProblem
+ * This method is not actually required by the TwoPTwoC model, but provided
+ * for the convenience of the TwoPTwoCLensProblem
*
* \param globalPos A global coordinate vector
*/
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition &globalPos) const
{
- return outerMaterialParams_;
+ if (isInLens_(globalPos))
+ return materialParamsLens_;
+ else
+ return materialParams_;
}
+ /*!
+ * \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 1480; /*specific heat capacity of granite [J / (kg K)]*/ }
+
+ /*!
+ * \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 1600; /*density of granite [kg/m^3]*/ }
+
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param globalPos The global position
+ */
+ Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
+ { return 2; }
+
private:
+ bool isInLens_(const GlobalPosition &globalPos) const
+ {
+ if(!hasLens_) return false;
+
+ for (int i = 0; i < dimWorld; ++i)
+ if (globalPos[i] < lensLowerLeft_[i] - eps_ || globalPos[i] > lensUpperRight_[i] + eps_)
+ return false;
+
+ return true;
+ }
+ static constexpr Scalar eps_ = 1e-7;
+
+ MaterialLawParams materialParams_;
+ MaterialLawParams materialParamsLens_;
+ Scalar permeability_, lensPermeability_;
+ Scalar porosity_, lensPorosity_;
- static constexpr Scalar eps_ = 1.5e-7;
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
- Scalar outerK_;
- MaterialLawParams outerMaterialParams_;
+ bool hasLens_;
};
} // end namespace Dumux
diff --git a/tutorial/capillaryfringe/richardsfringeproblem.hh b/tutorial/capillaryfringe/richardsfringeproblem.hh
index de03cba329..334bcc2fc9 100644
--- a/tutorial/capillaryfringe/richardsfringeproblem.hh
+++ b/tutorial/capillaryfringe/richardsfringeproblem.hh
@@ -130,7 +130,7 @@ public:
RichardsFringeProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name) + "_richards";
}
/*!
@@ -183,6 +183,9 @@ public:
bcTypes.setAllDirichlet();
else
bcTypes.setAllNeumann();
+ static const bool bottomDirichlet = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Problem, BottomDirichlet);
+ if (bottomDirichlet && onLowerBoundary_(globalPos))
+ bcTypes.setAllDirichlet();
return bcTypes;
}
@@ -196,20 +199,22 @@ public:
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- return initialAtPos(globalPos);
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
- *
- * \param globalPos The position for which the Dirichlet value is set
- *
- * For this method, the \a values parameter stores primary variables.
- */
- PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
- {
- return PrimaryVariables(0.0);
+ PrimaryVariables values(0.0);
+ if (onUpperBoundary_(globalPos))
+ {
+ static const Scalar sw = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TopSaturation);
+ const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
+ values[pressureIdx] = nonWettingReferencePressure() - pc;
+ values.setState(Indices::bothPhases);
+ }
+ else if(onLowerBoundary_(globalPos))
+ {
+ static const Scalar sw = 1.0;
+ const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
+ values[pressureIdx] = nonWettingReferencePressure() - pc;
+ values.setState(Indices::bothPhases);
+ }
+ return values;
}
/*!
@@ -227,21 +232,30 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
+ static const Scalar levelFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, WaterLevelFactor);
+ static const Scalar swBottom = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, SwBottomInitial);
+ static const Scalar swTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, SwTopInitial);
+
PrimaryVariables values(0.0);
- if (globalPos[1] > 0.5*(this->bBoxMax()[1] - this->bBoxMin()[1]))
+ if (globalPos[1] > levelFactor*(this->bBoxMax()[1] - this->bBoxMin()[1]))
{
- const Scalar sw = 0.05;
+ const Scalar sw = swTop;
const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
values[pressureIdx] = nonWettingReferencePressure() - pc;
+ values.setState(Indices::bothPhases);
}
else
{
- const Scalar sw = 1.0;
+ const Scalar sw = swBottom;
const Scalar pc = MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos(globalPos), sw);
values[pressureIdx] = nonWettingReferencePressure() - pc;
+ values.setState(Indices::bothPhases);
}
return values;
- };
+ }
+
+ bool shouldWriteRestartFile() const
+ { return false; }
// \}
@@ -252,6 +266,11 @@ private:
return globalPos[1] > this->bBoxMax()[1] - eps_;
}
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bBoxMin()[1] + eps_;
+ }
+
static constexpr Scalar eps_ = 1.5e-7;
GlobalPosition lensLowerLeft_;
diff --git a/tutorial/capillaryfringe/soil.input b/tutorial/capillaryfringe/soil.input
new file mode 100644
index 0000000000..106238a27f
--- /dev/null
+++ b/tutorial/capillaryfringe/soil.input
@@ -0,0 +1,42 @@
+[TimeManager]
+DtInitial = 600 # [s]
+MaxTimeStepSize = 600 # [s]
+TEnd = 36000 # [s]
+
+[Grid]
+UpperRight = 2 2
+Cells = 1 100
+
+[Problem]
+Name = fringe
+EnableGravity = true # enable gravity
+UsePrimaryVariableSwitch = false
+TopSaturation = 0.05
+SwBottomInitial = 1.0
+SwTopInitial = 0.05
+WaterLevelFactor = 0.5
+BottomDirichlet = true
+
+[SpatialParams]
+Permeability = 2.57e-12 # [m^2]
+Porosity = 0.3
+Swr = 0.05
+Snr = 0.0
+RegThresholdSw = 0.01
+VanGenAlpha = 5.8e-4
+VanGenN = 1.8
+
+[LensSpatialParams]
+HasLens = true
+LowerLeft = 0.1 1.4
+UpperRight = 1.8 1.7
+Permeability = 2.57e-16 # [m^2]
+Porosity = 0.03
+Swr = 0.1
+Snr = 0.0
+RegThresholdSw = 0.1
+VanGenAlpha = 0.008 #5.8e-4
+VanGenN = 1.68 #17.8
+
+[Newton]
+EnableChop = true # chop for better convergence
--
GitLab