Commit 09ad102d authored by Simon Scholz's avatar Simon Scholz Committed by Timo Koch
Browse files

[biomin] update biomin example, free from wrong perm matrix

parent 02ae47b3
...@@ -95,11 +95,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> { ...@@ -95,11 +95,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> {
}; };
template<class TypeTag> 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> 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> 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 } // end namespace properties
/*! /*!
...@@ -169,9 +169,6 @@ public: ...@@ -169,9 +169,6 @@ public:
ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(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"); name_ = getParam<std::string>("Problem.Name");
//initial values //initial values
densityW_ = getParam<Scalar>("Initial.InitDensityW"); densityW_ = getParam<Scalar>("Initial.InitDensityW");
......
...@@ -229,10 +229,14 @@ public: ...@@ -229,10 +229,14 @@ public:
const auto& dofPosition = scv.dofPosition(); const auto& dofPosition = scv.dofPosition();
const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition); const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition);
auto kInit = isInAquitardNotFaultZone ? aquitardPermeability_ : initialPermeability_; 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 refPoro = referencePorosity(element, scv);
const auto poro = porosity(element, scv, elemSol); const auto poro = porosity(element, scv, elemSol);
auto kEvaluated = permLaw_.evaluatePermeability(kInit, refPoro, poro); auto kEvaluated = permLaw_.evaluatePermeability(K, refPoro, poro);
referencePermeability_[eIdx][scv.indexInElement()] = calculateKRef(kInit, kEvaluated); K *= K[0][0]/kEvaluated[0][0];
referencePermeability_[eIdx][scv.indexInElement()] = K;
} }
} }
} }
...@@ -275,39 +279,8 @@ public: ...@@ -275,39 +279,8 @@ public:
{ return FluidSystem::liquidPhaseIdx; } { return FluidSystem::liquidPhaseIdx; }
private: private:
static constexpr Scalar eps_ = 1e-6; 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) Scalar calculatephiRef(Scalar phiInit, Scalar phiEvaluated)
{ return 2*phiInit - phiEvaluated;} { return 2*phiInit - phiEvaluated;}
......
...@@ -95,28 +95,23 @@ public: ...@@ -95,28 +95,23 @@ public:
Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx), Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
xwCa, xwCa,
volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O] volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O]
if (molalityUrea < 0)
molalityUrea = 0;
// TODO: dumux-course-task // TODO: dumux-course-task
// compute rate of ureolysis by implementing Z_urease,biofilm and r_urea // compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
// compute/set dissolution and precipitation rate of calcite // compute/set dissolution and precipitation rate of calcite
Scalar rprec = 0; Scalar rprec = 0.0;
rprec = rurea;
// regularize so that we do not get precipitation where there is no calcium
if (xwCa > 1e-8)
rprec = rurea;
// set source terms // set source terms
// TODO: dumux-course-task // TODO: dumux-course-task
// update terms according to stochiometry // update terms according to stochiometry
q[H2OIdx] += 0; q[H2OIdx] += 0.0;
q[CO2Idx] += 0; q[CO2Idx] += 0.0;
q[CaIdx] += 0; q[CaIdx] += 0.0;
q[UreaIdx] += 0; q[UreaIdx] += 0.0;
q[BiofilmIdx] += 0; q[BiofilmIdx] += 0.0;
q[CalciteIdx] += 0; q[CalciteIdx] += 0.0;
} }
private: private:
......
...@@ -11,9 +11,6 @@ Name = biomin ...@@ -11,9 +11,6 @@ Name = biomin
UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m] UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m]
Cells = 20 15 # x-/y-resolution of the grid Cells = 20 15 # x-/y-resolution of the grid
[Newton]
MaxRelativeShift = 1e-6 #
[Initial] [Initial]
InitDensityW = 997 # [kg/m³] initial wetting density InitDensityW = 997 # [kg/m³] initial wetting density
InitPressure = 1e5 # [Pa] initial pressure InitPressure = 1e5 # [Pa] initial pressure
......
...@@ -96,11 +96,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> { ...@@ -96,11 +96,11 @@ struct SpatialParams<TypeTag, TTag::ExerciseFourBioMin> {
}; };
template<class TypeTag> 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> 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> 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 } // end namespace properties
/*! /*!
...@@ -172,9 +172,6 @@ public: ...@@ -172,9 +172,6 @@ public:
ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(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"); name_ = getParam<std::string>("Problem.Name");
//initial values //initial values
densityW_ = getParam<Scalar>("Initial.InitDensityW"); densityW_ = getParam<Scalar>("Initial.InitDensityW");
......
...@@ -229,10 +229,14 @@ public: ...@@ -229,10 +229,14 @@ public:
const auto& dofPosition = scv.dofPosition(); const auto& dofPosition = scv.dofPosition();
const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition); const bool isInAquitardNotFaultZone = isInAquitard_(dofPosition) && !isFaultZone_(dofPosition);
auto kInit = isInAquitardNotFaultZone ? aquitardPermeability_ : initialPermeability_; 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 refPoro = referencePorosity(element, scv);
const auto poro = porosity(element, scv, elemSol); const auto poro = porosity(element, scv, elemSol);
auto kEvaluated = permLaw_.evaluatePermeability(kInit, refPoro, poro); auto kEvaluated = permLaw_.evaluatePermeability(K, refPoro, poro);
referencePermeability_[eIdx][scv.indexInElement()] = calculateKRef(kInit, kEvaluated); K *= K[0][0]/kEvaluated[0][0];
referencePermeability_[eIdx][scv.indexInElement()] = K;
} }
} }
} }
...@@ -275,39 +279,8 @@ public: ...@@ -275,39 +279,8 @@ public:
{ return FluidSystem::liquidPhaseIdx; } { return FluidSystem::liquidPhaseIdx; }
private: private:
static constexpr Scalar eps_ = 1e-6; 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) Scalar calculatephiRef(Scalar phiInit, Scalar phiEvaluated)
{ return 2*phiInit - phiEvaluated;} { return 2*phiInit - phiEvaluated;}
......
...@@ -96,8 +96,6 @@ public: ...@@ -96,8 +96,6 @@ public:
Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx), Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
xwCa, xwCa,
volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O] volVars.moleFraction(liquidPhaseIdx,CO2Idx)); // [mol_urea/kg_H2O]
if (molalityUrea < 0)
molalityUrea = 0;
// TODO: dumux-course-task // TODO: dumux-course-task
// compute rate of ureolysis by implementing Z_urease,biofilm and r_urea // compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
...@@ -105,20 +103,17 @@ public: ...@@ -105,20 +103,17 @@ public:
Scalar rurea = kurease_ * Zub * molalityUrea / (ku_ + molalityUrea); // [mol/m³s] Scalar rurea = kurease_ * Zub * molalityUrea / (ku_ + molalityUrea); // [mol/m³s]
// compute/set dissolution and precipitation rate of calcite // compute/set dissolution and precipitation rate of calcite
Scalar rprec = 0; Scalar rprec = 0.0;
rprec = rurea;
// regularize so that we do not get precipitation where there is no calcium
if (xwCa > 1e-8)
rprec = rurea;
// set source terms // set source terms
// TODO: dumux-course-task // TODO: dumux-course-task
// update terms according to stochiometry // update terms according to stochiometry
q[H2OIdx] += 0; q[H2OIdx] += 0.0;
q[CO2Idx] += rurea - rprec ; q[CO2Idx] += rurea - rprec ;
q[CaIdx] += - rprec; q[CaIdx] += - rprec;
q[UreaIdx] += - rurea; q[UreaIdx] += - rurea;
q[BiofilmIdx] += 0; q[BiofilmIdx] += 0.0;
q[CalciteIdx] += rprec; q[CalciteIdx] += rprec;
} }
......
...@@ -11,9 +11,6 @@ Name = biomin ...@@ -11,9 +11,6 @@ Name = biomin
UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m] UpperRight = 20 15 # x-/y-coordinates of the upper-right corner of the grid [m]
Cells = 20 15 # x-/y-resolution of the grid Cells = 20 15 # x-/y-resolution of the grid
[Newton]
MaxRelativeShift = 1e-6 #
[Initial] [Initial]
InitDensityW = 997 # [kg/m³] initial wetting density InitDensityW = 997 # [kg/m³] initial wetting density
InitPressure = 1e5 # [Pa] initial pressure 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