Commit 6b5aff33 authored by Timo Koch's avatar Timo Koch
Browse files

[2pncmin][cleanup] Beautifications, remove white space, unused typedefs, members, functions.

Introduce epsilon for floating point comparison in test.
parent 04324a52
......@@ -129,7 +129,7 @@ protected:
{
// calculate the mean intrinsic permeability
const SpatialParams &spatialParams = problem.spatialParams();
DimWorldMatrix K;
DimWorldMatrix K(0.0);
const VolumeVariables &volVarsI = elemVolVars[this->face().i];
const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
......
......@@ -35,7 +35,7 @@ namespace Dumux
* \ingroup TwoPNCMinModel
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the two-phase n-component mineralization fully implicit box model.
* using the two-phase n-component mineralization fully implicit model.
*
* This class is used to fill the gaps in ImplicitLocalResidual for the two-phase n-component flow.
*/
......@@ -93,17 +93,6 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
public:
/*!
* \brief Constructor. Sets the upwind weight.
*/
TwoPNCMinLocalResidual()
{
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
// it by the run-time parameter from the Dune::ParameterTree
this->massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
};
/*!
* \brief Evaluate the amount all conservation quantities
* (e.g. phase mass) within a sub-control volume.
......@@ -116,25 +105,25 @@ public:
* \param scvIdx The SCV (sub-control-volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
{
//call parenttype function
ParentType::computeStorage(storage, scvIdx, usePrevSol);
void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
{
//call parenttype function
ParentType::computeStorage(storage, scvIdx, usePrevSol);
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
: this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
const auto& elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
// Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
for (int phaseIdx = numPhases; phaseIdx < numPhases + numSPhases; ++phaseIdx)
{
int eqIdx = conti0EqIdx + numComponents-numPhases + phaseIdx;
storage[eqIdx] += volVars.precipitateVolumeFraction(phaseIdx)*volVars.molarDensity(phaseIdx);
}
// Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
for (int phaseIdx = numPhases; phaseIdx < numPhases + numSPhases; ++phaseIdx)
{
auto eqIdx = conti0EqIdx + numComponents-numPhases + phaseIdx;
storage[eqIdx] += volVars.precipitateVolumeFraction(phaseIdx)*volVars.molarDensity(phaseIdx);
}
Valgrind::CheckDefined(storage);
}
Valgrind::CheckDefined(storage);
}
};
} // end namespace
#endif
......@@ -393,10 +393,9 @@ public:
for (unsigned i = 0; i < this->staticDat_.size(); ++i)
this->staticDat_[i].visited = false;
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
for (const auto& element : elements(this->gridView_()))
{
FVElementGeometry fvGeometry;
fvGeometry.update(this->gridView_(), element);
for (int i = 0; i < fvGeometry.numScv; ++i)
{
......@@ -406,18 +405,20 @@ public:
continue;
this->staticDat_[globalIdx].visited = true;
VolumeVariables volVars;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
element,
fvGeometry,
i,
false);
const GlobalPosition &global = element.geometry().corner(i);
auto global = element.geometry().corner(i);
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
global))
{ wasSwitched = true;
{
wasSwitched = true;
}
}
}
......@@ -426,7 +427,7 @@ public:
// other partition we will also set the switch flag
// for our partition.
if (this->gridView_().comm().size() > 1)
wasSwitched = this->gridView_().comm().max(wasSwitched);
wasSwitched = this->gridView_().comm().max(wasSwitched);
setSwitched_(wasSwitched);
}
......@@ -488,39 +489,38 @@ protected:
newPhasePresence = wPhaseOnly;
//switch "Sl" to "xlN2"
globalSol[globalIdx][switchIdx]
= volVars.moleFraction(wPhaseIdx, nCompIdx /*N2*/);
globalSol[globalIdx][switchIdx] = volVars.moleFraction(wPhaseIdx, nCompIdx /*N2*/);
}
}
else if (phasePresence == nPhaseOnly)
{
Scalar sumxl = 0;
//Calculate sum of mole fractions (water and air) in the hypothetical liquid phase
for (int compIdx = 0; compIdx < numComponents; compIdx++)
Scalar sumxl = 0;
//Calculate sum of mole fractions (water and air) in the hypothetical liquid phase
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
sumxl += volVars.moleFraction(wPhaseIdx, compIdx);
}
Scalar xlmax = 1.0;
if (sumxl > xlmax)
wouldSwitch = true;
if (this->staticDat_[globalIdx].wasSwitched)
xlmax *=1.02;
//if the sum of the mole fractions would be larger than
//1, wetting phase appears
if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
{
// liquid phase appears
std::cout << "Liquid Phase appears at vertex " << globalIdx
<< ", coordinated: " << globalPos << ", sumxl: "
<< sumxl << std::endl;
newPhasePresence = bothPhases;
if (formulation == pgSl)
globalSol[globalIdx][switchIdx] = 0.0;
else if (formulation == plSg)
globalSol[globalIdx][switchIdx] = 1.0;
Scalar xlmax = 1.0;
if (sumxl > xlmax)
wouldSwitch = true;
if (this->staticDat_[globalIdx].wasSwitched)
xlmax *=1.02;
//if the sum of the mole fractions would be larger than
//1, wetting phase appears
if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
{
// liquid phase appears
std::cout << "Liquid Phase appears at vertex " << globalIdx
<< ", coordinated: " << globalPos << ", sumxl: "
<< sumxl << std::endl;
newPhasePresence = bothPhases;
if (formulation == pgSl)
globalSol[globalIdx][switchIdx] = 0.0;
else if (formulation == plSg)
globalSol[globalIdx][switchIdx] = 1.0;
//Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
}
}
}
else if (phasePresence == wPhaseOnly)
{
......@@ -539,8 +539,8 @@ protected:
if (sumxg > xgmax)
{
std::cout << "Gas Phase appears at vertex " << globalIdx
<< ", coordinated: " << globalPos << ", sumxg: "
<< sumxg << std::endl;
<< ", coordinated: " << globalPos << ", sumxg: "
<< sumxg << std::endl;
newPhasePresence = bothPhases;
//saturation of the liquid phase set to 0.9999 (if formulation pgSl and vice versa)
if (formulation == pgSl)
......
......@@ -115,7 +115,7 @@ class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
conti0EqIdx = Indices::conti0EqIdx,
contiTotalMassIdx = conti0EqIdx + FluidSystem::AirIdx,
precipNaClEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents,
contiWEqIdx = conti0EqIdx + FluidSystem::H2OIdx,
contiWEqIdx = conti0EqIdx + FluidSystem::H2OIdx,
// Phase State
wPhaseOnly = Indices::wPhaseOnly,
......@@ -133,7 +133,6 @@ class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::Intersection Intersection;
......@@ -141,17 +140,16 @@ class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
DissolutionProblem(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, GridCreator::grid().leafGridView())
DissolutionProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
outerSalinity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterSalinity);
temperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, Temperature);
reservoirPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ReservoirPressure);
initLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, LiquidSaturation);
initPrecipitatedSalt1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt1);
initPrecipitatedSalt2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt2);
initPrecipitatedSalt1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt1);
initPrecipitatedSalt2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt2);
outerLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterLiqSaturation);
innerLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerLiqSaturation);
......@@ -168,8 +166,8 @@ public:
temperatureHigh_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureHigh);
name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, OutputName);
freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
storageLastTimestep_ = Scalar(0);
lastMassOutputTime_ = Scalar(0);
storageLastTimestep_ = 0.0;
lastMassOutputTime_ = 0.0;
outfile.open("evaporation.out");
outfile << "time; evaporationRate" << std::endl;
......@@ -194,10 +192,9 @@ public:
bool shouldWriteOutput() const
{
return
this->timeManager().timeStepIndex() % 1 == 0 ||
this->timeManager().episodeWillBeOver() ||
this->timeManager().willBeFinished();
return this->timeManager().timeStepIndex() % 1 == 0 ||
this->timeManager().episodeWillBeOver() ||
this->timeManager().willBeFinished();
}
......@@ -233,20 +230,19 @@ public:
*/
void boundaryTypes(BoundaryTypes &bcTypes, const Vertex &vertex) const
{
const GlobalPosition &globalPos = vertex.geometry().center();
auto globalPos = vertex.geometry().center();
Scalar rmax = this->bBoxMax()[0]; // outerRadius_;
Scalar rmin = this->bBoxMin()[0]; // well radius equal to the first value of the dgf grid file
const Scalar rmax = this->bBoxMax()[0]; // outerRadius_;
const Scalar rmin = this->bBoxMin()[0]; // well radius equal to the first value of the dgf grid file
//Neumann
// default to Neumann
bcTypes.setAllNeumann();
//Constant pressure at reservoir boundary (Dirichlet condition)
// Constant pressure at reservoir boundary (Dirichlet condition)
if(globalPos[0] > rmax - eps_)
bcTypes.setAllDirichlet();
//Constant pressure at well (Dirichlet condition)
// Constant pressure at well (Dirichlet condition)
if(globalPos[0] < rmin + eps_)
bcTypes.setAllDirichlet();
}
......@@ -259,10 +255,10 @@ public:
*/
void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
const GlobalPosition &globalPos = vertex.geometry().center();
auto globalPos = vertex.geometry().center();
Scalar rmax = this->bBoxMax()[0];
Scalar rmin = this->bBoxMin()[0];
const Scalar rmax = this->bBoxMax()[0];
const Scalar rmin = this->bBoxMin()[0];
if(globalPos[0] > rmax - eps_)
{
......@@ -291,11 +287,11 @@ public:
* influx.
*/
void neumann(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
{
values = 0.0;
}
......@@ -307,16 +303,16 @@ public:
* variables.
*/
void initial(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
auto globalPos = element.geometry().corner(scvIdx);
values[pressureIdx] = reservoirPressure_;
values[switchIdx] = initLiqSaturation_; // Sl primary variable
values[xlNaClIdx] = massTomoleFrac_(outerSalinity_); // mole fraction
if(globalPos[0]>5.0 && globalPos[0]< 20.0)
if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 20.0 - eps_)
values[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
else
values[precipNaClIdx] = initPrecipitatedSalt1_; // [kg/m^3]
......@@ -339,10 +335,11 @@ public:
void solDependentSource(PrimaryVariables &source,
const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx, const ElementVolumeVariables &elemVolVars) const
int scvIdx,
const ElementVolumeVariables &elemVolVars) const
{
source = 0;
const VolumeVariables &volVars = elemVolVars[scvIdx];
const auto& volVars = elemVolVars[scvIdx];
Scalar moleFracNaCl_lPhase = volVars.moleFraction(wPhaseIdx, NaClIdx);
Scalar moleFracNaCl_gPhase = volVars.moleFraction(nPhaseIdx, NaClIdx);
Scalar massFracNaCl_Max_lPhase = this->spatialParams().SolubilityLimit();
......@@ -352,8 +349,8 @@ public:
// liquid phase
Scalar precipSalt = volVars.porosity() * volVars.molarDensity(wPhaseIdx)
* volVars.saturation(wPhaseIdx)
* pow(std::abs(moleFracNaCl_lPhase - moleFracNaCl_Max_lPhase), 1.0);
* volVars.saturation(wPhaseIdx)
* std::abs(moleFracNaCl_lPhase - moleFracNaCl_Max_lPhase);
if (moleFracNaCl_lPhase < moleFracNaCl_Max_lPhase)
precipSalt *= -1;
......@@ -361,8 +358,8 @@ public:
// gas phase
if (moleFracNaCl_gPhase > moleFracNaCl_Max_gPhase)
precipSalt += volVars.porosity() * volVars.molarDensity(nPhaseIdx)
* volVars.saturation(nPhaseIdx)
* pow(std::abs(moleFracNaCl_gPhase - moleFracNaCl_Max_gPhase), 1.0);
* volVars.saturation(nPhaseIdx)
* std::abs(moleFracNaCl_gPhase - moleFracNaCl_Max_gPhase);
// make sure we don't disolve more salt than previously precipitated
if (precipSalt*this->timeManager().timeStepSize() + volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)< 0)
......@@ -380,22 +377,12 @@ public:
* \brief Return the initial phase state inside a control volume.
*/
int initialPhasePresence(const Vertex &vert,
int &globalIdx,
int globalIdx,
const GlobalPosition &globalPos) const
{
return bothPhases;
}
void preTimeStep()
{}
void postTimeStep()
{}
void episodeEnd()
{}
private:
/*!
......@@ -410,7 +397,7 @@ private:
const Scalar X_NaCl = XlNaCl;
/* XlNaCl: conversion from mass fraction to mol fraction */
const Scalar xlNaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
auto xlNaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
return xlNaCl;
}
......
......@@ -43,13 +43,11 @@ SET_TYPE_PROP(DissolutionSpatialparams, SpatialParams, Dumux::DissolutionSpatial
SET_PROP(DissolutionSpatialparams, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
// define the material law which is parameterized by effective saturations
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
public:
// define the material law parameterized by absolute saturations
typedef EffToAbsLaw<EffMaterialLaw> type;
typedef Dumux::EffToAbsLaw<Dumux::RegularizedBrooksCorey<Scalar> > type;
};
}
......@@ -61,12 +59,11 @@ template<class TypeTag>
class DissolutionSpatialparams : public ImplicitSpatialParams<TypeTag>
{
typedef ImplicitSpatialParams<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename Grid::ctype CoordScalar;
typedef typename GridView::ctype CoordScalar;
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld,
......@@ -78,53 +75,29 @@ class DissolutionSpatialparams : public ImplicitSpatialParams<TypeTag>
nPhaseIdx = FluidSystem::nPhaseIdx,
};
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
typedef std::vector<Scalar> PermeabilityType;
typedef std::vector<MaterialLawParams> MaterialLawParamsVector;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
public:
DissolutionSpatialparams(const GridView &gridView)
: ParentType(gridView),
K_(0)
: ParentType(gridView), K_(0.0)
{
//set main diagonal entries of the permeability tensor to a value
//setting to one value means: isotropic, homogeneous
// set main diagonal entries of the permeability tensor to a value
// setting to one value means: isotropic, homogeneous
for (int i = 0; i < dim; i++)
K_[i][i] = 2.23e-14;
// residual saturations
materialParams_.setSwr(0.2);
materialParams_.setSnr(1E-3);
materialParams_.setSnr(1e-3);
//parameters of Brooks & Corey Law
// parameters of Brooks & Corey Law
materialParams_.setPe(500);
materialParams_.setLambda(2);
}
~DissolutionSpatialparams()
{}
/*!
* \brief Update the spatial parameters with the flow solution
* after a timestep.
*
* \param globalSolution The global solution vector
*/
void update(const SolutionVector &globalSolution)
{ };
/*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
* on the position in the domain
*
......@@ -136,9 +109,9 @@ public:
* could be defined, where globalPos is the vector including the global coordinates
* of the finite volume.
*/
const Dune::FieldMatrix<Scalar, dim, dim> &intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
const Tensor& intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return K_;
}
......@@ -151,7 +124,7 @@ public:
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
double porosityMin(const Element &element,
Scalar porosityMin(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
......@@ -166,30 +139,30 @@ public:
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
double porosity(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
Scalar porosity(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
return 0.11;
}
double solidity(const Element &element,
Scalar solidity(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
return 1 - 0.11;
return 1.0 - porosity(element, fvGeometry, scvIdx);
}
double SolubilityLimit() const
Scalar SolubilityLimit() const
{
return 0.26;
}
double theta(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
Scalar theta(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
return 10.0;
}
......@@ -197,14 +170,14 @@ public:
// return the brooks-corey context depending on the position
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
return materialParams_;
}
private:
Dune::FieldMatrix<Scalar, dim, dim> K_;
Tensor K_;
MaterialLawParams materialParams_;
};
......
......@@ -23,12 +23,12 @@ TemperatureLow = 273.15 # [Pa]low end for tabularization of fluid pr
TemperatureHigh = 400.00 # [Pa]high end for tabularization of fluid properties
[Problem]
OutputName = saltflushbox # [-] name for output files
OutputName = saltflushbox # [-] name for output files
ReservoirPressure = 11E6 # [Pa] Initial reservoir pressure
reservoirSaturation = 0.4 # [-] Initial saturation
solidity = 0.0 # [-] Initial solid salt precipitate
Temperature = 418.15 # [K] reservoir temperature
InnerPressure = 13E6 # [Pa]
InnerPressure = 13E6 # [Pa]
InnerLiqSaturation = 0.95 # [-] liquid saturation at inner boundary
InnerSalinity = 0.0001 # [-] salinity of inner liquid
OuterPressure = 11E6 # [Pa] reservoir boundary pressure
......@@ -41,7 +41,7 @@ InitPrecipitatedSalt2 = 0.05 # [-] initial precipitated salt in the bloc
[SpatialParams]