Commit 090083c8 authored by Matthias Schnebele's avatar Matthias Schnebele

co2plume should be working but the Boundary-Conditions have to be reviewed

convetivemixing compiles but the results aren't convincing
parent b6811e2e
Pipeline #483 failed with stages
in 30 seconds
......@@ -13,11 +13,12 @@ PressureLow = 1e5# [Pa] low end for tabularization of fluid properties
PressureHigh = 3e7# [Pa] high end for tabularization of fluid properties
TemperatureLow = 283.15 # [Pa] low end for tabularization of fluid properties
TemperatureHigh = 320.15 # [Pa] high end for tabularization of fluid properties
Salinity = 0.1
[SpatialParams]
Permeability = 1e-13 # [m^2] intrinsic permeability
DipAngle = 0.0 # [deg] dip angle for the domain
Porosity = 0.3 # porosity
Porosity = 0.35 # porosity
[MaterialLaw]
Swr = 0.2 # [-] residual wetting phase sat.
......
......@@ -43,8 +43,11 @@ namespace Dumux {
template <class TypeTag>
class PlumeShapeProblem;
namespace Properties {
NEW_TYPE_TAG(PlumeShapeTypeTag, INHERITS_FROM(BoxModel, TwoPTwoCCO2, PlumeShapeSpatialParamsTypeTag));
namespace Properties
{
NEW_TYPE_TAG(PlumeShapeTypeTag, INHERITS_FROM(BoxModel, TwoPTwoCCO2NI, PlumeShapeSpatialParams));
// Set grid to a two-dimensional structured quadrilateral grid
SET_TYPE_PROP(PlumeShapeTypeTag, Grid, Dune::YaspGrid<2>);
......@@ -72,16 +75,28 @@ class PlumeShapeProblem : public PorousMediumFlowProblem<TypeTag>
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
static constexpr int dimWorld = GridView::dimensionworld;
// copy some indices for convenience
static constexpr int pressureIdx = Indices::pressureIdx;
static constexpr int switchIdx = Indices::switchIdx;
static constexpr int wPhaseIdx = Indices::wPhaseIdx;
static constexpr int nPhaseIdx = Indices::nPhaseIdx;
static constexpr int wPhaseOnly = Indices::wPhaseOnly;
static constexpr int nCompIdx = FluidSystem::nCompIdx;
static constexpr int BrineIdx = FluidSystem::BrineIdx;
static constexpr int CO2Idx = FluidSystem::CO2Idx;
static constexpr int conti0EqIdx = Indices::conti0EqIdx;
static constexpr int contiCO2EqIdx = conti0EqIdx + CO2Idx;
enum {
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wPhaseOnly = Indices::wPhaseOnly,
nCompIdx = FluidSystem::nCompIdx,
BrineIdx = FluidSystem::BrineIdx,
CO2Idx = FluidSystem::CO2Idx,
#if !ISOTHERMAL
temperatureIdx = Indices::temperatureIdx,
energyEqIdx = Indices::energyEqIdx,
#endif
conti0EqIdx = Indices::conti0EqIdx,
contiCO2EqIdx = conti0EqIdx + CO2Idx
};
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
......@@ -93,7 +108,10 @@ class PlumeShapeProblem : public PorousMediumFlowProblem<TypeTag>
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CO2 = Dumux::CO2<Scalar, CO2TablesBenchmarkThree::CO2Tables>;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
//static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
static constexpr bool useMoles = true;
// the discretization method we are using
static constexpr auto discMethod = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod;
......@@ -110,6 +128,7 @@ public:
pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh");
temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow");
temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh");
salinity_ = getParam<Scalar>("FluidSystem.Salinity");
name_ = getParam<std::string>("Problem.Name");
pressure_ = getParam<Scalar>("InitialConditions.Pressure"); // hydrodynamic pressure at top layer
temperature_ = getParam<Scalar>("InitialConditions.Temperature");
......@@ -123,12 +142,13 @@ public:
gravity_[0] = -9.81 * sin(dipAngleRadians_);
// initialize the tables of the fluid system
FluidSystem::init(/*Tmin=*/temperatureLow_,
/*Tmax=*/temperatureHigh_,
/*nT=*/nTemperature_,
/*pmin=*/pressureLow_,
/*pmax=*/pressureHigh_,
/*np=*/nPressure_);
FluidSystem::init(/*Tmin*/temperatureLow_,
/*Tmax*/temperatureHigh_,
/*nT*/nTemperature_,
/*pmin*/pressureLow_,
/*pmax*/pressureHigh_,
/*np*/nPressure_,
/*salinity*/ salinity_);
}
/*!
......@@ -181,7 +201,15 @@ public:
{
values.setAllDirichlet();
}
/* TEST
// set part of West boundaries to Dirichlet & Neumann values
else if ( globalPos[1] < 30.0 + eps_ && globalPos[0] < eps_)
{
#if !ISOTHERMAL
values.setDirichlet(temperatureIdx, energyEqIdx);
#endif
}
*/
return values;
}
......@@ -194,8 +222,17 @@ public:
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
return initial_(globalPos);
{ PrimaryVariables values;
values = initial_(globalPos);
/*
#if !ISOTHERMAL
if (globalPos[1] < 30.0 + eps_ && globalPos[0] < eps_)
{
values[temperatureIdx] = injectionTemperature_; //in K
}
#endif
*/
return values;
}
/*!
......@@ -224,11 +261,13 @@ public:
{
NumEqVector fluxes(0.0);
// kg/(m^2*s) or mole/(m^2*s) depending on useMoles
const Scalar pressure = elemVolvars[scvf.insideScvIdx()].pressure(nPhaseIdx);
auto globalPos = scvf.ipGlobal();
if ( globalPos[1] < 30.0 + eps_ && globalPos[0] < eps_ )
{
Scalar massFlux = injectionRate_ /30.0; // [kg/(s*m^2)]
fluxes[contiCO2EqIdx] = -massFlux; // kg/(s*m^2)
fluxes[energyEqIdx] = -massFlux*CO2::gasEnthalpy(injectionTemperature_, pressure);
}
return fluxes;
......@@ -289,7 +328,9 @@ private:
densityW * 9.81 * distanceTopBoundaryToSurface_;
values[Indices::switchIdx] = massFracLiquidCO2;
#if !ISOTHERMAL
values[temperatureIdx] = temperature_;//290.0 + (depthBOR_ - globalPos[1])*0.03;
#endif
return values;
}
......@@ -308,6 +349,7 @@ private:
Scalar pressureLow_, pressureHigh_, pressure_;
Scalar temperatureLow_, temperatureHigh_;
Scalar salinity_;
};
} //end namespace Dumux
......
......@@ -140,6 +140,34 @@ public:
return materialParams_;
}
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param globalPos The global positio
*/
Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{ return 790; /*specific heat capacity of granite [J / (kg K)]*/ }
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param globalPos The global position
*/
Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
{ return 2700; /*density of granite [kg/m^3]*/ }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param globalPos The global position
*/
Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{ return 2.8; }
private:
Scalar permeability_;
Scalar porosity_;
......
......@@ -16,8 +16,12 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#include "config.h"
#include "convmixproblem.hh"
/*!
* \file
*
* \brief Test for the two-phase two-component CC model.
*/
#include <config.h>
#include <ctime>
#include <iostream>
......@@ -27,52 +31,33 @@
#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/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonmethod.hh>
#include <dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
// the problem definitions
#include "convmixproblem.hh"
//! 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-coupled: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.UpperRightX The x-coordinate value of the upper-right corner of the grid\n"
<< "\t-Grid.UpperRightY The y-coordinate value of the upper-right corner of the grid\n"
<< "\t-Grid.NumberOfCellsX The grid's x-resolution\n"
<< "\t-Grid.NumberOfCellsY The grid's y-resolution\n"
<< "\t-Problem.Name The name of the output files\n"
<< "\t-Problem.DepthBOR The depth of the bottom of reservoir\n"
<< "\t-Problem.Salinity The salinity of the brine phase [kg/kg]\n"
<< "\t-SpatialParams.Permeability The intrinsic permeability\n"
<< "\n";
}
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(ConvectiveMixingProblemTypeTag);
using TypeTag = TTAG(ConvectiveMixingTypeTag);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -82,18 +67,19 @@ int main(int argc, char** argv) try
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
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 = gridManager.grid().leafGridView();
const auto& leafGridView = GridCreator::grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -117,13 +103,14 @@ int main(int argc, char** argv) try
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const Scalar tEnd = getParam<Scalar>("TimeManager.TEnd");
const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
auto dt = getParam<Scalar>("TimeManager.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
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
// intialize the vtk output module
......@@ -133,8 +120,10 @@ int main(int argc, char** argv) try
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
timeLoop->setPeriodicCheckPoint(getParam<Scalar>("TimeLoop.EpisodeLength", std::numeric_limits<Scalar>::max()));
// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
......@@ -145,8 +134,11 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
using NewtonController = PriVarSwitchNewtonController<TypeTag>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
timeLoop->start(); do
......@@ -154,8 +146,26 @@ int main(int argc, char** argv) try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// 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;
......@@ -164,14 +174,17 @@ int main(int argc, char** argv) try
// 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()));
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
// write vtk output
// if episode length was specificied output only at the end of episodes
if (!haveParam("TimeLoop.EpisodeLength") || timeLoop->isCheckPoint() || timeLoop->finished() || timeLoop->timeStepIndex() == 1)
vtkWriter.write(timeLoop->time());
} while (!timeLoop->finished());
......@@ -189,7 +202,8 @@ int main(int argc, char** argv) try
}
return 0;
}
} // end main
catch (Dumux::ParameterException &e)
{
......@@ -199,10 +213,12 @@ catch (Dumux::ParameterException &e)
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;
"). 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)
......
tau = 0.706
[TimeLoop]
TEnd = 1.e9 # duration of the simulation [s]
[TimeManager]
TEnd = 1.e10 # duration of the simulation [s]
DtInitial = 10 # initial time step size [s]
EpisodeLength = 1.e7
EpisodeLength = 1.e8
[Grid]
UpperRight = 5 5 # upper right corner of the grid [m]
Cells = 10 10 # number of cells in (x, y) direction [-]
Cells = 70 70 # number of cells in (x, y) direction [-]
[FluidSystem]
NTemperature = 3# [-] number of tabularization entries
NPressure = 200# [-] number of tabularization entries
PressureLow = 1e5# [Pa] low end for tabularization of fluid properties
PressureHigh = 3e7# [Pa] high end for tabularization of fluid properties
TemperatureLow = 283.15 # [Pa] low end for tabularization of fluid properties
TemperatureHigh = 320.15 # [Pa] high end for tabularization of fluid properties
[Problem]
Name = convmix # [-] the name of the output files
DepthBOR = 1400 # [m] depth below ground surface
Salinity = 0.05 # [kg/kg] brine salinity
[MaterialLaw]
Swr = 0.2 # [-] residual wetting phase sat.
Snr = 0.2 # [-] residual non-wetting phase sat.
Pe = 5e3 # [Pa] capillary entry pressure
Lambda = 2 # [-] Brooks Corey parameter
[InitialConditions]
Temperature = 283.15 # [K] initial temperature of injected CO2
Pressure = 1.5e7 # [Pa] initial pressure
[SpatialParams]
Permeability = 3e-13 # [m^2] reservoir permeability
DipAngle = 0.0 # [deg] dip angle for the domain
Porosity = 0.2 # porosity
tau = 0.706
......@@ -23,6 +23,8 @@
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/porousmediumflow/co2/model.hh>
namespace Dumux {
......@@ -36,81 +38,122 @@ namespace Properties {
NEW_TYPE_TAG(ConvMixSpatialParamsTypeTag);
// Set the spatial parameters
SET_TYPE_PROP(ConvMixSpatialParamsTypeTag, SpatialParams, Dumux::ConvMixSpatialParamsTypeTag<TypeTag>);
} // end namespace Properties
SET_TYPE_PROP(ConvMixSpatialParams, SpatialParams, ConvMixSpatialParams<TypeTag>);
// Set the material Law
SET_PROP(ConvMixSpatialParams, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
// define the material law parameterized by absolute saturations
using type = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
};
}
/*!
* \ingroup OnePTwoCBoxModel
*
* \brief Definition of the spatial parameters for the
* convective mixing problem
*/
template<class TypeTag>
class ConvMixSpatialParams : public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry), typename GET_PROP_TYPE(TypeTag, Scalar), ConvMixSpatialParams<TypeTag>>
class ConvMixSpatialParams : public FVSpatialParams<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using ParentType = FVSpatialParams<TypeTag>;
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
enum { dimWorld = GridView::dimensionworld };
using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
public:
ConvMixSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
ConvMixSpatialParams(const Problem& problem)
: ParentType(problem)
{
// intrinsic permeability [m^2]
// intrinsic permeability
permeability_ = getParam<Scalar>("SpatialParams.Permeability");
}
~ConvMixSpatialParams()
{}
// porosity
//porosity_ = getParam<Scalar>("SpatialParams.Porosity");
// residual saturations
materialParams_.setSwr(getParam<Scalar>("MaterialLaw.Swr"));
materialParams_.setSnr(getParam<Scalar>("MaterialLaw.Snr"));
// parameters for the Brooks-Corey law
materialParams_.setPe(getParam<Scalar>("MaterialLaw.Pe"));
materialParams_.setLambda(getParam<Scalar>("MaterialLaw.Lambda"));
}
/*!
* \brief Define the intrinsic permeability \f$[m^2]\f$.
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$
* \note It is possibly solution dependent.
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return instrinsic permeability
*/
const Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
template<class ElementSolution>
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return permeability_;
}
/*!
* \brief Define the porosity \f$[-]\f$.
* \brief Returns the porosity \f$[-]\f$
*
* \param element The finite element
* \param fvElemGeom The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element
* \return porosity
*/
const Scalar porosity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
template<class ElementSolution>
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return 0.2;
}
/*!
* \brief Define the dispersivity \f$[?]\f$.
* \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
*
* \param element The finite element
* \param fvElemGeom The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* \return the material parameters object
* \param globalPos The position of the center of the element
*/
const Scalar dispersivity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
return materialParams_;
}
template<class ElementSolution>
const Scalar dispersivity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return 0.0;
}
private:
Scalar permeability_;
MaterialLawParams materialParams_;
};
} // end namespace Dumux
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment