Dear Dumux community,
I'm experiencing some difficulty with the attached simulation.
I borrowed it from 2p2c implicit examples and performed minimal
modifications which to me should not create any problem.
It seems I am wrong but I cannot figure out where!
May I ask you help to detect anything wrong, please?
Kind regards,
Lorenzo
// -*- 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 3 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
* \ingroup TwoPTwoCTests
* \brief Definition of the spatial parameters for the water-polymer-oil problem.
*/
#ifndef DUMUX_WATER_POLYMER_OIL_SPATIAL_PARAMS_HH
#define DUMUX_WATER_POLYMER_OIL_SPATIAL_PARAMS_HH
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/io/ploteffectivediffusivitymodel.hh>
#include <dumux/io/plotmateriallaw.hh>
#include <dumux/io/plotthermalconductivitymodel.hh>
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
namespace Dumux {
/*!
* \ingroup TwoPTwoCModel
* \brief Definition of the spatial parameters for the water-air problem.
*/
template<class GridGeometry, class Scalar>
class WaterPolymerOilSpatialParams
: public FVSpatialParams<GridGeometry, Scalar,
WaterPolymerOilSpatialParams<GridGeometry, Scalar>>
{
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
using ParentType = FVSpatialParams<GridGeometry, Scalar,
WaterPolymerOilSpatialParams<GridGeometry, Scalar>>;
static constexpr int dimWorld = GridView::dimensionworld;
using EffectiveLaw = RegularizedBrooksCorey<Scalar>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
//! Export the type used for the permeability
using PermeabilityType = Scalar;
//! Export the type used for the material law
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
WaterPolymerOilSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry)
{
layerBottom_ = 22.0;
// intrinsic permeabilities
fineK_ = 1e-13;
coarseK_ = 1e-12;
// porosities
finePorosity_ = 0.3;
coarsePorosity_ = 0.3;
// residual saturations
fineMaterialParams_.setSwr(0.2);
fineMaterialParams_.setSnr(0.0);
coarseMaterialParams_.setSwr(0.2);
coarseMaterialParams_.setSnr(0.0);
// parameters for the Brooks-Corey law
fineMaterialParams_.setPe(1e4);
coarseMaterialParams_.setPe(1e4);
fineMaterialParams_.setLambda(2.0);
coarseMaterialParams_.setLambda(2.0);
plotFluidMatrixInteractions_ = getParam<bool>("Output.PlotFluidMatrixInteractions");
}
/*!
* \brief This is called from the problem and creates a gnuplot output
* of e.g the pc-Sw curve
*/
void plotMaterialLaw()
{
PlotMaterialLaw<Scalar, MaterialLaw> plotMaterialLaw;
GnuplotInterface<Scalar> gnuplot(plotFluidMatrixInteractions_);
gnuplot.setOpenPlotWindow(plotFluidMatrixInteractions_);
plotMaterialLaw.addpcswcurve(gnuplot, fineMaterialParams_, 0.2, 1.0, "fine", "w lp");
plotMaterialLaw.addpcswcurve(gnuplot, coarseMaterialParams_, 0.2, 1.0, "coarse", "w l");
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, fineMaterialParams_, 0.2, 1.0, "fine");
plotMaterialLaw.addkrcurves(gnuplot, coarseMaterialParams_, 0.2, 1.0, "coarse");
gnuplot.plot("kr");
}
/*!
* \brief Applies the intrinsic permeability tensor to a pressure
* potential gradient.
*
* \param globalPos The global position
*/
Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
{
if (isFineMaterial_(globalPos))
return fineK_;
return coarseK_;
}
/*!
* \brief Defines the porosity \f$[-]\f$ of the spatial parameters
*
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
if (isFineMaterial_(globalPos))
return finePorosity_;
else
return coarsePorosity_;
}
/*!
* \brief Returns the parameter object for the Brooks-Corey material law
* which depends on the position
*
* \param globalPos The global position
*/
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
if (isFineMaterial_(globalPos))
return fineMaterialParams_;
else
return coarseMaterialParams_;
}
/*!
* \brief Function for defining which phase is to be considered as the wetting phase.
*
* \param globalPos The position of the center of the element
* \return The wetting phase index
*/
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{ return FluidSystem::H2OIdx; }
private:
bool isFineMaterial_(const GlobalPosition &globalPos) const
{ return globalPos[dimWorld-1] > layerBottom_; }
Scalar fineK_;
Scalar coarseK_;
Scalar layerBottom_;
Scalar finePorosity_;
Scalar coarsePorosity_;
MaterialLawParams fineMaterialParams_;
MaterialLawParams coarseMaterialParams_;
bool plotFluidMatrixInteractions_;
};
} // end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 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
* \ingroup TwoPTwoCTests
* \brief (Non)-isothermal gas injection problem where a gas (e.g. air)
* is injected into a fully water saturated medium.
*/
#ifndef DUMUX_WATER_POLYMER_OIL_PROBLEM_HH
#define DUMUX_WATER_POLYMER_OIL_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/method.hh>
#include <dumux/material/solidstates/inertsolidstate.hh>
#include <dumux/material/solidsystems/inertsolidphase.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/components/n2.hh>
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/h2oairxylene.hh>
#include <dumux/material/fluidsystems/h2oairmesitylene.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include <dumux/material/fluidstates/compositional.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/porousmediumflow/2p2c/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/io/gnuplotinterface.hh>
#include "spatialparams.hh"
#ifndef GRIDTYPE
#define GRIDTYPE Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
#endif
#ifndef ENABLECACHING
#define ENABLECACHING 0
#endif
namespace Dumux {
template <class TypeTag>
class WaterPolymerOilProblem;
namespace Properties {
// Create new type tags
namespace TTag {
// uncomment for non-isothermal simulations and also temperature and energy fields
//struct WaterPolymerOil { using InheritsFrom = std::tuple<TwoPTwoCNI>; };
struct WaterPolymerOil { using InheritsFrom = std::tuple<TwoPTwoC>; };
struct WaterPolymerOilBox { using InheritsFrom = std::tuple<WaterPolymerOil, BoxModel>; };
struct WaterPolymerOilCCTpfa { using InheritsFrom = std::tuple<WaterPolymerOil, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
//#if HAVE_UG
//template<class TypeTag>
//struct Grid<TypeTag, TTag::WaterPolymerOil> { using type = Dune::UGGrid<2>; };
//#else
template<class TypeTag>
//struct Grid<TypeTag, TTag::WaterPolymerOil> { using type = Dune::YaspGrid<1>; };
struct Grid<TypeTag, TTag::WaterPolymerOil> { using type = Dune::YaspGrid<2>; };
//#endif
//struct Grid<TypeTag, TTag::WaterPolymerOil> { using type = GRIDTYPE; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::WaterPolymerOil> { using type = WaterPolymerOilProblem<TypeTag>; };
// Set the wetting phase TODO: change with the appropriate fluidsystem
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::WaterPolymerOil>
{
//using type = FluidSystems::H2ON2<GetPropType<TypeTag, Properties::Scalar>>;
using type = FluidSystems::H2ON2<GetPropType<TypeTag, Properties::Scalar>,
FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/true>>;
//using type = FluidSystems::H2OAirMesitylene<GetPropType<TypeTag, Properties::Scalar>>;
//using type = FluidSystems::H2OAir<Scalar,
// Components::TabulatedComponent<Components::H2O<Scalar>>,
// FluidSystems::H2OAirDefaultPolicy</*fastButSimplifiedRelations=*/true>,
// true /*useKelvinEquation*/>;
};
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::WaterPolymerOil>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = WaterPolymerOilSpatialParams<GridGeometry, Scalar>;
};
// decide which type to use for floating values (double / quad)
template<class TypeTag>
struct Scalar<TypeTag, TTag::WaterPolymerOil>
{
using type = double;
};
template<class TypeTag>
struct Formulation<TypeTag, TTag::WaterPolymerOil>
{
public:
static const TwoPFormulation value = TwoPFormulation::p0s1;
};
// Define whether mole(true) or mass (false) fractions are used
template<class TypeTag>
struct UseMoles<TypeTag, TTag::WaterPolymerOil> { static constexpr bool value = true; };
} // end namespace Dumux
/*!
* \ingroup TwoPTwoCModel TODO: fix description
* \brief Non-isothermal gas injection problem where a gas (e.g. air)
* is injected into a fully water saturated medium.
*
* This problem uses the \ref TwoPTwoCModel and \ref NIModel model.
*
* To run the simulation execute the following line in shell:
* <tt>./test_box2p2cni</tt> or
* <tt>./test_cc2p2cni</tt>
* */
template <class TypeTag >
class WaterPolymerOilProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename ModelTraits::Indices;
// primary variable indices
enum
{
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx
//temperatureIdx = Indices::temperatureIdx,
//energyEqIdx = Indices::energyEqIdx
};
// equation indices
enum
{
contiH2OEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx, //!< Index of the mass conservation equation for the water component
contiN2EqIdx = Indices::conti0EqIdx + FluidSystem::N2Idx
};
// phase presence
enum { wPhaseOnly = Indices::firstPhaseOnly,
bothPhases = Indices::bothPhases
};
// component index
enum { N2Idx = FluidSystem::N2Idx };
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
//! Property that defines whether mole or mass fractions are used
static constexpr bool useMoles = ModelTraits::useMoles();
//! Property that defines the discretizatin method
static const bool isBox = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box;
//! Property that defines the dimension of the problem
static const int dimWorld = GridView::dimensionworld;
//static constexpr int dim = GridView::dimension;
//static constexpr int dimWorld = GridView::dimensionworld;
//static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
//enum { dofCodim = isBox ? dim : 0 };
public:
WaterPolymerOilProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
// nTemperature_ = getParam<int>("Problem.NTemperature");
// nPressure_ = getParam<int>("Problem.NPressure");
// pressureLow_ = getParam<Scalar>("Problem.PressureLow");
// pressureHigh_ = getParam<Scalar>("Problem.PressureHigh");
// temperatureLow_ = getParam<Scalar>("Problem.TemperatureLow");
// temperatureHigh_ = getParam<Scalar>("Problem.TemperatureHigh");
// temperature_ = getParam<Scalar>("Problem.InitialTemperature");
// depthBOR_ = getParam<Scalar>("Problem.DepthBOR");
name_ = getParam<std::string>("Problem.Name");
//
// // initialize the tables of the fluid system
// FluidSystem::init(/*Tmin=*/temperatureLow_,
// /*Tmax=*/temperatureHigh_,
// /*nT=*/nTemperature_,
// /*pmin=*/pressureLow_,
// /*pmax=*/pressureHigh_,
// /*np=*/nPressure_);
//
// not used
// maxDepth_ = 1000.0; // [m]
FluidSystem::init();
name_ = getParam<std::string>("Problem.Name");
// not used
//useDirichlet_ = name_.find("buoyancy") != std::string::npos;
//stating in the console whether mole or mass fractions are used
if(useMoles)
std::cout << "The problem uses mole-fractions" << std::endl;
else
std::cout << "The problem uses mass-fractions" << std::endl;
this->spatialParams().plotMaterialLaw();
// in case this BC is used ...
//InjRate_ = getParam<Scalar>("Problem.InjectionRate"); //[m/s]
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{ return name_; }
/*!
* \brief Returns the temperature within the domain [K].
*
* This problem assumes a temperature of 20 degrees Celsius.
*/
Scalar temperature() const
{ return 273.15 + 20; } // in [K]
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
// if(globalPos[0] > 40 - eps_ || globalPos[0] < eps_)
// bcTypes.setAllDirichlet();
// else
// bcTypes.setAllNeumann();
//
// if (useDirichlet_)
// {
// if (isInjectionArea_(globalPos))
// bcTypes.setAllDirichlet();
// }
if(globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
{
bcTypes.setAllDirichlet();
}
else
{
bcTypes.setAllNeumann();
}
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
*
* \param globalPos The position for which the Dirichlet condition should be evaluated
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables priVars = initial_(globalPos);
// if (useDirichlet_)
// {
// if (isInjectionArea_(globalPos))
// {
// priVars.setState(Indices::bothPhases);
// priVars[switchIdx] = 0.2;
// return priVars;
// }
// }
// condition for the polymer mole fraction at left boundary
if (globalPos[0] < eps_ )
{
//values[saturationOILIdx] = 0.25; // inject water
//values[pressureH2OIdx] = 34e6;
//values[PolymerIdx] = 0.01;
//values[pressureIdx] = 34e6;
// TODO: this may change with time to allow
// water+polymer+water slug scenario
priVars.setState(Indices::bothPhases);
priVars[switchIdx] = 0.1;
priVars[pressureIdx] = 0.3e5;
return priVars;
}
}
/*!
* \brief Evaluates the boundary conditions for a Neumann boundary segment.
*
* \param globalPos The position for which the Neumann condition should be evaluated
* \return the mole/mass flux of each component. Negative values mean influx.
*
* The units must be according to using either mole or mass fractions (mole/(m^2*s) or kg/(m^2*s)).
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
// we inject pure gasous nitrogen at the initial condition temperature and pressure from the bottom (negative values mean injection)
//if (isInjectionArea_(globalPos))
//{
// values[contiN2EqIdx] = useMoles ? -1e-3/FluidSystem::molarMass(N2Idx) : -1e-3; // kg/(m^2*s) or mole/(m^2*s)
// const auto initialValues = initial_(globalPos);
// const auto& mParams = this->spatialParams().materialLawParamsAtPos(globalPos);
// using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
// const auto pn = initialValues[pressureIdx] + MaterialLaw::endPointPc(mParams);
// //const auto t = initialValues[temperatureIdx];
// // note: energy equation is always formulated in terms of mass specific quantities, not per mole
// //values[energyEqIdx] = -1e-3*Components::N2<Scalar>::gasEnthalpy(t, pn);
//}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluates the initial value for a control volume.
*
* \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
auto priVars = initial_(globalPos);
// initially there is a heated lens in the domain
//if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] < 30 + eps_)
// priVars[temperatureIdx] = 380.0;
return priVars;
}
/*!
* \brief Evaluates the source term for all phases within a given
* sub-control volume.
*
* For this method, the \a priVars parameter stores the rate mass
* of a component is generated or annihilated per volume
* unit. Positive values mean that mass is created, negative ones
* mean that it vanishes.
*
* The units must be according to either using mole or mass fractions (mole/(m^3*s) or kg/(m^3*s)).
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
private:
// internal method for the initial condition (reused for the
// dirichlet conditions!)
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
const auto xMax = this->gridGeometry().bBoxMax()[0];
PrimaryVariables priVars(0.0);
priVars.setState(wPhaseOnly);
//priVars.setState(bothPhases);
//Scalar densityW = 1000.0;
//Scalar densityW = FluidSystem::H2O::liquidDensity(temperature_, 1e5);
//Scalar pl = 1e5 - densityW*this->spatialParams().gravity(globalPos)[1]*(depthBOR_ - globalPos[1]);
//Scalar moleFracLiquidN2 = pl*0.95/BinaryCoeff::H2O_N2::henry(temperature_);
//Scalar moleFracLiquidH2O = 1.0 - moleFracLiquidN2;
//Scalar meanM =
// FluidSystem::molarMass(H2OIdx)*moleFracLiquidH2O +
// FluidSystem::molarMass(N2Idx)*moleFracLiquidN2;
//if(useMoles)
//{
// //mole-fraction formulation
// priVars[switchIdx] = moleFracLiquidN2;
//}
//else
//{
// //mass fraction formulation
// Scalar massFracLiquidN2 = moleFracLiquidN2*FluidSystem::molarMass(N2Idx)/meanM;
// priVars[switchIdx] = massFracLiquidN2;
//}
//priVars[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
//priVars[pressureIdx] = 1e5;
priVars[switchIdx] = 0.99; // should be full of oil (?)
//priVars[temperatureIdx] = initialTemperatureProfile_(globalPos);
if (globalPos[0] > xMax - eps_ )
priVars[pressureIdx] = -0.3e5;
//if (globalPos[0] < 0.1 )
// priVars[switchIdx] = 0.01;
return priVars;
}
// not used
//Scalar initialTemperatureProfile_(const GlobalPosition &globalPos) const
//{ return 283.0 + (maxDepth_ - globalPos[1])*0.03; }
// not used
//bool isInjectionArea_(const GlobalPosition& globalPos) const
//{
// return globalPos[0] > 14.8 - eps_ && globalPos[0] < 25.2 + eps_ && globalPos[1] < eps_;
//}
Scalar maxDepth_;
static constexpr Scalar eps_ = 1e-6;
std::string name_;
bool useDirichlet_;
};
} // end namespace Dumux
#endif
[TimeLoop]
DtInitial = 1 # [s]
TEnd = 432000 # [s]
MaxTimeStepSize = 3600 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 10 10
Cells = 100 100
[Problem]
Name = waterpolymeroil
[SpatialParams]
RandomField = false
Porosity = 0.17
Permeability = 3.62E-14
[Gstat]
ControlFile = control.gstat
InputFile = gstatInput.txt
OutputFilePrefix = permeability
[Output]
PlotFluidMatrixInteractions = false
[LinearSolver]
ResidualReduction = 1e-10
[Newton]
MaxRelativeShift = 1e-13
EnableShiftCriterion = true # (relative shift convergence
criterion)
MinSteps = 2
MaxSteps = 18
TargetSteps = 10 # set the "desirable" number of
Newton iterations of a time step
RetryTimeStepReductionFactor = 5.000e-01
MaxTimeStepDivisions = 10
#EnableChop = true # chop for better convergence
[Component]
SolidDensity = 2700
SolidThermalConductivity = 2.8
SolidHeatCapacity = 790
Density = 2246 # Inert solid component density
(kg/m3)
AdsorptionGamma1 = 1000
AdsorptionGamma2 = 0
AdsorptionN = 1
adsMult = 5000
[Adaptive]
RefineAtDirichletBC = 0
RefineAtFluxBC = 1
MinLevel = 0
MaxLevel = 2
CoarsenTolerance = 1e-4
RefineTolerance = 1e-4
// -*- 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 3 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
* \ingroup TwoPTwoCTests
* \brief Test for the two-phase two-component CC model.
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/loadsolution.hh>
// the problem definitions
#include "problem.hh"
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::TYPETAG;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
if (restartTime > 0)
{
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
const auto fileName = getParam<std::string>("Restart.File");
const auto pvName = createPVNameFunction<IOFields, PrimaryVariables, ModelTraits, FluidSystem>();
loadSolution(x, fileName, pvName, *gridGeometry);
}
else
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.write(restartTime);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
//using LinearSolver = AMGBackend<TypeTag>;
//auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
// write vtk output
vtkWriter.write(timeLoop->time());
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
catch (Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
add_input_file_links(FILES params.input)
dumux_add_test(NAME test_2p2c_waterpolymeroil_box
LABELS porousmediumflow 2p2c
SOURCES main.cc
TIMEOUT 1500
COMPILE_DEFINITIONS TYPETAG=WaterPolymerOilBox ENABLECACHING=0
CMD_ARGS --script fuzzy
--command
"${CMAKE_CURRENT_BINARY_DIR}/test_2p2c_waterpolymeroil_box params.input
-Problem.Name test_2p2c_waterpolymeroil_box")
dumux_add_test(NAME test_2p2c_waterpolymeroil_tpfa
LABELS porousmediumflow 2p2c
SOURCES main.cc
TIMEOUT 1500
COMPILE_DEFINITIONS TYPETAG=WaterPolymerOilCCTpfa ENABLECACHING=0
CMD_ARGS --script fuzzy
--command
"${CMAKE_CURRENT_BINARY_DIR}/test_2p2c_waterpolymeroil_tpfa params.input
-Problem.Name test_2p2c_waterpolymeroil_tpfa")
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
