Dear developers, maintainers and users,
I'm encountering some difficulties in applying a new Material law which
is not regularized and doesn't use EffToAbsLaw.
I did as suggested previously, but when I comment the following lines I
get an error:
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
//typedef EffToAbsLaw<RawMaterialLaw> type;
In addition, when I use:
#include <dumux/io/plotmateriallaw.hh>
to plot Pc-Sw and kr profiles, I get constant behaviors.
I start to suspect that I didn't implement the Material Law correctly
and these effects are correlated.
To be more clear I add my files.
Can you give me some hint about mistakes, please?
Thanks for help,
Lorenzo
On 17/10/2018 16:53, Timo Koch wrote:
Hi Lorenzo,
On 17.10.2018 15:21, lc wrote:
Dear developers and users,
I have some other questions.
1) I successfully ran the tutorial_implicit example from the
handbook 2.12. Now, I need to modify such test case in order to be
consistent with my requirements.
In particular, I need to use a different relative permeability and
capillarity pressure laws. In practice:
for the capillarity:
pc(s) = constant1 * (s-s_min)^alpha
for the permeabilities:
kr_w(s) = constant2 * (s-s_min)^beta
kr_n(s) = constant3 * (s_max-s)^beta
with s_min and s_max other two given constants.
So, which is the best way to introduce such laws? And what about the
"regularization"? It should be applied to them also or I use them "as
they are"?
You would create a new class "MyLaw" (or whatever fits best) and a
class "MyLawParams" that has the same interface as BrooksCorey and
BrooksCoreyParams. In particular you need to implement, the
pc-function and the kr-functions. Then you change the property
"MaterialLaw" (the property is usually either set in problem.hh or
spatialparams.hh) to use your class instead of the BrooksCorey class.
You don't need to use regularization, if you don't want to, the same
goes for the generalized EffToAbsLaw which takes care of converting
absolute to effective saturations and the other way around. So you can
just replace the whole thing by your new class.
The parameters go into your "MyLawParams" class and you need an
interface to set them there.
2) About boundary conditions.
Where can I find all the BC that can I impose? Can I impose an
injection velocity and a constant saturation on the left and a
constant pressure and a zero saturation flux on the right? As well as
normal velocity and zero saturation flux up and down (of a 2d
rectangle domain)?
You can impose whatever boundary conditions you want to impose, given
that they are valid for your type of problem. As boundary types Dumux
differs between Dirichlet boundaries and Neumann boundaries. The
Neumann boundary type includes Robin boundary conditions types (which
are called solution dependent Neumann in the code sometimes). At the
boundaries marked as Dirichlet boundary the
"dirichlet"/"dirichletAtPos" function is called, on boundaries marked
as Neumann, the "neumann"/"neumannAtPos" function is called. You can
specify where to use which type of boundary condition using the global
position argument.
How you impose boundary conditions also depends on the discretization
scheme you are using. Are you using the Box method or the
cell-centered TPFA method? You cannot impose velocities directly as
velocity is not a unknown in your equations (for the discretization
schemes available in Dumux), you can however impose fluxes on Neumann
boundaries, and compute the flux with your given velocity. You can
enforce the pressure and saturation as Dirichlet boundary conditions.
Also note that Dumux is short before a major version update (from 2.12
to 3.0) with a lot of new features / changes. There will be continued
bugfix support for 2.12, as far as we can handle it. But you might
want to have a look at the new version (master branch). The simulation
steps are easier to understand in the new version. Everything I
explained above also applies for the new version.
There are exercises in the dumux-course module I told you about in my
last mail.
Best wishes
Timo
Thank you,
Lorenzo
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
_______________________________________________________________
Timo Koch phone: +49 711 685 64676
IWS, Universität Stuttgart fax: +49 711 685 60430
Pfaffenwaldring 61 email:[email protected]
D-70569 Stuttgart url:www.hydrosys.uni-stuttgart.de
_______________________________________________________________
// -*- 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 spatial parameters for the fully coupled tutorial problem
* which uses the twophase box model.
*/
#ifndef DUMUX_TUTORIAL_SPATIAL_PARAMS_IMPLICIT_HH
#define DUMUX_TUTORIAL_SPATIAL_PARAMS_IMPLICIT_HH
#include <iostream>
// include parent spatialparameters
#include <dumux/material/spatialparams/implicit.hh>
// include material laws
//#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
//#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
//#include <dumux/material/fluidmatrixinteractions/2p/chebyshev.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
//#include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
//#include <dumux/material/fluidmatrixinteractions/2p/mybrookscorey.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
//#include <dumux/material/fluidmatrixinteractions/2p/mod_regularizedbrookscorey.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
#include <dumux/material/fluidmatrixinteractions/2p/mod_brookscorey.hh> /*@\label{tutorial-implicit:rawLawInclude}@*/
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/io/plotmateriallaw.hh>
namespace Dumux {
//forward declaration
template<class TypeTag>
class TutorialSpatialParamsImplicit;
namespace Properties
{
// The spatial parameters TypeTag
NEW_TYPE_TAG(TutorialSpatialParamsImplicit);/*@\label{tutorial-implicit:define-spatialparameters-typetag}@*/
// Set the spatial parameters
SET_TYPE_PROP(TutorialSpatialParamsImplicit, SpatialParams,
TutorialSpatialParamsImplicit<TypeTag>); /*@\label{tutorial-implicit:set-spatialparameters}@*/
// Set the material law
SET_PROP(TutorialSpatialParamsImplicit, MaterialLaw)
{
private:
// material law typedefs
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
// select material law to be used
// typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
// typedef LinearMaterial<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
// typedef Chebyshev<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
// typedef BrooksCorey<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
// typedef ModRegularizedBrooksCorey<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
typedef ModBrooksCorey<Scalar> RawMaterialLaw; /*@\label{tutorial-implicit:rawlaw}@*/
public:
// adapter for absolute law
typedef EffToAbsLaw<RawMaterialLaw> type; /*@\label{tutorial-implicit:eff2abs}@*/
};
}
/*!
* \ingroup TwoPBoxModel
*
* \brief The spatial parameters for the fully coupled tutorial problem
* which uses the twophase box model.
*/
template<class TypeTag>
class TutorialSpatialParamsImplicit: public ImplicitSpatialParams<TypeTag> /*@\label{tutorial-implicit:tutorialSpatialParameters}@*/
{
// Get informations for current implementation via property system
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
enum
{
dim = Grid::dimension
};
// Get object types for function arguments
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename Grid::Traits::template Codim<0>::Entity Element;
public:
// get material law from property system
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
// determine appropriate parameters depending on selected materialLaw
typedef typename MaterialLaw::Params MaterialLawParams; /*@\label{tutorial-implicit:matLawObjectType}@*/
/*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
* on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function intrinsicPermeabilityAtPos(const GlobalPosition& globalPos)
* could be defined, where globalPos is the vector including the global coordinates
* of the finite volume.
*/
const Dune::FieldMatrix<Scalar, dim, dim> &intrinsicPermeability(const Element &element, /*@\label{tutorial-implicit:permeability}@*/
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
//std::cout << "Permeability = "<< K_ << std::endl;
return K_;
}
/*! Defines the porosity \f$[-]\f$ of the porous medium depending
* on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function porosityAtPos(const GlobalPosition& globalPos)
* could be defined, where globalPos is the vector including the global coordinates
* of the finite volume.
*/
Scalar porosity(const Element &element, /*@\label{tutorial-implicit:porosity}@*/
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{ return 0.255; }
/*! Returns the parameter object for the material law (i.e. Brooks-Corey)
* depending on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function materialLawParamsAtPos(const GlobalPosition& globalPos)
* could be defined, where globalPos is the vector including the global coordinates
* of the finite volume.
*/
const MaterialLawParams& materialLawParams(const Element &element, /*@\label{tutorial-implicit:matLawParams}@*/
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return materialParams_;
}
/*!
* \brief The constructor
*
* \param gridView The grid view
* \param K_(0) The permeability tensor
*/
TutorialSpatialParamsImplicit(const GridView& gridView) :
ImplicitSpatialParams<TypeTag>(gridView),
K_(0)
{
//set main diagonal entries of the permeability tensor to a value
//setting to one value means: isotropic, homogeneous
for (int i = 0; i < dim; i++)
K_[i][i] = 2e-13;
//set residual saturations
materialParams_.setSwr(0.2121); /*@\label{tutorial-implicit:setLawParams}@*/
materialParams_.setSnr(1.-0.2121);
//parameters of Mod Brooks & Corey Law
//materialParams_.setPe(500.0);
//materialParams_.setPe(80000.);
//materialParams_.setLambda(2);
//materialParams_.setLambda(-0.476190476190476);
//parameters of Linear Law
//materialParams_.setMaxPc(2000.0);
//materialParams_.setEntryPc(0.);
//parameters of Chebyshev Law
//materialParams_.setPe(80000.);
//materialParams_.setAlpha(2.1);
//materialParams_.setSmin(0.2121);
//materialParams_.setSmax(0.7856);
//materialParams_.setLambda(-0.476190476190476);
plotFluidMatrixInteractions_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Output, PlotFluidMatrixInteractions);
}
/*!
* \brief This is called from the problem and creates a gnuplot output
* of e.g the pc-Sw curve
*/
void plotMaterialLaw()
{
// some problem in retrieving the correct law
// or some problem in the law itself!
PlotMaterialLaw<TypeTag> plotMaterialLaw;
GnuplotInterface<Scalar> gnuplot(plotFluidMatrixInteractions_);
gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_);
plotMaterialLaw.addpcswcurve(gnuplot, materialParams_, 0.212101, 0.7856, "Pc(Sw)", "w lp");
gnuplot.setOption("set xrange [0:1]");
gnuplot.setOption("set label \"residual\\nsaturation\" at 0.1,100000 center");
gnuplot.plot("pc-Sw");
gnuplot.resetAll();
plotMaterialLaw.addkrcurves(gnuplot, materialParams_, 0.212101, 0.7856, "kr(Sw)");
gnuplot.plot("kr");
}
private:
Dune::FieldMatrix<Scalar, dim, dim> K_;
// Object that holds the values/parameters of the selected material law.
MaterialLawParams materialParams_; /*@\label{tutorial-implicit:matParamsObject}@*/
bool plotFluidMatrixInteractions_;
};
} // end namespace
#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 Tutorial problem for a fully coupled twophase box model.
*/
#ifndef DUMUX_TUTORIAL_PROBLEM_IMPLICIT_HH // guardian macro /*@\label{tutorial-implicit:guardian1}@*/
#define DUMUX_TUTORIAL_PROBLEM_IMPLICIT_HH // guardian macro /*@\label{tutorial-implicit:guardian2}@*/
// The numerical model
#include <dumux/porousmediumflow/2p/implicit/model.hh>
// The base porous media box problem
#include <dumux/porousmediumflow/implicit/problem.hh>
// Spatially dependent parameters
#include "tutorialspatialparams_implicit.hh"
// The components that are used
//#include <dumux/material/components/chebyshev_h2o.hh>
//#include <dumux/material/components/chebyshev_oil.hh>
//#include <dumux/material/components/h2o.hh>
//#include <dumux/material/components/lnapl.hh>
#include <dumux/material/components/pseudooil.hh>
#include <dumux/material/components/pseudoh2o.hh>
namespace Dumux{
// Forward declaration of the problem class
template <class TypeTag>
class TutorialProblemImplicit;
namespace Properties {
// Create a new type tag for the problem
NEW_TYPE_TAG(TutorialProblemImplicit, INHERITS_FROM(BoxTwoP, TutorialSpatialParamsImplicit)); /*@\label{tutorial-implicit:create-type-tag}@*/
// Set the "Problem" property
SET_PROP(TutorialProblemImplicit, Problem) /*@\label{tutorial-implicit:set-problem}@*/
{ typedef TutorialProblemImplicit<TypeTag> type;};
// Set grid and the grid creator to be used
#if HAVE_DUNE_ALUGRID /*@\label{tutorial-implicit:set-grid}@*/
//SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>); /*@\label{tutorial-implicit:set-grid-ALU}@*/
// The following line is needed in order to use gmsh grids
SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::simplex, Dune::conforming>); /*@\label{tutorial-implicit:set-grid-ALU}@*/
#elif HAVE_UG
SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::YaspGrid<2>);
#endif // HAVE_DUNE_ALUGRID
// Set the wetting phase
SET_PROP(TutorialProblemImplicit, WettingPhase) /*@\label{tutorial-implicit:2p-system-start}@*/
{
private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
//public: typedef FluidSystems::LiquidPhase<Scalar, ChebyshevH2O<Scalar> > type; /*@\label{tutorial-implicit:wettingPhase}@*/
public: typedef FluidSystems::LiquidPhase<Scalar, PseudoH2O<Scalar> > type;
};
// Set the non-wetting phase
SET_PROP(TutorialProblemImplicit, NonwettingPhase)
{
private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
//public: typedef FluidSystems::LiquidPhase<Scalar, ChebyshevOil<Scalar> > type; /*@\label{tutorial-implicit:nonwettingPhase}@*/
//public: typedef FluidSystems::LiquidPhase<Scalar, LNAPL<Scalar> > type;
public: typedef FluidSystems::LiquidPhase<Scalar, PseudoOil<Scalar> > type; /*@\label{tutorial-implicit:nonwettingPhase}@*/
}; /*@\label{tutorial-implicit:2p-system-end}@*/
SET_TYPE_PROP(TutorialProblemImplicit, FluidSystem, TwoPImmiscibleFluidSystem<TypeTag>);/*@\label{tutorial-implicit:set-fluidsystem}@*/
// Disable gravity
SET_BOOL_PROP(TutorialProblemImplicit, ProblemEnableGravity, false); /*@\label{tutorial-implicit:gravity}@*/
}
/*!
* \ingroup TwoPBoxModel
*
* \brief Tutorial problem for a fully coupled twophase box model.
*/
template <class TypeTag>
class TutorialProblemImplicit : public ImplicitPorousMediaProblem<TypeTag> /*@\label{tutorial-implicit:def-problem}@*/
{
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
// Grid dimension
enum { dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
// Types from DUNE-Grid
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
// Dumux specific types
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
public:
TutorialProblemImplicit(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, gridView)
, eps_(1e-6)
{
#if !(HAVE_DUNE_ALUGRID || HAVE_UG)
std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl;
#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG)
this->spatialParams().plotMaterialLaw();
}
//! Specifies the problem name. This is used as a prefix for files
//! generated by the simulation.
std::string name() const
{ return "tutorial_implicit"; }
//! Returns true if a restart file should be written.
bool shouldWriteRestartFile() const /*@\label{tutorial-implicit:restart}@*/
{ return false; }
//! Returns true if the current solution should be written to disk
//! as a VTK file
bool shouldWriteOutput() const /*@\label{tutorial-implicit:output}@*/
{
return
(this->timeManager().timeStepIndex() > 0
&& (this->timeManager().timeStepIndex() % 1 == 0));
}
//! Returns the temperature within a finite volume. We use constant
//! 10 degrees Celsius.
Scalar temperature() const
{ return 273.15 + 10.; }
//! Specifies which kind of boundary condition should be used for
//! which equation for a finite volume on the boundary.
void boundaryTypes(BoundaryTypes &bcTypes, const Vertex &vertex) const
{
const GlobalPosition &globalPos = vertex.geometry().center();
if (globalPos[0] < eps_) // Dirichlet conditions on left boundary
bcTypes.setAllDirichlet();
else // neuman for the remaining boundaries
bcTypes.setAllNeumann();
}
//! Evaluates the Dirichlet boundary conditions for a finite volume
//! on the grid boundary. Here, the 'values' parameter stores
//! primary variables.
void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
values[Indices::pwIdx] = 200.0e3; // 200 kPa = 2 bar
values[Indices::snIdx] = 0.0; // 0 % oil saturation on left boundary
}
//! Evaluates the boundary conditions for a Neumann boundary
//! segment. Here, the 'values' parameter stores the mass flux in
//! [kg/(m^2 * s)] in normal direction of each phase. Negative
//! values mean influx.
void neumann(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &intersection,
int scvIdx,
int boundaryFaceIdx) const
{
const GlobalPosition &globalPos =
fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar right = this->bBoxMax()[0];
// extraction of oil on the right boundary for approx. 1.e6 seconds
if (globalPos[0] > right - eps_) {
// oil outflux of 30 g/(m * s) on the right boundary.
values[Indices::contiWEqIdx] = 0;
values[Indices::contiNEqIdx] = 3e-2;
} else {
// no-flow on the remaining Neumann-boundaries.
values[Indices::contiWEqIdx] = 0;
values[Indices::contiNEqIdx] = 0;
}
}
//! Evaluates the initial value for a control volume. For this
//! method, the 'values' parameter stores primary variables.
void initial(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
values[Indices::pwIdx] = 200.0e3; // 200 kPa = 2 bar
//values[Indices::snIdx] = 1.0;
// To retrieve x- and y- local coordinates
const auto pos = fvGeometry.subContVol[scvIdx].global;
double x_0 = 0;
double smin = 0.2121;
double smax = 0.7856;
double exp_arg = 0.5 - pow((pos[0]-x_0)/(2.0+cos(3.0*sqrt(std::abs(pos[1]-75.0)))), 2);
double exp_value = exp_arg >= -100 ? exp(exp_arg) : 0;
//double small_number = 1e-4;
double small_number = 1e-1;
values[Indices::snIdx] = std::min(smax, smin + exp_value + small_number);
}
//! Evaluates the source term for all phases within a given
//! sub-control-volume. In this case, the 'values' parameter
//! stores the rate mass generated or annihilated per volume unit
//! in [kg / (m^3 * s)]. Positive values mean that mass is created.
void source(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
values[Indices::contiWEqIdx] = 0.0;
values[Indices::contiNEqIdx] = 0.0;
}
private:
// small epsilon value
Scalar eps_;
};
}
#endif
[TimeManager]
TEnd = 1000000 # duration of the simulation [s]
DtInitial = 10 # initial time step size [s]
[Grid]
File = ./chebyshev_unstruct.msh
Verbosity = 1 # verbose grid file parsing
BoundarySegments = 1 # read parametrized boundary segments
DomainMarkers = 1 # read domain markers (gmsh physical entities)
#UpperRight = 150 150 # x-/y-coordinates of the upper-right corner of the grid
[m]
#Cells = 150 15 # x-/y-resolution of the grid
#CellType = Simplex
[Problem]
OnlyPlotMaterialLaws = false
EnableGravity = false
// -*- 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 Main file of the tutorial for a fully coupled twophase box model.
*/
#include <config.h> /*@\label{tutorial-implicit:include-begin}@*/
#include "tutorialproblem_implicit.hh" /*@\label{tutorial-implicit:include-problem-header}@*/
#include <dumux/common/start.hh> /*@\label{tutorial-implicit:include-end}@*/
//! Prints a usage/help message if something goes wrong or the user asks for help
void usage(const char *progName, const std::string &errorMsg) /*@\label{tutorial-implicit:usage-function}@*/
{
std::cout
<< "\nUsage: " << progName << " [options]\n";
if (errorMsg.size() > 0)
std::cout << errorMsg << "\n";
std::cout
<< "\n"
<< "The list of mandatory arguments for this program is:\n"
<< "\t-TEnd The end of the simulation [s]\n"
<< "\t-DtInitial The initial timestep size [s]\n"
<< "\t-Grid.UpperRight The x-/y-coordinates of the grid's upper-right corner [m]\n"
<< "\t-Grid.Cells The grid's x-/y-resolution\n"
<< "\n";
}
int main(int argc, char** argv)
{
typedef TTAG(TutorialProblemImplicit) TypeTag; /*@\label{tutorial-implicit:set-type-tag}@*/
return Dumux::start<TypeTag>(argc, argv, usage); /*@\label{tutorial-implicit:call-start}@*/
}
// -*- 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 Implementation of the capillary pressure and
* relative permeability <-> saturation relations according to the Chebyshev model.
*
*/
#ifndef DUMUX_MOD_BROOKS_COREY_HH
#define DUMUX_MOD_BROOKS_COREY_HH
#include "mod_brookscoreyparams.hh"
#include <algorithm>
#include <cmath>
#include <iostream>
namespace Dumux
{
/*!
* \ingroup fluidmatrixinteractionslaws
*
* \brief Implementation of the Brooks-Corey capillary pressure <->
* saturation relation. This class bundles the "raw" curves
* as static members and doesn't concern itself converting
* absolute to effective saturations and vice versa.
*
* For general info: EffToAbsLaw
*
*\see BrooksCoreyParams
*/
template <class ScalarT, class ParamsT = ModBrooksCoreyParams<ScalarT> >
class ModBrooksCorey
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
/*!
* \brief The capillary pressure-saturation curve according to the Chebyshev law.
*
* The Chebyshev law is given by:
*
* \f$\mathrm{
p_C = p_e\overline{S}_w^{-1/\lambda}
* }\f$
*
* \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Capillary pressure calculated by Brooks & Corey constitutive relation.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
// static Scalar pc(const Params ¶ms, Scalar swe)
// {
// using std::pow;
// using std::min;
// using std::max;
//
// //swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
// swe = min(max(swe, 0.2121), 0.7856); // the equation below is only defined for 0.0 <= sw <= 1.0
//
// // return params.pe()*pow(swe, -1.0/params.lambda());
// // return params.pe()*pow(-log(swe - 0.2121), -1.0/params.lambda());
// // return params.pe()*pow(-log(swe - 0.2121), 2.1);
// // return 0.8*1e5*pow(-log(swe - 0.2121), 2.1);
// return 0.8*1e5*pow(-log(swe), 2.1);
// }
static Scalar pc(const Params ¶ms, Scalar Sw)
{
using std::pow;
using std::min;
using std::max;
// eps_smin introduced to avoid pc(0.)
// which would give -inf from the log.
const Scalar eps_smin = 1e-6;
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
if (Sw >= smax)
return smax;
else if (Sw <= smin)
return smin;
// TODO: add unit test here ...
//std::cout << Sw << " " << 80000*pow(-log(0.24 + eps_smin - smin), 2.1) << std::endl;
return 80000*(pow(-log(Sw + eps_smin - smin), 2.1));
}
/*!
* \brief The saturation-capillary pressure curve according to Brooks & Corey.
*
* This is the inverse of the capillary pressure-saturation curve:
* \f$\mathrm{
\overline{S}_w = (\frac{p_C}{p_e})^{-\lambda}}\f$
*
* \param pc Capillary pressure \f$\mathrm{[p_C]}\f$ in \f$\mathrm{[Pa]}\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Effective wetting phase saturation calculated as inverse of BrooksCorey constitutive relation.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar sw(const Params ¶ms, Scalar pc)
{
using std::pow;
using std::max;
pc = max(pc, 0.0); // the equation below is undefined for negative pcs
const Scalar smin = 0.2121;
//return pow(pc/params.pe(), -params.lambda());
return smin + exp(-(pow(pc/(80000)), 1./2.1));
//DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::Sw");
}
/*!
* \brief The partial derivative of the capillary
* pressure w.r.t. the effective saturation according to Brooks & Corey.
*
* This is equivalent to
* \f$\mathrm{
\frac{\partial p_C}{\partial \overline{S}_w} =
-\frac{p_e}{\lambda} \overline{S}_w^{-1/\lambda - 1}
}\f$
*
* \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of \f$\mathrm{[p_c]}\f$ w.r.t. effective saturation according to Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar dpc_dswe(const Params ¶ms, Scalar swe)
{
using std::pow;
using std::min;
using std::max;
swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
//return - params.pe()/params.lambda() * pow(swe, -1/params.lambda() - 1);
//return 0.;
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dpc_dsw");
}
/*!
* \brief The partial derivative of the effective
* saturation w.r.t. the capillary pressure according to Brooks & Corey.
*
* \param pc Capillary pressure \f$\mathrm{[p_c]}\f$ in \f$\mathrm{[Pa]}\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of effective saturation w.r.t. \f$\mathrm{[p_c]}\f$ according to Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar dswe_dpc(const Params ¶ms, Scalar pc)
{
using std::pow;
using std::max;
pc = max(pc, 0.0); // the equation below is undefined for negative pcs
//return -params.lambda()/params.pe() * pow(pc/params.pe(), - params.lambda() - 1);
//return 0.;
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dswe_dpc");
}
/*!
* \brief The relative permeability for the wetting phase of
* the medium implied by the Brooks-Corey
* parameterization.
*
* \param swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
* and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar krw(const Params ¶ms, Scalar swe)
{
using std::pow;
using std::min;
using std::max;
//swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
if (swe >= smax) {
//std::cout << " swe = " << swe << std::endl;
return smax;
} else if (swe <= smin) {
//std::cout << " swe = " << swe << std::endl;
return smin;
}
//std::cout << swe << " " << 5.0*pow((swe - smin), 5.) << std::endl;
// return pow(swe, 2.0/params.lambda() + 3);
return 5.0*pow((swe - smin), 5.);
}
// static Scalar krw(const Params ¶ms, Scalar Sw)
// {
// return kr_(Sw);
// }
/*!
* \brief The derivative of the relative permeability for the
* wetting phase with regard to the wetting saturation of the
* medium implied by the Brooks-Corey parameterization.
*
* \param swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
* and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the wetting phase w.r.t. effective wetting phase
* saturation calculated as implied by Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar dkrw_dswe(const Params ¶ms, Scalar swe)
{
using std::pow;
using std::min;
using std::max;
swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
//return (2.0/params.lambda() + 3)*pow(swe, 2.0/params.lambda() + 2);
//return 25.0*pow((swe - 0.2121), 4.);
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dkrw_dswe");
}
/*!
* \brief The relative permeability for the non-wetting phase of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar krn(const Params ¶ms, Scalar swe)
{
using std::pow;
using std::min;
using std::max;
//swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
if (swe >= smax) {
//std::cout << " swe = " << swe << std::endl;
return smax;
} else if (swe <= smin) {
//std::cout << " swe = " << swe << std::endl;
return smin;
}
//std::cout << swe << " " << 16.0*pow((smax - swe), 5.) << std::endl;
//const Scalar exponent = 2.0/params.lambda() + 1;
//const Scalar tmp = 1.0 - swe;
//return tmp*tmp*(1.0 - pow(swe, exponent));
return 16.0*pow((smax - swe), 5.);
}
// static Scalar krn(const Params ¶ms, Scalar Sw)
// {
// Scalar Sn = 1 - Sw;
// return kr_(Sn);
// }
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
* and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the non-wetting phase w.r.t. effective wetting phase
* saturation calculated as implied by Brooks & Corey.
*
* \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
* by clamping the input.
*/
static Scalar dkrn_dswe(const Params ¶ms, Scalar swe)
{
using std::pow;
using std::min;
using std::max;
swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
//return 2.0*(swe - 1)*(1 + pow(swe, 2.0/params.lambda())*(1.0/params.lambda() + 1.0/2
// - swe*(1.0/params.lambda() + 1.0/2)
// )
// );
//return 80.0*pow((0.7856 - swe), 4.);
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dkrn_dswe");
}
// private:
//
// Scalar kr_(Scalar S)
// {
// const Scalar smin = 0.2121;
// const Scalar smax = 0.7856;
//
// if (S >= smax)
// return smax;
// else if (S <= smin)
// return smin;
//
// return 16.0*pow((0.7856 - S), 5.);
// }
//
// Scalar kr_(Scalar S)
// {
// const Scalar smin = 0.2121;
// const Scalar smax = 0.7856;
//
// if (S >= smax)
// return smax;
// else if (S <= smin)
// return smin;
//
// return 5.0*pow((S-smin), 5.);
// }
};
} // end namespace Dumux
#endif // MOD_BROOKS_COREY_HH
// -*- 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 Specification of the material parameters
* for the Brooks Corey constitutive relations.
*/
#ifndef DUMUX_MOD_BROOKS_COREY_PARAMS_HH
#define DUMUX_MOD_BROOKS_COREY_PARAMS_HH
#include <dumux/common/valgrind.hh>
namespace Dumux
{
/*!
* \brief Specification of the material parameters
* for the Brooks Corey constitutive relations.
*
* \ingroup fluidmatrixinteractionsparams
*
*\see BrooksCorey
*/
template <class ScalarT>
class ModBrooksCoreyParams
{
public:
typedef ScalarT Scalar;
ModBrooksCoreyParams()
{
Valgrind::SetUndefined(*this);
}
ModBrooksCoreyParams(Scalar pe, Scalar lambda)
: pe_(pe), lambda_(lambda)
{
}
/*!
* \brief Returns the entry pressure in \f$\mathrm{[Pa]}\f$
*/
Scalar pe() const
{ return pe_; }
/*!
* \brief Set the entry pressure in \f$\mathrm{[Pa]}\f$]
*/
void setPe(Scalar v)
{ pe_ = v; }
/*!
* \brief Returns the lambda shape parameter \f$\mathrm{[-]}\f$
*/
Scalar lambda() const
{ return lambda_; }
/*!
* \brief Set the lambda shape parameter \f$\mathrm{[-]}\f$
*/
void setLambda(Scalar v)
{ lambda_ = v; }
private:
Scalar pe_;
Scalar lambda_;
};
} // namespace Dumux
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
