From 09ad102d408790421492905009a637ae0148ac47 Mon Sep 17 00:00:00 2001
From: Simon Scholz <simon.scholz@iws.uni-stuttgart.de>
Date: Fri, 21 Dec 2018 11:50:11 +0100
Subject: [PATCH] [biomin] update biomin example, free from wrong perm matrix
---
.../biominproblem.hh | 9 ++---
.../biominspatialparams.hh | 39 +++----------------
.../chemistry/simplebiominreactions.hh | 21 ++++------
.../exercisebiomin.input | 3 --
.../biominproblem.hh | 9 ++---
.../biominspatialparams.hh | 39 +++----------------
.../chemistry/simplebiominreactions.hh | 13 ++-----
.../exercisebiomin_solution.input | 3 --
8 files changed, 30 insertions(+), 106 deletions(-)
diff --git a/exercises/exercise-biomineralization/biominproblem.hh b/exercises/exercise-biomineralization/biominproblem.hh
index c72dff46..20d0790b 100644
--- a/exercises/exercise-biomineralization/biominproblem.hh
+++ b/exercises/exercise-biomineralization/biominproblem.hh
@@ -95,11 +95,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> {
};
template<class TypeTag>
-struct EnableFVGridGeometryCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableFVGridGeometryCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
} // end namespace properties
/*!
@@ -169,9 +169,6 @@ public:
ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
-#if !DUNE_VERSION_NEWER(DUNE_COMMON,2,7)
- Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
-#endif
name_ = getParam<std::string>("Problem.Name");
//initial values
densityW_ = getParam<Scalar>("Initial.InitDensityW");
diff --git a/exercises/exercise-biomineralization/biominspatialparams.hh b/exercises/exercise-biomineralization/biominspatialparams.hh
index 393bd57d..3815f373 100644
--- a/exercises/exercise-biomineralization/biominspatialparams.hh
+++ b/exercises/exercise-biomineralization/biominspatialparams.hh
@@ -229,10 +229,14 @@ public:
const auto& dofPosition = scv.dofPosition();
const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition);
auto kInit = isInAquitardNotFaultZone ? aquitardPermeability_ : initialPermeability_;
+ PermeabilityType K (0.0);
+ for (int i = 0; i < dimWorld; ++i)
+ K[i][i] = kInit;
const auto refPoro = referencePorosity(element, scv);
const auto poro = porosity(element, scv, elemSol);
- auto kEvaluated = permLaw_.evaluatePermeability(kInit, refPoro, poro);
- referencePermeability_[eIdx][scv.indexInElement()] = calculateKRef(kInit, kEvaluated);
+ auto kEvaluated = permLaw_.evaluatePermeability(K, refPoro, poro);
+ K *= K[0][0]/kEvaluated[0][0];
+ referencePermeability_[eIdx][scv.indexInElement()] = K;
}
}
}
@@ -275,39 +279,8 @@ public:
{ return FluidSystem::liquidPhaseIdx; }
private:
-
static constexpr Scalar eps_ = 1e-6;
- //! calculateKRef for tensorial permeabilities
- PermeabilityType calculateKRef(PermeabilityType K, PermeabilityType kEvaluated)
- {
- //TODO this assumes the change in permeability is homogeneous as currently also assumed in PermeabilityKozenyCarman.
- // there also might be nicer/safer ways to calculate this.
- Scalar factor = 0.0;
- int counter = 0;
- for (int i = 0; i < dimWorld; ++i)
- for (int j = 0; i < dimWorld; ++i)
- if(kEvaluated[i][j]!=0)
- {
- factor += K[i][j] / kEvaluated[i][j];
- ++counter;
- }
- factor /= counter;
- K *= factor;
-
- return K;
- }
-
- //! calculateKRef for scalar permeabilities
- template<class PT = PermeabilityType,
- std::enable_if_t<Dune::IsNumber< PT >::value, int> = 0>
- PermeabilityType calculateKRef(PermeabilityType K, PermeabilityType kEvaluated)
- {
- Scalar factor = K / kEvaluated;
- K *= factor;
- return K;
- }
-
Scalar calculatephiRef(Scalar phiInit, Scalar phiEvaluated)
{ return 2*phiInit - phiEvaluated;}
diff --git a/exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh b/exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh
index 8a18560d..f7a0f62a 100644
--- a/exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh
+++ b/exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh
@@ -95,28 +95,23 @@ public:
Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
xwCa,
volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O]
- if (molalityUrea < 0)
- molalityUrea = 0;
// TODO: dumux-course-task
// compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
// compute/set dissolution and precipitation rate of calcite
- Scalar rprec = 0;
-
- // regularize so that we do not get precipitation where there is no calcium
- if (xwCa > 1e-8)
- rprec = rurea;
+ Scalar rprec = 0.0;
+ rprec = rurea;
// set source terms
// TODO: dumux-course-task
// update terms according to stochiometry
- q[H2OIdx] += 0;
- q[CO2Idx] += 0;
- q[CaIdx] += 0;
- q[UreaIdx] += 0;
- q[BiofilmIdx] += 0;
- q[CalciteIdx] += 0;
+ q[H2OIdx] += 0.0;
+ q[CO2Idx] += 0.0;
+ q[CaIdx] += 0.0;
+ q[UreaIdx] += 0.0;
+ q[BiofilmIdx] += 0.0;
+ q[CalciteIdx] += 0.0;
}
private:
diff --git a/exercises/exercise-biomineralization/exercisebiomin.input b/exercises/exercise-biomineralization/exercisebiomin.input
index 0f3d3e41..fb2df919 100644
--- a/exercises/exercise-biomineralization/exercisebiomin.input
+++ b/exercises/exercise-biomineralization/exercisebiomin.input
@@ -11,9 +11,6 @@ Name = biomin
UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m]
Cells = 20 15 # x-/y-resolution of the grid
-[Newton]
-MaxRelativeShift = 1e-6 #
-
[Initial]
InitDensityW = 997 # [kg/m³] initial wetting density
InitPressure = 1e5 # [Pa] initial pressure
diff --git a/exercises/solution/exercise-biomineralization/biominproblem.hh b/exercises/solution/exercise-biomineralization/biominproblem.hh
index 2c41e437..756b7537 100644
--- a/exercises/solution/exercise-biomineralization/biominproblem.hh
+++ b/exercises/solution/exercise-biomineralization/biominproblem.hh
@@ -96,11 +96,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> {
};
template<class TypeTag>
-struct EnableFVGridGeometryCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableFVGridGeometryCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = false; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::ExerciseFourBioMin> { static constexpr bool value = true; };
} // end namespace properties
/*!
@@ -172,9 +172,6 @@ public:
ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
-#if !DUNE_VERSION_NEWER(DUNE_COMMON,2,7)
- Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
-#endif
name_ = getParam<std::string>("Problem.Name");
//initial values
densityW_ = getParam<Scalar>("Initial.InitDensityW");
diff --git a/exercises/solution/exercise-biomineralization/biominspatialparams.hh b/exercises/solution/exercise-biomineralization/biominspatialparams.hh
index 393bd57d..3815f373 100644
--- a/exercises/solution/exercise-biomineralization/biominspatialparams.hh
+++ b/exercises/solution/exercise-biomineralization/biominspatialparams.hh
@@ -229,10 +229,14 @@ public:
const auto& dofPosition = scv.dofPosition();
const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition);
auto kInit = isInAquitardNotFaultZone ? aquitardPermeability_ : initialPermeability_;
+ PermeabilityType K (0.0);
+ for (int i = 0; i < dimWorld; ++i)
+ K[i][i] = kInit;
const auto refPoro = referencePorosity(element, scv);
const auto poro = porosity(element, scv, elemSol);
- auto kEvaluated = permLaw_.evaluatePermeability(kInit, refPoro, poro);
- referencePermeability_[eIdx][scv.indexInElement()] = calculateKRef(kInit, kEvaluated);
+ auto kEvaluated = permLaw_.evaluatePermeability(K, refPoro, poro);
+ K *= K[0][0]/kEvaluated[0][0];
+ referencePermeability_[eIdx][scv.indexInElement()] = K;
}
}
}
@@ -275,39 +279,8 @@ public:
{ return FluidSystem::liquidPhaseIdx; }
private:
-
static constexpr Scalar eps_ = 1e-6;
- //! calculateKRef for tensorial permeabilities
- PermeabilityType calculateKRef(PermeabilityType K, PermeabilityType kEvaluated)
- {
- //TODO this assumes the change in permeability is homogeneous as currently also assumed in PermeabilityKozenyCarman.
- // there also might be nicer/safer ways to calculate this.
- Scalar factor = 0.0;
- int counter = 0;
- for (int i = 0; i < dimWorld; ++i)
- for (int j = 0; i < dimWorld; ++i)
- if(kEvaluated[i][j]!=0)
- {
- factor += K[i][j] / kEvaluated[i][j];
- ++counter;
- }
- factor /= counter;
- K *= factor;
-
- return K;
- }
-
- //! calculateKRef for scalar permeabilities
- template<class PT = PermeabilityType,
- std::enable_if_t<Dune::IsNumber< PT >::value, int> = 0>
- PermeabilityType calculateKRef(PermeabilityType K, PermeabilityType kEvaluated)
- {
- Scalar factor = K / kEvaluated;
- K *= factor;
- return K;
- }
-
Scalar calculatephiRef(Scalar phiInit, Scalar phiEvaluated)
{ return 2*phiInit - phiEvaluated;}
diff --git a/exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh b/exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh
index 66331ec4..250333fa 100644
--- a/exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh
+++ b/exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh
@@ -96,8 +96,6 @@ public:
Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
xwCa,
volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O]
- if (molalityUrea < 0)
- molalityUrea = 0;
// TODO: dumux-course-task
// compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
@@ -105,20 +103,17 @@ public:
Scalar rurea = kurease_ * Zub * molalityUrea / (ku_ + molalityUrea); // [mol/m³s]
// compute/set dissolution and precipitation rate of calcite
- Scalar rprec = 0;
-
- // regularize so that we do not get precipitation where there is no calcium
- if (xwCa > 1e-8)
- rprec = rurea;
+ Scalar rprec = 0.0;
+ rprec = rurea;
// set source terms
// TODO: dumux-course-task
// update terms according to stochiometry
- q[H2OIdx] += 0;
+ q[H2OIdx] += 0.0;
q[CO2Idx] += rurea - rprec ;
q[CaIdx] += - rprec;
q[UreaIdx] += - rurea;
- q[BiofilmIdx] += 0;
+ q[BiofilmIdx] += 0.0;
q[CalciteIdx] += rprec;
}
diff --git a/exercises/solution/exercise-biomineralization/exercisebiomin_solution.input b/exercises/solution/exercise-biomineralization/exercisebiomin_solution.input
index 0f3d3e41..fb2df919 100644
--- a/exercises/solution/exercise-biomineralization/exercisebiomin_solution.input
+++ b/exercises/solution/exercise-biomineralization/exercisebiomin_solution.input
@@ -11,9 +11,6 @@ Name = biomin
UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m]
Cells = 20 15 # x-/y-resolution of the grid
-[Newton]
-MaxRelativeShift = 1e-6 #
-
[Initial]
InitDensityW = 997 # [kg/m³] initial wetting density
InitPressure = 1e5 # [Pa] initial pressure
--
GitLab