Commit d87a7704 authored by Gabi Seitz's avatar Gabi Seitz Committed by Dennis Gläser
Browse files

[1pncmin] delete folder implicit and inherit from 1pnc and mineralization

parent 5eb5cb17
add_subdirectory("implicit")
#install headers
install(FILES
model.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1pncmin)
#install headers
install(FILES
indices.hh
localresidual.hh
model.hh
properties.hh
propertydefaults.hh
volumevariables.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1pncmin/implicit)
// -*- 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/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
// -*- 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 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 CompositionalLocalResidual<TypeTag>
{
protected:
using ParentType = CompositionalLocalResidual<TypeTag>;
using ThisType = OnePNCMinLocalResidual<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);
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,
};
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
// -*- 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/porousmediumflow/1pnc/model.hh>
namespace Dumux
{
/*!
* \ingroup OnePNCMinModel
* \brief Adaption of the fully implicit scheme to the
* one-phase n-component fully implicit model with additional solid/mineral phases.
*
* This model implements one-phase n-component flow of a compressible fluid composed of
* the n components \f$\kappa \f$ in combination with mineral precipitation and dissolution
* of the solid phases. The standard multiphase Darcy
* approach is used as the equation for the conservation of momentum:
* \f[
v = - \frac{k_{r}}{\mu} \mbox{\bf K}
\left(\text{grad}\, p - \varrho_{f} \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 ( \varrho_f X^\kappa \phi )}
{\partial t} - \text{div} \left\{ \varrho_f X^\kappa
\frac{k_{r}}{\mu} \mbox{\bf K}
(\text{grad}\, p - \varrho_{f} \mbox{\bf g}) \right\}
\nonumber \\ \nonumber \\
&-& \text{div} \left\{{\bf D_{pm}^\kappa} \varrho_{f} \text{grad}\, X^\kappa \right\}
- q_\kappa = 0 \qquad \kappa \in \{w, a,\cdots \}
\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 as spatial and the implicit Euler method as time
* discretization.
*
* The primary variables are the pressure \f$p\f$ and the mole fractions of the
* dissolved components \f$x^k\f$. 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>
{
using ThisType = Dumux::OnePNCMinModel<TypeTag>;
using ParentType = Dumux::OnePNCModel<TypeTag>;
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);
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 = FluidSystem::numSComponents,
pressureIdx = Indices::pressureIdx,
firstMoleFracIdx = Indices::firstMoleFracIdx,
phaseIdx = Indices::phaseIdx,
};
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("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]; });
vtkOutputModule.addSecondaryVariable("permeabilityFactor", [](const VolumeVariables& v)
{ return v.permeabilityFactor(); });
for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
vtkOutputModule.addSecondaryVariable("precipitateVolumeFraction_" + FluidSystem::phaseName(numPhases + sPhaseIdx),
[sPhaseIdx](const VolumeVariables& v)
{ return v.precipitateVolumeFraction(numPhases + sPhaseIdx); });
}
/*!
* \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);
}
private:
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"
#endif
// -**- 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/model.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
// -**- 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/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 SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
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),
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);
// calculate the remaining quantities
auto&& priVars = isBox ? elemSol[scv.index()] : elemSol[0];
// porosity evaluation
initialPorosity_ = problem.spatialParams().initialPorosity(element, scv);
sumPrecipitates_ = 0.0;
for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
{
precipitateVolumeFraction_[sPhaseIdx] = priVars[numComponents + sPhaseIdx];
sumPrecipitates_+= precipitateVolumeFraction_[sPhaseIdx];
}
// 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 = 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;