Skip to content

GitLab

Commit b3f68777 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'cleanup/free-biominreactions-typetag' into 'master' 

Cleanup/free biominreactions typetag

Closes #7 and #9

See merge request !57
parents 2a81724d 09ad102d
Pipeline #755 passed with stage
in 30 seconds
Hide whitespace changes
Inline Side-by-side
exercises/exercise-biomineralization/CMakeLists.txt
View file @ b3f68777
# executables for exercisebiomin
dune_add_test(NAME exercisebiomin
SOURCES exercisebiomin.cc)
SOURCES exercisebiomin.cc
CMD_ARGS -TimeLoop.TEnd 1e5)
# add tutorial to the common target
add_dependencies(test_exercises exercisebiomin)
......
exercises/exercise-biomineralization/README.md
View file @ b3f68777
......@@ -75,7 +75,7 @@ To enable your newly created chemical equation, the chemistry file has to be inc
Additionally the TypeTag of your chemistry file needs to be set in the problem file (line 125):
```c++
using Chemistry = typename Dumux::SimpleBiominReactions<TypeTag>;
using Chemistry = typename Dumux::SimpleBiominReactions<NumEqVector, VolumeVariables>;
```
__Task__
......
exercises/exercise-biomineralization/biominproblem.hh
View file @ b3f68777
......@@ -50,8 +50,7 @@ namespace Dumux {
template <class TypeTag>
class ExerciseFourBioMinProblem;
namespace Properties
{
namespace Properties {
//! Create new type tag for the problem
// Create new type tags
namespace TTag {
......@@ -96,15 +95,14 @@ 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
/*!
*
* \brief Problem biomineralization (MICP) in an experimental setup.
*/
template <class TypeTag>
......@@ -171,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");
......
exercises/exercise-biomineralization/biominspatialparams.hh
View file @ b3f68777
......@@ -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;}
......
exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh
View file @ b3f68777
......@@ -30,15 +30,12 @@ namespace Dumux {
* \ingroup Chemistry
* \brief The source and sink terms due to reactions are calculated in this class.
*/
template <class TypeTag>
template< class NumEqVector, class VolumeVariables >
class SimpleBiominReactions
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using ThisType = SimpleBiominReactions<TypeTag>;
using Sources = GetPropType<TypeTag, Properties::NumEqVector>;
using SolidSystem = typename VolumeVariables::SolidSystem;
using FluidSystem = typename VolumeVariables::FluidSystem;
using Scalar = typename FluidSystem::Scalar;
public:
......@@ -81,7 +78,7 @@ public:
* \param Source The source
* \param volVars The volume variables
*/
void reactionSource(Sources &q, const VolumeVariables &volVars)
void reactionSource(NumEqVector &q, const VolumeVariables &volVars)
{
// define and compute some parameters for convenience:
Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
......@@ -91,36 +88,32 @@ public:
if (volFracBiofilm < 0)
volFracBiofilm = 0;
// TODO: dumux-course-task
// implement mass of biofilm
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:
// urease parameters
Scalar kub_;
......
exercises/exercise-biomineralization/exercisebiomin.input
View file @ b3f68777
......@@ -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
......
exercises/solution/exercise-biomineralization/CMakeLists.txt
View file @ b3f68777
# executables for exercisebiomin
dune_add_test(NAME exercisebiomin_solution
SOURCES exercisebiomin.cc)
SOURCES exercisebiomin.cc
CMD_ARGS -TimeLoop.TEnd 1e5)
# add tutorial to the common target
add_dependencies(test_exercises exercisebiomin_solution)
......
exercises/solution/exercise-biomineralization/biominproblem.hh
View file @ b3f68777
......@@ -36,7 +36,7 @@
#include "fluidsystems/biomin.hh" // The biomineralization fluid system
// TODO: dumux-course-task
// include chemistry file here
#include "chemistry/simplebiominreactions.hh" // chemical reactions
#include "chemistry/simplebiominreactions.hh" // chemical reactions
#include "biominspatialparams.hh" // Spatially dependent parameters
namespace Dumux {
......@@ -51,8 +51,7 @@ namespace Dumux {
template <class TypeTag>
class ExerciseFourBioMinProblem;
namespace Properties
{
namespace Properties {
//! Create new type tag for the problem
// Create new type tags
namespace TTag {
......@@ -97,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
/*!
......@@ -139,7 +138,7 @@ class ExerciseFourBioMinProblem : public PorousMediumFlowProblem<TypeTag>
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
// TODO: dumux-course-task
// set the chemistry TypeTag
using Chemistry = typename Dumux::SimpleBiominReactions<TypeTag>;
using Chemistry = typename Dumux::SimpleBiominReactions<NumEqVector, VolumeVariables>;
enum
{
......@@ -173,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");
......@@ -381,7 +377,6 @@ public:
return source;
}
const std::vector<Scalar>& getKxx()
{
return Kxx_;
......
exercises/solution/exercise-biomineralization/biominspatialparams.hh
View file @ b3f68777
......@@ -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;}
......
exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh
View file @ b3f68777
......@@ -30,15 +30,12 @@ namespace Dumux {
* \ingroup Chemistry
* \brief The source and sink terms due to reactions are calculated in this class.
*/
template <class TypeTag>
template< class NumEqVector, class VolumeVariables >
class SimpleBiominReactions
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using ThisType = SimpleBiominReactions<TypeTag>;
using Sources = GetPropType<TypeTag, Properties::NumEqVector>;
using SolidSystem = typename VolumeVariables::SolidSystem;
using FluidSystem = typename VolumeVariables::FluidSystem;
using Scalar = typename FluidSystem::Scalar;
public:
......@@ -81,7 +78,7 @@ public:
* \param Source The source
* \param volVars The volume variables
*/
void reactionSource(Sources &q, const VolumeVariables &volVars)
void reactionSource(NumEqVector &q, const VolumeVariables &volVars)
{
// define and compute some parameters for convenience:
Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
......@@ -91,7 +88,7 @@ public:
if (volFracBiofilm < 0)
volFracBiofilm = 0;
// TODO: dumux-course-task
// implement mass of biofilm
Scalar massBiofilm = densityBiofilm * volFracBiofilm;
......@@ -99,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
......@@ -108,23 +103,19 @@ 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;
}
private:
// urease parameters
......
exercises/solution/exercise-biomineralization/exercisebiomin_solution.input
View file @ b3f68777
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment