Commit 1b6b0611 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

use opm-material in the fully-implicit 2p2c test (first try)

Many of the required changes are at the top level of problem and
spatial parameters. Some are required in 2p2c's VolumeVariables.
parent d6e4323e
......@@ -29,6 +29,8 @@
#include <dumux/material/fluidstates/compositional.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
#include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
#include "properties.hh"
#include "indices.hh"
......@@ -93,10 +95,10 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem, useKelvinEquation> MiscibleMultiPhaseComposition;
typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
static_assert(useMoles || (!useMoles && useConstraintSolver),
"if UseMoles is set false, UseConstraintSolver has to be set to true");
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
......@@ -136,27 +138,26 @@ public:
// Second instance of a parameter cache.
// Could be avoided if diffusion coefficients also
// became part of the fluid state.
typename FluidSystem::ParameterCache paramCache;
typename FluidSystem::template ParameterCache<Scalar> paramCache;
paramCache.updateAll(fluidState_);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// relative permeabilities
Scalar kr;
if (phaseIdx == wPhaseIdx)
kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
kr = MaterialLaw::krw(materialParams, fluidState_);
else // ATTENTION: krn requires the wetting phase saturation
// as parameter!
kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
kr = MaterialLaw::krn(materialParams, fluidState_);
relativePermeability_[phaseIdx] = kr;
Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
// binary diffusion coefficients
diffCoeff_[phaseIdx] =
FluidSystem::binaryDiffusionCoefficient(fluidState_,
FluidSystem::diffusionCoefficient(fluidState_,
paramCache,
phaseIdx,
wCompIdx,
nCompIdx);
wCompIdx);
Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
}
......@@ -215,7 +216,7 @@ public:
// calculate capillary pressure
const auto& materialParams =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
Scalar pc = MaterialLaw::pc(materialParams, 1 - sn);
Scalar pc = MaterialLaw::pcnw(materialParams, fluidState);
if (formulation == pwsn) {
fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
......@@ -230,7 +231,7 @@ public:
/////////////
// calculate the phase compositions
/////////////
typename FluidSystem::ParameterCache paramCache;
typename FluidSystem::template ParameterCache<Scalar> paramCache;
//get the phase pressures and set the fugacity coefficients here if constraintsolver is not used
Scalar pn = 0;
......@@ -273,11 +274,12 @@ public:
//get the partial pressure of the main component of the the wetting phase ("H20") within the nonwetting (gas) phase == vapor pressure due to equilibrium
//note that in this case the fugacityCoefficient * pw is the vapor pressure (see implementation in respective fluidsystem)
Scalar partPressH2O = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
wPhaseIdx,
wCompIdx) * pw;
if (useKelvinEquation)
partPressH2O = FluidSystem::kelvinVaporPressure(fluidState, wPhaseIdx, wCompIdx);
//if (useKelvinEquation)
// partPressH2O = FluidSystem::kelvinVaporPressure(fluidState, wPhaseIdx, wCompIdx);
// get the partial pressure of the main component of the the nonwetting (gas) phase ("Air")
Scalar partPressAir = pn - partPressH2O;
......@@ -289,7 +291,7 @@ public:
// calculate the mole fractions of the components within the wetting phase
//note that in this case the fugacityCoefficient * pw is the Henry Coefficient (see implementation in respective fluidsystem)
Scalar xwn = partPressAir
/ (FluidSystem::fugacityCoefficient(fluidState,
/ (FluidSystem::fugacityCoefficient(fluidState, paramCache,
wPhaseIdx,nCompIdx)
* pw);
......@@ -323,7 +325,7 @@ public:
else
{
// setMassFraction() has only to be called 1-numComponents times
fluidState.setMassFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
//fluidState.setMassFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
}
// calculate the composition of the remaining phases (as
......@@ -351,7 +353,7 @@ public:
// If xww > 1 : gas is over-saturated with water vapor,
// condensation takes place (see switch criterion in model)
Scalar xww = xnw * pn
/ (FluidSystem::fugacityCoefficient(fluidState,
/ (FluidSystem::fugacityCoefficient(fluidState, paramCache,
wPhaseIdx,wCompIdx)
* pw);
......@@ -360,7 +362,7 @@ public:
//partialPressure = xnn * pn
//xwn = xnn * pn / Henry
// Henry = fugacityCoefficient * pw
Scalar xwn = xnn * pn / (FluidSystem::fugacityCoefficient(fluidState,
Scalar xwn = xnn * pn / (FluidSystem::fugacityCoefficient(fluidState, paramCache,
wPhaseIdx,nCompIdx)
* pw);
......@@ -388,7 +390,7 @@ public:
else // mass-fraction formulation
{
// setMassFraction() has only to be called 1-numComponents times
fluidState.setMassFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
//fluidState.setMassFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
}
// calculate the composition of the remaining phases (as
......@@ -410,7 +412,7 @@ public:
//first, xnw:
//psteam = xnw * pn = partial pressure of water in gas phase
//psteam = fugacityCoefficient * pw
Scalar xnw = (FluidSystem::fugacityCoefficient(fluidState,
Scalar xnw = (FluidSystem::fugacityCoefficient(fluidState, paramCache,
wPhaseIdx,wCompIdx)
* pw) / pn ;
......@@ -420,7 +422,7 @@ public:
// xwn = pn * xnn / Henry
// xnn = xwn * Henry / pn
// Henry = fugacityCoefficient * pw
Scalar xnn = xwn * (FluidSystem::fugacityCoefficient(fluidState,
Scalar xnn = xwn * (FluidSystem::fugacityCoefficient(fluidState, paramCache,
wPhaseIdx,nCompIdx)
* pw) / pn ;
......
......@@ -2,5 +2,5 @@ Module: dumux
Version: 2.11-git
Maintainer: dumux@listserv.uni-stuttgart.de
Depends: dune-common (>= 2.4.1) dune-grid (>= 2.4) dune-localfunctions (>= 2.4) dune-istl (>= 2.4)
Suggests: dune-alugrid (>=2.4) dune-pdelab (>=2.0) dune-multidomain dune-foamgrid (>=2.4) opm-grid
Suggests: dune-alugrid (>=2.4) dune-pdelab (>=2.0) dune-multidomain dune-foamgrid (>=2.4) opm-material
Whitespace-Hook: Yes
......@@ -27,6 +27,8 @@
#include <dumux/porousmediumflow/2p2c/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include "injectionspatialparams.hh"
......@@ -51,7 +53,11 @@ SET_TYPE_PROP(InjectionProblem, Problem, InjectionProblem<TypeTag>);
// Set fluid configuration
SET_TYPE_PROP(InjectionProblem,
FluidSystem,
FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false /*useComplexRelations*/>);
Opm::FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar)>);
SET_TYPE_PROP(InjectionProblem,
FluidState,
Opm::CompositionalFluidState<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, FluidSystem)>);
// Define whether mole(true) or mass (false) fractions are used
SET_BOOL_PROP(InjectionProblem, UseMoles, true);
......
......@@ -30,6 +30,8 @@
#include <dumux/material/spatialparams/implicit.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
#include <dumux/porousmediumflow/2p2c/implicit/properties.hh>
......@@ -40,6 +42,13 @@ namespace Dumux
template<class TypeTag>
class InjectionSpatialParams;
struct BCTraits {
using Scalar = double;
static const int numPhases = 2;
static const int wettingPhaseIdx = 0;
static const int nonWettingPhaseIdx = 1;
};
namespace Properties
{
// The spatial parameters TypeTag
......@@ -51,7 +60,7 @@ SET_TYPE_PROP(InjectionSpatialParams, SpatialParams, InjectionSpatialParams<Type
// Set the material law parameterized by absolute saturations
SET_TYPE_PROP(InjectionSpatialParams,
MaterialLaw,
EffToAbsLaw<RegularizedBrooksCorey<typename GET_PROP_TYPE(TypeTag, Scalar)> >);
Opm::EffToAbsLaw<Opm::RegularizedBrooksCorey<BCTraits> >);
}
/*!
......@@ -107,14 +116,14 @@ public:
lambdaSolid_ = 2.8;
// residual saturations
fineMaterialParams_.setSwr(0.2);
fineMaterialParams_.setSnr(0.0);
coarseMaterialParams_.setSwr(0.2);
coarseMaterialParams_.setSnr(0.0);
fineMaterialParams_.setResidualSaturation(0, 0.2);
fineMaterialParams_.setResidualSaturation(1, 0.0);
coarseMaterialParams_.setResidualSaturation(0, 0.2);
coarseMaterialParams_.setResidualSaturation(1, 0.0);
// parameters for the Brooks-Corey law
fineMaterialParams_.setPe(1e4);
coarseMaterialParams_.setPe(1e4);
fineMaterialParams_.setEntryPressure(1e4);
coarseMaterialParams_.setEntryPressure(1e4);
fineMaterialParams_.setLambda(2.0);
coarseMaterialParams_.setLambda(2.0);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment