From 5f3b6fdd088a812f8e224eefe474f0a05709378e Mon Sep 17 00:00:00 2001
From: Gabriele Seitz <gabriele.seitz@iws.uni-stuttgart.de>
Date: Fri, 23 Jun 2017 17:12:13 +0200
Subject: [PATCH] [1pncmin] add model to next
not yet working
---
dumux/porousmediumflow/1pncmin/CMakeLists.txt | 1 +
.../1pncmin/implicit/CMakeLists.txt | 10 +
.../1pncmin/implicit/indices.hh | 48 ++
.../1pncmin/implicit/localresidual.hh | 115 +++++
.../1pncmin/implicit/model.hh | 478 ++++++++++++++++++
.../1pncmin/implicit/properties.hh | 69 +++
.../1pncmin/implicit/propertydefaults.hh | 180 +++++++
.../1pncmin/implicit/volumevariables.hh | 390 ++++++++++++++
dumux/porousmediumflow/CMakeLists.txt | 1 +
9 files changed, 1292 insertions(+)
create mode 100644 dumux/porousmediumflow/1pncmin/CMakeLists.txt
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/CMakeLists.txt
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/indices.hh
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/localresidual.hh
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/model.hh
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/properties.hh
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
create mode 100644 dumux/porousmediumflow/1pncmin/implicit/volumevariables.hh
diff --git a/dumux/porousmediumflow/1pncmin/CMakeLists.txt b/dumux/porousmediumflow/1pncmin/CMakeLists.txt
new file mode 100644
index 0000000000..ba8341c614
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/dumux/porousmediumflow/1pncmin/implicit/CMakeLists.txt b/dumux/porousmediumflow/1pncmin/implicit/CMakeLists.txt
new file mode 100644
index 0000000000..92b891fb7b
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+localresidual.hh
+model.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1pncmin/implicit)
diff --git a/dumux/porousmediumflow/1pncmin/implicit/indices.hh b/dumux/porousmediumflow/1pncmin/implicit/indices.hh
new file mode 100644
index 0000000000..4c5a2085c1
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/indices.hh
@@ -0,0 +1,48 @@
+// -*- 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 Defines the indices required for the one-phase n-component mineralization
+ * fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_INDICES_HH
+#define DUMUX_1PNCMIN_INDICES_HH
+
+#include <dumux/porousmediumflow/1pnc/implicit/indices.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCMinModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal one-phase n-component mineralization model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+ class OnePNCMinIndices: public OnePNCIndices<TypeTag, PVOffset>
+{
+};
+
+// \}
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pncmin/implicit/localresidual.hh b/dumux/porousmediumflow/1pncmin/implicit/localresidual.hh
new file mode 100644
index 0000000000..0a9cdae818
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/localresidual.hh
@@ -0,0 +1,115 @@
+// -*- 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the one-phase n-component mineralisation box model.
+ */
+
+#ifndef DUMUX_1PNCMIN_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_1PNCMIN_LOCAL_RESIDUAL_BASE_HH
+
+#include "properties.hh"
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCMinModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the one-phase n-component mineralization fully implicit model.
+ *
+ * This class is used to fill the gaps in ImplicitLocalResidual for the one-phase n-component flow.
+ */
+template<class TypeTag>
+class OnePNCMinLocalResidual: public OnePNCLocalResidual<TypeTag>
+{
+protected:
+// typedef OnePNCLocalResidual<TypeTag> ParentType;
+ using ParentType = CompositionalLocalResidual<TypeTag>;
+ using ThisType = TwoPNCMinLocalResidual<TypeTag>;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables)
+// 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, BoundaryTypes) BoundaryTypes;
+// typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+
+ 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),
+
+ pressureIdx = Indices::pressureIdx,
+ firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+ phaseIdx = Indices::phaseIdx,
+ sPhaseIdx = 1, // don't use the sPhaseIdx of the fluidsystem
+
+ conti0EqIdx = Indices::conti0EqIdx,
+ };
+
+// typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+// typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+// typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+// typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume).
+ * In contrast to the 1pnc model, here, the storage of solid phases is included too.
+ *
+ * \param scv the SCV (sub-control-volume)
+ * \param volVars The volume variables of the right time step
+ */
+ PrimaryVariables computeStorage(const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ //call parenttype function
+ auto storage = ParentType::computeStorage(scv, volVars);
+
+ // Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
+ for (int Idx = sPhaseIdx; Idx < numPhases + numSPhases; ++Idx)
+ {
+ auto eqIdx = conti0EqIdx + numComponents-numPhases + Idx;
+ storage[eqIdx] += volVars.precipitateVolumeFraction(Idx)*volVars.molarDensity(Idx);
+ }
+
+ return storage;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/1pncmin/implicit/model.hh b/dumux/porousmediumflow/1pncmin/implicit/model.hh
new file mode 100644
index 0000000000..419fbec093
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/model.hh
@@ -0,0 +1,478 @@
+// -*- 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 Adaption of the fully implicit box scheme to the two-phase n-component flow model.
+*/
+
+#ifndef DUMUX_1PNCMIN_MODEL_HH
+#define DUMUX_1PNCMIN_MODEL_HH
+
+#include "properties.hh"
+#include "indices.hh"
+#include "localresidual.hh"
+
+// #include <dumux/material/constants.hh>
+#include <dumux/porousmediumflow/1pnc/implicit/model.hh>
+#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCMinModel
+ * \brief Adaption of the fully implicit scheme to the
+ * two-phase n-component fully implicit model with additional solid/mineral phases.
+ *
+ * This model implements two-phase n-component flow of two compressible and
+ * partially miscible fluids \f$\alpha \in \{ w, n \}\f$ composed of the n components
+ * \f$\kappa \in \{ w, n,\cdots \}\f$ in combination with mineral precipitation and dissolution.
+ * The solid phases. The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum:
+ * \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ * \f]
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, one gets one transport equation for each component
+ * \f{eqnarray}
+ && \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa \phi S_\alpha )}
+ {\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
+ \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ (\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}
+ \nonumber \\ \nonumber \\
+ &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a,\cdots \} \, ,
+ \alpha \in \{w, g\}
+ \f}
+ *
+ * The solid or mineral phases are assumed to consist of a single component.
+ * Their mass balance consist only of a storage and a source term:
+ * \f$\frac{\partial \varrho_\lambda \phi_\lambda )} {\partial t}
+ * = q_\lambda\f$
+ *
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme (this is not done for 2pnc approach yet, however possible) as
+ * spatial and the implicit Euler method as time discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
+ * unknowns can be reduced to number of components.
+ *
+ * The used primary variables are, like in the two-phase model, either \f$p_w\f$ and \f$S_n\f$
+ * or \f$p_n\f$ and \f$S_w\f$. The formulation which ought to be used can be
+ * specified by setting the <tt>Formulation</tt> property to either
+ * TwoPTwoCIndices::pWsN or TwoPTwoCIndices::pNsW. By
+ * default, the model uses \f$p_w\f$ and \f$S_n\f$.
+ *
+ * Moreover, the second primary variable depends on the phase state, since a
+ * primary variable switch is included. The phase state is stored for all nodes
+ * of the system. The model is uses mole fractions.
+ *Following cases can be distinguished:
+ * <ul>
+ * <li> Both phases are present: The saturation is used (either \f$S_n\f$ or \f$S_w\f$, dependent on the chosen <tt>Formulation</tt>),
+ * as long as \f$ 0 < S_\alpha < 1\f$</li>.
+ * <li> Only wetting phase is present: The mole fraction of, e.g., air in the wetting phase \f$x^a_w\f$ is used,
+ * as long as the maximum mole fraction is not exceeded (\f$x^a_w<x^a_{w,max}\f$)</li>
+ * <li> Only non-wetting phase is present: The mole fraction of, e.g., water in the non-wetting phase, \f$x^w_n\f$, is used,
+ * as long as the maximum mole fraction is not exceeded (\f$x^w_n<x^w_{n,max}\f$)</li>
+ * </ul>
+ *
+ * For the other components, the mole fraction \f$x^\kappa_w\f$ is the primary variable.
+ * The primary variable of the solid phases is the volume fraction \f$\phi_\lambda = \frac{V_\lambda}{V_{total}}\f$.
+ *
+ * The source an sink terms link the mass balances of the n-transported component to the solid phases.
+ * The porosity \f$\phi\f$ is updated according to the reduction of the initial (or solid-phase-free porous medium) porosity \f$\phi_0\f$
+ * by the accumulated volume fractions of the solid phases:
+ * \f$ \phi = \phi_0 - \sum (\phi_\lambda)\f$
+ * Additionally, the permeability is updated depending on the current porosity.
+ */
+
+template<class TypeTag>
+class OnePNCMinModel: public OnePNCModel<TypeTag>
+{
+ typedef Dumux::OnePNCMinModel<TypeTag> ThisType;
+ typedef Dumux::OnePNCModel<TypeTag> ParentType;
+
+ 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 Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ //old
+// typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+// typedef Dumux::Constants<Scalar> Constant;
+// typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ 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>;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+ numSComponents = FluidSystem::numSComponents,
+
+ pressureIdx = Indices::pressureIdx,
+ firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+ phaseIdx = Indices::phaseIdx,
+ };
+
+ //old
+// typedef typename GridView::template Codim<dim>::Entity Vertex;
+// typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+
+ /*!
+ * \brief Apply the initial conditions to the model.
+ *
+ * \param problem The object representing the problem which needs to
+ * be simulated.
+ */
+ void init(Problem& problem)
+ {
+ ParentType::init(problem);
+
+ // register standardized vtk output fields
+ auto& vtkOutputModule = problem.vtkOutputModule();
+ vtkOutputModule.addSecondaryVariable("p", [](const VolumeVariables& v){ return v.pressure(phaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(phaseIdx); });
+ vtkOutputModule.addSecondaryVariable("rho", [](const VolumeVariables& v){ return v.density(phaseIdx); });
+ vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+ vtkOutputModule.addSecondaryVariable("permeabilityFactor", [](const VolumeVariables& v){ return v.permeabilityFactor(); });
+
+ vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+
+ for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+ vtkOutputModule.addSecondaryVariable("precipitateVolumeFraction_" + FluidSystem::phaseName(numPhases + sPhaseIdx),
+ [sPhaseIdx](const VolumeVariables& v)
+ { return v.precipitateVolumeFraction(numPhases + sPhaseIdx); });
+
+ vtkOutputModule.addSecondaryVariable("Kxx",
+ [this](const VolumeVariables& v){ return this->perm_(v.permeability())[0][0]; });
+ if (dim >= 2)
+ vtkOutputModule.addSecondaryVariable("Kyy",
+ [this](const VolumeVariables& v){ return this->perm_(v.permeability())[1][1]; });
+ if (dim >= 3)
+ vtkOutputModule.addSecondaryVariable("Kzz",
+ [this](const VolumeVariables& v){ return this->perm_(v.permeability())[2][2]; });
+
+ for (int i = 0; i < numComponents; ++i)
+ vtkOutputModule.addSecondaryVariable("x_" + FluidSystem::componentName(i),[i](const VolumeVariables& v){ return v.moleFraction(phaseIdx,i); });
+
+ for (int i = 0; i < numComponents; ++i)
+ vtkOutputModule.addSecondaryVariable("m^w_" + FluidSystem::componentName(i),
+ [i](const VolumeVariables& v){ return v.molarity(phaseIdx,i); });
+ }
+
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ *
+ * \param sol The solution vector
+ * \param writer The writer for multi-file VTK datasets
+ */
+// template<class MultiWriter>
+// //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;
+//
+// // get the number of degrees of freedom
+// auto numDofs = this->numDofs();
+// ScalarField *p = writer.allocateManagedBuffer (numDofs);
+// ScalarField *rho = writer.allocateManagedBuffer (numDofs);
+// ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
+// ScalarField *poro = writer.allocateManagedBuffer (numDofs);
+// ScalarField *permeabilityFactor = writer.allocateManagedBuffer (numDofs);
+// ScalarField *precipitateVolumeFraction[numSPhases];
+//
+// for (int i = 0; i < numSPhases; ++i)
+// precipitateVolumeFraction[i] = writer.allocateManagedBuffer(numDofs);
+//
+// // 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];
+// for (int j = 0; j < numComponents ; ++j)
+// molarity[j] = writer.allocateManagedBuffer(numDofs);
+//
+// ScalarField *moleFraction[numComponents];
+// for (int j = 0; j < numComponents; ++j)
+// moleFraction[j] = writer.allocateManagedBuffer(numDofs);
+//
+// ScalarField *moleFractionSolid[numSComponents];
+// for (int j = numComponents; j < numSComponents; ++j)
+// moleFractionSolid[j] = writer.allocateManagedBuffer(numDofs);
+//
+// ScalarField *Perm[dim];
+// for (int j = 0; j < dim; ++j) //Permeability only in main directions xx and yy
+// Perm[j] = writer.allocateManagedBuffer(numDofs);
+//
+// VectorField *velocity = 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)
+// {
+// (*velocity)[i] = Scalar(0);
+// }
+// }
+//
+// auto numElements = this->gridView_().size(0);
+// ScalarField *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? */);
+//
+// for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+// {
+// auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+//
+// (*p)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
+// (*rho)[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
+// (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+//
+// for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx){
+// (*precipitateVolumeFraction[sPhaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].precipitateVolumeFraction(sPhaseIdx + numPhases);
+// }
+//
+// (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+// (*permeabilityFactor)[dofIdxGlobal] = elemVolVars[scvIdx].permeabilityFactor();
+//
+// // for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+// // (*massFraction[phaseIdx][compIdx])[dofIdxGlobal]= elemVolVars[scvIdx].massFraction(phaseIdx,compIdx);
+//
+// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+// (*molarity[compIdx])[dofIdxGlobal] = (elemVolVars[scvIdx].molarity(compIdx));
+//
+// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+// (*moleFraction[compIdx])[dofIdxGlobal]= elemVolVars[scvIdx].moleFraction(compIdx);
+//
+// for (int compIdx = numComponents; compIdx < numSComponents; ++compIdx)
+// (*moleFractionSolid[compIdx])[dofIdxGlobal]= elemVolVars[scvIdx].moleFraction(compIdx);
+//
+// Tensor K = this->perm_(this->problem_().spatialParams().intrinsicPermeability(element, fvGeometry, scvIdx));
+//
+// for (int j = 0; j<dim; ++j)
+// (*Perm[j])[dofIdxGlobal] = K[j][j] * elemVolVars[scvIdx].permeabilityFactor();
+// };
+//
+// // velocity output
+// if(velocityOutput.enableOutput()){
+// velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, phaseIdx);
+// }
+// } // loop over element
+//
+// writer.attachDofData(*p, "pressure", isBox);
+// writer.attachDofData(*rho, "rho", isBox);
+// writer.attachDofData(*poro, "porosity", isBox);
+// writer.attachDofData(*permeabilityFactor, "permeabilityFactor", isBox);
+// writer.attachDofData(*temperature, "temperature", isBox);
+//
+// for (int i = 0; i < numSPhases; ++i)
+// {
+// std::cout << "test1" << "\n";
+// std::ostringstream oss;
+// oss << "precipitateVolumeFraction_" << FluidSystem::componentName(numComponents + i);
+// writer.attachDofData(*precipitateVolumeFraction[i], oss.str(), isBox);
+// }
+//
+// writer.attachDofData(*Perm[0], "Kxx", isBox);
+// if (dim >= 2)
+// writer.attachDofData(*Perm[1], "Kyy", isBox);
+// if (dim == 3)
+// writer.attachDofData(*Perm[2], "Kzz", isBox);
+//
+//
+// // for (int j = 0; j < numComponents; ++j)
+// // {
+// // std::ostringstream oss;
+// // oss << "X^" << FluidSystem::phaseName(i) << "_" << FluidSystem::componentName(j);
+// // writer.attachDofData(*massFraction[i][j], oss.str(), isBox);
+// // }
+//
+// for (int j = 0; j < numComponents; ++j)
+// {
+// std::ostringstream oss;
+// oss << "m^w_" << FluidSystem::componentName(j);
+// writer.attachDofData(*molarity[j], oss.str(), isBox);
+// }
+//
+// for (int j = 0; j < numComponents; ++j)
+// {
+// std::ostringstream oss;
+// oss << "x"
+// << FluidSystem::componentName(j);
+// writer.attachDofData(*moleFraction[j], oss.str(), isBox);
+// }
+//
+// for (int j = numComponents; j < numSComponents; ++j)
+// {
+// std::ostringstream oss;
+// oss << "m^w_" << FluidSystem::componentName(j);
+// writer.attachDofData(*molarity[j], oss.str(), isBox);
+// }
+//
+// for (int j = numComponents; j < numSComponents; ++j)
+// {
+// std::ostringstream oss;
+// oss << "x"
+// << FluidSystem::componentName(j);
+// writer.attachDofData(*moleFraction[j], oss.str(), isBox);
+// }
+//
+// if (velocityOutput.enableOutput()) // check if velocity output is demanded
+// {
+// writer.attachDofData(*velocity, "velocity", isBox, dim);
+// }
+//
+// writer.attachCellData(*rank, "process rank");
+// }
+
+ /*!
+ * \brief Update the static data of all vertices in the grid.
+ *
+ * \param curGlobalSol The current global solution
+ * \param oldGlobalSol The previous global solution
+ */
+ void updateStaticData(SolutionVector &curGlobalSol,
+ const SolutionVector &oldGlobalSol)
+ {
+ for (unsigned i = 0; i < this->staticDat_.size(); ++i)
+ this->staticDat_[i].visited = false;
+
+ 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);
+ }
+ }
+ }
+
+ /*!
+ * \brief Write the current solution to a restart file.
+ *
+ * \param outStream The output stream of one entity for the restart file
+ * \param entity The entity, either a vertex or an element
+ */
+ template<class Entity>
+ void serializeEntity(std::ostream &outStream, const Entity &entity)
+ {
+ // write primary variables
+ ParentType::serializeEntity(outStream, entity);
+
+ int dofIdxGlobal = this->dofMapper().index(entity);
+
+ if (!outStream.good())
+ DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal);
+
+ }
+
+ /*!
+ * \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);
+ }
+
+
+ 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;
+ }
+};
+
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/porousmediumflow/1pncmin/implicit/properties.hh b/dumux/porousmediumflow/1pncmin/implicit/properties.hh
new file mode 100644
index 0000000000..bfe7ee8998
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/properties.hh
@@ -0,0 +1,69 @@
+// -**- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePNCMinModel
+ *
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase n-component mineralization
+ * fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_PROPERTIES_HH
+#define DUMUX_1PNCMIN_PROPERTIES_HH
+
+#include <dumux/porousmediumflow/1pnc/implicit/properties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal two phase n component mineralisation problems
+NEW_TYPE_TAG(OnePNCMin, INHERITS_FROM(OnePNC));
+NEW_TYPE_TAG(BoxOnePNCMin, INHERITS_FROM(BoxModel, OnePNCMin));
+NEW_TYPE_TAG(CCOnePNCMin, INHERITS_FROM(CCModel, OnePNCMin));
+
+//! The type tags for the corresponding non-isothermal problems
+NEW_TYPE_TAG(OnePNCMinNI, INHERITS_FROM(OnePNCMin, NonIsothermal));
+NEW_TYPE_TAG(BoxOnePNCMinNI, INHERITS_FROM(BoxModel, OnePNCMinNI));
+NEW_TYPE_TAG(CCOnePNCMinNI, INHERITS_FROM(CCModel, OnePNCMinNI));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumSPhases); //!< Number of solid phases in the system
+NEW_PROP_TAG(NumFSPhases); //!< Number of fluid and solid phases in the system
+NEW_PROP_TAG(NumSComponents); //!< Number of solid components in the system
+NEW_PROP_TAG(NumPSComponents); //!< Number of fluid and solid components in the system
+NEW_PROP_TAG(NumTraceComponents); //!< Number of trace fluid components which are not considered in the calculation of the phase density
+NEW_PROP_TAG(NumSecComponents); //!< Number of secondary components which are not primary variables
+NEW_PROP_TAG(OnePNCMinIndices); //!< Enumerations for the 2pncMin models
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+
+}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh b/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
new file mode 100644
index 0000000000..a17763e6f1
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
@@ -0,0 +1,180 @@
+// -**- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePNCMinModel
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * two-phase n-component mineralization fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_PROPERTY_DEFAULTS_HH
+#define DUMUX_1PNCMIN_PROPERTY_DEFAULTS_HH
+
+#include "model.hh"
+#include "indices.hh"
+#include "volumevariables.hh"
+#include "properties.hh"
+
+#include <dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh>
+#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
+#include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of secondary components.
+ * Secondary components are components calculated from
+ * primary components by equilibrium relations and
+ * do not have mass balance equation on their own.
+ * These components are important in the context of bio-mineralization applications.
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(OnePNCMin, NumSComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+public:
+ static const int value = FluidSystem::numSComponents;
+};
+/*!
+ * \brief Set the property for the number of solid phases, excluding the non-reactive matrix.
+ *
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(OnePNCMin, NumSPhases)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numSPhases;
+};
+
+/*!
+ * \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(OnePNCMin, NumEq)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+// static const int value = FluidSystem::numComponents + FluidSystem::numSPhases;
+ static const int value = FluidSystem::numComponents + FluidSystem::numSComponents; //steamaircao2h2 has 2 components in the fluidphase
+};
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ * store the thermodynamic state. This should be chosen
+ * appropriately for the model ((non-)isothermal, equilibrium, ...).
+ * This can be done in the problem.
+ */
+SET_PROP(OnePNCMin, FluidState){
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ public:
+ typedef CompositionalFluidState<Scalar, FluidSystem> type;
+};
+
+//! Use the 2pncmin local residual operator
+SET_TYPE_PROP(OnePNCMin,
+ LocalResidual,
+ OnePNCMinLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(OnePNCMin, Model, OnePNCMinModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(OnePNCMin, VolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+// SET_TYPE_PROP(OnePNCMin, FluxVariables, OnePNCMinFluxVariables<TypeTag>);
+
+//! The indices required by the isothermal 2pNcMin model
+SET_TYPE_PROP(OnePNCMin, Indices, OnePNCMinIndices <TypeTag, /*PVOffset=*/0>);
+
+//! default value for the forchheimer coefficient
+// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
+// Actually the Forchheimer coefficient is also a function of the dimensions of the
+// porous medium. Taking it as a constant is only a first approximation
+// (Nield, Bejan, Convection in porous media, 2006, p. 10)
+SET_SCALAR_PROP(OnePNCMin, SpatialParamsForchCoeff, 0.55);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNCMin, IsothermalVolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(OnePNCMinNI, ThermalConductivityModel)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+public:
+ typedef ThermalConductivityAverage<Scalar> type;
+};
+
+
+SET_BOOL_PROP(OnePNCMinNI, NiOutputLevel, 0);
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(OnePNCMinNI, IsothermalModel, OnePNCMinModel<TypeTag>);
+
+// set isothermal FluxVariables
+//SET_TYPE_PROP(OnePNCMinNI, IsothermalFluxVariables, OnePNCMinFluxVariables<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNCMinNI, IsothermalVolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(OnePNCMinNI, IsothermalLocalResidual, OnePNCMinLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(OnePNCMinNI, IsothermalIndices, OnePNCMinIndices<TypeTag, /*PVOffset=*/0>);
+
+//set isothermal NumEq
+SET_PROP(OnePNCMinNI, IsothermalNumEq)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents +FluidSystem::numSComponents;// in NonIsothermal 1 is substracted by default
+};
+}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pncmin/implicit/volumevariables.hh b/dumux/porousmediumflow/1pncmin/implicit/volumevariables.hh
new file mode 100644
index 0000000000..70fafae134
--- /dev/null
+++ b/dumux/porousmediumflow/1pncmin/implicit/volumevariables.hh
@@ -0,0 +1,390 @@
+// -**- 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 quantities which are constant within a
+ * finite volume in the two-phase, n-component mineralization model.
+ */
+#ifndef DUMUX_1PNCMIN_VOLUME_VARIABLES_HH
+#define DUMUX_1PNCMIN_VOLUME_VARIABLES_HH
+
+
+#include <dumux/common/math.hh>
+#include <dumux/implicit/model.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+#include <dumux/porousmediumflow/1pnc/implicit/volumevariables.hh>
+
+#include "properties.hh"
+#include "indices.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup OnePNCMinModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase, n-component model.
+ */
+template <class TypeTag>
+class OnePNCMinVolumeVariables : public OnePNCVolumeVariables<TypeTag>
+{
+ // base type is used for energy related quantites
+ using BaseType = ImplicitVolumeVariables<TypeTag>;
+
+ using ParentType = OnePNCVolumeVariables<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 ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+// typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+ numSComponents = GET_PROP_VALUE(TypeTag, NumSComponents), //if there is more than 1 component in the solid phase
+
+ phaseCompIdx = Indices::phaseCompIdx,
+
+ // phase indices
+ phaseIdx = FluidSystem::gPhaseIdx,
+ cPhaseIdx = Indices::phaseIdx +1,
+ hPhaseIdx = Indices::phaseIdx +2,
+
+ // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using CoordScalar = typename Grid::ctype;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const ElementSolutionVector &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume& scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+
+// ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
+// completeFluidState(priVars, problem, element, fvGeometry, scvIdx, this->fluidState_, isOldSol);
+
+ /////////////
+ // calculate the remaining quantities
+ /////////////
+
+ auto&& priVars = isBox ? elemSol[scv.index()] : elemSol[0];
+
+ // porosity evaluation
+ initialPorosity_ = problem.spatialParams().porosity(element, scv);
+ minimumPorosity_ = problem.spatialParams().porosityMin(element, scv);
+ maximumPorosity_ = problem.spatialParams().porosityMax(element, scv);
+
+ sumPrecipitates_ = 0.0;
+ for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+ {
+ precipitateVolumeFraction_[sPhaseIdx] = priVars[numComponents + sPhaseIdx];
+ sumPrecipitates_+= precipitateVolumeFraction_[sPhaseIdx];
+
+ }
+
+ // TODO/FIXME: The salt crust porosity is not clearly defined. However form literature review it is
+ // found that the salt crust have porosity of approx. 10 %. Thus we restrict the decrease in porosity
+ // to this limit. Moreover in the Problem files the precipitation should also be made dependent on local
+ // porosity value, as the porous media media properties change related to salt precipitation will not be
+ // accounted otherwise.
+
+ this->porosity_ = 1 - sumPrecipitates_;
+
+ 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_(elemSol, problem, element, scv);
+ }
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ * \param isOldSol Specifies whether this is the previous solution or the current one
+ */
+ static void completeFluidState(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ FluidState& fluidState)
+
+ {
+// Scalar t = ParentType::temperature(elemSol, problem, element, scv); //old
+ Scalar t = BaseType::temperature(elemSol, problem, element, scv);
+ fluidState.setTemperature(t);
+
+ /////////////
+ // set the saturations
+ /////////////
+
+ fluidState.setSaturation(phaseIdx, 1.0 );
+
+ /////////////
+ // set the pressures of the fluid phase
+ /////////////
+ const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+ fluidState.setPressure(phaseIdx, priVars[pressureIdx]);
+
+ /////////////
+ // calculate the phase compositions
+ /////////////
+
+ typename FluidSystem::ParameterCache paramCache;
+
+ Dune::FieldVector<Scalar, numComponents> moleFrac;
+
+ Scalar sumMoleFracNotWater = 0;
+ for (int compIdx=firstMoleFracIdx; compIdx<numComponents; ++compIdx)
+ {
+ moleFrac[compIdx] = priVars[compIdx];
+ sumMoleFracNotWater+=moleFrac[compIdx];
+ }
+
+ moleFrac[0] = 1 -sumMoleFracNotWater;
+
+// //mole fractions for the solid phase
+// moleFrac[firstMoleFracIdx+1]= priVars[firstMoleFracIdx+1];
+// moleFrac[firstMoleFracIdx +2] = 1- moleFrac[firstMoleFracIdx+1];
+
+ // convert mass to mole fractions and set the fluid state
+ for (int compIdx=0; compIdx<numComponents; ++compIdx)
+ {
+ fluidState.setMoleFraction(phaseIdx, compIdx, moleFrac[compIdx]);
+ }
+
+// // set mole fractions for the solid phase
+// fluidState.setMoleFraction(cPhaseIdx, firstMoleFracIdx+1, moleFrac[firstMoleFracIdx+1]);
+// fluidState.setMoleFraction(hPhaseIdx, firstMoleFracIdx+2, moleFrac[firstMoleFracIdx+2]);
+
+ paramCache.updateAll(fluidState);
+
+// Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+ Scalar h = BaseType::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, h);
+ }
+ /*!
+ * \brief Returns the volume fraction of the precipitate (solid phase)
+ * for the given phaseIdx
+ *
+ * \param phaseIdx the index of the solid phase
+ */
+ Scalar precipitateVolumeFraction(int phaseIdx) const
+ {
+ return precipitateVolumeFraction_[phaseIdx - numPhases];
+ }
+
+ /*!
+ * \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 density of the phase for all fluid and solid phases
+ *
+ * \param phaseIdx the index of the fluid phase
+ */
+ Scalar density(int phaseIdx) const
+ {
+ if (phaseIdx <= numPhases)
+ return this->fluidState_.density(phaseIdx);
+ else if (phaseIdx > numPhases)
+ return FluidSystem::precipitateDensity(phaseIdx);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(int phaseIdx) const
+ {
+ if (phaseIdx < 1)
+ return this->fluidState_.molarDensity(phaseIdx);
+ else if (phaseIdx >= 1){
+ /*Attention: sPhaseIdx of the fluidsystem and the model can be different.*/
+// std::cout << "FluidSystem::precipitateMolarDensity("<<phaseIdx<<") = " << FluidSystem::precipitateMolarDensity(phaseIdx) << "\n";
+// for (int phaseIdx=1; phaseIdx<numSPhases; ++phaseIdx)
+ return FluidSystem::precipitateMolarDensity(phaseIdx);
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Returns the molality of a component in the phase
+ *
+ * \param phaseIdx the index of the fluid phase
+ * \param compIdx the index of the component
+ * \f$\mathrm{molality}=\frac{n_\mathrm{component}}{m_\mathrm{solvent}}
+ * =\frac{n_\mathrm{component}}{n_\mathrm{solvent}*M_\mathrm{solvent}}\f$
+ * compIdx of the main component (solvent) in the
+ * phase is equal to the phaseIdx
+ */
+ Scalar molality(int phaseIdx, int compIdx) const // [moles/Kg]
+ { return this->fluidState_.moleFraction(phaseIdx, compIdx)
+ /(this->fluidState_.moleFraction(phaseIdx, phaseIdx)
+ * FluidSystem::molarMass(phaseIdx));}
+
+ /*!
+ * Circumvents the inheritance architecture of the nonisothermal model
+ */
+ static Scalar callProtectedTemperature(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+// return Implementation::temperature_(priVars, problem,element, fvGeometry, scvIdx);
+ return BaseType::temperature(elemSol, problem, element, scv);
+ }
+
+ /*!
+ * Circumvents the inheritance architecture of the ninisothermal model
+ */
+ void callProtectedUpdateEnergy(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ asImp_().updateEnergy_(elemSol, problem, element, scv);
+ };
+
+protected:
+ friend class OnePNCVolumeVariables<TypeTag>;
+ static Scalar temperature_(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ {
+ return problem.temperatureAtPos(scv);
+ }
+
+ 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 ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv)
+ { };
+
+ Scalar precipitateVolumeFraction_[numSPhases];
+ Scalar permeabilityFactor_;
+ Scalar initialPorosity_;
+ Scalar initialPermeability_;
+ Scalar minimumPorosity_;
+ Scalar maximumPorosity_;
+ Scalar sumPrecipitates_;
+
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index 6bc244e01f..013eb52802 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory("1p")
add_subdirectory("1p2c")
add_subdirectory("1pnc")
+add_subdirectory("1pncmin")
add_subdirectory("2p")
add_subdirectory("2p1c")
add_subdirectory("2p2c")
--
GitLab