diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index 4f9898ae54280cc58545869af74459254a2ed5f6..f75fd9d1f47f4396ee7632efe25e7b0bc99b8142 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -745,10 +745,10 @@ private:
*/
void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::false_type) const
{
- for (const auto& element : elements(fvGridGeometry_->gridView()))
+ for (const auto& vertex : vertices(fvGridGeometry_->gridView()))
{
- const auto dofIdxGlobal = fvGridGeometry_->elementMapper().index(element);
- sol[dofIdxGlobal] = asImp_().initial(element);
+ const auto dofIdxGlobal = fvGridGeometry_->vertexMapper().index(vertex);
+ sol[dofIdxGlobal] = asImp_().initial(vertex);
}
}
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 1c334bf8ca0d9761051d0e8172514c0e34092e92..a281388ec0d805470062ec495e2aa7ea3fe6db21 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -136,6 +136,11 @@ NEW_PROP_TAG(WettingPhase); //! The wetting phase for two
NEW_PROP_TAG(NonwettingPhase); //! The non-wetting phase for two-phase models
NEW_PROP_TAG(Formulation); //! The formulation of the model
+/////////////////////////////////////////////////////////////
+// Properties used by models involving mineralization:
+/////////////////////////////////////////////////////////////
+NEW_PROP_TAG(NumSPhases);
+
/////////////////////////////////////////////////////////////
// non-isothermal porous medium flow models
/////////////////////////////////////////////////////////////
diff --git a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
index b1fe9cecd4eafb234014490aab2c621a1ffd30bb..6f20c5f7024a538af05a781e6ef5e35870a1e457 100644
--- a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
+++ b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
@@ -34,7 +34,7 @@ namespace Dumux
*/
/**
- * \brief Calculates the porosity depeding on the volume fractions of precipitated minerals.
+ * \brief Calculates the porosity depending on the volume fractions of precipitated minerals.
*/
template<class TypeTag>
class PorosityPrecipitation
diff --git a/dumux/porousmediumflow/2pncmin/implicit/model.hh b/dumux/porousmediumflow/2pncmin/implicit/model.hh
index b5741a3f2dcaa1591c33a13f9d15db06d2212090..cdd99c28ac7789417d02fe82b02cd9f2d6a415b7 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/model.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh
@@ -25,384 +25,8 @@
#ifndef DUMUX_2PNCMIN_MODEL_HH
#define DUMUX_2PNCMIN_MODEL_HH
-#include "properties.hh"
-#include "indices.hh"
-#include "primaryvariableswitch.hh"
-#include "localresidual.hh"
-
-#include <dumux/porousmediumflow/2pnc/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
-namespace Dumux
-{
-/*!
- * \ingroup TwoPNCMinModel
- * \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 TwoPNCMinModel: public GET_PROP_TYPE(TypeTag, BaseModel)
-{
- // 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 ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
- 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),
-
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx
- };
-
- 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 { 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("Sn", [](const VolumeVariables& v){ return v.saturation(nPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(wPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
- vtkOutputModule.addSecondaryVariable("rhoW", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(nPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(nPhaseIdx); });
- vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
- 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 < numPhases; ++i)
- for (int j = 0; j < numComponents; ++j)
- vtkOutputModule.addSecondaryVariable("x_" + FluidSystem::phaseName(i) + "^" + FluidSystem::componentName(j),
- [i,j](const VolumeVariables& v){ return v.moleFraction(i,j); });
-
- for (int j = 0; j < numComponents; ++j)
- vtkOutputModule.addSecondaryVariable("m_" + FluidSystem::phaseName(wPhaseIdx) + "^" + FluidSystem::componentName(j),
- [j](const VolumeVariables& v){ return v.molarity(wPhaseIdx,j); });
- }
-
- /*!
- * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
- */
- template<class VtkOutputModule>
- void addVtkOutputFields(VtkOutputModule& outputModule) const
- {
- auto& phasePresence = outputModule.createScalarField("phase presence", dofCodim);
- for (std::size_t i = 0; i < phasePresence.size(); ++i)
- phasePresence[i] = this->curSol()[i].state();
- }
-
- /*!
- * \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
- for (const auto& element : elements(this->problem_().gridView()))
- {
- // make sure FVElementGeometry & vol vars are bound to the element
- auto fvGeometry = localView(this->fvGridGeometry());
- fvGeometry.bindElement(element);
-
- if (switchFlag_)
- {
- for (auto&& scv : scvs(fvGeometry))
- {
- auto dofIdxGlobal = scv.dofIndex();
- if (priVarSwitch_().wasSwitched(dofIdxGlobal))
- {
- const auto eIdx = this->problem_().elementMapper().index(element);
- const auto elemSol = this->elementSolution(element, this->curSol());
- this->nonConstCurGlobalVolVars().volVars(eIdx, scv.indexInElement()).update(elemSol,
- this->problem_(),
- element,
- scv);
- }
- }
- }
-
- // handle the boundary volume variables for cell-centered models
- if(!isBox)
- {
- for (auto&& scvf : scvfs(fvGeometry))
- {
- // if we are not on a boundary, skip the rest
- if (!scvf.boundary())
- continue;
-
- // check if boundary is a pure dirichlet boundary
- const auto bcTypes = this->problem_().boundaryTypes(element, scvf);
- if (bcTypes.hasOnlyDirichlet())
- {
- const auto insideScvIdx = scvf.insideScvIdx();
- const auto& insideScv = fvGeometry.scv(insideScvIdx);
- const auto elemSol = ElementSolutionVector{this->problem_().dirichlet(element, scvf)};
-
- this->nonConstCurGlobalVolVars().volVars(scvf.outsideScvIdx(), 0/*indexInElement*/).update(elemSol, this->problem_(), element, insideScv);
- }
- }
- }
- }
- }
-
- /*!
- * \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();
- // reset privar switch flag
- switchFlag_ = false;
- }
-
- /*!
- * \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();
- // reset privar switch flag
- 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 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);
-
- outStream << this->curSol()[dofIdxGlobal].state() << " ";
- }
-
- /*!
- * \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;
-
- this->curSol()[dofIdxGlobal].setState(phasePresence);
- this->prevSol()[dofIdxGlobal].setState(phasePresence);
- }
-
- const Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag>& priVarSwitch() const
- { return switch_; }
-
-private:
-
- Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag>& priVarSwitch_()
- { return switch_; }
-
- /*!
- * \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...
- */
- void applyInitialSolution_()
- {
- ParentType::applyInitialSolution_();
-
- // initialize the primary variable switch
- priVarSwitch_().init(this->problem_());
- }
-
- //! the class handling the primary variable switch
- Dumux::TwoPNCMinPrimaryVariableSwitch<TypeTag> switch_;
- bool switchFlag_;
-
- Tensor perm_(Scalar perm) const
- {
- Tensor K(0.0);
-
- for(int i=0; i<dimWorld; i++)
- K[i][i] = perm;
-
- return K;
- }
-
- const Tensor& perm_(const Tensor& perm) const
- {
- return perm;
- }
-};
-
-} // end namespace Dumux
-
-#include "propertydefaults.hh"
+#include "properties.hh"
#endif
diff --git a/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh b/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh
deleted file mode 100644
index 6d20ae0309bd08ad11695ae2560db63042ae9027..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2pncmin/implicit/primaryvariableswitch.hh
+++ /dev/null
@@ -1,190 +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 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 {
- pwsn = TwoPNCFormulation::pwsn,
- pnsw = TwoPNCFormulation::pnsw,
- 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 = priVars.state();
- 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 == pnsw)
- priVars[switchIdx] = 0.0;
- else if (formulation == pwsn)
- 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 pnsw and vice versa)
- if (formulation == pnsw)
- priVars[switchIdx] = 0.999;
- else if (formulation == pwsn)
- priVars[switchIdx] = 0.001;
-
- }
- }
- priVars.setState(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 d024b672cba96849357931400f8aceba7d64b05e..84510c7c1717464e38e711157306585de07ed6ea 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/properties.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/properties.hh
@@ -31,6 +31,9 @@
#include <dumux/porousmediumflow/2pnc/implicit/properties.hh>
+#include "volumevariables.hh"
+#include "vtkoutputfields.hh"
+
namespace Dumux
{
@@ -39,30 +42,124 @@ namespace Properties
//////////////////////////////////////////////////////////////////
// Type tags
//////////////////////////////////////////////////////////////////
-
-//! The type tag for the isothermal two phase n component mineralisation problems
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
+// Property tags for the isothermal 2pncmin model
//////////////////////////////////////////////////////////////////
-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(TwoPNCMinIndices); //!< Enumerations for the 2pncMin models
-NEW_PROP_TAG(useSalinity); //!< Determines if salinity is used
-NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+SET_TYPE_PROP(TwoPNCMin, PrimaryVariables, SwitchablePrimaryVariables<TypeTag, int>); //! The primary variables vector for the 2pncmin model
+SET_TYPE_PROP(TwoPNCMin, PrimaryVariableSwitch, TwoPNCPrimaryVariableSwitch<TypeTag>); //! Using the 2pnc primary variable switch for the 2pncmin model
+SET_TYPE_PROP(TwoPNCMin, VolumeVariables, TwoPNCMinVolumeVariables<TypeTag>); //! the VolumeVariables property
+SET_TYPE_PROP(TwoPNCMin, Indices, TwoPNCIndices <TypeTag, /*PVOffset=*/0>); //! Using the 2pnc indices required by the isothermal 2pncmin model
+SET_TYPE_PROP(TwoPNCMin, SpatialParams, ImplicitSpatialParams<TypeTag>); //! Use the ImplicitSpatialParams by default
+SET_TYPE_PROP(TwoPNCMin, VtkOutputFields, TwoPNCMinVtkOutputFields<TypeTag>); //! Set the vtk output fields specific to the TwoPNCMin model
+SET_TYPE_PROP(TwoPNCMin, LocalResidual, CompositionalLocalResidual<TypeTag>); //! Use the compositional local residual
+
+SET_INT_PROP(TwoPNCMin, NumComponents, GET_PROP_TYPE(TypeTag, FluidSystem)::numComponents); //! Use the number of components of the fluid system
+SET_INT_PROP(TwoPNCMin, ReplaceCompEqIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::numComponents); //! Per default, no component mass balance is replaced
+SET_INT_PROP(TwoPNCMin, Formulation, TwoPNCFormulation::pwsn); //! Using the default 2pnc formulation: pw-Sn, overwrite if necessary
+
+SET_BOOL_PROP(TwoPNCMin, SetMoleFractionsForWettingPhase, true); //! Set the primary variables mole fractions for the wetting or non-wetting phase
+SET_BOOL_PROP(TwoPNCMin, EnableAdvection, true); //! Enable advection
+SET_BOOL_PROP(TwoPNCMin, EnableMolecularDiffusion, true); //! Enable molecular diffusion
+SET_BOOL_PROP(TwoPNCMin, EnableEnergyBalance, false); //! This is the isothermal variant of the model
+SET_BOOL_PROP(TwoPNCMin, UseMoles, true); //! Use mole fractions in the balance equations by default
+
+//! Use the model after Millington (1961) for the effective diffusivity
+SET_TYPE_PROP(TwoPNCMin, EffectiveDiffusivityModel, DiffusivityMillingtonQuirk<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+//! The major components belonging to the existing phases, e.g. 2 for water and air being the major components in a liquid-gas-phase system
+SET_PROP(TwoPNCMin, NumMajorComponents)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 2, "The model is restricted to two phases, thus number of major components must also be two.");
+};
+
+//! We use the number of phases of the fluid system. Make sure it is 2!
+SET_PROP(TwoPNCMin, NumPhases)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 2, "Only fluid systems with 2 fluid phases are supported by the 2p-nc-min model!");
+};
+
+//! This model uses the compositional fluid state
+SET_PROP(TwoPNCMin, FluidState)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+public:
+ using type = CompositionalFluidState<Scalar, FluidSystem>;
+};
+
+//! Set the property for the material parameters by extracting it from the material law.
+SET_PROP(TwoPNCMin, MaterialLawParams)
+{
+private:
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw));
+
+public:
+ using type = typename MaterialLaw::Params;
+};
+
+//! Set the property for the number of solid phases, excluding the non-reactive matrix.
+SET_PROP(TwoPNCMin, NumSPhases)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem));
+
+public:
+ static const int value = FluidSystem::numSPhases;
+};
+
+//! 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(TwoPNCMin, NumEq)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem));
+
+public:
+ static const int value = FluidSystem::numComponents + FluidSystem::numSPhases;
+};
+/////////////////////////////////////////////////
+// Properties for the non-isothermal 2pncmin model
+/////////////////////////////////////////////////
+SET_TYPE_PROP(TwoPNCMinNI, IsothermalVolumeVariables, TwoPNCMinVolumeVariables<TypeTag>); //! set isothermal VolumeVariables
+SET_TYPE_PROP(TwoPNCMinNI, IsothermalLocalResidual, CompositionalLocalResidual<TypeTag>); //! set isothermal LocalResidual
+SET_TYPE_PROP(TwoPNCMinNI, IsothermalIndices, TwoPNCMinIndices<TypeTag, /*PVOffset=*/0>); //! set isothermal Indices
+SET_TYPE_PROP(TwoPNCMinNI, IsothermalVtkOutputFields, TwoPNCMinVtkOutputFields<TypeTag>); //! set isothermal output fields
+
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(TwoPNCMinNI, ThermalConductivityModel)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+public:
+ typedef ThermalConductivitySomerton<Scalar, Indices> type;
+};
+
+//! set isothermal NumEq
+SET_PROP(TwoPNCMinNI, IsothermalNumEq)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+};
}
}
diff --git a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
index 7efe9cd602f74284b7c8430b2883d79216b729e5..63966c1016303b5a85d310bdb6e469b7ae7e8195 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
@@ -48,7 +48,7 @@ namespace Dumux
template <class TypeTag>
class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
{
- // base type is used for energy related quantites
+ // base type is used for energy related quantities
using BaseType = ImplicitVolumeVariables<TypeTag>;
using ParentType = TwoPNCVolumeVariables<TypeTag>;
@@ -97,7 +97,6 @@ class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
- useSalinity = GET_PROP_VALUE(TypeTag, useSalinity)
};
using Element = typename GridView::template Codim<0>::Entity;
@@ -106,9 +105,6 @@ class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
using Miscible2pNCComposition = Dumux::Miscible2pNCComposition<Scalar, FluidSystem>;
using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
-
public:
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
@@ -129,7 +125,7 @@ public:
/////////////
// calculate the remaining quantities
/////////////
- auto&& priVars = isBox ? elemSol[scv.indexInElement()] : elemSol[0];
+ auto&& priVars = elemSol[scv.indexInElement()];
sumPrecipitates_ = 0.0;
for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
@@ -195,6 +191,7 @@ protected:
Scalar precipitateVolumeFraction_[numSPhases];
Scalar sumPrecipitates_;
+ FluidState fluidState_;
private:
Implementation &asImp_()
diff --git a/dumux/porousmediumflow/2pncmin/implicit/vtkoutputfields.hh b/dumux/porousmediumflow/2pncmin/implicit/vtkoutputfields.hh
index 6ea93e7f3aae64e610109a9705143c90cc2c6d9f..b4b61ba5762263b50611457372f2eeba822c6478 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/vtkoutputfields.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/vtkoutputfields.hh
@@ -41,7 +41,7 @@ class TwoPNCMinVtkOutputFields
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static constexpr int numSPhases = GET_PROP_VALE(TypteTag, NumSPhases);
+ static constexpr int numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases);
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
static constexpr int dim = GridView::dimension;
@@ -52,14 +52,11 @@ public:
// use default fields from the 2pnc model
TwoPNCVtkOutputFields<TypeTag>::init(vtk);
- //output additional to TwoPNC output:
+ //output additional to TwoPNCMin output:
for (int i = 0; i < numSPhases; ++i)
+ {
vtk.addVolumeVariable([i](const VolumeVariables& v){ return v.precipitateVolumeFraction(numPhases + i); },"precipVolFrac_"+ FluidSystem::phaseName(numPhases + i));
- vtk.addVolumeVariable([](const VolumeVariables& v){ return this->perm_(v.permeability())[0][0]; }, "Kxx"); //TODO: get correct permeability from where? add perm_ function in private?
- if (dim >= 2)
- vtk.addVolumeVariable([](const VolumeVariables& v){ return this->perm_(v.permeability())[1][1]; }, "Kyy"); //TODO: get correct permeability from where? add perm_ function in private?
- if (dim >= 3)
- vtk.addVolumeVariable([](const VolumeVariables& v){ return this->perm_(v.permeability())[2][2]; }, "Kzz"); //TODO: get correct permeability from where? add perm_ function in private?
+ }
}
};
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
index 6fa7d6bc5ffcf12bb9adbcabb6789aa21b4ee053..fc44675de35c1f58dd8a5a345af9606b5d4a619a 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
@@ -60,9 +60,6 @@ SET_PROP(DissolutionProblem, FluidSystem)
// Set the spatial parameters
SET_TYPE_PROP(DissolutionProblem, SpatialParams, DissolutionSpatialparams<TypeTag>);
-// Enable gravity
-SET_BOOL_PROP(DissolutionProblem, ProblemEnableGravity, true);
-
//Set properties here to override the default property settings in the model.
SET_INT_PROP(DissolutionProblem, ReplaceCompEqIdx, 1); //! Replace gas balance by total mass balance
SET_INT_PROP(DissolutionProblem, Formulation, TwoPNCFormulation::pnsw);
@@ -133,9 +130,11 @@ class DissolutionProblem : public PorousMediumFlowProblem<TypeTag>
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
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 FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using Vertex = typename GridView::template Codim<dim>::Entity;
using Intersection = typename GridView::Intersection;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
@@ -143,31 +142,34 @@ class DissolutionProblem : public PorousMediumFlowProblem<TypeTag>
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
- DissolutionProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ DissolutionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
- outerSalinity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterSalinity);
- temperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, Temperature);
- reservoirPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, ReservoirPressure);
- initLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, LiquidSaturation);
- initPrecipitatedSalt1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt1);
- initPrecipitatedSalt2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InitPrecipitatedSalt2);
-
- outerLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterLiqSaturation);
- innerLiqSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerLiqSaturation);
- innerSalinity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerSalinity);
- innerPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InnerPressure);
- outerPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OuterPressure);
- reservoirSaturation_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, reservoirSaturation);
-
- nTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NTemperature);
- nPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NPressure);
- pressureLow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureLow);
- pressureHigh_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureHigh);
- 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);
+ outerSalinity_ = getParam<Scalar>("Problem.OuterSalinity");
+ temperature_ = getParam<Scalar>("Problem.Temperature");
+ reservoirPressure_ = getParam<Scalar>("Problem.ReservoirPressure");
+ initLiqSaturation_ = getParam<Scalar>("Problem.LiquidSaturation");
+ initPrecipitatedSalt1_ = getParam<Scalar>("Problem.InitPrecipitatedSalt1");
+ initPrecipitatedSalt2_ = getParam<Scalar>("Problem.InitPrecipitatedSalt2");
+
+ outerLiqSaturation_ = getParam<Scalar>("Problem.OuterLiqSaturation");
+ innerLiqSaturation_ = getParam<Scalar>("Problem.InnerLiqSaturation");
+ innerSalinity_ = getParam<Scalar>("Problem.InnerSalinity");
+ innerPressure_ = getParam<Scalar>("Problem.InnerPressure");
+ outerPressure_ = getParam<Scalar>("Problem.OuterPressure");
+ reservoirSaturation_ = getParam<Scalar>("Problem.reservoirSaturation");
+
+ nTemperature_ = getParam<int>("FluidSystem.NTemperature");
+ nPressure_ = getParam<int>("FluidSystem.NPressure");
+ pressureLow_ = getParam<Scalar>("FluidSystem.PressureLow");
+ pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh");
+ temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow");
+ temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh");
+ name_ = getParam<std::string>("Problem.Name");
+
+ Kxx_.resize(fvGridGeometry->gridView().size(0));
+ Kyy_.resize(fvGridGeometry->gridView().size(0));
outfile.open("evaporation.out");
outfile << "time; evaporationRate" << std::endl;
@@ -185,13 +187,22 @@ public:
outfile.close();
}
- bool shouldWriteOutput() const
+ void setTime( Scalar time )
{
- return this->timeManager().timeStepIndex() % 1 == 0 ||
- this->timeManager().episodeWillBeFinished() ||
- this->timeManager().willBeFinished();
+ time_ = time;
}
+ void setTimeStepSize( Scalar timeStepSize )
+ {
+ timeStepSize_ = timeStepSize;
+ }
+// bool shouldWriteOutput() const
+// {
+// return this->timeManager().timeStepIndex() % 1 == 0 ||
+// this->timeManager().episodeWillBeFinished() ||
+// this->timeManager().willBeFinished();
+// }
+
/*!
* \name Problem parameters
@@ -227,8 +238,8 @@ public:
{
BoundaryTypes bcTypes;
- const Scalar rmax = this->bBoxMax()[0];
- const Scalar rmin = this->bBoxMin()[0];
+ const Scalar rmax = this->fvGridGeometry().bBoxMax()[0];
+ const Scalar rmin = this->fvGridGeometry().bBoxMin()[0];
// default to Neumann
bcTypes.setAllNeumann();
@@ -253,8 +264,8 @@ public:
PrimaryVariables priVars(0.0);
priVars.setState(bothPhases);
- const Scalar rmax = this->bBoxMax()[0];
- const Scalar rmin = this->bBoxMin()[0];
+ const Scalar rmax = this->fvGridGeometry().bBoxMax()[0];
+ const Scalar rmin = this->fvGridGeometry().bBoxMin()[0];
if(globalPos[0] > rmax - eps_)
{
@@ -279,23 +290,20 @@ public:
/*!
* \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
+ * \param Vertex The finite element
*
* For this method, the \a values parameter stores primary
* variables.
*/
- PrimaryVariables initial(const Element &element) const
+ PrimaryVariables initial(const Vertex &vertex) const
{
PrimaryVariables priVars(0.0);
priVars.setState(bothPhases);
- const auto& globalPos = element.geometry().center();
+ const auto& globalPos = vertex.geometry().corner(0);
priVars[pressureIdx] = reservoirPressure_;
- priVars[switchIdx] = initLiqSaturation_; // Sl primary variable
+ priVars[switchIdx] = initLiqSaturation_; // Sw primary variable
priVars[xlNaClIdx] = massToMoleFrac_(outerSalinity_); // mole fraction
if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 19.0 + eps_)
priVars[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
@@ -357,12 +365,11 @@ public:
* 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();
+ // make sure we don't dissolve more salt than previously precipitated
+ if (precipSalt*time_ + volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)< 0)
+ precipSalt = -volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)/timeStepSize_;
- auto curElemSol = this->model().elementSolution(element, this->model().curSol());
- if (volVars.precipitateVolumeFraction(sPhaseIdx) >= this->spatialParams().porosity(element, scv, curElemSol) - saltPorosity && precipSalt > 0)
+ if (volVars.precipitateVolumeFraction(sPhaseIdx) >= volVars.porosity() - saltPorosity && precipSalt > 0)
precipSalt = 0;
source[conti0EqIdx + NaClIdx] += -precipSalt;
@@ -371,6 +378,40 @@ public:
return source;
}
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+
+ const std::vector<Scalar>& getKxx()
+ {
+ return Kxx_;
+ }
+
+ const std::vector<Scalar>& getKyy()
+ {
+ return Kyy_;
+ }
+
+ void updateVtkOutput(const SolutionVector& curSol)
+ {
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ ElementSolutionVector elemSol(element, curSol, this->fvGridGeometry());
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(elemSol, *this, element, scv);
+ const auto dofIdxGlobal = scv.dofIndex();
+ Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
+ Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
+ }
+ }
+ }
+
private:
/*!
@@ -406,9 +447,13 @@ private:
Scalar initPrecipitatedSalt1_;
Scalar initPrecipitatedSalt2_;
Scalar innerSalinity_;
+ Scalar time_ = 0.0;
+ Scalar timeStepSize_ = 0.0;
static constexpr Scalar eps_ = 1e-6;
Scalar reservoirSaturation_;
std::ofstream outfile;
+ std::vector<double> Kxx_;
+ std::vector<double> Kyy_;
};
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
index 26901b9144a94f352ac3aa64e4bc2fcc6c4bf983..c7cb9a4e3dbcf8fd129fd648b56e1c2beabce91a 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
@@ -92,10 +92,10 @@ class DissolutionSpatialparams : public ImplicitSpatialParams<TypeTag>
public:
// type used for the permeability (i.e. tensor or scalar)
- using PermeabilityType = Scalar;
+ using PermeabilityType = Tensor;
- DissolutionSpatialparams(const Problem& problem, const GridView &gridView)
- : ParentType(problem, gridView)
+ DissolutionSpatialparams(const Problem& problem)
+ : ParentType(problem)
{
// residual saturations
materialParams_.setSwr(0.2);
@@ -124,7 +124,7 @@ public:
*
* Solution dependent permeability function
*/
- Scalar permeability(const Element& element,
+ PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{ return permLaw_.evaluatePermeability(element, scv, elemSol); }
diff --git a/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input b/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
index 02470b48386c7409c891fe21d01057b8cc3ecd83..265623fff89c0e833252d053dfc98e03ab998feb 100644
--- a/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
+++ b/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
@@ -2,7 +2,7 @@
# Everything behind a '#' is a comment.
# Type "./test_2pncmin --help" for more information.
-[TimeManager]
+[TimeLoop]
TEnd = 1e6 # duration of the simulation [s]
DtInitial = 10 # initial time step size [s]
MaxTimeStepSize = 50000 # maximum time step size
diff --git a/test/porousmediumflow/2pncmin/implicit/test_box2pncmin.cc b/test/porousmediumflow/2pncmin/implicit/test_box2pncmin.cc
index b4975cb1781eb5ff29be5cce18eacd38bcf5d5e3..41e1a022fd4233e2fd51555bb5c7d21ca02bde83 100644
--- a/test/porousmediumflow/2pncmin/implicit/test_box2pncmin.cc
+++ b/test/porousmediumflow/2pncmin/implicit/test_box2pncmin.cc
@@ -74,7 +74,7 @@ void usage(const char *progName, const std::string &errorMsg)
}
}
-int main(int argc, char** argv)
+int main(int argc, char** argv) try
{
using namespace Dumux;
@@ -141,9 +141,6 @@ int main(int argc, char** argv)
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
//add specific output
- vtkWriter.addField(problem->getCurrentDensity(), "currentDensity [A/cm^2]");
- vtkWriter.addField(problem->getReactionSourceH2O(), "reactionSourceH2O [mol/(sm^2)]");
- vtkWriter.addField(problem->getReactionSourceO2(), "reactionSourceO2 [mol/(sm^2)]");
vtkWriter.addField(problem->getKxx(), "Kxx");
vtkWriter.addField(problem->getKyy(), "Kyy");
vtkWriter.write(0.0);
@@ -169,6 +166,10 @@ int main(int argc, char** argv)
// time loop
timeLoop->start(); do
{
+ // set time for problem for implicit Euler scheme
+ problem->setTime( timeLoop->time() + timeLoop->timeStepSize() );
+ problem->setTimeStepSize( timeLoop->timeStepSize() );
+
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
diff --git a/test/porousmediumflow/2pncmin/implicit/test_cc2pncmin.cc b/test/porousmediumflow/2pncmin/implicit/test_cc2pncmin.cc
index ac4d872988b6af0af44004aa13a19f1ce6512d83..fdbafc6593fcbacb5fd13afc2fb4f358de0b0319 100644
--- a/test/porousmediumflow/2pncmin/implicit/test_cc2pncmin.cc
+++ b/test/porousmediumflow/2pncmin/implicit/test_cc2pncmin.cc
@@ -113,8 +113,9 @@ int main(int argc, char** argv) try
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- SolutionVector x(leafGridView.size(0));
+ SolutionVector x(leafGridView.size(GridView::dimension));
problem->applyInitialSolution(x);
auto xOld = x;
@@ -140,9 +141,6 @@ int main(int argc, char** argv) try
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
//add specific output
- vtkWriter.addField(problem->getCurrentDensity(), "currentDensity [A/cm^2]");
- vtkWriter.addField(problem->getReactionSourceH2O(), "reactionSourceH2O [mol/(sm^2)]");
- vtkWriter.addField(problem->getReactionSourceO2(), "reactionSourceO2 [mol/(sm^2)]");
vtkWriter.addField(problem->getKxx(), "Kxx");
vtkWriter.addField(problem->getKyy(), "Kyy");
vtkWriter.write(0.0);
@@ -168,6 +166,10 @@ int main(int argc, char** argv) try
// time loop
timeLoop->start(); do
{
+ // set time for problem for implicit Euler scheme
+ problem->setTime( timeLoop->time() + timeLoop->timeStepSize() );
+ problem->setTimeStepSize( timeLoop->timeStepSize() );
+
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);