Commit fa076953 authored by Katherina Baber's avatar Katherina Baber
Browse files

more interesting and more reasonable test-problems for the 1p and 1p2c

model

-1ptest_problem is now a lens problem, the  lens can be defined in the  parameter
file, documentation is corrected

-1p2ctest is not a bio-problem anymore, but uses a new 1phase
fluidsystem for component transport of N2 in water. The problem includes
an outflow condition, the documentation is corrected (for example:
compressible fluids are also possible)

-1p2cmodel: use phaseIdx (set in 1p2cindices) instead of hard coded
index 0 for the phase, so that it is in general possible to use i.e. a
gas phase (see stokes2c as example)


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7686 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 3e6b86d2
......@@ -66,7 +66,7 @@ class OnePFluxVariables
public:
/*
* \brief The constructor
* \brief The constructor.
*
* \param problem The problem
* \param element The finite element
......@@ -87,27 +87,28 @@ public:
calculateGradients_(problem, element, elemDat);
};
/*!
* \brief The face of the current sub-control volume.
*/
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
/*!
* \brief Return the intrinsic permeability.
* \brief Return the intrinsic permeability \f$\mathrm{[m^2]}\f$.
*/
const Tensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the pressure potential gradient.
* \brief Return the pressure potential gradient \f$\mathrm{[Pa/m]}\f$.
*/
const Vector &potentialGrad() const
{ return potentialGrad_; }
/*!
* \brief Given the intrinisc permeability times the pressure
* potential gradient and SCV face normal for a phase,
* return the local index of the upstream control volume
* for a given phase.
* \brief Return the local index of the upstream control volume
* for a given phase, given the intrinisc permeability times the pressure
* potential gradient and SCV face normal for a phase.
*
* \param normalFlux The normal flux i.e. the given intrinsic permeability
* times the pressure potential gradient and SCV face normal.
......
......@@ -76,7 +76,7 @@ public:
/*!
* \brief Evaluate the rate of change of all conservation
* quantites (e.g. phase mass) within a sub control
* quantites (e.g. phase mass) within a sub-control
* volume of a finite volume element for the OneP
* model.
*
......@@ -101,7 +101,7 @@ public:
/*!
* \brief Evaluates the mass flux over a face of a subcontrol
* \brief Evaluate the mass flux over a face of a sub-control
* volume.
*
* \param flux The flux over the SCV (sub-control-volume) face
......@@ -131,7 +131,7 @@ public:
}
/*!
* \brief Calculate the source term of the equation
* \brief Calculate the source term of the equation.
*
* \param q The source/sink in the SCV
* \param localVertexIdx The index of the SCV
......
......@@ -24,7 +24,7 @@
* \file
*
* \brief Base class for all models which use the one-phase,
* box model
* box model.
* Adaption of the BOX scheme to the one-phase flow model.
*/
......@@ -50,13 +50,14 @@ namespace Dumux
* centered finite volume (box) scheme as spatial and
* the implicit Euler method as time discretization.
* Of course, the model can also be used for incompressible
* single phase flow modeling, if in the problem file a fluid with constant density is chosen.
* single phase flow modeling, if a fluid with constant density is chosen in the problem file.
*/
template<class TypeTag >
class OnePBoxModel : public BoxModel<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
......@@ -82,6 +83,7 @@ public:
// create the required scalar fields
unsigned numVertices = this->problem_().gridView().size(dim);
ScalarField *p = writer.allocateManagedBuffer(numVertices);
ScalarField *K = writer.allocateManagedBuffer(numVertices);
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
......@@ -110,12 +112,17 @@ public:
fvElemGeom,
i,
false);
const SpatialParameters &spatialParams = this->problem_().spatialParameters();
(*p)[globalIdx] = volVars.pressure();
(*K)[globalIdx] = spatialParams.intrinsicPermeability(*elemIt,
fvElemGeom,
i);
};
}
writer.attachVertexData(*p, "p");
writer.attachVertexData(*K, "K");
writer.attachCellData(*rank, "process rank");
}
};
......
......@@ -23,7 +23,7 @@
* \file
*
* \brief Base class for all problems which use the one-phase box
* model
* model.
*/
#ifndef DUMUX_1P_PROBLEM_HH
#define DUMUX_1P_PROBLEM_HH
......@@ -36,7 +36,7 @@ namespace Dumux
/*!
* \ingroup OnePBoxModel
* \ingroup BoxBaseProblems
* \brief Base class for all problems which use the single-phase box model
* \brief Base class for all problems which use the single-phase box model.
*
*/
template<class TypeTag>
......@@ -58,7 +58,7 @@ class OnePBoxProblem : public BoxProblem<TypeTag>
public:
/*!
* \brief The constructor
* \brief The constructor.
*
* \param timeManager The time manager
* \param gridView The grid view
......@@ -79,14 +79,14 @@ public:
// \{
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume.
* \brief Return the temperature \f$\mathrm{[K]}\f$ within a control volume.
*
* This is the discretization specific interface for the box
* method. By default it just calls temperature(pos).
*
* \param element The DUNE Codim<0> enitiy which intersects with
* the finite volume.
* \param fvGeom The finite volume geometry of the element.
* the finite volume
* \param fvGeom The finite volume geometry of the element
* \param scvIdx The local index of the sub control volume inside the element
*/
Scalar boxTemperature(const Element &element,
......@@ -95,18 +95,18 @@ public:
{ return asImp_().temperatureAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at a given global position.
* \brief Return the temperature \f$\mathrm{[K]}\f$ at a given global position.
*
* This is not specific to the discretization. By default it just
* calls temperature().
*
* \param pos The position in global coordinates where the temperature should be specified.
* \param pos The position in global coordinates where the temperature should be specified
*/
Scalar temperatureAtPos(const GlobalPosition &pos) const
{ return asImp_().temperature(); }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
* \brief Return the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
*
* This is not specific to the discretization. By default it just
* throws an exception so it must be overloaded by the problem if
......@@ -116,10 +116,15 @@ public:
{ DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem"); };
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
* \brief Return the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* This is the box discretization specific interface. By default
* it just calls gravityAtPos().
*
* \param element The DUNE Codim<0> enitiy which intersects with
* the finite volume
* \param fvGeom The finite volume geometry of the element
* \param scvIdx The local index of the sub control volume inside the element
*/
const GlobalPosition &boxGravity(const Element &element,
const FVElementGeometry &fvGeom,
......@@ -127,16 +132,18 @@ public:
{ return asImp_().gravityAtPos(fvGeom.subContVol[scvIdx].global); }
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
* \brief Return the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* This is discretization independent interface. By default it
* just calls gravity().
*
* \param pos The position in global coordinates
*/
const GlobalPosition &gravityAtPos(const GlobalPosition &pos) const
{ return asImp_().gravity(); }
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
* \brief Return the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* This method is used for problems where the gravitational
* acceleration does not depend on the spatial position. The
......@@ -148,7 +155,7 @@ public:
{ return gravity_; }
/*!
* \brief Returns the spatial parameters object.
* \brief Return the spatial parameters object.
*/
SpatialParameters &spatialParameters()
{ return *spatialParameters_; }
......@@ -160,7 +167,7 @@ public:
{ return *spatialParameters_; }
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
//! Return the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
......@@ -26,7 +26,7 @@
* \ingroup OnePBoxModel
* \file
*
* \brief Defines the properties required for the onephase BOX model.
* \brief Defines the properties required for the one-phase BOX model.
*/
#ifndef DUMUX_1P_PROPERTIES_DATA_HH
#define DUMUX_1P_PROPERTIES_DATA_HH
......
......@@ -39,7 +39,7 @@ namespace Dumux
/*!
* \ingroup OnePBoxModel
* \ingroup BoxVolumeVariables
* \brief Contains the quantities which are are constant within a
* \brief Contains the quantities which are constant within a
* finite volume in the one-phase model.
*/
template <class TypeTag>
......@@ -123,7 +123,7 @@ public:
}
/*!
* \brief Returns temperature inside the sub-control volume.
* \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
*
* Note that we assume thermodynamic equilibrium, i.e. the
* temperature of the rock matrix and of all fluid phases are
......@@ -133,34 +133,35 @@ public:
{ return fluidState_.temperature(); }
/*!
* \brief Returns the effective pressure of a given phase within
* \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
* the control volume.
*
*/
Scalar pressure() const
{ return fluidState_.pressure(/*phaseIdx=*/0); }
/*!
* \brief Returns the mass density of a given phase within the
* \brief Return the mass density \f$\mathrm{[kg/m13]}\f$ of a given phase within the
* control volume.
*/
Scalar density() const
{ return fluidState_.density(/*phaseIdx=*/0); }
/*!
* \brief Returns the dynamic viscosity of the fluid within the
* \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
* control volume.
*/
Scalar viscosity() const
{ return fluidState_.viscosity(/*phaseIdx=*/0); }
/*!
* \brief Returns the average porosity within the control volume.
* \brief Return the average porosity \f$\mathrm{[-]}\f within the control volume.
*/
Scalar porosity() const
{ return porosity_; }
/*!
* \brief Returns the fluid state of the control volume.
* \brief Return the fluid state of the control volume.
*/
const FluidState &fluidState() const
{ return fluidState_; }
......@@ -176,7 +177,7 @@ protected:
}
/*!
* \brief Called by update() to compute the energy related quantities
* \brief Called by update() to compute the energy related quantities.
*/
void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
......
......@@ -27,7 +27,7 @@
* \brief This file contains the data which is required to calculate
* all fluxes of fluid phases over a face of a finite volume.
*
* This means pressure and temperature gradients, phase densities at
* This means pressure and mole-/mass-fraction gradients, phase densities at
* the integration point, etc.
*
*/
......@@ -49,7 +49,7 @@ namespace Dumux
* calculate the fluxes of the fluid phases over a face of a
* finite volume for the one-phase, two-component model.
*
* This means pressure and concentration gradients, phase densities at
* This means pressure and mole-/mass-fraction gradients, phase densities at
* the intergration point, etc.
*/
template <class TypeTag>
......@@ -90,6 +90,8 @@ public:
* \param elemGeom The finite-volume geometry in the box scheme
* \param scvfIdx The local index of the SCV (sub-control-volume) face
* \param elemDat The volume variables of the current element
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
OnePTwoCFluxVariables(const Problem &problem,
const Element &element,
......@@ -121,7 +123,7 @@ public:
* vertex j.
*
* Note that the length of the face's normal is the area of the
* phase, so this is not the actual velocity by the integral of
* phase, so this is not the actual velocity but the integral of
* the velocity over the face's area. Also note that the phase
* mobility is not yet included here since this would require a
* decision on the upwinding approach (which is done in the
......@@ -150,25 +152,25 @@ public:
}
/*!
* \brief Return the intrinsic permeability tensor.
* \brief Return the intrinsic permeability tensor \f$\mathrm{[m^2]}\f$.
*/
const Tensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the dispersion tensor.
* \brief Return the dispersion tensor \f$\mathrm{[m^2/s]}\f$.
*/
const Tensor &dispersionTensor() const
{ return dispersionTensor_; }
/*!
* \brief Return the pressure potential gradient.
* \brief Return the pressure potential gradient \f$\mathrm{[Pa/m]}\f$.
*/
const Vector &potentialGrad() const
{ return potentialGrad_; }
/*!
* \brief Return the concentration gradient.
* \brief Return the concentration gradient \f$\mathrm{[mol/m^3/m]}\f$.
*
* \param compIdx The index of the considered component
*/
......@@ -183,7 +185,7 @@ public:
};
/*!
* \brief The molar concentration gradient of a component in a phase.
* \brief Return the mole-fraction gradient of a component in a phase \f$\mathrm{[mol/mol/m)]}\f$.
*
* \param compIdx The index of the considered component
*/
......@@ -198,7 +200,7 @@ public:
};
/*!
* \brief The mass fraction gradient of a component in a phase.
* \brief Return the mass-fraction gradient of a component in a phase \f$\mathrm{[kg/kg/m)]}\f$.
*
* \param compIdx The index of the considered component
*/
......@@ -213,7 +215,7 @@ public:
};
/*!
* \brief The binary diffusion coefficient for each fluid phase in the porous medium.
* \brief The binary diffusion coefficient for each fluid phase in the porous medium \f$\mathrm{[m^2/s]}\f$.
*/
Scalar porousDiffCoeff() const
{
......@@ -243,23 +245,23 @@ public:
{ return densityAtIP_; }
/*!
* \brief Given the intrinisc permeability times the pressure
* \brief Given the intrinsic permeability times the pressure
* potential gradient and SCV face normal for a phase,
* return the local index of the upstream control volume
* for a given phase.
*
* \param normalFlux The flux over a face of the sub control volume
* \param normalFlux The flux over a face of the sub-control volume
*/
int upstreamIdx(Scalar normalFlux) const
{ return (normalFlux >= 0)?face().i:face().j; }
/*!
* \brief Given the intrinisc permeability times the pressure
* \brief Given the intrinsic permeability times the pressure
* potential gradient and SCV face normal for a phase,
* return the local index of the downstream control volume
* for a given phase.
*
* \param normalFlux The flux over a face of the sub control volume
* \param normalFlux The flux over a face of the sub-control volume
*/
int downstreamIdx(Scalar normalFlux) const
{ return (normalFlux > 0)?face().j:face().i; }
......@@ -281,7 +283,7 @@ public:
protected:
/*!
* \brief Calculation of the pressure and concentration gradients
* \brief Calculation of the pressure and mole-/mass-fraction gradients.
*
* \param problem The considered problem file
* \param element The considered element of the grid
......@@ -317,17 +319,20 @@ protected:
tmp *= elemVolVars[idx].molarity(comp1Idx);
concentrationGrad_ += tmp;
// the mole-fraction gradient
tmp = feGrad;
tmp *= elemVolVars[idx].moleFraction(comp1Idx);
moleFracGrad_ += tmp;
// the mass-fraction gradient
tmp = feGrad;
tmp *= elemVolVars[idx].massFraction(comp1Idx);
massFracGrad_ += tmp;
// phase viscosity
viscosityAtIP_ += elemVolVars[idx].viscosity()*face().shapeValue[idx];
//phase moledensity
//phase molar density
molarDensityAtIP_ += elemVolVars[idx].molarDensity()*face().shapeValue[idx];
//phase density
......@@ -403,8 +408,8 @@ protected:
/*!
* \brief Calculation of the velocity normal to face using Darcy's law.
* Tensorial permeability is multiplied with the potential Gradient and the face normal.
* Identify upstream node of face
* Tensorial permeability is multiplied with the potential gradient and the face normal.
* Identify upstream node of face.
*
* \param problem The considered problem file
* \param element The considered element of the grid
......@@ -428,7 +433,7 @@ protected:
}
}
/*!
* \brief Calculation of the effective diffusion coefficient
* \brief Calculation of the effective diffusion coefficient.
*
* \param problem The considered problem file
* \param element The considered element of the grid
......@@ -450,7 +455,7 @@ protected:
}
/*!
* \brief Calculation of the dispersion
* \brief Calculation of the dispersion.
*
* \param problem The considered problem file
* \param element The considered element of the grid
......@@ -504,8 +509,9 @@ protected:
Vector potentialGrad_;
//! concentratrion gradient
Vector concentrationGrad_;
//! molar concentratrion gradient
//! mole-fraction gradient
Vector moleFracGrad_;
//! mass-fraction gradient
Vector massFracGrad_;
//! the effective diffusion coefficent in the porous medium
Scalar diffCoeffPM_;
......
......@@ -59,7 +59,6 @@ struct OnePTwoCIndices
// primary variable indices
static const int pressureIdx = PVOffset + 0; //!< pressure
static const int x1Idx = PVOffset + 1; //!< mole fraction of the second component
//in my case the therapeutic agent
};
// \}
......
......@@ -52,7 +52,7 @@ namespace Dumux
* \brief Calculate the local Jacobian for the single-phase,
* two-component model in the BOX scheme.
*
* This class is used to fill the gaps in BoxLocalResidual for the 1P-2C flow.
* This class is used to fill the gaps in BoxLocalResidual for the 1p2c flow.
*/
template<class TypeTag>
class OnePTwoCLocalResidual : public BoxLocalResidual<TypeTag>
......@@ -76,8 +76,11 @@ protected:
{
numEq = GET_PROP_VALUE(TypeTag, NumEq),
// indices of the primary variables
//phase index
phaseIdx = Indices::phaseIdx,
// indices of the primary variables
pressuerIdx = Indices::pressureIdx,
comp1Idx = Indices::comp1Idx,
// indices of the equations
......@@ -85,6 +88,7 @@ protected:
transEqIdx = Indices::transEqIdx
};
//property that defines whether mole or mass fractions are used
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
......@@ -106,7 +110,7 @@ public:
* (e.g. phase mass) within a finite volume.
*
* \param result The mass of the component within the sub-control volume
* \param scvIdx The index of the considered face of the sub control volume
* \param scvIdx The index of the considered face of the sub-control volume
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
......@@ -124,10 +128,10 @@ public:
{
// storage term of continuity equation - massfractions
result[contiEqIdx] +=
volVars.fluidState().density(/*phaseIdx*/0)*volVars.porosity();
volVars.fluidState().density(phaseIdx)*volVars.porosity();
//storage term of the transport equation - massfractions
result[transEqIdx] +=
volVars.fluidState().density(/*phaseIdx*/0) * volVars.fluidState().massFraction(/*phaseIdx*/0, comp1Idx) * volVars.porosity();
volVars.fluidState().density(phaseIdx) * volVars.fluidState().massFraction(phaseIdx, comp1Idx) * volVars.porosity();
}
else
{
......@@ -136,19 +140,20 @@ public:
result[contiEqIdx] += volVars.molarDensity()*volVars.porosity();
// storage term of the transport equation - molefractions
result[transEqIdx] +=
volVars.fluidState().molarDensity(/*phaseIdx*/0)*volVars.fluidState().moleFraction(/*phaseIdx*/0, comp1Idx) *
volVars.fluidState().molarDensity(phaseIdx)*volVars.fluidState().moleFraction(phaseIdx, comp1Idx) *
volVars.porosity();
}
}
/*!
* \brief Evaluates the mass flux over a face of a subcontrol
* \brief Evaluate the mass flux over a face of a sub-control
* volume.
*
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param faceId The index of the considered face of the sub control volume
* \param onBoundary If the considered face exists at the boundary
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
void computeFlux(PrimaryVariables &flux, int faceId, bool onBoundary=false) const
{
......@@ -165,8 +170,8 @@ public:
}
/*!
* \brief Evaluates the advective mass flux of all components over
* a face of a subcontrol volume.
* \brief Evaluate the advective mass flux of all components over
* a face of a sub-control volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current SCV
......@@ -197,9 +202,9 @@ public:
// advective flux of the second component - massfraction
flux[transEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)*up.fluidState().density(/*phaseIdx=*/0) * up.fluidState().massFraction(/*phaseIdx=*/0, comp1Idx)/up.viscosity()
(( upwindWeight_)*up.fluidState().density(phaseIdx) * up.fluidState().massFraction(phaseIdx, comp1Idx)/up.viscosity()
+
(1 - upwindWeight_)*dn.fluidState().density(/*phaseIdx=*/0)*dn.fluidState().massFraction(/*phaseIdx=*/0, comp1Idx)/dn.viscosity());
(1 - upwindWeight_)*dn.fluidState().density(phaseIdx)*dn.fluidState().massFraction(phaseIdx, comp1Idx)/dn.viscosity());
}
else
{
......@@ -214,16 +219,16 @@ public:
// advective flux of the second component -molefraction
flux[transEqIdx] +=
fluxVars.KmvpNormal() *
(( upwindWeight_)*up.molarDensity() * up.fluidState().moleFraction(/*phaseIdx=*/0, comp1Idx)/up.viscosity()
(( upwindWeight_)*up.molarDensity() * up.fluidState().moleFraction(phaseIdx, comp1Idx)/up.viscosity()
+
(1 - upwindWeight_)*dn.molarDensity() * dn.fluidState().moleFraction(/*phaseIdx=*/0, comp1Idx)/dn.viscosity());
(1 - upwindWeight_)*dn.molarDensity() * dn.fluidState().moleFraction(phaseIdx, comp1Idx)/dn.viscosity());
}
}
/*!
* \brief Adds the diffusive mass flux of all components over
* a face of a subcontrol volume.
* a face of a sub-control volume.