Hello,
I enclose all the files that I use.
Up to now in tutorialproblem_implicit.hh I introduce a double
small_number = 1e-6; which works ok when I initialize the oil saturation
profile (snIdx) but it is not effective on swIdx.
Also as I general question I'm wondering which is the correct model to
approach this problem. I see that in this case it is used a different model.
NEW_TYPE_TAG(McWhorterProblem, INHERITS_FROM(FVPressureTwoP,
FVTransportTwoP, IMPESTwoP, BuckleyLeverettSpatialParams));
And in general, in terms of accuracy and speed should I prefer
sequential or implicit setup?
Thanks a lot,
Lorenzo
On 02/11/2018 13:56, Timo Koch wrote:
Hi Lorenzo,
it all depends on your setup. Probably something wrong with your
custom material law.
What have you tried so far to fix your issue?
Timo
On 2. Nov 2018, at 11:53, lc <lorenzo.camp...@uniroma1.it
<mailto:lorenzo.camp...@uniroma1.it>> wrote:
Hello Dennis,
thank you. I did as you suggested but as I do this, the Newton solver
does not converge anymore on any grid level.
Do you have any hints?
Thank you,
Lorenzo
On 02/11/2018 10:45, Dennis Gläser wrote:
Good morning Lorenzo,
per default, the formulation for the 2p model is pw-sn. That means
your primary variables are the water pressure and the non-wetting
phase saturation (in your case oil I assume). Therefore,
Indices::swIdx does not exists as it is not part of your primary
variables.
Since you are considering a two-phase system, just set the oil
saturation by doing
*values[Indices::snIdx] = 1.0 - std::min(smax, smin + exp_value +
small_number);*
This is what you want I guess.
Best wishes,
Dennis
On 02.11.2018 08:37, lc wrote:
Good morning,
I'd like to ask how to impose an initial condition on the water
saturation and not on the oil saturation.
Actually, what I implemented is the following:
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.+cos(3.*sqrt(std::abs(pos[1]-75.0)))), 2.) );
double exp_value = exp_arg >= -100 ? exp(exp_arg) : 0;
double small_number = 1e-6;
* values[Indices::snIdx] = std::min(smax, smin + exp_value +
small_number);*
but this last line, should be imposed on values[Indices::swIdx] but
if I write it I get an error, not declared variable.
Thank you,
Lorenzo
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de <mailto:Dumux@listserv.uni-stuttgart.de>
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
// -*- 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 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);
// Set the spatial parameters
SET_TYPE_PROP(TutorialSpatialParamsImplicit, SpatialParams, TutorialSpatialParamsImplicit<TypeTag>);
// 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;
// typedef LinearMaterial<Scalar> RawMaterialLaw;
// typedef Chebyshev<Scalar> RawMaterialLaw;
// typedef BrooksCorey<Scalar> RawMaterialLaw;
// typedef ModRegularizedBrooksCorey<Scalar> RawMaterialLaw;
// typedef ModBrooksCorey<Scalar> RawMaterialLaw;
// typedef ModBrooksCorey<Scalar> type;
public:
// adapter for absolute law
//typedef EffToAbsLaw<RawMaterialLaw> type;
// select material law to be used (used absolute)
typedef ModBrooksCorey<Scalar> type;
};
}
/*!
* \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.7856);
//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 the pc-Sw and Kr curves. It is a good way to check the laws!
*/
void plotMaterialLaw()
{
PlotMaterialLaw<TypeTag> plotMaterialLaw;
GnuplotInterface<Scalar> gnuplot(plotFluidMatrixInteractions_);
gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_);
plotMaterialLaw.addpcswcurve(gnuplot, materialParams_, 0.2121, 0.7856, "Pc(Sw)", "w lp");
//plotMaterialLaw.addpcswcurve(gnuplot, materialParams_, 0.24, 0.79, "Pc(Sw)", "w lp");
//plotMaterialLaw.addpcswcurve(gnuplot, materialParams_, 0., 1., "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.2121, 0.7856, "kr(Sw)");
//plotMaterialLaw.addkrcurves(gnuplot, materialParams_, 0.22, 0.79, "kr(Sw)");
plotMaterialLaw.addkrcurves(gnuplot, materialParams_, 0., 1., "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
//SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
// Uncomment the line below to use gmsh grids (and comment the previous one)
SET_TYPE_PROP(TutorialProblemImplicit, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::simplex, Dune::conforming>);
#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)
{
private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
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, PseudoOil<Scalar> > type;
};
SET_TYPE_PROP(TutorialProblemImplicit, FluidSystem, TwoPImmiscibleFluidSystem<TypeTag>);
// Disable gravity
SET_BOOL_PROP(TutorialProblemImplicit, ProblemEnableGravity, false);
}
/*!
* \ingroup TwoPBoxModel
*
* \brief Tutorial problem for a fully coupled twophase box model.
*/
template <class TypeTag>
class TutorialProblemImplicit : public ImplicitPorousMediaProblem<TypeTag>
{
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
};
// imported from buckleyleverettproblem.hh -->
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
//typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
// <--
// 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)
// points to the plotting class for
// capillarity and permeability laws
this->spatialParams().plotMaterialLaw();
// here we can set run-time options from input file
//upperRight_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, GlobalPosition, Grid, UpperRight);
//WettingPhase::Component::setViscosity(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Fluid, ViscosityW));
//NonwettingPhase::Component::setViscosity(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Fluid, ViscosityNW));
//WettingPhase::Component::setDensity(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Fluid, DensityW));
//NonwettingPhase::Component::setDensity(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Fluid, DensityNW));
//densityNonWetting_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Fluid, DensityNW);
}
//! 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. TODO: which temperature for us?
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();
// Dirichlet condition on left boundary
if (globalPos[0] < eps_)
bcTypes.setAllDirichlet();
// Neumann condition on the other boundaries
else
bcTypes.setAllNeumann();
//bcTypes.setNeumann(Indices::pressEqIdx);
}
//! 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
{
const Scalar smax = 0.7856;
const Scalar smin = 0.2121;
const GlobalPosition &globalPos = vertex.geometry().center();
if (globalPos[0] < eps_) { // inlet, left boundary
values[Indices::pwIdx] = 30*1e5; //60.0*1e5; // pressure
values[Indices::snIdx] = (1.-smax); // oil saturation
} else { // outlet, right boundary
// It should be irrelevant which constant
// value of pressure we impose at outlet ...
values[Indices::pwIdx] = -30*1e5;
values[Indices::snIdx] = (1.-smin);
}
}
//! 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;
// define boundaries of the domain (maybe not necessary, just faster)
Scalar right = this->bBoxMax()[0]; //std::cout << right << std::endl;
Scalar left = this->bBoxMin()[0]; //std::cout << left << std::endl;
Scalar up = this->bBoxMax()[1]; //std::cout << up << std::endl;
Scalar down = this->bBoxMin()[1]; //std::cout << down << std::endl;
// everywhere imposes zero-flux
// then overwrites it
values = 0.;
// oil outflow on the right boundary
if (globalPos[0] > right - eps_) {
// oil outflux of 30 g/(m * s) on the right boundary.
//values[Indices::contiNEqIdx] = 3e-2;
// no water outflux on the right boundary.
values[Indices::contiWEqIdx] = 0;
} else if (globalPos[0] < eps_) {
const Scalar waterDensity = 1000.0;
const Scalar oilDensity = 850.0;
//values[Indices::wPhaseIdx] = 283.68; // N/m^3
values[Indices::wPhaseIdx] = (0.0034722*0.001/(2.*5.1*1e-3*15*800));
// ... other possible fields
//Indices::pressureEqIdx
//Indices::satEqIdx
//Indices::wPhaseIdx
//Indices::nPhaseIdx
//Indices::swIdx
}
else {
// no-flow on the remaining Neumann-boundaries.
// actually this is redundant because we have
// already set value = 0. previously
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
{
// Assign p=0, makes the matrix singular, and why should be vacuum?
// Since it shouldn't make any difference just chooce p!=0
// TODO: check that it makes no difference
values[Indices::pwIdx] = 30.*1e5;
//values[Indices::snIdx] = (1.0-0.7856);
//values[Indices::swIdx] = (0.7856);
// 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.+cos(3.*sqrt(std::abs(pos[1]-75.0)))), 2.));
//double exp_arg = 1. - (0.5 - pow((pos[0]-x_0)/(2.+cos(3.*sqrt(std::abs(pos[1]-75.0)))), 2.));
double exp_value = exp_arg >= -100 ? exp(exp_arg) : 0;
double small_number = 1e-6; // help the Newton solver convergence
// Warning: the initial distribution should be defined on sw and not on sn!
values[Indices::snIdx] = std::min(smax, smin + exp_value + small_number);
// with this law the Newton solver does not converge easily
//values[Indices::snIdx] = 1.0 - std::min(smax, smin + exp_value + small_number);
// values[Indices::snIdx] = std::min(smax, smin + exp_value + small_number);
//std::cout << std::min(smax, smin + exp_value + small_number) << " " <<
// 1.0 - std::min(smax, smin + exp_value + small_number) << std::endl;
// from COMSOL
// min(smax, exp(-((x+wellDistance/2-wellRadius/2)/1/(2+ cos(3*sqrt(max(-y, y)))))^2+1/2+0*1500/(1+max(y, -y)^4))+sMin)
// values[Indices::snIdx] = std::min(smax, exp(-pow(((pos[0]+150/2-0.1/2)/1/(2 + cos(3*sqrt(std::max(-pos[1], pos[1]))))),2)
// +1/2+0*1500/pow(1+std::max(pos[1],-pos[1]),4))+smin);
// values[Indices::snIdx] = 1.0 -
// std::min(smax, exp(-pow(((pos[0]+150/2-0.1/2)/1/(2 + cos(3*sqrt(std::max(-pos[1], pos[1]))))),2)
// +1/2+0*1500/pow(1+std::max(pos[1],-pos[1]),4))+smin);
}
//! 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
{
// no source term
values[Indices::contiWEqIdx] = 0.0;
values[Indices::contiNEqIdx] = 0.0;
}
private:
Scalar eps_;
GlobalPosition upperRight_;
Scalar swr_, snr_;
Scalar pLeftBc_;
Scalar densityNonWetting_;
bool paraviewOutput_;
};
}
#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 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 Chebyshev 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 ModBrooksCoreyParams
*/
template <class ScalarT, class ParamsT = ModBrooksCoreyParams<ScalarT> >
class ModBrooksCorey
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
// Capillarity-saturation law according to the Chebyshev parametrization.
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-5;
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
sw = min(max(sw, 0.2121), 0.7856);
//if (sw > smax)
// return 80000*(pow(-log(smax - smin), 2.1));
//else if (sw <= smin)
// return 80000*(pow(-log(smin - smin + eps_smin), 2.1));
return 80000*(pow(-log(sw - smin), 2.1));
}
// Saturation-capillarity law according to the Chebyshev parametrization.
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 smin + exp(-(pow(pc/(80000)), 1./2.1));
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::Sw");
}
// The partial derivative of the capillary pressure
// w.r.t. the effective saturation according to Chebyshev law.
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
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dpc_dsw");
}
// The partial derivative of the effective saturation
// w.r.t. the capillary pressure according to Chebyshev law.
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
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dswe_dpc");
}
// The relative permeability for the wetting phase of
// the medium implied by the Chebyshev parameterization
// given in absolute terms (not effective).
static Scalar krw(const Params ¶ms, Scalar sw)
{
using std::pow;
using std::min;
using std::max;
sw = min(max(sw, 0.2121), 0.7856);
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
//if (sw >= smax)
// return 5.0*pow((smax-smin), 5.);
//else if (sw <= smin)
// return 0.;
return 5.0*(sw - smin)*(sw - smin)*(sw - smin)*(sw - smin)*(sw - smin);
}
// The derivative of the relative permeability for the
// wetting phase with regard to the wetting saturation of the
// medium implied by the Chebyshev parameterization.
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 25.0*pow((swe - 0.2121), 4.);
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dkrw_dswe");
}
// The relative permeability for the non-wetting phase of
// the medium implied by the Chebyshev parameterization
// given in absolute (not effective) terms.
static Scalar krn(const Params ¶ms, Scalar sw)
{
using std::pow;
using std::min;
using std::max;
sw = min(max(sw, 0.2121), 0.7856);
const Scalar smin = 0.2121;
const Scalar smax = 0.7856;
//if (sw >= smax)
// return 0.;
//else if (sw <= smin)
// return 1.;
return 16.0*(smax - sw)*(smax - sw)*(smax - sw)*(smax - sw)*(smax - sw);
}
// The derivative of the relative permeability for the
// non-wetting phase in regard to the wetting saturation of
// the medium as implied by the Chebyshev parameterization.
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 80.0*pow((0.7856 - swe), 4.);
DUNE_THROW(Dune::NotImplemented, "mod_brookscoreyparams::dkrn_dswe");
}
};
} // 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
// -*- 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 Properties of pure water \f$H_2O\f$.
*/
#ifndef DUMUX_PSEUDOH2O_HH
#define DUMUX_PSEUDOH2O_HH
#include <dumux/material/components/component.hh>
namespace Dumux
{
/*!
* \brief Rough estimate for testing purposes of some water.
*/
template <class ScalarT>
class PseudoH2O : public Component<ScalarT, PseudoH2O<ScalarT> >
{
typedef Component<ScalarT, PseudoH2O<ScalarT> > ParentType;
typedef ScalarT Scalar;
public:
/*!
* \brief A human readable name for the water.
*/
static std::string name()
{ return "H2O"; }
/*!
* \brief Rough estimate of the density of oil [kg/m^3].
* \param temperature DOC ME!
* \param pressure DOC ME!
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
return density_;
}
/*!
* \brief Rough estimate of the viscosity of oil kg/(ms).
* \param temperature DOC ME!
* \param pressure DOC ME!
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{
return viscosity_;
};
/*!
* \brief DOC ME!
* \param viscosity DOC ME!
*/
static void setViscosity(Scalar viscosity)
{
viscosity_ = viscosity;
}
/*!
* \brief DOC ME!
* \param density DOC ME!
*/
static void setDensity(Scalar density)
{
density_ = density;
}
private:
static Scalar viscosity_;
static Scalar density_;
};
template <class ScalarT>
typename PseudoH2O<ScalarT>::Scalar PseudoH2O<ScalarT>::viscosity_ = 0.001;
template <class ScalarT>
typename PseudoH2O<ScalarT>::Scalar PseudoH2O<ScalarT>::density_ = 1000;
} // end namepace
#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 Properties of pure water \f$H_2O\f$.
*/
#ifndef DUMUX_PSEUDOOIL_HH
#define DUMUX_PSEUDOOIL_HH
#include <dumux/material/components/component.hh>
//#include "interface_BL.xml"
namespace Dumux
{
/*!
* \brief Rough estimate for testing purposes of some oil.
*/
template <class ScalarT>
class PseudoOil : public Component<ScalarT, PseudoOil<ScalarT> >
{
typedef Component<ScalarT, PseudoOil<ScalarT> > ParentType;
typedef ScalarT Scalar;
public:
/*!
* \brief A human readable name for the water.
*/
static std::string name()
{ return "Oil"; }
/*!
* \brief Rough estimate of the density of oil [kg/m^3].
* \param temperature DOC ME!
* \param pressure DOC ME!
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
return density_;
}
// double densityOil = 1000;
// double viscosityWater = interfaceFluidProps.IFP_ViscosityWettingFluid;
// double viscosityOil = interfaceFluidProps.IFP_ViscosityNonWettingFluid;
/*!
* \brief Rough estimate of the viscosity of oil kg/(ms).
* \param temperature DOC ME!
* \param pressure DOC ME!
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{
return viscosity_;
};
/*!
* \brief DOC ME!
* \param viscosity DOC ME!
*/
static void setViscosity(Scalar viscosity)
{
viscosity_ = viscosity;
}
/*!
* \brief DOC ME!
* \param density DOC ME!
*/
static void setDensity(Scalar density)
{
density_ = density;
}
private:
static Scalar viscosity_;
static Scalar density_;
};
template <class ScalarT>
typename PseudoOil<ScalarT>::Scalar PseudoOil<ScalarT>::viscosity_ = 0.104;
template <class ScalarT>
typename PseudoOil<ScalarT>::Scalar PseudoOil<ScalarT>::density_ = 850;
} // end namepace
#endif
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
