Commit d555b16a authored by Timo Koch's avatar Timo Koch
Browse files

[mpnc] Replace preprocessor var USE_PCMAX by runtime param

parent 318373bb
......@@ -180,41 +180,40 @@ public:
//obtain parameters for interfacial area constitutive relations
const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
const Scalar Sw = fluidState.saturation(phase0Idx);
Scalar awn;
using AwnSurface = typename Problem::SpatialParams::AwnSurface;
awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
interfacialArea_[phase0Idx][phase1Idx] = awn;
interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
interfacialArea_[phase0Idx][phase0Idx] = 0.;
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc);
const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter // That value is obtained by regularization of the pc(Sw) function.
#if USE_PCMAX
const Scalar pcMax = problem.spatialParams().pcMax(element, scv, elemSol);
//I know the solid surface from the pore network. But it is more consistent to use the fitvalue.
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter
// That value is obtained by regularization of the pc(Sw) function.
static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
if (computeAwsFromAnsAndPcMax)
{
// I know the solid surface from the pore network. But it is more consistent to use the fit value.
const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
}
else
{
using AwsSurface = typename Problem::SpatialParams::AwsSurface;
const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
}
const Scalar aws = solidSurface_ - ans;
interfacialArea_[phase0Idx][sPhaseIdx] = aws;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
#else
using AwsSurface = typename Problem::SpatialParams::AwsSurface;
const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
const auto aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc);
interfacialArea_[phase0Idx][sPhaseIdx] = aws ;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
#endif
interfacialArea_[phase1Idx][sPhaseIdx] = ans;
interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
interfacialArea_[phase1Idx][phase1Idx] = 0.;
......@@ -259,7 +258,6 @@ private:
Scalar characteristicLength_;
Scalar factorEnergyTransfer_;
Scalar solidSurface_ ;
InterfacialAreasArray interfacialArea_;
};
......@@ -720,46 +718,43 @@ public:
const Scv& scv)
{
// obtain (standard) material parameters (needed for the residual saturations)
const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
//obtain parameters for interfacial area constitutive relations
const auto& aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
const Scalar Sw = fluidState.saturation(phase0Idx);
Scalar awn;
using AwnSurface = typename Problem::SpatialParams::AwnSurface;
awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
const auto awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc);
interfacialArea_[phase0Idx][phase1Idx] = awn;
interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
interfacialArea_[phase0Idx][phase0Idx] = 0.;
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc);
const auto& aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
const auto ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, Sw, pc);
// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter.
// That value is obtained by regularization of the pc(Sw) function.
#if USE_PCMAX
const Scalar pcMax = problem.spatialParams().pcMax(element, scv, elemSol);
// I know the solid surface from the pore network. But it is more consistent to use the fit value.
using AnsSurface = typename Problem::SpatialParams::AnsSurface;
solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
if (computeAwsFromAnsAndPcMax)
{
// I know the solid surface from the pore network. But it is more consistent to use the fit value.
const Scalar pcMax = aWettingNonWettingSurfaceParams.pcMax();
const auto solidSurface = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - ans;
}
else
{
using AwsSurface = typename Problem::SpatialParams::AwsSurface;
const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
interfacialArea_[phase0Idx][sPhaseIdx] = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, materialParams, Sw, pc);
}
const Scalar aws = solidSurface_ - ans;
interfacialArea_[phase0Idx][sPhaseIdx] = aws;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
#else
using AwsSurface = typename Problem::SpatialParams::AwsSurface;
const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams();
const auto aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc );
interfacialArea_[phase0Idx][sPhaseIdx] = aws ;
interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
#endif
interfacialArea_[phase1Idx][sPhaseIdx] = ans;
interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
......@@ -819,7 +814,6 @@ private:
Scalar characteristicLength_;
Scalar factorEnergyTransfer_;
Scalar factorMassTransfer_;
Scalar solidSurface_ ;
InterfacialAreasArray interfacialArea_;
};
......
......@@ -35,12 +35,6 @@
#ifndef DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
#define DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
// this sets that the relation using pc_max is used.
// i.e. - only parameters for awn, ans are given,
// - the fit for ans involves the maximum value for pc, where Sw, awn are zero.
// setting it here, because it impacts volume variables and spatialparameters
#define USE_PCMAX 1
#include <dune/grid/yaspgrid.hh>
#include <dumux/common/properties.hh>
......
......@@ -128,9 +128,17 @@ public:
materialParamsFF_.setPe(0.);
// determine maximum capillary pressure for wetting-nonwetting surface
/* Of course physically there is no such thing as a maximum capillary pressure.
* The parametrization (VG/BC) goes to infinity and physically there is only one pressure
* for single phase conditions.
* Here, this is used for fitting the interfacial area surface: the capillary pressure,
* where the interfacial area is zero.
* Technically this value is obtained as the capillary pressure of saturation zero.
* This value of course only exists for the case of a regularized pc-Sw relation.
*/
using TwoPLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
pcMax_ = TwoPLaw::pc(materialParamsPM_, /*sw = */0.0);
aWettingNonWettingSurfaceParams_.setPcMax(pcMax_);
const auto pcMax = TwoPLaw::pc(materialParamsPM_, /*sw = */0.0);
aWettingNonWettingSurfaceParams_.setPcMax(pcMax);
// wetting-non wetting: surface which goes to zero on the edges, but is a polynomial
aWettingNonWettingSurfaceParams_.setA1(aWettingNonWettingA1_);
......@@ -228,7 +236,7 @@ public:
else if (inPM_(globalPos))
return aWettingNonWettingSurfaceParams_ ;
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
}
/*!\brief Returns a reference to the container object for the
* parametrization of the surface between non-Wetting and solid phase.
......@@ -250,26 +258,24 @@ public:
else if (inPM_(globalPos))
return aNonWettingSolidSurfaceParams_ ;
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
}
/*!\brief Returns the maximum capillary pressure for the given pc-Sw curve
/*!\brief Returns a reference to the container object for the
* parametrization of the surface between wetting and solid phase.
*
* Of course physically there is no such thing as a maximum capillary pressure.
* The parametrization (VG/BC) goes to infinity and physically there is only one pressure
* for single phase conditions.
* Here, this is used for fitting the interfacial area surface: the capillary pressure,
* where the interfacial area is zero.
* Technically this value is obtained as the capillary pressure of saturation zero.
* This value of course only exists for the case of a regularized pc-Sw relation.
* The position is determined based on the coordinate of
* the vertex belonging to the considered sub-control volume.
* \param element The finite element
* \param scv The sub-control volume
* \param elemSol The element solution
*/
template<class ElementSolution>
const Scalar pcMax(const Element &element,
const SubControlVolume &scv,
const ElementSolution &elemSol) const
{ return aWettingNonWettingSurfaceParams_.pcMax() ; }
const AwsSurfaceParams& aWettingSolidSurfaceParams(const Element &element,
const SubControlVolume &scv,
const ElementSolution &elemSol) const
{
DUNE_THROW(Dune::NotImplemented, "wetting-solid-interface surface params");
}
/*!
* \brief Returns the characteristic length for the mass transfer.
......@@ -367,8 +373,6 @@ private:
AwnSurfaceParams aWettingNonWettingSurfaceParamsFreeFlow_;
AnsSurfaceParams aNonWettingSolidSurfaceParamsFreeFlow_ ;
Scalar pcMax_ ;
// Porous Medium Domain
Scalar intrinsicPermeabilityPM_ ;
Scalar porosityPM_ ;
......@@ -401,6 +405,6 @@ private:
std::vector<Scalar> gridVector_;
};
}
} // end namespace Dumux
#endif // GUARDIAN
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