diff --git a/dumux/discretization/cellcentered/globalvolumevariables.hh b/dumux/discretization/cellcentered/globalvolumevariables.hh
index 371791f1e494ed2eaf7fe57fb2b7ddcb93f81848..b4d691d2c0f01150db1a057d6951644a0f0f7d85 100644
--- a/dumux/discretization/cellcentered/globalvolumevariables.hh
+++ b/dumux/discretization/cellcentered/globalvolumevariables.hh
@@ -45,6 +45,9 @@ class CCGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
friend CCElementVolumeVariables<TypeTag, true>;
// The local jacobian needs to access and change volVars for derivative calculation
friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
+ // as does the primary variable switch
+ friend class PrimaryVariableSwitch<TypeTag>;
+ friend typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -103,13 +106,12 @@ public:
friend inline ElementVolumeVariables localView(const CCGlobalVolumeVariables& global)
{ return ElementVolumeVariables(global); }
-private:
const VolumeVariables& volVars(const IndexType scvIdx) const
{ return volumeVariables_[scvIdx]; }
VolumeVariables& volVars(const IndexType scvIdx)
{ return volumeVariables_[scvIdx]; }
-
+private:
const Problem& problem_() const
{ return *problemPtr_; }
diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index cdbe9ca41bfc887977c35b165cda8443dccd792e..24499d68bcb11f0eca2a53035f99c12e2676c3be 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -155,7 +155,7 @@ public:
const auto insideScvIdx = scvFace.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = elemVolVars[insideScvIdx];
- const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv, insideVolVars);
+ const auto insideK = problem.spatialParams().solDependentIntrinsicPermeability(insideScv, insideVolVars);
Scalar ti = calculateOmega_(problem, scvFace, insideK, element, insideScv);
if (!scvFace.boundary())
@@ -166,7 +166,7 @@ public:
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto outsideElement = fvGeometry.globalFvGeometry().element(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
- const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv, outsideVolVars);
+ const auto outsideK = problem.spatialParams().solDependentIntrinsicPermeability(outsideScv, outsideVolVars);
Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideElement, outsideScv);
tij = scvFace.area()*(ti * tj)/(ti + tj);
diff --git a/dumux/porousmediumflow/2pncmin/implicit/fluxvariables.hh b/dumux/porousmediumflow/2pncmin/implicit/fluxvariables.hh
deleted file mode 100644
index 71c26b96cb3f8fd5978c8164ce1d3144cbdad6bc..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2pncmin/implicit/fluxvariables.hh
+++ /dev/null
@@ -1,179 +0,0 @@
-// -**- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Contains the data which is required to calculate
- * all fluxes of components over a face of a finite volume for
- * the two-phase two-component mineralization model fully implicit model.
- */
-#ifndef DUMUX_2PNCMIN_FLUX_VARIABLES_HH
-#define DUMUX_2PNCMIN_FLUX_VARIABLES_HH
-
-#include <dumux/common/math.hh>
-#include <dumux/common/spline.hh>
-#include <dumux/porousmediumflow/2pnc/implicit/fluxvariables.hh>
-#include "properties.hh"
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPNCMinModel
- * \ingroup ImplicitFluxVariables
- * \brief Contains the data which is required to calculate
- * all fluxes of components over a face of a finite volume for
- * the two-phase n-component mineralization fully implicit model.
- *
- * This means pressure and concentration gradients, phase densities at
- * the integration point, etc.
- */
-
-template <class TypeTag>
-class TwoPNCMinFluxVariables : public TwoPNCFluxVariables<TypeTag>
-{
- friend typename GET_PROP_TYPE(TypeTag, BaseFluxVariables); // be friends with base class
- friend class TwoPNCFluxVariables<TypeTag>; // be friends with parent class
-
- typedef TwoPNCFluxVariables<TypeTag> ParentType;
- typedef TwoPNCMinFluxVariables<TypeTag> ThisType;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GridView::ctype CoordScalar;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
- enum {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
- numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> DimWorldMatrix;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- enum {
- wPhaseIdx = FluidSystem::wPhaseIdx,
- nPhaseIdx = FluidSystem::nPhaseIdx,
- wCompIdx = FluidSystem::wCompIdx,
- };
-
-protected:
- /*!
- * \brief Actual calculation of the normal Darcy velocities.
- * \note this overloads the darcy flux variables velocity calculation
- * This is only necessary because of the permeability factor.
- * \todo Remove this once we have solDependent spatialParams!
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateNormalVelocity_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the mean intrinsic permeability
- const auto& spatialParams = problem.spatialParams();
- DimWorldMatrix K(0.0);
- const auto& volVarsI = elemVolVars[this->face().i];
- const auto& volVarsJ = elemVolVars[this->face().j];
-
- if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
- {
- auto Ki = spatialParams.intrinsicPermeability(element, this->fvGeometry_(), this->face().i);
- Ki *= volVarsI.permeabilityFactor();
-
- auto Kj = spatialParams.intrinsicPermeability(element, this->fvGeometry_(), this->face().j);
- Kj *= volVarsJ.permeabilityFactor();
-
- spatialParams.meanK(K, Ki, Kj);
- }
- else
- {
- const Element& elementI = this->fvGeometry_().neighbors[this->face().i];
- FVElementGeometry fvGeometryI;
- fvGeometryI.subContVol[0].global = elementI.geometry().center();
-
- const Element& elementJ = this->fvGeometry_().neighbors[this->face().j];
- FVElementGeometry fvGeometryJ;
- fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
-
- auto Ki = spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0);
- Ki *= volVarsI.permeabilityFactor();
-
- auto Kj = spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0);
- Kj *= volVarsJ.permeabilityFactor();
-
- spatialParams.meanK(K, Ki, Kj);
- }
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // calculate the flux in the normal direction of the
- // current sub control volume face:
- //
- // v = - (K_f grad phi) * n
- // with K_f = rho g / mu K
- //
- // Mind, that the normal has the length of it's area.
- // This means that we are actually calculating
- // Q = - (K grad phi) dot n /|n| * A
-
-
- K.mv(this->potentialGrad_[phaseIdx], this->kGradP_[phaseIdx]);
- this->kGradPNormal_[phaseIdx] = this->kGradP_[phaseIdx]*this->face().normal;
-
- // determine the upwind direction
- if (this->kGradPNormal_[phaseIdx] < 0)
- {
- this->upstreamIdx_[phaseIdx] = this->face().i;
- this->downstreamIdx_[phaseIdx] = this->face().j;
- }
- else
- {
- this->upstreamIdx_[phaseIdx] = this->face().j;
- this->downstreamIdx_[phaseIdx] = this->face().i;
- }
-
- // obtain the upwind volume variables
- const auto& upVolVars = elemVolVars[ this->upstreamIdx(phaseIdx) ];
- const auto& downVolVars = elemVolVars[ this->downstreamIdx(phaseIdx) ];
-
- // the minus comes from the Darcy relation which states that
- // the flux is from high to low potentials.
- // set the velocity
- this->velocity_[phaseIdx] = this->kGradP_[phaseIdx];
- this->velocity_[phaseIdx] *= - ( this->mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
- + (1.0 - this->mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx)) ;
-
- // set the volume flux
- this->volumeFlux_[phaseIdx] = this->velocity_[phaseIdx] * this->face().normal;
- }// loop all phases
- }
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/2pncmin/implicit/localresidual.hh b/dumux/porousmediumflow/2pncmin/implicit/localresidual.hh
index 1f53e0eabf48d1b83bb47cb7032824f53ec277c6..113877310a801e6d3f097a5d7636cc64ad295409 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/localresidual.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/localresidual.hh
@@ -19,79 +19,45 @@
/*!
* \file
*
- * \brief Element-wise calculation of the Jacobian matrix for problems
+ * \brief Element-wise calculation of the local residual for problems
* using the two-phase n-component mineralisation box model.
*/
-#ifndef DUMUX_2PNCMIN_LOCAL_RESIDUAL_BASE_HH
-#define DUMUX_2PNCMIN_LOCAL_RESIDUAL_BASE_HH
+#ifndef DUMUX_2PNCMIN_LOCAL_RESIDUAL_HH
+#define DUMUX_2PNCMIN_LOCAL_RESIDUAL_HH
#include "properties.hh"
-#include <dumux/porousmediumflow/2pnc/implicit/localresidual.hh>
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
namespace Dumux
{
/*!
* \ingroup TwoPNCMinModel
* \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
+ * \brief Element-wise calculation of the local residual for problems
* 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.
*/
template<class TypeTag>
-class TwoPNCMinLocalResidual: public TwoPNCLocalResidual<TypeTag>
+class TwoPNCMinLocalResidual: public CompositionalLocalResidual<TypeTag>
{
-protected:
- typedef TwoPNCLocalResidual<TypeTag> ParentType;
- typedef TwoPNCMinLocalResidual<TypeTag> ThisType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
+ using ParentType = CompositionalLocalResidual<TypeTag>;
+ using ThisType = TwoPNCMinLocalResidual<TypeTag>;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum
{
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
- replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
-
- pressureIdx = Indices::pressureIdx,
- switchIdx = Indices::switchIdx,
-
- wPhaseIdx = FluidSystem::wPhaseIdx,
- nPhaseIdx = FluidSystem::nPhaseIdx,
-
- wCompIdx = FluidSystem::wCompIdx,
- nCompIdx = FluidSystem::nCompIdx,
-
- conti0EqIdx = Indices::conti0EqIdx,
-
- wPhaseOnly = Indices::wPhaseOnly,
- nPhaseOnly = Indices::nPhaseOnly,
- bothPhases = Indices::bothPhases,
-
- plSg = TwoPNCFormulation::plSg,
- pgSl = TwoPNCFormulation::pgSl,
- formulation = GET_PROP_VALUE(TypeTag, Formulation)
+ conti0EqIdx = Indices::conti0EqIdx
};
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
public:
/*!
* \brief Evaluate the amount all conservation quantities
@@ -101,17 +67,14 @@ public:
* inside a sub control volume divided by the volume).
* In contrast to the 2pnc model, here, the storage of solid phases is included too.
*
- * \param storage the mass of the component within the sub-control volume
- * \param scvIdx The SCV (sub-control-volume) index
- * \param usePrevSol Evaluate function with solution of current or previous time step
+ * \param scv the SCV (sub-control-volume)
+ * \param volVars The volume variables of the right time step
*/
- void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
+ PrimaryVariables computeStorage(const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
{
- //call parenttype function
- ParentType::computeStorage(storage, scvIdx, usePrevSol);
-
- const auto& elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
- const VolumeVariables &volVars = elemVolVars[scvIdx];
+ // call parenttype function
+ auto storage = ParentType::computeStorage(scv, volVars);
// Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
for (int phaseIdx = numPhases; phaseIdx < numPhases + numSPhases; ++phaseIdx)
@@ -120,10 +83,10 @@ public:
storage[eqIdx] += volVars.precipitateVolumeFraction(phaseIdx)*volVars.molarDensity(phaseIdx);
}
- Valgrind::CheckDefined(storage);
+ return storage;
}
};
-} // end namespace
+} // end namespace Dumux
#endif
diff --git a/dumux/porousmediumflow/2pncmin/implicit/model.hh b/dumux/porousmediumflow/2pncmin/implicit/model.hh
index 55758e6d4ac41508d2ffe2ed3aba2faf9b06e54c..d5a83a4f15e2b34a1c0c107fe0c4e56998726cd4 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/model.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh
@@ -27,10 +27,11 @@
#include "properties.hh"
#include "indices.hh"
+#include "primaryvariableswitch.hh"
+#include "localresidual.hh"
#include <dumux/material/constants.hh>
#include <dumux/porousmediumflow/2pnc/implicit/model.hh>
-#include "localresidual.hh"
#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
namespace Dumux
@@ -108,22 +109,23 @@ namespace Dumux
*/
template<class TypeTag>
-class TwoPNCMinModel: public TwoPNCModel<TypeTag>
+class TwoPNCMinModel: public GET_PROP_TYPE(TypeTag, BaseModel)
{
- typedef TwoPNCMinModel<TypeTag> ThisType;
- typedef TwoPNCModel<TypeTag> ParentType;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- typedef Constants<Scalar> Constant;
+ // the parent class needs to access the variable switch
+ friend typename GET_PROP_TYPE(TypeTag, BaseModel);
+
+ using ThisType = Dumux::TwoPNCMinModel<TypeTag>;
+ using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using Constant = Constants<Scalar>;
enum {
dim = GridView::dimension,
@@ -154,18 +156,105 @@ class TwoPNCMinModel: public TwoPNCModel<TypeTag>
formulation = GET_PROP_VALUE(TypeTag, Formulation)
};
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename GridView::ctype CoordScalar;
- typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
- typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ using Vertex = typename GridView::template Codim<dim>::Entity;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using CoordScalar = typename GridView::ctype;
+ using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
+ using PhasesVector = Dune::FieldVector<Scalar, numPhases>;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
public:
+ /*!
+ * \brief One Newton iteration was finished.
+ * \param uCurrent The solution after the current Newton iteration
+ */
+ template<typename T = TypeTag>
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), void>::type
+ newtonEndStep()
+ {
+ // \todo resize volvars vector if grid was adapted
+
+ // update the variable switch
+ switchFlag_ = priVarSwitch_().update(this->problem_(), this->curSol());
+
+ // update the secondary variables if global caching is enabled
+ // \note we only updated if phase presence changed as the volume variables
+ // are already updated once by the switch
+ if (switchFlag_)
+ {
+ for (const auto& element : elements(this->problem_().gridView()))
+ {
+ // make sure FVElementGeometry & vol vars are bound to the element
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto dofIdxGlobal = scv.dofIndex();
+ if (priVarSwitch_().wasSwitched(dofIdxGlobal))
+ {
+ this->nonConstCurGlobalVolVars().volVars(dofIdxGlobal).update(this->curSol()[dofIdxGlobal],
+ this->problem_(),
+ element,
+ scv);
+ }
+ }
+
+ }
+ }
+ }
+
+ /*!
+ * \brief One Newton iteration was finished.
+ * \param uCurrent The solution after the current Newton iteration
+ */
+ template<typename T = TypeTag>
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), void>::type
+ newtonEndStep()
+ {
+ // update the variable switch
+ switchFlag_ = priVarSwitch_().update(this->problem_(), this->curSol());
+ }
+
+ /*!
+ * \brief Called by the update() method if applying the Newton
+ * method was unsuccessful.
+ */
+ void updateFailed()
+ {
+ ParentType::updateFailed();
+
+ switchFlag_ = false;
+ priVarSwitch_().resetPhasePresence();
+ }
+
+ /*!
+ * \brief Called by the problem if a time integration was
+ * successful, post processing of the solution is done and the
+ * result has been written to disk.
+ *
+ * This should prepare the model for the next time integration.
+ */
+ void advanceTimeLevel()
+ {
+ ParentType::advanceTimeLevel();
+
+ // update the phase state
+ priVarSwitch_().updateOldPhasePresence();
+ switchFlag_ = false;
+ }
+
+ /*!
+ * \brief Returns true if the primary variables were switched for
+ * at least one dof after the last timestep.
+ */
+ bool switched() const
+ {
+ return switchFlag_;
+ }
/*!
* \brief Append all quantities of interest which can be derived
@@ -179,114 +268,112 @@ public:
//additional output of the permeability and the precipitate volume fractions
void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
{
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
- typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
+ using ScalarField = Dune::BlockVector<Dune::FieldVector<double, 1> >;
// get the number of degrees of freedom
auto numDofs = this->numDofs();
// create the required scalar fields
- ScalarField *Sg = writer.allocateManagedBuffer(numDofs);
- ScalarField *Sl = writer.allocateManagedBuffer(numDofs);
- ScalarField *pg = writer.allocateManagedBuffer (numDofs);
- ScalarField *pl = writer.allocateManagedBuffer (numDofs);
- ScalarField *pc = writer.allocateManagedBuffer (numDofs);
- ScalarField *rhoL = writer.allocateManagedBuffer (numDofs);
- ScalarField *rhoG = writer.allocateManagedBuffer (numDofs);
- ScalarField *mobL = writer.allocateManagedBuffer (numDofs);
- ScalarField *mobG = writer.allocateManagedBuffer (numDofs);
- ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
- ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
- ScalarField *poro = writer.allocateManagedBuffer (numDofs);
- ScalarField *permeabilityFactor = writer.allocateManagedBuffer (numDofs);
- ScalarField *precipitateVolumeFraction[numSPhases];
+ auto* Sg = writer.allocateManagedBuffer(numDofs);
+ auto* Sl = writer.allocateManagedBuffer(numDofs);
+ auto* pg = writer.allocateManagedBuffer (numDofs);
+ auto* pl = writer.allocateManagedBuffer (numDofs);
+ auto* pc = writer.allocateManagedBuffer (numDofs);
+ auto* rhoL = writer.allocateManagedBuffer (numDofs);
+ auto* rhoG = writer.allocateManagedBuffer (numDofs);
+ auto* mobL = writer.allocateManagedBuffer (numDofs);
+ auto* mobG = writer.allocateManagedBuffer (numDofs);
+ auto* phasePresence = writer.allocateManagedBuffer (numDofs);
+ auto* temperature = writer.allocateManagedBuffer (numDofs);
+ auto* poro = writer.allocateManagedBuffer (numDofs);
+ auto* permeabilityFactor = writer.allocateManagedBuffer (numDofs);
+ ScalarField* precipitateVolumeFraction[numSPhases];
for (int i = 0; i < numSPhases; ++i)
precipitateVolumeFraction[i] = writer.allocateManagedBuffer(numDofs);
- ScalarField *massFraction[numPhases][numComponents];
+ ScalarField* massFraction[numPhases][numComponents];
for (int i = 0; i < numPhases; ++i)
for (int j = 0; j < numComponents; ++j)
massFraction[i][j] = writer.allocateManagedBuffer(numDofs);
- ScalarField *molarity[numComponents];
+ ScalarField* molarity[numComponents];
for (int j = 0; j < numComponents ; ++j)
molarity[j] = writer.allocateManagedBuffer(numDofs);
- ScalarField *Perm[dim];
+ ScalarField* Perm[dim];
for (int j = 0; j < dim; ++j) //Permeability only in main directions xx and yy
Perm[j] = writer.allocateManagedBuffer(numDofs);
- VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
- VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
- ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+ // auto* velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ // auto* velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ // ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
- if (velocityOutput.enableOutput()) // check if velocity output is demanded
- {
- // initialize velocity fields
- for (unsigned int i = 0; i < numDofs; ++i)
- {
- (*velocityN)[i] = Scalar(0);
- (*velocityW)[i] = Scalar(0);
- }
- }
+ // if (velocityOutput.enableOutput()) // check if velocity output is demanded
+ // {
+ // // initialize velocity fields
+ // for (unsigned int i = 0; i < numDofs; ++i)
+ // {
+ // (*velocityN)[i] = Scalar(0);
+ // (*velocityW)[i] = Scalar(0);
+ // }
+ // }
auto numElements = this->gridView_().size(0);
- ScalarField *rank = writer.allocateManagedBuffer(numElements);
+ auto* rank = writer.allocateManagedBuffer(numElements);
for (const auto& element : elements(this->gridView_()))
{
auto eIdxGlobal = this->problem_().elementMapper().index(element);
(*rank)[eIdxGlobal] = this->gridView_().comm().rank();
- FVElementGeometry fvGeometry;
- fvGeometry.update(this->gridView_(), element);
- ElementVolumeVariables elemVolVars;
- elemVolVars.update(this->problem_(),
- element,
- fvGeometry,
- false /* oldSol? */);
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ auto elemVolVars = localView(this->curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, this->curSol());
+
+ for (auto&& scv : scvs(fvGeometry))
{
- auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
- (*Sg)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
- (*Sl)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
- (*pg)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
- (*pl)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
- (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
- (*rhoL)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
- (*rhoG)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
- (*mobL)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
- (*mobG)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
- (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+ auto dofIdxGlobal = scv.dofIndex();
+ const auto& volVars = elemVolVars[scv];
+
+ (*Sg)[dofIdxGlobal] = volVars.saturation(nPhaseIdx);
+ (*Sl)[dofIdxGlobal] = volVars.saturation(wPhaseIdx);
+ (*pg)[dofIdxGlobal] = volVars.pressure(nPhaseIdx);
+ (*pl)[dofIdxGlobal] = volVars.pressure(wPhaseIdx);
+ (*pc)[dofIdxGlobal] = volVars.capillaryPressure();
+ (*rhoL)[dofIdxGlobal] = volVars.density(wPhaseIdx);
+ (*rhoG)[dofIdxGlobal] = volVars.density(nPhaseIdx);
+ (*mobL)[dofIdxGlobal] = volVars.mobility(wPhaseIdx);
+ (*mobG)[dofIdxGlobal] = volVars.mobility(nPhaseIdx);
+ (*poro)[dofIdxGlobal] = volVars.porosity();
for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
- (*precipitateVolumeFraction[sPhaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].precipitateVolumeFraction(sPhaseIdx + numPhases);
+ (*precipitateVolumeFraction[sPhaseIdx])[dofIdxGlobal] = volVars.precipitateVolumeFraction(sPhaseIdx + numPhases);
- (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
- (*permeabilityFactor)[dofIdxGlobal] = elemVolVars[scvIdx].permeabilityFactor();
- (*phasePresence)[dofIdxGlobal] = this->staticDat_[dofIdxGlobal].phasePresence;
+ (*temperature)[dofIdxGlobal] = volVars.temperature();
+ // (*permeabilityFactor)[dofIdxGlobal] = volVars.permeabilityFactor();
+ (*phasePresence)[dofIdxGlobal] = priVarSwitch().phasePresence(dofIdxGlobal);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- (*massFraction[phaseIdx][compIdx])[dofIdxGlobal]= elemVolVars[scvIdx].massFraction(phaseIdx,compIdx);
+ (*massFraction[phaseIdx][compIdx])[dofIdxGlobal]= volVars.massFraction(phaseIdx,compIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- (*molarity[compIdx])[dofIdxGlobal] = (elemVolVars[scvIdx].molarity(wPhaseIdx, compIdx));
+ (*molarity[compIdx])[dofIdxGlobal] = (volVars.molarity(wPhaseIdx, compIdx));
- Tensor K = this->perm_(this->problem_().spatialParams().intrinsicPermeability(element, fvGeometry, scvIdx));
+ auto K = this->perm_(this->problem_().spatialParams().solDependentIntrinsicPermeability(scv, volVars));
for (int j = 0; j<dim; ++j)
- (*Perm[j])[dofIdxGlobal] = K[j][j] * elemVolVars[scvIdx].permeabilityFactor();
+ (*Perm[j])[dofIdxGlobal] = K[j][j];
};
- // velocity output
- if(velocityOutput.enableOutput()){
- velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
- velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
- }
+ // // velocity output
+ // if(velocityOutput.enableOutput()){
+ // velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+ // velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+ // }
} // loop over element
writer.attachDofData(*Sg, "Sg", isBox);
@@ -333,180 +420,103 @@ public:
writer.attachDofData(*molarity[j], oss.str(), isBox);
}
- if (velocityOutput.enableOutput()) // check if velocity output is demanded
- {
- writer.attachDofData(*velocityW, "velocityW", isBox, dim);
- writer.attachDofData(*velocityN, "velocityN", isBox, dim);
- }
+ // if (velocityOutput.enableOutput()) // check if velocity output is demanded
+ // {
+ // writer.attachDofData(*velocityW, "velocityW", isBox, dim);
+ // writer.attachDofData(*velocityN, "velocityN", isBox, dim);
+ // }
writer.attachCellData(*rank, "process rank");
}
/*!
- * \brief Update the static data of all vertices in the grid.
+ * \brief Write the current solution to a restart file.
*
- * \param curGlobalSol The current global solution
- * \param oldGlobalSol The previous global solution
+ * \param outStream The output stream of one entity for the restart file
+ * \param entity The entity, either a vertex or an element
*/
- void updateStaticData(SolutionVector &curGlobalSol,
- const SolutionVector &oldGlobalSol)
+ template<class Entity>
+ void serializeEntity(std::ostream &outStream, const Entity &entity)
{
- bool wasSwitched = false;
+ // write primary variables
+ ParentType::serializeEntity(outStream, entity);
- for (unsigned i = 0; i < this->staticDat_.size(); ++i)
- this->staticDat_[i].visited = false;
+ int dofIdxGlobal = this->dofMapper().index(entity);
- for (const auto& element : elements(this->gridView_()))
- {
- FVElementGeometry fvGeometry;
- fvGeometry.update(this->gridView_(), element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
- if (this->staticDat_[dofIdxGlobal].visited)
- continue;
-
- this->staticDat_[dofIdxGlobal].visited = true;
- VolumeVariables volVars;
- volVars.update(curGlobalSol[dofIdxGlobal],
- this->problem_(),
- element,
- fvGeometry,
- scvIdx,
- false);
- auto global = element.geometry().corner(scvIdx);
- if (primaryVarSwitch_(curGlobalSol, volVars, dofIdxGlobal, global))
- wasSwitched = true;
- }
- }
+ if (!outStream.good())
+ DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal);
- // make sure that if there was a variable switch in an
- // other partition we will also set the switch flag
- // for our partition.
- if (this->gridView_().comm().size() > 1)
- wasSwitched = this->gridView_().comm().max(wasSwitched);
+ outStream << priVarSwitch().phasePresence(dofIdxGlobal) << " ";
+ }
- this->setSwitched_(wasSwitched);
+ /*!
+ * \brief Reads the current solution from a restart file.
+ *
+ * \param inStream The input stream of one entity from the restart file
+ * \param entity The entity, either a vertex or an element
+ */
+ template<class Entity>
+ void deserializeEntity(std::istream &inStream, const Entity &entity)
+ {
+ // read primary variables
+ ParentType::deserializeEntity(inStream, entity);
+
+ // read phase presence
+ int dofIdxGlobal = this->dofMapper().index(entity);
+
+ if (!inStream.good())
+ DUNE_THROW(Dune::IOError, "Could not deserialize entity " << dofIdxGlobal);
+
+ int phasePresence;
+ inStream >> phasePresence;
+
+ priVarSwitch_().setPhasePresence(dofIdxGlobal, phasePresence);
+ priVarSwitch_().setOldPhasePresence(dofIdxGlobal, phasePresence);
}
-protected:
+
+ const Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag>& priVarSwitch() const
+ { return switch_; }
+
+private:
+
+ Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag>& priVarSwitch_()
+ { return switch_; }
/*!
- * \brief Set whether there was a primary variable switch after in
- * the last timestep.
+ * \brief Applies the initial solution for all vertices of the grid.
+ *
+ * \todo the initial condition needs to be unique for
+ * each vertex. we should think about the API...
*/
- bool primaryVarSwitch_(SolutionVector &globalSol,
- const VolumeVariables &volVars,
- int dofIdxGlobal,
- const GlobalPosition &globalPos)
+ void applyInitialSolution_()
{
- // evaluate primary variable switch
- bool wouldSwitch = false;
- int phasePresence = this->staticDat_[dofIdxGlobal].phasePresence;
- int newPhasePresence = phasePresence;
+ ParentType::applyInitialSolution_();
- //check if a primary variable switch is necessary
- if (phasePresence == bothPhases)
- {
- Scalar Smin = 0.0; //saturation threshold
- if (this->staticDat_[dofIdxGlobal].wasSwitched)
- Smin = -0.01;
+ // initialize the primary variable switch
+ priVarSwitch_().init(this->problem_());
+ }
- //if saturation of liquid phase is smaller 0 switch
- if (volVars.saturation(wPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- //liquid phase has to disappear
- std::cout << "Liquid Phase disappears at vertex " << dofIdxGlobal
- << ", coordinated: " << globalPos << ", Sl: "
- << volVars.saturation(wPhaseIdx) << std::endl;
- newPhasePresence = nPhaseOnly;
-
- //switch not depending on formulation
- //switch "Sl" to "xgH20"
- globalSol[dofIdxGlobal][switchIdx]
- = volVars.moleFraction(nPhaseIdx, wCompIdx /*H2O*/);
- //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
- }
- //if saturation of gas phase is smaller than 0 switch
- else if (volVars.saturation(nPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- //gas phase has to disappear
- std::cout << "Gas Phase disappears at vertex " << dofIdxGlobal
- << ", coordinated: " << globalPos << ", Sg: "
- << volVars.saturation(nPhaseIdx) << std::endl;
- newPhasePresence = wPhaseOnly;
-
- //switch "Sl" to "xlN2"
- globalSol[dofIdxGlobal][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++)
- {
- sumxl += volVars.moleFraction(wPhaseIdx, compIdx);
- }
- Scalar xlmax = 1.0;
- if (sumxl > xlmax)
- wouldSwitch = true;
- if (this->staticDat_[dofIdxGlobal].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 " << dofIdxGlobal
- << ", coordinated: " << globalPos << ", sumxl: "
- << sumxl << std::endl;
- newPhasePresence = bothPhases;
- if (formulation == pgSl)
- globalSol[dofIdxGlobal][switchIdx] = 0.0;
- else if (formulation == plSg)
- globalSol[dofIdxGlobal][switchIdx] = 1.0;
- //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
- }
- }
- else if (phasePresence == wPhaseOnly)
- {
- Scalar xgmax = 1;
- Scalar sumxg = 0;
- //Calculate sum of mole fractions in the hypothetical gas phase
- for (int compIdx = 0; compIdx < numComponents; compIdx++)
- {
- sumxg += volVars.moleFraction(nPhaseIdx, compIdx);
- }
- if (sumxg > xgmax)
- wouldSwitch = true;
- if (this->staticDat_[dofIdxGlobal].wasSwitched)
- xgmax *=1.02;
- //liquid phase appears if sum is larger than one
- if (sumxg > xgmax)
- {
- std::cout << "Gas Phase appears at vertex " << dofIdxGlobal
- << ", 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)
- globalSol[dofIdxGlobal][switchIdx] = 0.999;
- else if (formulation == plSg)
- globalSol[dofIdxGlobal][switchIdx] = 0.001;
+ //! the class handling the primary variable switch
+ Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag> switch_;
+ bool switchFlag_;
- }
- }
- this->staticDat_[dofIdxGlobal].phasePresence = newPhasePresence;
- this->staticDat_[dofIdxGlobal].wasSwitched = wouldSwitch;
- return phasePresence != newPhasePresence;
+ Tensor perm_(Scalar perm) const
+ {
+ Tensor K(0.0);
+
+ for(int i=0; i<dim; i++)
+ K[i][i] = perm;
+
+ return K;
+ }
+
+ const Tensor& perm_(const Tensor& perm) const
+ {
+ return perm;
}
};
-}
+} // end namespace Dumux
#include "propertydefaults.hh"
diff --git a/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh b/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0f3ab6d825fdb70a2dcdc89b846b00632503fbdc
--- /dev/null
+++ b/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh
@@ -0,0 +1,190 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief The primary variable switch for the 2pncmin model
+ */
+#ifndef DUMUX_2PNCMIN_PRIMARY_VARIABLE_SWITCH_HH
+#define DUMUX_2PNCMIN_PRIMARY_VARIABLE_SWITCH_HH
+
+#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCMinModel
+ * \brief The primary variable switch controlling the phase presence state variable
+ */
+template<class TypeTag>
+class TwoPNCMinPrimaryVariableSwitch : public Dumux::PrimaryVariableSwitch<TypeTag>
+{
+ friend typename Dumux::PrimaryVariableSwitch<TypeTag>;
+ using ParentType = Dumux::PrimaryVariableSwitch<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+ static const int numMajorComponents = GET_PROP_VALUE(TypeTag, NumMajorComponents);
+
+ enum {
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ wCompIdx = FluidSystem::wCompIdx,
+ nCompIdx = FluidSystem::nCompIdx,
+
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+
+ enum {
+ plSg = TwoPNCFormulation::plSg,
+ pgSl = TwoPNCFormulation::pgSl,
+ formulation = GET_PROP_VALUE(TypeTag, Formulation)
+ };
+
+protected:
+
+ // perform variable switch at a degree of freedom location
+ bool update_(PrimaryVariables& priVars,
+ const VolumeVariables& volVars,
+ IndexType dofIdxGlobal,
+ const GlobalPosition& globalPos)
+ {
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = this->phasePresence_[dofIdxGlobal];
+ int newPhasePresence = phasePresence;
+
+ //check if a primary variable switch is necessary
+ if (phasePresence == bothPhases)
+ {
+ Scalar Smin = 0.0; //saturation threshold
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ //if saturation of liquid phase is smaller 0 switch
+ if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ //liquid phase has to disappear
+ std::cout << "Liquid Phase disappears at vertex " << dofIdxGlobal
+ << ", coordinated: " << globalPos << ", Sl: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ newPhasePresence = nPhaseOnly;
+
+ //switch not depending on formulation
+ //switch "Sl" to "xgH20"
+ priVars[switchIdx] = volVars.moleFraction(nPhaseIdx, wCompIdx /*H2O*/);
+ //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+ }
+ //if saturation of gas phase is smaller than 0 switch
+ else if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ //gas phase has to disappear
+ std::cout << "Gas Phase disappears at vertex " << dofIdxGlobal
+ << ", coordinated: " << globalPos << ", Sg: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ newPhasePresence = wPhaseOnly;
+
+ //switch "Sl" to "xlN2"
+ priVars[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++)
+ {
+ sumxl += volVars.moleFraction(wPhaseIdx, compIdx);
+ }
+ Scalar xlmax = 1.0;
+ if (sumxl > xlmax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ 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 " << dofIdxGlobal
+ << ", coordinated: " << globalPos << ", sumxl: "
+ << sumxl << std::endl;
+ newPhasePresence = bothPhases;
+ if (formulation == pgSl)
+ priVars[switchIdx] = 0.0;
+ else if (formulation == plSg)
+ priVars[switchIdx] = 1.0;
+ //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+ }
+ }
+ else if (phasePresence == wPhaseOnly)
+ {
+ Scalar xgmax = 1;
+ Scalar sumxg = 0;
+ //Calculate sum of mole fractions in the hypothetical gas phase
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ sumxg += volVars.moleFraction(nPhaseIdx, compIdx);
+ }
+ if (sumxg > xgmax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xgmax *=1.02;
+ //liquid phase appears if sum is larger than one
+ if (sumxg > xgmax)
+ {
+ std::cout << "Gas Phase appears at vertex " << dofIdxGlobal
+ << ", 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)
+ priVars[switchIdx] = 0.999;
+ else if (formulation == plSg)
+ priVars[switchIdx] = 0.001;
+
+ }
+ }
+ this->phasePresence_[dofIdxGlobal] = newPhasePresence;
+ this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+};
+
+} // end namespace dumux
+
+#endif
diff --git a/dumux/porousmediumflow/2pncmin/implicit/properties.hh b/dumux/porousmediumflow/2pncmin/implicit/properties.hh
index 2171d0bc83d893c6433d16a41e84b7a8428bc503..d024b672cba96849357931400f8aceba7d64b05e 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/properties.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/properties.hh
@@ -45,6 +45,11 @@ NEW_TYPE_TAG(TwoPNCMin, INHERITS_FROM(TwoPNC));
NEW_TYPE_TAG(BoxTwoPNCMin, INHERITS_FROM(BoxModel, TwoPNCMin));
NEW_TYPE_TAG(CCTwoPNCMin, INHERITS_FROM(CCModel, TwoPNCMin));
+//! The type tag for the non-isothermal two phase n component mineralisation problems
+NEW_TYPE_TAG(TwoPNCMinNI, INHERITS_FROM(TwoPNCMin, NonIsothermal));
+NEW_TYPE_TAG(BoxTwoPNCMinNI, INHERITS_FROM(BoxModel, TwoPNCMinNI));
+NEW_TYPE_TAG(CCTwoPNCMinNI, INHERITS_FROM(CCModel, TwoPNCMinNI));
+
//////////////////////////////////////////////////////////////////
// Property tags
//////////////////////////////////////////////////////////////////
diff --git a/dumux/porousmediumflow/2pncmin/implicit/propertydefaults.hh b/dumux/porousmediumflow/2pncmin/implicit/propertydefaults.hh
index b5228ad6100f5b92ec3f5a69485ef2921de76264..6f9ef2fab548587006fdac99158f3c53c0d09250 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/propertydefaults.hh
@@ -31,9 +31,9 @@
#include "indices.hh"
#include "model.hh"
#include "indices.hh"
-#include "fluxvariables.hh"
#include "volumevariables.hh"
#include "properties.hh"
+#include "primaryvariableswitch.hh"
#include <dumux/porousmediumflow/2pnc/implicit/newtoncontroller.hh>
#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
@@ -119,8 +119,8 @@ SET_TYPE_PROP(TwoPNCMin, Model, TwoPNCMinModel<TypeTag>);
//! the VolumeVariables property
SET_TYPE_PROP(TwoPNCMin, VolumeVariables, TwoPNCMinVolumeVariables<TypeTag>);
-//! the FluxVariables property
-SET_TYPE_PROP(TwoPNCMin, FluxVariables, TwoPNCMinFluxVariables<TypeTag>);
+//! the primary variable switch
+SET_TYPE_PROP(TwoPNCMin, PrimaryVariableSwitch, TwoPNCMinPrimaryVariableSwitch<TypeTag>);
//! The indices required by the isothermal 2pNcMin model
SET_TYPE_PROP(TwoPNCMin, Indices, TwoPNCMinIndices <TypeTag, /*PVOffset=*/0>);
@@ -136,6 +136,33 @@ SET_BOOL_PROP(TwoPNCMin, useSalinity, false);
// (Nield, Bejan, Convection in porous media, 2006, p. 10)
SET_SCALAR_PROP(TwoPNCMin, SpatialParamsForchCoeff, 0.55);
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(TwoPNCMinNI, ThermalConductivityModel)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+public:
+ using type = ThermalConductivitySomerton<Scalar, Indices>;
+};
+
+/*!
+ * \brief Set the property for the number of equations.
+ * For each component and each precipitated mineral/solid phase one equation has to
+ * be solved.
+ */
+SET_PROP(TwoPNCMinNI, IsothermalNumEq)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+public:
+ static const int value = FluidSystem::numComponents + FluidSystem::numSPhases;
+};
+
+//! The indices required by the isothermal 2pNcMin model
+SET_TYPE_PROP(TwoPNCMinNI, IsothermalIndices, TwoPNCMinIndices <TypeTag, /*PVOffset=*/0>);
+
}
}
diff --git a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
index e775435e054bbf4d34b6b088d151d734cde7a686..f681c635fbe0b23f4fe7b70a3de140f335836825 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
@@ -25,9 +25,6 @@
#ifndef DUMUX_2PNCMIN_VOLUME_VARIABLES_HH
#define DUMUX_2PNCMIN_VOLUME_VARIABLES_HH
-#include <vector>
-#include <iostream>
-
#include <dumux/common/math.hh>
#include <dumux/implicit/model.hh>
#include <dumux/material/fluidstates/compositional.hh>
@@ -50,19 +47,22 @@ namespace Dumux
template <class TypeTag>
class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
{
- typedef TwoPNCVolumeVariables<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ // base type is used for energy related quantites
+ using BaseType = ImplicitVolumeVariables<TypeTag>;
+
+ using ParentType = TwoPNCVolumeVariables<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLawParams);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum
{
@@ -99,17 +99,18 @@ class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
useSalinity = GET_PROP_VALUE(TypeTag, useSalinity)
};
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename Grid::ctype CoordScalar;
- typedef Dumux::Miscible2pNCComposition<Scalar, FluidSystem> Miscible2pNCComposition;
- typedef Dumux::ComputeFromReferencePhase2pNCMin<Scalar, FluidSystem> ComputeFromReferencePhase2pNCMin;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using CoordScalar = typename Grid::ctype;
+ using Miscible2pNCComposition = Dumux::Miscible2pNCComposition<Scalar, FluidSystem>;
+ using ComputeFromReferencePhase2pNCMin = Dumux::ComputeFromReferencePhase2pNCMin<Scalar, FluidSystem>;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
+
public:
- typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
/*!
* \copydoc ImplicitVolumeVariables::update
@@ -117,22 +118,15 @@ public:
void update(const PrimaryVariables &priVars,
const Problem &problem,
const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx,
- bool isOldSol)
+ const SubControlVolume& scv)
{
- ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
- completeFluidState(priVars, problem, element, fvGeometry, scvIdx, this->fluidState_, isOldSol);
+ // Update parent type (also completes the fluid state)
+ ParentType::update(priVars, problem, element, scv);
/////////////
// calculate the remaining quantities
/////////////
- // porosity evaluation
- initialPorosity_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx);
- minimumPorosity_ = problem.spatialParams().porosityMin(element, fvGeometry, scvIdx);
-
-
sumPrecipitates_ = 0.0;
for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
{
@@ -151,73 +145,39 @@ public:
// porosity value, as the porous media media properties change related to salt precipitation will not be
// accounted otherwise.
- // this->porosity_ = initialPorosity_ - sumPrecipitates_;
-
- this->porosity_ = std::max(minimumPorosity_, std::max(0.0, initialPorosity_ - sumPrecipitates_));
+ // overwrite porosity
+ // porosity evaluation
+ auto initialPorosity = problem.spatialParams().porosity(scv);
+ auto minimumPorosity = problem.spatialParams().porosityMin(scv);
+ this->porosity_ = std::max(minimumPorosity, std::max(0.0, initialPorosity - sumPrecipitates_));
salinity_= 0.0;
moleFractionSalinity_ = 0.0;
for (int compIdx = numMajorComponents; compIdx< numComponents; compIdx++) //sum of the mass fraction of the components
{
- if(this->fluidState_.moleFraction(wPhaseIdx, compIdx)> 0)
+ if(this->fluidState_.moleFraction(wPhaseIdx, compIdx) > 0)
{
salinity_+= this->fluidState_.massFraction(wPhaseIdx, compIdx);
moleFractionSalinity_ += this->fluidState_.moleFraction(wPhaseIdx, compIdx);
}
}
-
- // TODO/FIXME: Different relations for the porosoty-permeability changes are given here. We have to fins a way
- // so that one can select the relation form the input file.
-
- // kozeny-Carman relation
- permeabilityFactor_ = std::pow(((1-initialPorosity_)/(1-this->porosity_)), 2)
- * std::pow((this->porosity_/initialPorosity_), 3);
-
- // Verma-Pruess relation
- // permeabilityFactor_ = 100 * std::pow(((this->porosity_/initialPorosity_)-0.9),2);
-
- // Modified Fair-Hatch relation with final porosity set to 0.2 and E1=1
- // permeabilityFactor_ = std::pow((this->porosity_/initialPorosity_),3)
- // * std::pow((std::pow((1 - initialPorosity_),2/3))+(std::pow((0.2 - initialPorosity_),2/3)),2)
- // / std::pow((std::pow((1 -this->porosity_),2/3))+(std::pow((0.2 -this->porosity_),2/3)),2);
-
- //Timur relation with residual water saturation set to 0.001
- // permeabilityFactor_ = 0.136 * (std::pow(this->porosity_,4.4)) / (2000 * (std::pow(0.001,2)));
-
- //Timur relation1 with residual water saturation set to 0.001
- // permeabilityFactor_ = 0.136 * (std::pow(this->porosity_,4.4)) / (200000 * (std::pow(0.001,2)));
-
- // Bern. relation
- // permeabilityFactor_ = std::pow((this->porosity_/initialPorosity_),8);
-
- //Tixier relation with residual water saturation set to 0.001
- // permeabilityFactor_ = (std::pow((250 * (std::pow(this->porosity_,3)) / 0.001),2)) / initialPermeability_;
-
- //Coates relation with residual water saturation set to 0.001
- // permeabilityFactor_ = (std::pow((100 * (std::pow(this->porosity_,2)) * (1-0.001) / 0.001,2))) / initialPermeability_ ;
-
- // energy related quantities not contained in the fluid state
- //asImp_().updateEnergy_(priVars, problem,element, fvGeometry, scvIdx, isOldSol);
}
- /*!
- * \copydoc ImplicitModel::completeFluidState
- * \param isOldSol Specifies whether this is the previous solution or the current one
- */
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ * \param isOldSol Specifies whether this is the previous solution or the current one
+ */
static void completeFluidState(const PrimaryVariables& priVars,
const Problem& problem,
const Element& element,
- const FVElementGeometry& fvGeometry,
- int scvIdx,
- FluidState& fluidState,
- bool isOldSol = false)
+ const SubControlVolume& scv,
+ FluidState& fluidState)
{
- Scalar t = Implementation::temperature_(priVars, problem, element,fvGeometry, scvIdx);
+ Scalar t = BaseType::temperature(priVars, problem, element, scv);
fluidState.setTemperature(t);
- int dofIdxGlobal = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
- int phasePresence = problem.model().phasePresence(dofIdxGlobal, isOldSol);
+ auto phasePresence = problem.model().priVarSwitch().phasePresence(scv.dofIndex());
/////////////
// set the saturations
@@ -242,7 +202,10 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
}
else
+ {
DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ }
+
fluidState.setSaturation(nPhaseIdx, Sg);
fluidState.setSaturation(wPhaseIdx, 1.0 - Sg);
@@ -251,19 +214,32 @@ public:
/////////////
// calculate capillary pressure
- const MaterialLawParams &materialParams = problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
- Scalar pc = MaterialLaw::pc(materialParams, 1 - Sg);
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv);
+ auto pc = MaterialLaw::pc(materialParams, 1 - Sg);
// extract the pressures
- if (formulation == plSg) {
+ if (formulation == plSg)
+ {
fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
+ if (priVars[pressureIdx] + pc < 0.0)
+ DUNE_THROW(Dumux::NumericalProblem, "Capillary pressure is too low");
fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
}
- else if (formulation == pgSl) {
+ else if (formulation == pgSl)
+ {
fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
+ // Here we check for (p_g - pc) in order to ensure that (p_l > 0)
+ if (priVars[pressureIdx] - pc < 0.0)
+ {
+ std::cout<< " The gas pressure is p_g = "<< priVars[pressureIdx]<<", the capillary pressure p_c = "<< pc << std::endl;
+ DUNE_THROW(Dumux::NumericalProblem, "Capillary pressure is too high");
+ }
fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
}
- else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+ }
/////////////
// calculate the phase compositions
@@ -287,10 +263,10 @@ public:
}
Miscible2pNCComposition::solve(fluidState,
- paramCache,
- wPhaseIdx, //known phaseIdx
- /*setViscosity=*/true,
- /* =*/false);
+ paramCache,
+ wPhaseIdx, //known phaseIdx
+ /*setViscosity=*/true,
+ /*setInternalEnergy=*/false);
}
else if (phasePresence == nPhaseOnly)
{
@@ -301,24 +277,23 @@ public:
for (int compIdx=0; compIdx<numComponents; ++compIdx)
{
fugCoeffL[compIdx] = FluidSystem::fugacityCoefficient(fluidState,
- paramCache,
- wPhaseIdx,
- compIdx);
+ paramCache,
+ wPhaseIdx,
+ compIdx);
fugCoeffG[compIdx] = FluidSystem::fugacityCoefficient(fluidState,
- paramCache,
- nPhaseIdx,
- compIdx);
+ paramCache,
+ nPhaseIdx,
+ compIdx);
}
for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
moleFrac[compIdx] = (priVars[compIdx]*fugCoeffL[compIdx]*fluidState.pressure(wPhaseIdx))
- /(fugCoeffG[compIdx]*fluidState.pressure(nPhaseIdx));
+ /(fugCoeffG[compIdx]*fluidState.pressure(nPhaseIdx));
moleFrac[wCompIdx] = priVars[switchIdx];
Scalar sumMoleFracNotGas = 0;
for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
- {
sumMoleFracNotGas+=moleFrac[compIdx];
- }
+
sumMoleFracNotGas += moleFrac[wCompIdx];
moleFrac[nCompIdx] = 1 - sumMoleFracNotGas;
@@ -338,10 +313,9 @@ public:
/*setViscosity=*/true,
/*setEnthalpy=*/false);
- }
+ }
else if (phasePresence == wPhaseOnly)
{
-
// only the liquid phase is present, i.e. liquid phase
// composition is stored explicitly.
// extract _mass_ fractions in the gas phase
@@ -379,7 +353,12 @@ public:
paramCache.updateAll(fluidState);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+ Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+ Scalar h = BaseType::enthalpy(fluidState, paramCache, phaseIdx);
+
+ fluidState.setDensity(phaseIdx, rho);
+ fluidState.setViscosity(phaseIdx, mu);
fluidState.setEnthalpy(phaseIdx, h);
}
}
@@ -396,38 +375,24 @@ public:
* \brief Returns the inital porosity of the
* pure, precipitate-free porous medium
*/
- Scalar initialPorosity() const
- { return initialPorosity_;}
-
- /*!
- * \brief Returns the inital permeability of the
- * pure, precipitate-free porous medium
- */
- Scalar initialPermeability() const
- { return initialPermeability_;}
-
- /*!
- * \brief Returns the factor for the reduction of the initial permeability
- * due precipitates in the porous medium
- */
- Scalar permeabilityFactor() const
- { return permeabilityFactor_; }
-
- /*!
- * \brief Returns the mole fraction of the salinity in the liquid phase
- */
- Scalar moleFracSalinity() const
- {
- return moleFractionSalinity_;
- }
-
- /*!
- * \brief Returns the salinity (mass fraction) in the liquid phase
- */
- Scalar salinity() const
- {
- return salinity_;
- }
+ // Scalar initialPorosity() const
+ // { return initialPorosity_;}
+
+ // /*!
+ // * \brief Returns the mole fraction of the salinity in the liquid phase
+ // */
+ // Scalar moleFracSalinity() const
+ // {
+ // return moleFractionSalinity_;
+ // }
+
+ // /*!
+ // * \brief Returns the salinity (mass fraction) in the liquid phase
+ // */
+ // Scalar salinity() const
+ // {
+ // return salinity_;
+ // }
/*!
* \brief Returns the density of the phase for all fluid and solid phases
@@ -470,56 +435,18 @@ public:
* phase is equal to the phaseIdx
*/
Scalar molality(int phaseIdx, int compIdx) const // [moles/Kg]
- { return this->fluidState_.moleFraction(phaseIdx, compIdx)
- /(fluidState_.moleFraction(phaseIdx, phaseIdx)
- * FluidSystem::molarMass(phaseIdx));}
-
-protected:
- friend class TwoPNCVolumeVariables<TypeTag>;
- static Scalar temperature_(const PrimaryVariables &priVars,
- const Problem& problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx)
{
- return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
+ return this->fluidState_.moleFraction(phaseIdx, compIdx)
+ /(this->fluidState_.moleFraction(phaseIdx, phaseIdx)
+ * FluidSystem::molarMass(phaseIdx));
}
- template<class ParameterCache>
- static Scalar enthalpy_(const FluidState& fluidState,
- const ParameterCache& paramCache,
- int phaseIdx)
- {
- return 0;
- }
-
- /*!
- * \brief Update all quantities for a given control volume.
- *
- * \param priVars The solution primary variables
- * \param problem The problem
- * \param element The element
- * \param fvGeometry Evaluate function with solution of current or previous time step
- * \param scvIdx The local index of the SCV (sub-control volume)
- * \param isOldSol Evaluate function with solution of current or previous time step
- */
- void updateEnergy_(const PrimaryVariables &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx,
- bool isOldSol)
- {};
+protected:
Scalar precipitateVolumeFraction_[numSPhases];
- Scalar permeabilityFactor_;
- Scalar initialPorosity_;
- Scalar initialPermeability_;
- Scalar minimumPorosity_;
Scalar sumPrecipitates_;
Scalar salinity_;
Scalar moleFractionSalinity_;
- FluidState fluidState_;
private:
Implementation &asImp_()
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
index 2c3bb297725c19db78bc2b5b3619aa9b32d44ce0..c08f541cc8c7bb36c6e80ff52a129e23a0d38973 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
@@ -24,6 +24,7 @@
#ifndef DUMUX_DISSOLUTION_PROBLEM_HH
#define DUMUX_DISSOLUTION_PROBLEM_HH
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/2pncmin/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/brineair.hh>
@@ -40,7 +41,7 @@ namespace Properties
{
NEW_TYPE_TAG(DissolutionProblem, INHERITS_FROM(TwoPNCMin, DissolutionSpatialparams));
NEW_TYPE_TAG(DissolutionBoxProblem, INHERITS_FROM(BoxModel, DissolutionProblem));
-NEW_TYPE_TAG(DissolutionCCProblem, INHERITS_FROM(CCModel, DissolutionProblem));
+NEW_TYPE_TAG(DissolutionCCProblem, INHERITS_FROM(CCTpfaModel, DissolutionProblem));
// Set the grid type
SET_TYPE_PROP(DissolutionProblem, Grid, Dune::YaspGrid<2>);
@@ -51,8 +52,8 @@ SET_TYPE_PROP(DissolutionProblem, Problem, DissolutionProblem<TypeTag>);
// Set fluid configuration
SET_PROP(DissolutionProblem, FluidSystem)
{
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef FluidSystems::BrineAir<Scalar, H2O<Scalar>, true/*useComplexrelations=*/> type;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::BrineAir<Scalar, H2O<Scalar>, true/*useComplexrelations=*/>;
};
// Set the spatial parameters
@@ -62,7 +63,7 @@ SET_TYPE_PROP(DissolutionProblem, SpatialParams, DissolutionSpatialparams<TypeTa
SET_BOOL_PROP(DissolutionProblem, ProblemEnableGravity, true);
//Set properties here to override the default property settings in the model.
-SET_INT_PROP(DissolutionProblem, ReplaceCompEqIdx, 1);
+SET_INT_PROP(DissolutionProblem, ReplaceCompEqIdx, 1); //! Replace gas balance by total mass balance
SET_INT_PROP(DissolutionProblem, Formulation, TwoPNCFormulation::pgSl);
}
@@ -87,12 +88,12 @@ SET_INT_PROP(DissolutionProblem, Formulation, TwoPNCFormulation::pgSl);
template <class TypeTag>
class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
{
- typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using ParentType = ImplicitPorousMediaProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
@@ -128,20 +129,20 @@ class DissolutionProblem : public ImplicitPorousMediaProblem<TypeTag>
};
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::Intersection Intersection;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<dim>::Entity;
+ using Intersection = typename GridView::Intersection;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
DissolutionProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ : ParentType(timeManager, gridView)
{
outerSalinity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterSalinity);
@@ -165,9 +166,6 @@ public:
temperatureLow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureLow);
temperatureHigh_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureHigh);
name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
- freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
- storageLastTimestep_ = 0.0;
- lastMassOutputTime_ = 0.0;
outfile.open("evaporation.out");
outfile << "time; evaporationRate" << std::endl;
@@ -203,7 +201,7 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const std::string name() const
+ const std::string& name() const
{ return name_; }
/*!
@@ -223,9 +221,11 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
- void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
- const Scalar rmax = this->bBoxMax()[0]; // outerRadius_;
+ BoundaryTypes bcTypes;
+
+ const Scalar rmax = this->bBoxMax()[0];
const Scalar rmin = this->bBoxMin()[0];
// default to Neumann
@@ -238,70 +238,80 @@ public:
// Constant pressure at well (Dirichlet condition)
if(globalPos[0] < rmin + eps_)
bcTypes.setAllDirichlet();
+
+ return bcTypes;
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* boundary segment.
- *
- * For this method, the \a values parameter stores primary variables.
*/
- void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
+ PrimaryVariables priVars(0.0);
+
const Scalar rmax = this->bBoxMax()[0];
const Scalar rmin = this->bBoxMin()[0];
if(globalPos[0] > rmax - eps_)
{
- values[pressureIdx] = outerPressure_ ; // Outer boundary pressure bar
- values[switchIdx] = outerLiqSaturation_; // Saturation outer boundary
- values[xlNaClIdx] = massTomoleFrac_(outerSalinity_);// mole fraction salt
- values[precipNaClIdx] = 0.0;// precipitated salt
+ priVars[pressureIdx] = outerPressure_ ; // Outer boundary pressure bar
+ priVars[switchIdx] = outerLiqSaturation_; // Saturation outer boundary
+ priVars[xlNaClIdx] = massToMoleFrac_(outerSalinity_);// mole fraction salt
+ priVars[precipNaClIdx] = 0.0;// precipitated salt
}
if(globalPos[0] < rmin + eps_)
{
- values[pressureIdx] = innerPressure_ ; // Inner boundary pressure bar
- values[switchIdx] = innerLiqSaturation_; // Saturation inner boundary
- values[xlNaClIdx] = massTomoleFrac_(innerSalinity_);// mole fraction salt
- values[precipNaClIdx] = 0.0;// precipitated salt
+ priVars[pressureIdx] = innerPressure_ ; // Inner boundary pressure bar
+ priVars[switchIdx] = innerLiqSaturation_; // Saturation inner boundary
+ priVars[xlNaClIdx] = massToMoleFrac_(innerSalinity_);// mole fraction salt
+ priVars[precipNaClIdx] = 0.0;// precipitated salt
}
+
+ return priVars;
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
- * For this method, the \a values parameter stores the mass flux
+ * For this method, the \a priVars parameter stores the mass flux
* in normal direction of each component. Negative values mean
* influx.
*/
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &is,
- int scvIdx,
- int boundaryFaceIdx) const
+ PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
{
- values = 0.0;
+ return PrimaryVariables(0.0);
}
/*!
* \brief Evaluate the initial value for a control volume.
*
+ * \param values The initial values for the primary variables
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scvIdx The local subcontrolvolume index
+ *
* For this method, the \a values parameter stores primary
* variables.
*/
- void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables initial(const SubControlVolume &scv) const
{
- values[pressureIdx] = reservoirPressure_;
- values[switchIdx] = initLiqSaturation_; // Sl primary variable
- values[xlNaClIdx] = massTomoleFrac_(outerSalinity_); // mole fraction
+ PrimaryVariables priVars(0.0);
+
+ const auto& globalPos = scv.dofPosition();
+
+ priVars[pressureIdx] = reservoirPressure_;
+ priVars[switchIdx] = initLiqSaturation_; // Sl primary variable
+ priVars[xlNaClIdx] = massToMoleFrac_(outerSalinity_); // mole fraction
if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 20.0 - eps_)
- values[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
+ priVars[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
else
- values[precipNaClIdx] = initPrecipitatedSalt1_; // [kg/m^3]
+ priVars[precipNaClIdx] = initPrecipitatedSalt1_; // [kg/m^3]
+
+ return priVars;
}
/*!
@@ -313,58 +323,67 @@ public:
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
- * For this method, the \a values parameter stores the rate mass
- * of a component is generated or annihilate per volume
- * unit. Positive values mean that mass is created, negative ones
- * mean that it vanishes.
+ * This is the method for the case where the source term is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param values The source and sink values for the conservation equations in units of
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scv The subcontrolvolume
+ *
+ * For this method, the \a values parameter stores the conserved quantity rate
+ * generated or annihilate per volume unit. Positive values mean
+ * that the conserved quantity is created, negative ones mean that it vanishes.
+ * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
- void solDependentSource(PrimaryVariables &source,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx,
- const ElementVolumeVariables &elemVolVars) const
+ PrimaryVariables source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
{
- source = 0;
- const auto& volVars = elemVolVars[scvIdx];
+ PrimaryVariables source(0.0);
+
+ const auto& volVars = elemVolVars[scv];
+
Scalar moleFracNaCl_lPhase = volVars.moleFraction(wPhaseIdx, NaClIdx);
Scalar moleFracNaCl_gPhase = volVars.moleFraction(nPhaseIdx, NaClIdx);
- Scalar massFracNaCl_Max_lPhase = this->spatialParams().SolubilityLimit();
- Scalar moleFracNaCl_Max_lPhase = massTomoleFrac_(massFracNaCl_Max_lPhase);
+ Scalar massFracNaCl_Max_lPhase = this->spatialParams().solubilityLimit();
+ Scalar moleFracNaCl_Max_lPhase = massToMoleFrac_(massFracNaCl_Max_lPhase);
Scalar moleFracNaCl_Max_gPhase = moleFracNaCl_Max_lPhase / volVars.pressure(nPhaseIdx);
- Scalar saltPorosity = this->spatialParams().porosityMin(element, fvGeometry, scvIdx);
+ Scalar saltPorosity = this->spatialParams().porosityMin(scv);
// liquid phase
Scalar precipSalt = volVars.porosity() * volVars.molarDensity(wPhaseIdx)
* volVars.saturation(wPhaseIdx)
- * std::abs(moleFracNaCl_lPhase - moleFracNaCl_Max_lPhase);
-
- if (moleFracNaCl_lPhase < moleFracNaCl_Max_lPhase)
- precipSalt *= -1;
+ * (moleFracNaCl_lPhase - moleFracNaCl_Max_lPhase);
// gas phase
- if (moleFracNaCl_gPhase > moleFracNaCl_Max_gPhase)
- precipSalt += volVars.porosity() * volVars.molarDensity(nPhaseIdx)
- * volVars.saturation(nPhaseIdx)
- * std::abs(moleFracNaCl_gPhase - moleFracNaCl_Max_gPhase);
+ precipSalt += volVars.porosity() * volVars.molarDensity(nPhaseIdx)
+ * volVars.saturation(nPhaseIdx)
+ * (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)
- precipSalt = - volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)/this->timeManager().timeStepSize();
+ precipSalt = -volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)/this->timeManager().timeStepSize();
- if (volVars.precipitateVolumeFraction(sPhaseIdx) >= volVars.initialPorosity() - saltPorosity && precipSalt > 0)
+ if (volVars.precipitateVolumeFraction(sPhaseIdx) >= this->spatialParams().porosity(scv) - saltPorosity && precipSalt > 0)
precipSalt = 0;
source[conti0EqIdx + NaClIdx] += -precipSalt;
source[precipNaClEqIdx] += precipSalt;
- Valgrind::CheckDefined(source);
+ return source;
}
+
/*!
* \brief Return the initial phase state inside a control volume.
+ *
+ * \param scv The sub control volume
*/
- int initialPhasePresence(const Element& element,
- const FVElementGeometry& fvGeometry,
- int scvIdx) const
+ int initialPhasePresence(const SubControlVolume& scv) const
{
return bothPhases;
}
@@ -376,7 +395,7 @@ private:
*
* \param XlNaCl the XlNaCl [kg NaCl / kg solution]
*/
- static Scalar massTomoleFrac_(Scalar XlNaCl)
+ static Scalar massToMoleFrac_(Scalar XlNaCl)
{
const Scalar Mw = 18.015e-3; /* molecular weight of water [kg/mol] */
const Scalar Ms = 58.44e-3; /* molecular weight of NaCl [kg/mol] */
@@ -389,9 +408,6 @@ private:
int nTemperature_;
int nPressure_;
- int freqMassOutput_;
- PrimaryVariables storageLastTimestep_;
- Scalar lastMassOutputTime_;
std::string name_;
Scalar pressureLow_, pressureHigh_;
@@ -412,6 +428,7 @@ private:
std::ofstream outfile;
};
-} //end namespace
+
+} //end namespace Dumux
#endif
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
index e0e993dea4aae8870f520a319b9fdcb4b0c2866c..6f7683be3392f307a5f0cc9f7ac8b5cbf204786b 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
@@ -44,12 +44,13 @@ SET_PROP(DissolutionSpatialparams, MaterialLaw)
{
private:
// define the material law which is parameterized by effective saturations
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
// define the material law parameterized by absolute saturations
- typedef EffToAbsLaw<RegularizedBrooksCorey<Scalar> > type;
+ using type = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
};
-}
+
+} // end namespace Properties
/**
* \brief Definition of the spatial parameters for the brine-co2 problem
@@ -58,36 +59,39 @@ public:
template<class TypeTag>
class DissolutionSpatialparams : public ImplicitSpatialParams<TypeTag>
{
- typedef ImplicitSpatialParams<TypeTag> ParentType;
- 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 GridView::ctype CoordScalar;
+ using ParentType = ImplicitSpatialParams<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLawParams);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using CoordScalar = typename GridView::ctype;
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld,
};
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
wPhaseIdx = FluidSystem::wPhaseIdx,
nPhaseIdx = FluidSystem::nPhaseIdx,
};
- typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
- typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GridView::template Codim<0>::Entity Element;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+ using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using Element = typename GridView::template Codim<0>::Entity;
public:
- DissolutionSpatialparams(const GridView &gridView)
- : ParentType(gridView), K_(0.0)
+ DissolutionSpatialparams(const Problem& problem, const GridView &gridView)
+ : ParentType(problem, gridView), initialPermeability_(0.0)
{
// 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;
+ initialPermeability_[i][i] = 2.23e-14;
// residual saturations
materialParams_.setSwr(0.2);
@@ -105,15 +109,22 @@ public:
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
- * Alternatively, the function intrinsicPermeabilityAtPos(const GlobalPosition& globalPos)
- * could be defined, where globalPos is the vector including the global coordinates
- * of the finite volume.
+ * Solution dependent permeability function
*/
- const Tensor& intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Tensor solDependentIntrinsicPermeability(const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
{
- return K_;
+ // Kozeny-Carman relation
+ auto initialPorosity = porosity(scv);
+ auto factor = std::pow((1-initialPorosity)/(1-volVars.porosity()), 2)
+ * std::pow(volVars.porosity()/initialPorosity, 3);
+
+ auto perm = initialPermeability_;
+
+ for (int i = 0; i < dim; i++)
+ perm[i][i] *= factor;
+
+ return perm;
}
/*!
@@ -124,12 +135,10 @@ public:
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
- Scalar porosityMin(const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
- {
+ Scalar porosityMin(const SubControlVolume &scv) const
+ {
return 1e-5;
- }
+ }
/*!
* \brief Define the minimum porosity \f$[-]\f$ after clogging caused by mineralization
@@ -139,30 +148,24 @@ public:
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
- Scalar porosity(const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
- {
+ Scalar porosity(const SubControlVolume &scv) const
+ {
return 0.11;
- }
+ }
- Scalar solidity(const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
+ Scalar solidity(const SubControlVolume &scv) const
{
- return 1.0 - porosity(element, fvGeometry, scvIdx);
+ return 1.0 - porosity(scv);
}
- Scalar SolubilityLimit() const
+ Scalar solubilityLimit() const
{
return 0.26;
}
- Scalar theta(const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
+ Scalar theta(const SubControlVolume &scv) const
{
return 10.0;
}
@@ -170,17 +173,16 @@ public:
// return the brooks-corey context depending on the position
const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
+ const SubControlVolume &scv) const
{
return materialParams_;
}
private:
- Tensor K_;
+ Tensor initialPermeability_;
MaterialLawParams materialParams_;
};
-}
+} // end namespace Dumux
#endif