From ed054eac030d9df55ef03779e0e197fe4c963e1a Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Mon, 13 Feb 2017 20:26:17 +0100
Subject: [PATCH] [test][richardsnc] Add output of tracer mass
---
.../1p2c_richards2c/richardstestproblem.hh | 35 ++++--
.../1p2c_richards2c/rootsystemtestproblem.hh | 28 +++--
.../1p2c_richards2c/rositestproblem.hh | 110 +++++++++++++++++-
.../1p2c_richards2c/test_rosi2c.input | 2 +-
4 files changed, 149 insertions(+), 26 deletions(-)
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
index d71331baba..e2892d6bd1 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/richardstestproblem.hh
@@ -29,6 +29,9 @@
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/porousmediumflow/richardsnc/implicit/model.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/liquidphase2c.hh>
//! get the properties needed for subproblems
#include <dumux/mixeddimension/subproblemproperties.hh>
@@ -59,6 +62,13 @@ SET_BOOL_PROP(RichardsTestProblem, SolutionDependentHeatConduction, false);
// Set the problem property
SET_TYPE_PROP(RichardsTestProblem, Problem, RichardsTestProblem<TypeTag>);
+// Set the fluid system
+SET_PROP(RichardsTestProblem, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::LiquidPhaseTwoC<TypeTag, SimpleH2O<Scalar>, Constant<TypeTag,Scalar>>;
+};
+
// Set the spatial parameters
SET_TYPE_PROP(RichardsTestProblem, SpatialParams, RichardsTestSpatialParams<TypeTag>);
@@ -70,6 +80,9 @@ SET_BOOL_PROP(RichardsTestProblem, VtkAddVelocity, true);
// Set the grid parameter group
SET_STRING_PROP(RichardsTestProblem, GridParameterGroup, "SoilGrid");
+
+// Use mass fractions
+SET_BOOL_PROP(RichardsTestProblem, UseMoles, false);
}
/*!
@@ -132,7 +145,7 @@ public:
: ParentType(timeManager, gridView)
{
name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name) + "-soil";
- contaminantMoleFraction_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ContaminantMoleFraction);
+ contaminantMassFraction_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ContaminantMassFraction);
// for initial conditions
const Scalar sw = 0.2; // start with 20% saturation on top
@@ -225,12 +238,14 @@ public:
// compute correct upwind concentration
if (sourceValues[conti0EqIdx] > 0)
sourceValues[transportEqIdx] = sourceValues[conti0EqIdx]
- * lowDimVolVars.moleFraction(wPhaseIdx, transportCompIdx)
- * lowDimVolVars.molarDensity(wPhaseIdx);
+ * lowDimVolVars.massFraction(wPhaseIdx, transportCompIdx);
else
sourceValues[transportEqIdx] = sourceValues[conti0EqIdx]
- * bulkVolVars.moleFraction(wPhaseIdx, transportCompIdx)
- * bulkVolVars.molarDensity(wPhaseIdx);
+ * bulkVolVars.massFraction(wPhaseIdx, transportCompIdx);
+
+ sourceValues[transportEqIdx] += 2* M_PI *rootRadius * 1.0e-8
+ *(lowDimVolVars.massFraction(wPhaseIdx, transportCompIdx) - bulkVolVars.massFraction(wPhaseIdx, transportCompIdx))
+ *0.5*bulkVolVars.density(wPhaseIdx)*lowDimVolVars.density(wPhaseIdx);
sourceValues *= source.quadratureWeight()*source.integrationElement();
source = sourceValues;
@@ -282,11 +297,7 @@ public:
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
-
- // if (globalPos[dimWorld-1] > this->bBoxMax()[dimWorld-1] - eps_)
- // bcTypes.setAllDirichlet();
- // else
- bcTypes.setAllNeumann();
+ bcTypes.setAllNeumann();
return bcTypes;
}
@@ -327,7 +338,7 @@ public:
contaminationPos += this->bBoxMin();
if ((globalPos - contaminationPos).two_norm() < 0.1*(this->bBoxMax()-this->bBoxMin()).two_norm() + eps_)
- return contaminantMoleFraction_;
+ return contaminantMassFraction_;
else
return 0.0;
}();
@@ -357,7 +368,7 @@ public:
private:
std::string name_;
std::shared_ptr<CouplingManager> couplingManager_;
- Scalar contaminantMoleFraction_;
+ Scalar contaminantMassFraction_;
Scalar pcTop_;
static constexpr Scalar eps_ = 1e-7;
};
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
index 64b0c402a9..3846bbf6e9 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/rootsystemtestproblem.hh
@@ -76,6 +76,9 @@ SET_BOOL_PROP(RootsystemTestProblem, ProblemEnableGravity, true);
// Enable velocity output
SET_BOOL_PROP(RootsystemTestProblem, VtkAddVelocity, true);
+
+// Use mass fractions
+SET_BOOL_PROP(RootsystemTestProblem, UseMoles, false);
}
/*!
@@ -107,7 +110,7 @@ class RootsystemTestProblem : public ImplicitPorousMediaProblem<TypeTag>
transportCompIdx = Indices::transportCompIdx
};
- static const int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static const int wPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
@@ -226,13 +229,12 @@ public:
PrimaryVariables values(0.0);
if (scvf.center()[2] + eps_ > this->bBoxMax()[2])
{
- const auto r = this->spatialParams().radius(this->gridView().indexSet().index(element));
- values[conti0EqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, BoundaryConditions.TranspirationRate)
- /(M_PI*r*r)/scvf.area();
const auto& volVars = elemVolvars[scvf.insideScvIdx()];
+ values[conti0EqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, BoundaryConditions.TranspirationRate)
+ /volVars.extrusionFactor()/scvf.area();
+
values[transportEqIdx] = values[conti0EqIdx]
- * volVars.moleFraction(phaseIdx, transportCompIdx)
- * volVars.molarDensity(phaseIdx);
+ * volVars.massFraction(wPhaseIdx, transportCompIdx);
}
return values;
@@ -308,18 +310,20 @@ public:
// sink defined as radial flow Jr * density [m^2 s-1]* [kg m-3]
PrimaryVariables sourceValues(0.0);
- sourceValues[conti0EqIdx] = 2* M_PI *rootRadius * Kr *(bulkVolVars.pressure(phaseIdx) - lowDimVolVars.pressure(phaseIdx))
- *bulkVolVars.density(phaseIdx);
+ sourceValues[conti0EqIdx] = 2* M_PI *rootRadius * Kr *(bulkVolVars.pressure(wPhaseIdx) - lowDimVolVars.pressure(wPhaseIdx))
+ *bulkVolVars.density(wPhaseIdx);
// compute correct upwind concentration
if (sourceValues[conti0EqIdx] < 0)
sourceValues[transportEqIdx] = sourceValues[conti0EqIdx]
- * lowDimVolVars.moleFraction(phaseIdx, transportCompIdx)
- * lowDimVolVars.molarDensity(phaseIdx);
+ * lowDimVolVars.massFraction(wPhaseIdx, transportCompIdx);
else
sourceValues[transportEqIdx] = sourceValues[conti0EqIdx]
- * bulkVolVars.moleFraction(phaseIdx, transportCompIdx)
- * bulkVolVars.molarDensity(phaseIdx);
+ * bulkVolVars.massFraction(wPhaseIdx, transportCompIdx);
+
+ sourceValues[transportEqIdx] += 2* M_PI *rootRadius * 1.0e-8
+ *(bulkVolVars.massFraction(wPhaseIdx, transportCompIdx) - lowDimVolVars.massFraction(wPhaseIdx, transportCompIdx))
+ *0.5*bulkVolVars.density(wPhaseIdx)*lowDimVolVars.density(wPhaseIdx);
sourceValues *= source.quadratureWeight()*source.integrationElement();
source = sourceValues;
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh b/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
index d96257a60c..1263c1ed09 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
+++ b/test/mixeddimension/embedded/1p2c_richards2c/rositestproblem.hh
@@ -88,6 +88,16 @@ class RosiTestProblem : public MixedDimensionProblem<TypeTag>
using BulkGridView = typename GET_PROP_TYPE(BulkProblemTypeTag, GridView);
using LowDimGridView = typename GET_PROP_TYPE(LowDimProblemTypeTag, GridView);
+ using BulkIndices = typename GET_PROP_TYPE(BulkProblemTypeTag, Indices);
+ using LowDimIndices = typename GET_PROP_TYPE(LowDimProblemTypeTag, Indices);
+
+ enum {
+ // indices of the primary variables
+ wPhaseIdx = BulkIndices::wPhaseIdx,
+ transportCompIdx = BulkIndices::compMainIdx + 1,
+ transportEqIdx = BulkIndices::conti0EqIdx + 1
+ };
+
public:
RosiTestProblem(TimeManager &timeManager, const BulkGridView &bulkGridView, const LowDimGridView &lowDimgridView)
: ParentType(timeManager, bulkGridView, lowDimgridView)
@@ -96,7 +106,105 @@ public:
void preTimeStep()
{
ParentType::preTimeStep();
- std::cout << '\n' << "Simulation time in hours: " << this->timeManager().time()/3600 << '\n';
+ std::cout << '\n' << "\033[1m" << "Simulation time in hours: " << this->timeManager().time()/3600 << "\033[0m" << '\n';
+ }
+
+ void init()
+ {
+ ParentType::init();
+
+ // compute the mass in the entire domain to make sure the tracer is conserved
+ Scalar mass = 0.0;
+
+ // low dim elements
+ for (const auto& element : elements(this->lowDimProblem().gridView()))
+ {
+ auto fvGeometry = localView(this->lowDimProblem().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->lowDimProblem().model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->lowDimProblem().model().curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ mass += volVars.massFraction(wPhaseIdx, transportCompIdx)*volVars.density(wPhaseIdx)
+ *scv.volume() * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+
+ // bulk elements
+ for (const auto& element : elements(this->bulkProblem().gridView()))
+ {
+ auto fvGeometry = localView(this->bulkProblem().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->bulkProblem().model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->bulkProblem().model().curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ mass += volVars.massFraction(wPhaseIdx, transportCompIdx)*volVars.density(wPhaseIdx)
+ * scv.volume() * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+
+ std::cout << "\033[1;36m" << "The domain initially contains " << mass*1e9 << " ug tracer.\033[0m" << '\n';
+ }
+
+ void postTimeStep()
+ {
+ ParentType::postTimeStep();
+
+ // compute the mass in the entire domain to make sure the tracer is conserved
+ Scalar mass = 0.0;
+ Scalar boundaryMass = 0.0;
+
+ // low dim elements
+ for (const auto& element : elements(this->lowDimProblem().gridView()))
+ {
+ auto fvGeometry = localView(this->lowDimProblem().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->lowDimProblem().model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->lowDimProblem().model().curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ mass += volVars.massFraction(wPhaseIdx, transportCompIdx)*volVars.density(wPhaseIdx)
+ *scv.volume() * volVars.porosity() * volVars.extrusionFactor();
+ }
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ if (scvf.boundary())
+ boundaryMass += this->lowDimProblem().neumann(element, fvGeometry, elemVolVars, scvf)[transportEqIdx]
+ * scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor()
+ * this->timeManager().timeStepSize();
+
+ }
+
+ // bulk elements
+ for (const auto& element : elements(this->bulkProblem().gridView()))
+ {
+ auto fvGeometry = localView(this->bulkProblem().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(this->bulkProblem().model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->bulkProblem().model().curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ mass += volVars.massFraction(wPhaseIdx, transportCompIdx)*volVars.density(wPhaseIdx)
+ *scv.volume() * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+
+ std::cout << "\033[1;36m" << "The domain contains " << mass*1e9 << " ug, "
+ << -boundaryMass*1e9 << " (over the boundary, ug) = "
+ << (mass - boundaryMass)*1e9 << " ug balanced.\033[0m" << '\n';
}
};
diff --git a/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input b/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
index 7c2bbf4ad2..9db0aeb9bd 100644
--- a/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
+++ b/test/mixeddimension/embedded/1p2c_richards2c/test_rosi2c.input
@@ -18,7 +18,7 @@ Cells = 20 20 20
[Problem]
Name = rosi2c
-ContaminantMoleFraction = 1e-9
+ContaminantMassFraction = 1e-6
[Component]
LiquidDiffusionCoefficient = 1e-9
--
GitLab