Commit 13ae37f7 authored by Kai Wendel's avatar Kai Wendel

[lectureEFM1p2c] updated

parent c31225e7
[TimeManager]
[TimeLoop]
MaxTimeStepSize = 5.0e6 # maximal time step size [s]
TEnd = 5.0e7 # end time of the simulation [s]
DtInitial = 1e0 # initial time step for the simulation [s]
......@@ -26,3 +26,7 @@ CoarseResidualSaturationNonWetting = 0.1 # residual saturation of the non-wet
[Boundary]
LowerPressure = 2.0e5 # Dirichlet pressure value for the boundary condition at the lower boundary [Pa]
UpperPressure = 4.0e5 # Dirichlet pressure value for the boundary condition at the upper boundary [Pa]
[Problem]
EnableGravity = false
Name = lens-1p2c
......@@ -19,7 +19,29 @@
#include "config.h"
#include "lens1p2cproblem.hh"
#include <dumux/common/start.hh>
#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/discretization/methods.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/io/vtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
......@@ -29,6 +51,10 @@
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
......@@ -58,8 +84,172 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv)
int main(int argc, char** argv) try
{
// typedef TTAG(LensOnePTwoCProblem) TypeTag;
// return Dumux::start<TypeTag>(argc, argv, usage);
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(LensOnePTwoCProblem);
// 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, usage);
// try to create a grid (from the given grid file or the input file)
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(fvGridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions"); // Q: What is this? Do we need this here?
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 = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) // Fehlt ebenso
restartTime = getParam<Scalar>("TimeLoop.Restart");
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
// 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, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
using LinearSolver = AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// try solving the non-linear system
for (int i = 0; i < maxDivisions; ++i)
{
// linearize & solve
auto converged = nonLinearSolver.solve(x);
if (converged)
break;
if (!converged && i == maxDivisions-1)
DUNE_THROW(Dune::MathError,
"Newton solver didn't converge after "
<< maxDivisions
<< " time-step divisions. dt="
<< timeLoop->timeStepSize()
<< ".\nThe solutions of the current and the previous time steps "
<< "have been saved to restart files.");
}
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the Newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} 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);
}
}
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 (...)
{
typedef TTAG(LensOnePTwoCProblem) TypeTag;
return Dumux::start<TypeTag>(argc, argv, usage);
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
This diff is collapsed.
......@@ -16,11 +16,16 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
// kaiw: file has been changed!
#ifndef DUMUX_1P2C_SPATIALPARAMS_HH
#define DUMUX_1P2C_SPATIALPARAMS_HH
#include <dumux/material/spatialparams/implicit1p.hh>
#include <dumux/porousmediumflow/1p2c/implicit/model.hh>
#include <dumux/material/spatialparams/fv1p.hh>
#include <dumux/material/fluidmatrixinteractions/porosityreactivebed.hh>
#include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
#include <dumux/material/fluidmatrixinteractions/mineralization/effectivesoliddensity.hh>
#include <dumux/material/fluidmatrixinteractions/mineralization/effectivesolidheatcapacity.hh>
/**
* @file
......@@ -34,41 +39,38 @@ namespace Dumux
/** \todo Please doc me! */
template<class TypeTag>
class Lens1p2cSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
class Lens1p2cSpatialParams : public FVSpatialParamsOneP<TypeTag>
{
typedef ImplicitSpatialParamsOneP<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename Grid::ctype CoordScalar;
// typedef typename GET_PROP_TYPE(TypeTag, OnePTwoCIndices) Indices;
using ParentType = FVSpatialParamsOneP<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using CoordScalar = typename GridView::ctype;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
//
// // indices of the primary variables
// conti0EqIdx = Indices::conti0EqIdx,
// transportEqIdx = Indices::transportEqIdx
};
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dim,dim> Tensor;
};
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP(TypeTag, ParameterTree) Params;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
using PorosityLaw = PorosityReactiveBed<TypeTag>;
using PermeabilityLaw = PermeabilityKozenyCarman<TypeTag>;
using EffectiveSolidRho = EffectiveSolidDensity<TypeTag>;
using EffectiveSolidCp = EffectiveSolidHeatCapacity<TypeTag>;
public:
// export permeability type
using PermeabilityType = Scalar;
/*!
* \brief DOC ME!
* \param gridView DOC ME!
*/
Lens1p2cSpatialParams(const GridView& gridView)
: ParentType(gridView),
lensK_(0),
outerK_(0)
Lens1p2cSpatialParams(const Problem& problem)
: ParentType(problem)
{
lensLowerLeft_[0] = 0.0;
......@@ -76,18 +78,15 @@ public:
lensUpperRight_[0]= 4.1;
lensUpperRight_[1]= 1.0;
lensPorosity_ = GET_RUNTIME_PARAM(TypeTag,double, SpatialParams.FinePorosity);
outerPorosity_ = GET_RUNTIME_PARAM(TypeTag, double, SpatialParams.CoarsePorosity);
lensPorosity_ = getParam<double>("SpatialParams.FinePorosity");
outerPorosity_ = getParam<double>("SpatialParams.CoarsePorosity");
longitudinalDispersivity_ = 1.0e-5;
transverseDispersivity_ = 1.0e-6;
//TODO: The example is very bad! The numerical diffusion is very high, so that the dispersion/diffusion coefficients nearly do not have any influence!
for(int i = 0; i < dim; i++)
{
lensK_[i][i] = GET_RUNTIME_PARAM(TypeTag, double, SpatialParams.FinePermeability);
outerK_[i][i] = GET_RUNTIME_PARAM(TypeTag, double, SpatialParams.CoarsePermeability);
}
lensK_ = getParam<double>("SpatialParams.FinePermeability");
outerK_ = getParam<double>("SpatialParams.CoarsePermeability");
}
/*!
......@@ -99,21 +98,15 @@ public:
* \param scvIdx The index sub-control volume face where the
* intrinsic velocity ought to be calculated.
*/
const Tensor &intrinsicPermeability(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
{
const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
if (isInLens_(globalPos))
return lensK_;
return outerK_;
}
Scalar porosity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
if (isInLens_(globalPos))
return lensPorosity_;
return outerPorosity_;
......@@ -125,7 +118,6 @@ public:
*/
//! Set the bounding box of the fine-sand lens
void setLensCoords(const GlobalPosition& lensLowerLeft,
const GlobalPosition& lensUpperRight)
{
......@@ -181,8 +173,8 @@ private:
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
Tensor lensK_;
Tensor outerK_;
Scalar lensK_;
Scalar outerK_;
Scalar lensPorosity_;
Scalar outerPorosity_;
Scalar longitudinalDispersivity_;
......
......@@ -124,6 +124,7 @@ public:
* \param timeManager The time manager
* \param gridView The grid view
*/
// QUESTION: Is this part to be actualized?
LensTwoPTwoCProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
......@@ -139,11 +140,11 @@ public:
/*pmax=*/3e7,
/*np=*/200);
episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
episodeLength_ = getParam<Scalar>("TimeManager.EpisodeLength"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
this->timeManager().startNextEpisode(episodeLength_);
lowerPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = getParam<Scalar>("Boundary.UpperPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
}
/*!
......@@ -206,13 +207,15 @@ public:
* \param values The boundary types for the conservation equations
* \param globalPos The position of the center of the finite volume
*/
void boundaryTypesAtPos(BoundaryTypes &values,
const GlobalPosition &globalPos) const
BoundaryTypes boundaryTypesAtPos( const GlobalPosition &globalPos) const
{
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
......@@ -224,9 +227,10 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
if (onUpperBoundary_(globalPos))
{
values[Indices::pressureIdx] = upperPressure_;
......@@ -237,6 +241,8 @@ public:
values[Indices::pressureIdx] = lowerPressure_;
values[Indices::switchIdx] = 0.0;
}
return values;
}
/*!
......@@ -253,10 +259,9 @@ public:
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables neumannAtPos( const GlobalPosition &globalPos) const
{
values = 0.0;
return PrimaryVariables(0.0);
}
// \}
......@@ -284,7 +289,8 @@ public:
{
values = 0;
}
PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
{ return PrimaryVariables(0.0); }
/*!
* \brief Evaluate the initial value for a control volume.
*
......@@ -296,9 +302,10 @@ public:
* For this method, the \a values parameter stores primary
* variables.
*/
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
const Scalar depth = this->bBoxMax()[1] - globalPos[1];
const Scalar height = this->bBoxMax()[1] - this->bBoxMin()[1];
......@@ -308,6 +315,8 @@ public:
//values[Indices::plIdx] = 1e5 - densityW*this->gravity()[1]*(depthBOR_ - globalPos[1]);
//values[Indices::plIdx] = 1.0e5;
//values[Indices::SgOrXIdx] = values[Indices::plIdx]*0.95/ BinaryCoeff::H2O_N2::henry(temperature_);
return values;
}
/*!
......
......@@ -168,12 +168,12 @@ public:
{
eps_ = 3e-6;
episodeLength_ = GET_RUNTIME_PARAM(TypeTag, Scalar, TimeManager.EpisodeLength);
episodeLength_ = getParam<Scalar>("TimeManager.EpisodeLength"); // GET_RUNTIME_PARAM(TypeTag, Scalar, TimeManager.EpisodeLength);
this->timeManager().startNextEpisode(episodeLength_);
// the boundary condition data
lowerPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = getParam<Scalar>("Boundary.UpperPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
}
/*!
......@@ -213,13 +213,15 @@ public:
* \param values The boundary types for the conservation equations
* \param globalPos The position of the center of the finite volume
*/
void boundaryTypesAtPos(BoundaryTypes &values,
const GlobalPosition &globalPos) const
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
......@@ -231,9 +233,9 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables dirichletAtPos( const GlobalPosition &globalPos) const
{
PrimaryVariables values;
FluidState fluidState;
fluidState.setTemperature(this->temperature());
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
......@@ -286,10 +288,9 @@ public:
* generated or annihilate per volume unit. Positive values mean
* that mass is created, negative ones mean that it vanishes.
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
numEqVector sourceAtPos( const GlobalPosition &globalPos) const
{
values = 0;
return numEqVector(0.0);
}
/*!
......@@ -301,9 +302,9 @@ public:
* For this method, the \a values parameter stores primary
* variables.
*/
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
FluidState fluidState;
fluidState.setTemperature(this->temperature());
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1.0e5);
......@@ -315,6 +316,8 @@ public:
values[pressureIdx] = upperPressure_ - depth/height*(upperPressure_-lowerPressure_);
values[snIdx] = (isInitial_(globalPos)) ? 0.9 :0.0;
return values;
}
// \}
......
......@@ -101,23 +101,23 @@ public:
lensLowerLeft_ = {0.0, 0.0};
lensUpperRight_= {4.0, 1.0};
lensPorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FinePorosity);
outerPorosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarsePorosity);
lensPorosity_ = getParam<Scalar> ("SpatialParams.FinePorosity"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FinePorosity);
outerPorosity_ = getParam<Scalar>("SpatialParams.CoarsePorosity");// GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarsePorosity);
lensK_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FinePermeability);
outerK_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarsePermeability);
lensK_ = getParam<Scalar>("SpatialParams.FinePermeability"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FinePermeability);
outerK_ = getParam<Scalar>("SpatialParams.CoarsePermeability"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarsePermeability);
// residual saturations
lensMaterialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FineResidualSaturationWetting));
lensMaterialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, FineResidualSaturationNonWetting));
outerMaterialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarseResidualSaturationWetting));
outerMaterialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, CoarseResidualSaturationNonWetting));
lensMaterialParams_.setSwr( getParam<Scalar>("SpatialParams.FineResidualSaturationWetting") );
lensMaterialParams_.setSnr( getParam<Scalar>("SpatialParams.FineResidualSaturationNonWetting") );
outerMaterialParams_.setSwr( getParam<Scalar>("SpatialParams.CoarseResidualSaturationWetting") );
outerMaterialParams_.setSnr( getParam<Scalar>("SpatialParams.CoarseResidualSaturationNonWetting") );
// parameters for the Brooks-Corey law
lensMaterialParams_.setPe(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,Scalar, SpatialParams, FineBrooksCoreyEntryPressure));
lensMaterialParams_.setLambda(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,Scalar, SpatialParams, FineBrooksCoreyLambda));
outerMaterialParams_.setPe(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,Scalar, SpatialParams, CoarseBrooksCoreyEntryPressure));
outerMaterialParams_.setLambda(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,Scalar, SpatialParams, CoarseBrooksCoreyLambda));
lensMaterialParams_.setPe( getParam<Scalar>("SpatialParams.FineBrooksCoreyEntryPressure") );
lensMaterialParams_.setLambda( getParam<Scalar>("SpatialParams.FineBrooksCoreyLambda") );
outerMaterialParams_.setPe( getParam<Scalar>("SpatialParams.CoarseBrooksCoreyEntryPressure") );
outerMaterialParams_.setLambda( getParam<Scalar>("SpatialParams,CoarseBrooksCoreyLambda") );
}
~Lens2pSpatialParams()
......
......@@ -127,6 +127,9 @@ class LensOnePTwoCProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
public:
/*!
* \brief The constructor
......@@ -142,22 +145,22 @@ public:
//This overwrites the lens settings made in the spatialparameter file
this->spatialParams().setLensCoords({0.8, 2.0}, {4.0, 3.0});
infiltrationRate_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, InfiltrationRate);
infiltrationRate_ = getParam<Scalar>("Boundary.InfiltrationRate"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, InfiltrationRate);
infiltrationStartTime_= 1.0e-9; //The infiltrations starts always after the first time step!
infiltrationEndTime_= GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, InfiltrationEndTime);
infiltrationEndTime_= getParam<Scalar>("Boundary.InfiltrationEndTime"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, InfiltrationEndTime);
episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
episodeLength_ = getParam<Scalar>("TimeManager.EpisodeLength"); // (TypeTag, Scalar, TimeManager, EpisodeLength);
this->timeManager().startNextEpisode(episodeLength_);
// the boundary condition data
lowerPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure"); //GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = getParam<Scalar>("Boundary.UpperPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
paraviewOutput_ = ParameterTree::tree().template get<bool>("Output.paraviewOutput", true);
if (!paraviewOutput_)
{
// the number of cells in x and y direction
numCells_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::vector<int>, Grid, Cells);
numCells_ = getParam<std::vector<int> >("Grid.Cells"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::vector<int>, Grid, Cells);
// Write header for output file
std::ofstream dataFile;
......@@ -215,8 +218,9 @@ public:
* \param globalPos DOC ME!
*
*/
void boundaryTypesAtPos(BoundaryTypes &values, const GlobalPosition& globalPos) const
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
......@@ -224,6 +228,8 @@ public:
if (onInlet_(globalPos))
values.setNeumann(transportEqIdx);
return BoundaryTypes;
}
/*!
......@@ -235,8 +241,10 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
PrimaryVariables dirichletAtPos( const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
if (onUpperBoundary_(globalPos))
{
values[pressureIdx] = upperPressure_;
......@@ -247,7 +255,9 @@ public:
values[pressureIdx] = lowerPressure_;
values[massOrMoleFracIdx] = 0.0;
}
else {values = 0.0;}
// else {values = 0.0;}
return values;
}
/*!
......@@ -259,10 +269,10 @@ public:
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
values = 0.0;
PrimaryVaribles values(0.0);
// TODO: Apply the new kind of time Management implementation. Or is it maybe already done????
const Scalar time = this->timeManager().time();
if (time <= infiltrationEndTime_ && infiltrationStartTime_ <= time && onInlet_(globalPos))
values[transportEqIdx] = -infiltrationRate_;
......@@ -282,10 +292,14 @@ public:
* generated or annihilate per volume unit. Positive values mean
* that mass is created, negative ones mean that it vanishes.
*/
/*
void sourceAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
values = Scalar(0.0);
}
}*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
/*!
* \brief Evaluate the initial value for a control volume.
......@@ -296,14 +310,17 @@ public:
* For this method, the \a values parameter stores primary