diff --git a/dumux/boxmodels/richards/richardsfluidstate.hh b/dumux/boxmodels/richards/richardsfluidstate.hh index a69d6b45e195ba6177c5d70729c65cd7bf834fd5..b4e6a4bef49f73e98da6cd813895e413d5680749 100644 --- a/dumux/boxmodels/richards/richardsfluidstate.hh +++ b/dumux/boxmodels/richards/richardsfluidstate.hh @@ -52,14 +52,23 @@ class RichardsFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTA }; public: + //! The number of fluid phases used by the fluid state enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) }; - void update(const Problem &problem, + /*! + * \brief Updates the fluid quantities from the primary variables + * of the Richards model. + * + * \param pnRef The reference pressure of the non-wetting fluid phase [Pa] + * \param matParams The parameters for the capillary pressure/relative permeability material law + * \param priVars The primary variables for which the fluid state ought to be calculated + * \param temperature The temperature which should be used + */ + void update(Scalar pnRef, const MaterialLawParams &matParams, const PrimaryVariables &priVars, Scalar temperature) { - Scalar pnRef = problem.pNreference(); Scalar minPc = MaterialLaw::pC(matParams, 1.0); phasePressure_[wPhaseIdx] = priVars[pwIdx]; phasePressure_[nPhaseIdx] = std::max(pnRef, priVars[pwIdx] + minPc); @@ -77,7 +86,9 @@ public: } /*! - * \brief Returns the saturation of a phase. + * \brief Returns the saturation [] of a phase. + * + * \param phaseIdx The index of the fluid phase */ Scalar saturation(int phaseIdx) const { @@ -88,7 +99,10 @@ public: }; /*! - * \brief Returns the mass fraction of a component in a phase. + * \brief Returns the mass fraction [] of a component in a phase. + * + * \param phaseIdx The index of the fluid phase + * \param compIdx The index of the chemical species */ Scalar massFrac(int phaseIdx, int compIdx) const { @@ -98,7 +112,10 @@ public: } /*! - * \brief Returns the molar fraction of a component in a fluid phase. + * \brief Returns the molar fraction [] of a component in a fluid phase. + * + * \param phaseIdx The index of the fluid phase + * \param compIdx The index of the chemical species */ Scalar moleFrac(int phaseIdx, int compIdx) const { @@ -109,6 +126,8 @@ public: * \brief Returns the total concentration of a phase [mol / m^3]. * * This is equivalent to the sum of all component concentrations. + * + * \param phaseIdx The index of the fluid phase */ Scalar phaseConcentration(int phaseIdx) const { @@ -117,6 +136,9 @@ public: /*! * \brief Returns the concentration of a component in a phase [mol / m^3]. + * + * \param phaseIdx The index of the fluid phase + * \param compIdx The index of the chemical species */ Scalar concentration(int phaseIdx, int compIdx) const { @@ -127,6 +149,8 @@ public: /*! * \brief Returns the density of a phase [kg / m^3]. + * + * \param phaseIdx The index of the fluid phase */ Scalar density(int phaseIdx) const { return density_[phaseIdx]; } @@ -136,12 +160,16 @@ public: * * This is equivalent to the sum of all component molar masses * weighted by their respective mole fraction. + * + * \param phaseIdx The index of the fluid phase */ Scalar averageMolarMass(int phaseIdx) const { return FluidSystem::molarMass(phaseIdx); }; /*! * \brief Returns the partial pressure of a component in the gas phase [Pa]. + * + * \param compIdx The index of the chemical species */ Scalar partialPressure(int compIdx) const { @@ -152,6 +180,8 @@ public: /*! * \brief Returns the pressure of a fluid phase [Pa]. + * + * \param phaseIdx The index of the fluid phase */ Scalar phasePressure(int phaseIdx) const { return phasePressure_[phaseIdx]; } diff --git a/dumux/boxmodels/richards/richardsfluxvariables.hh b/dumux/boxmodels/richards/richardsfluxvariables.hh index 1970551bd13cc5ba8cebb392f207a447364a28f8..7917be9b21cd44e48a24965c7e59944a34eb4f84 100644 --- a/dumux/boxmodels/richards/richardsfluxvariables.hh +++ b/dumux/boxmodels/richards/richardsfluxvariables.hh @@ -17,8 +17,8 @@ /*! * \file * - * \brief This file contains the data which is required to calculate - * the flux of fluid over a face of a finite volume. + * \brief Data which is required to calculate the flux of fluid over a + * face of a finite volume */ #ifndef DUMUX_RICHARDS_FLUX_VARIABLES_HH #define DUMUX_RICHARDS_FLUX_VARIABLES_HH @@ -30,7 +30,7 @@ namespace Dumux /*! * \ingroup RichardsModel - * \brief This template class contains the data which is required to + * \brief Calculates and stores the data which is required to * calculate the flux of fluid over a face of a finite volume. */ template <class TypeTag> @@ -61,27 +61,39 @@ class RichardsFluxVariables typedef typename GridView::template Codim<0>::Entity Element; public: + /*! + * \brief Constructor + * + * \param problem The representation of the physical problem + * \param element The DUNE Codim<0> entity which contains the face of + * the finite volume + * \param fvElemGeom The finite volume geometry of the element + * \param scvfIdx The local index of the sub-control volume face in the + * element's finite volume geometry. + * \param elemVolVars An array containing the volume variables for all + * sub-control volumes of the element. + */ RichardsFluxVariables(const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int faceIdx, + const FVElementGeometry &fvElemGeom, + int scfvIdx, const ElementVolumeVariables &elemVolVars) - : fvElemGeom_(elemGeom) + : fvElemGeom_(fvElemGeom) { - scvfIdx_ = faceIdx; + scvfIdx_ = scfvIdx; calculateGradients_(problem, element, elemVolVars); calculateK_(problem, element, elemVolVars); }; /* - * \brief Return the intrinsic permeability. + * \brief Return the intrinsic permeability [m^2]. */ const Tensor &intrinsicPermeability() const { return K_; } /*! - * \brief Return the pressure potential gradient. + * \brief Return the pressure potential gradient [Pa/m] */ const Vector &potentialGradW() const { return potentialGrad_; } @@ -91,6 +103,9 @@ public: * 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 mass flux over the face multiplied with + * the face's normal. */ int downstreamIdx(Scalar normalFlux) const { return (normalFlux >= 0)?face().j:face().i; } @@ -100,10 +115,17 @@ public: * 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 mass flux over the face multiplied with + * the face's normal. */ int upstreamIdx(Scalar normalFlux) const { return (normalFlux > 0)?face().i:face().j; } + /*! + * \brief Returns the face of the element's finite volume geometry + * which the flux variables object looks at + */ const SCVFace &face() const { return fvElemGeom_.subContVolFace[scvfIdx_]; } diff --git a/dumux/boxmodels/richards/richardsindices.hh b/dumux/boxmodels/richards/richardsindices.hh index f022f624ade0780a852ba73441ee4700cf4ec93c..363c2bfca62fe28cb0fd30a0fd9973a96374389a 100644 --- a/dumux/boxmodels/richards/richardsindices.hh +++ b/dumux/boxmodels/richards/richardsindices.hh @@ -16,7 +16,7 @@ /*! * \file * - * \brief Indices for the Richards model. + * \brief Index names for the Richards model. */ #ifndef DUMUX_RICHARDS_INDICES_HH #define DUMUX_RICHARDS_INDICES_HH @@ -29,7 +29,7 @@ namespace Dumux // \{ /*! - * \brief Indices for the Richards model. + * \brief Index names for the Richards model. */ struct RichardsIndices { @@ -37,19 +37,21 @@ struct RichardsIndices // primary variable indices ////////// - //! wetting phase pressure + //! Primary variable index for the wetting phase pressure static const int pwIdx = 0; ////////// // equation indices ////////// - //! continuity + //! Equation index for the mass conservation of the wetting phase static const int contiEqIdx = 0; ////////// // phase indices ////////// + //! Phase index for the wetting phase static const int wPhaseIdx = 0; + //! Phase index for the wetting phase static const int nPhaseIdx = 1; }; // \} diff --git a/dumux/boxmodels/richards/richardslocalresidual.hh b/dumux/boxmodels/richards/richardslocalresidual.hh index 5f8f80d44f81ece2f4016ae028dff771132f2cd7..d4c84f7dbc37a931eeb2a66ef27c83af9578813d 100644 --- a/dumux/boxmodels/richards/richardslocalresidual.hh +++ b/dumux/boxmodels/richards/richardslocalresidual.hh @@ -37,8 +37,7 @@ namespace Dumux { /*! * \ingroup RichardsModel - * \brief Element-wise calculation of the Jacobian matrix for problems - * using the Richards box model. + * \brief Element-wise calculation of the residual for the Richards box model. */ template<class TypeTag> class RichardsLocalResidual : public BoxLocalResidual<TypeTag> @@ -70,6 +69,11 @@ public: * model. * * This function should not include the source and sink terms. + * + * \param result Stores the average mass per unit volume for each phase [kg/m^3] + * \param scvIdx The sub control volume index of the current element + * \param usePrevSol Calculate the storage term of the previous solution + * instead of the model's current solution. */ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const { @@ -94,13 +98,19 @@ public: /*! * \brief Evaluates the mass flux over a face of a subcontrol * volume. + * + * + * \param flux Stores the total mass fluxes over a sub-control volume face + * of the current element [kg/s] + * \param scvfIdx The sub control volume face index inside the current + * element */ - void computeFlux(PrimaryVariables &flux, int faceId) const + void computeFlux(PrimaryVariables &flux, int scvfIdx) const { FluxVariables fluxVars(this->problem_(), this->elem_(), this->fvElemGeom_(), - faceId, + scvfIdx, this->curVolVars_()); // calculate the flux in the normal direction of the @@ -128,13 +138,18 @@ public: /*! * \brief Calculate the source term of the equation + * + * \param flux Stores the average source term of all phases inside a + * sub-control volume of the current element [kg/(m^3 s)] + * \param scvIdx The sub control volume index inside the current + * element */ - void computeSource(PrimaryVariables &q, int localVertexIdx) + void computeSource(PrimaryVariables &q, int scvIdx) { this->problem_().source(q, this->elem_(), this->fvElemGeom_(), - localVertexIdx); + scvIdx); } }; diff --git a/dumux/boxmodels/richards/richardsmodel.hh b/dumux/boxmodels/richards/richardsmodel.hh index 39175522b7a399c35f456e868075da08cf7b06e9..32e6b1ab288156efa759f19a7507b1d154fdc800 100644 --- a/dumux/boxmodels/richards/richardsmodel.hh +++ b/dumux/boxmodels/richards/richardsmodel.hh @@ -36,11 +36,6 @@ namespace Dumux { -/*! - * \ingroup BoxProblems - * \defgroup RichardsBoxProblems Richards box problems - */ - /*! * \ingroup BoxModels * \defgroup RichardsModel Richards box model @@ -48,26 +43,43 @@ namespace Dumux /*! * \ingroup RichardsModel - * \brief Adaption of the BOX scheme to the isothermal Richards model. - * - * - * In the unsaturated zone, Richards' equation can be used. - * Gas has resistance against the water flow in porous media. - * However, viscosity of air is about 1\% of the viscosity of water, - * which makes it highly mobile compared to the water phase. - * Therefore, in Richards` equation only water phase with capillary effects are considered, - * where pressure of the gas phase is set to a reference pressure (\f${p_n}_{ref}\f$). - * - * \f{align*} - * \varrho \hspace{1mm} \phi \hspace{1mm} \frac{\partial S_w}{\partial p_c} \frac{\partial p_c}{\partial t} - \nabla \cdot (\frac{kr_w}{\mu_w} \hspace{1mm} \varrho_w \hspace{1mm} K \hspace{1mm} - * (\nabla p_w - \varrho_w \hspace{1mm} \vec{g})) \hspace{1mm} = \hspace{1mm} q, - * \f} - * where \f$p_w = {p_n}_{ref} - p_c\f$. - * Here \f$ p_w \f$, \f$ p_c \f$, and \f$ {p_n}_{ref} \f$ - * denote water pressure, capillary pressure, and non-wetting phase reference pressure, repectively. - * - * To overcome convergence problems, \f$ \frac{\partial S_w}{\partial p_c} \f$ is taken from the old iteration step. + * \brief Implements the Richards model for quasi twophase flow. * + * In the unsaturated zone, Richards' equation can be used. Fundamentally, the Richards-equation + * is equivalent to the twophase model, i.e. + \f[ + \frac{\partial\;\phi S_\alpha \rho_\alpha}{\partial t} + - + \mathbf{div} \left\{ + \frac{k_{r\alpha}}{\mu_\alpha}\;K + \mathbf{grad}\left[ + p_alpha - g\rho_\alpha + \right] + \right\} + = + q_\alpha, + \f] + * where \f$\alpha \in \{w, n\}\f$ is the fluid phase, \f$\rho_\alpha\f$ is the fluid + * density, \f$S_\alpha\f$ is the fluid saturation, \f$\phi\f$ is the porosity, + * \f$k_{r\alpha}\f$ is the relative permeability of the fluid, \f$\mu_\alpha\f$ is + * the fluid's dynamic viscosity, \f$K\f$ is the intrinsic permeability, \f$p_\alpha\f$ + * is the fluid pressure and \f$g\f$ is the potential of the gravity. + * + * However, the Richards model assumes that the non-wetting is gas + * which typically exhibits a much low viscosity than the liquid + * wetting phase. (For example air has about \f$1\%\f$ of the viscosity of + * liquid water.) As a consequence, the + * \f$\frac{k_{r\alpha}}{\mu_\alpha}\f$ term is typically much larger for + * the non-wetting phase than for the wetting phase. In the Richards + * model it is now assumed that \f$\frac{k_{rn}}{\mu_n} \to \infty\f$ + * which means that the pressure of the non-wetting phase is + * equivalent to hydrostatic pressure of the gas or can be externally + * specified. Therefore, in Richards' equation mass conservation only + * needs to be considered for the wetting phase. The model thus choses + * \f$p_w\f$ as its only primary variable and calculates the wetting phase + * saturation using the inverse of the capilary pressure, i.e. + \f[ S_w = p_c^{-1}(p_n - p_w)\, \f] + * where \f$p_n\f$ is an externally given reference pressure. */ template<class TypeTag > class RichardsModel : public BoxModel<TypeTag> @@ -100,6 +112,9 @@ public: /*! * \brief Returns the relative weight of a primary variable for * calculating relative errors. + * + * \param vertIdx The global index of the vertex in question + * \param pvIdx The index of the primary variable */ Scalar primaryVarWeight(int vertIdx, int pvIdx) const { @@ -111,6 +126,9 @@ public: /*! * \brief All relevant primary and secondary of a given * solution to an ouput writer. + * + * \param sol The current solution which ought to be written to disk + * \param writer The Dumux::VtkMultiWriter which is be used to write the data */ template <class MultiWriter> void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) diff --git a/dumux/boxmodels/richards/richardsnewtoncontroller.hh b/dumux/boxmodels/richards/richardsnewtoncontroller.hh index b69ca41c856bcdcacac0888da8ba0a6dcfec6924..f3f6f2a1882a33d52e05e797ebd2127858e08103 100644 --- a/dumux/boxmodels/richards/richardsnewtoncontroller.hh +++ b/dumux/boxmodels/richards/richardsnewtoncontroller.hh @@ -26,11 +26,10 @@ namespace Dumux { /*! - * \brief A 2pNc specific controller for the newton solver. + * \brief A Richards model specific controller for the newton solver. * * This controller 'knows' what a 'physically meaningful' solution is - * which allows the newton method to abort quicker if the solution is - * way out of bounds. + * and can thus do update smarter than the plain Newton controller. */ template <class TypeTag> class RichardsNewtonController : public NewtonController<TypeTag> @@ -59,9 +58,25 @@ class RichardsNewtonController : public NewtonController<TypeTag> typedef Dune::FieldVector<Scalar, dim> GlobalPosition; public: + /*! + * \brief Constructor + */ RichardsNewtonController() { }; + /*! + * \brief Update the current solution of the newton method + * + * This is basically the step + * \f[ u^{k+1} = u^k - \Delta u^k \f] + * + * \param deltaU When the method is called, this contains the + * vector of differences between the current + * iterative solution and the next \f$\Delta + * u\f$. After the method is finished it should + * contain the next iterative solution \f$ u^{k+1} \f$. + * \param uOld The current iterative solution \f$ u^k \f$ + */ void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld) { this->writeConvergence_(uOld, deltaU); @@ -89,7 +104,8 @@ public: const MaterialLawParams &mp = sp.materialLawParams(*eIt, fvElemGeom, i); Scalar pcMin = MaterialLaw::pC(mp, 1.0); Scalar pW = uOld[globI][pwIdx]; - Scalar pN = std::max(this->problem_().pNreference(), pW + pcMin); + Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvElemGeom, i), + pW + pcMin); Scalar pcOld = pN - pW; Scalar SwOld = std::max(0.0, MaterialLaw::Sw(mp, pcOld)); diff --git a/dumux/boxmodels/richards/richardsproblem.hh b/dumux/boxmodels/richards/richardsproblem.hh index b89d757373ef9e6343f61d3615ec7d275acedf0a..731cc7d1db74685ffda9c093183912728d69ca9b 100644 --- a/dumux/boxmodels/richards/richardsproblem.hh +++ b/dumux/boxmodels/richards/richardsproblem.hh @@ -28,7 +28,7 @@ namespace Dumux * \ingroup RichardsProblems * \brief Base class for all problems which use the two-phase box model * - * \todo Please doc me more! + * For a description of the Richards model, see Dumux::RichardsModel */ template<class TypeTag> class RichardsBoxProblem : public BoxProblem<TypeTag> @@ -38,6 +38,7 @@ class RichardsBoxProblem : public BoxProblem<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry; // material properties typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters; @@ -46,10 +47,25 @@ class RichardsBoxProblem : public BoxProblem<TypeTag> dim = GridView::dimension, dimWorld = GridView::dimensionworld }; + + typedef typename GridView::template Codim<0>::Entity Element; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; public: + /*! + * \brief The constructor. + * + * The overloaded class must allocate all data structures + * required, but _must not_ do any calls to the model, the + * jacobian assembler, etc inside the constructor. + * + * If the problem requires information from these, the + * BoxProblem::init() method be overloaded. + * + * \param timeManager The TimeManager which keeps track of time + * \param gridView The GridView used by the problem. + */ RichardsBoxProblem(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView), gravity_(0), spatialParams_(gridView) @@ -65,23 +81,38 @@ public: // \{ /*! - * \brief Returns the temperature within the domain. + * \brief Returns the temperature [K] within a control volume. * * This method MUST be overwritten by the actual problem. + * + * \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 */ - Scalar temperature() const - { return asImp_()->temperature(); }; + Scalar temperature(const Element &element, + const FVElementGeometry fvGeom, + int scvIdx) const + { DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the actual problem"); }; /*! - * \brief Returns the reference pressure of the non-wetting phase within the domain. + * \brief Returns the reference pressure [Pa] of the non-wetting + * phase within a control volume. * * This method MUST be overwritten by the actual problem. + * + * \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 */ - Scalar pNreference() const - { return asImp_()->pNreference(); }; + Scalar referencePressure(const Element &element, + const FVElementGeometry fvGeom, + int scvIdx) const + { DUNE_THROW(Dune::NotImplemented, "referencePressure() method not implemented by the actual problem"); }; /*! - * \brief Returns the acceleration due to gravity. + * \brief Returns the acceleration due to gravity [m/s^2]. * * If the <tt>EnableGravity</tt> property is true, this means * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$ @@ -104,14 +135,7 @@ public: // \} private: - //! Returns the implementation of the problem (i.e. static polymorphism) - Implementation *asImp_() - { return static_cast<Implementation *>(this); } - - //! \copydoc asImp_() - const Implementation *asImp_() const - { return static_cast<const Implementation *>(this); } - + // the gravity vector GlobalPosition gravity_; // material properties diff --git a/dumux/boxmodels/richards/richardsproperties.hh b/dumux/boxmodels/richards/richardsproperties.hh index f3f5a3e47573920499b0b0a2cbb982aaf12d12ec..52f64170aa0d8e9cb103b5ac1ccc123d0c30009e 100644 --- a/dumux/boxmodels/richards/richardsproperties.hh +++ b/dumux/boxmodels/richards/richardsproperties.hh @@ -16,7 +16,8 @@ /*! * \file * - * \brief Contains the properties for the Richards BOX model. + * \brief Contains the property declarations for the Richards box + * model. */ #ifndef DUMUX_RICHARDS_PROPERTIES_HH #define DUMUX_RICHARDS_PROPERTIES_HH @@ -38,8 +39,7 @@ namespace Properties { // Type tags ////////////////////////////////////////////////////////////////// -//! The type tag for problems discretized using the isothermal -//! richards model +//! The type tag for problems discretized using the Richards model NEW_TYPE_TAG(BoxRichards, INHERITS_FROM(BoxModel)); ////////////////////////////////////////////////////////////////// @@ -47,12 +47,12 @@ NEW_TYPE_TAG(BoxRichards, INHERITS_FROM(BoxModel)); ////////////////////////////////////////////////////////////////// NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system -NEW_PROP_TAG(RichardsIndices); //!< Enumerations for the richards models +NEW_PROP_TAG(RichardsIndices); //!< Enumerations used by the Richards models NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object -NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters) -NEW_PROP_TAG(MaterialLawParams); //!< The context material law (extracted from the spatial parameters) -NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the richards model -NEW_PROP_TAG(FluidState); //!< The fluid state to be used for the richards model +NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (by default extracted from the spatial parameters) +NEW_PROP_TAG(MaterialLawParams); //!< The context material law (by default extracted from the spatial parameters) +NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the Richards model +NEW_PROP_TAG(FluidState); //!< The fluid state to be used for the Richards model NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility // \} diff --git a/dumux/boxmodels/richards/richardspropertydefaults.hh b/dumux/boxmodels/richards/richardspropertydefaults.hh index 9a725421a6fd847dbc331fac71235f635221a037..01bdabf8f8558bd9ab23af81d98e6fb338d09af8 100644 --- a/dumux/boxmodels/richards/richardspropertydefaults.hh +++ b/dumux/boxmodels/richards/richardspropertydefaults.hh @@ -16,8 +16,8 @@ /*! * \file * - * \brief Contains the default definitions for the properties for the - * Richards BOX model. + * \brief Contains the default definitions for the properties required + * by the Richards box model. */ #ifndef DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH #define DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH @@ -45,37 +45,41 @@ namespace Properties { ////////////////////////////////////////////////////////////////// // Properties values ////////////////////////////////////////////////////////////////// +//! Number of equations required by the model SET_INT_PROP(BoxRichards, NumEq, 1); +//! Number of fluid phases considered SET_INT_PROP(BoxRichards, NumPhases, 2); -//! Use the 2p local jacobian operator for the 2p model +//! The local residual operator SET_TYPE_PROP(BoxRichards, LocalResidual, RichardsLocalResidual<TypeTag>); -//! the Model property +//! The global model used SET_TYPE_PROP(BoxRichards, Model, RichardsModel<TypeTag>); -//! the VolumeVariables property +//! The class for the volume averaged quantities SET_TYPE_PROP(BoxRichards, VolumeVariables, RichardsVolumeVariables<TypeTag>); -//! the FluxVariables property +//! The class for the quantities required for the flux calculation SET_TYPE_PROP(BoxRichards, FluxVariables, RichardsFluxVariables<TypeTag>); -//! the NewtonController property +//! The class of the newton controller SET_TYPE_PROP(BoxRichards, NewtonController, RichardsNewtonController<TypeTag>); -//! the weight of the upwind vertex for the mobility +//! The weight of the upwind vertex when looking at the mobility SET_SCALAR_PROP(BoxRichards, MobilityUpwindAlpha, 1.0); -//! The indices required by the isothermal single-phase model +//! The class with all index definitions for the model SET_TYPE_PROP(BoxRichards, RichardsIndices, Dumux::RichardsIndices); /*! - * \brief Set the property for the material law by retrieving it from - * the spatial parameters. + * \brief The material law for capillary pressure and relative permeability + * + * By default this type is determined by retrieving it from the + * spatial parameters. */ SET_PROP(BoxRichards, MaterialLaw) { @@ -87,8 +91,9 @@ public: }; /*! - * \brief Set the property for the material parameters by extracting - * it from the material law. + * \brief Set type of the parameter objects for the material law + * + * By default this is just retrieved from the material law. */ SET_PROP(BoxRichards, MaterialLawParams) { @@ -99,13 +104,33 @@ public: typedef typename MaterialLaw::Params type; }; +/*! + *\brief The fluid system used by the model. + * + * By default this uses the immiscible twophase fluid system. The + * actual fluids used are specified using in the problem definition by + * the WettingPhase and NonwettingPhase properties. Be aware that + * using different fluid systems in conjunction with the Richards + * model only makes very limited sense. + */ SET_TYPE_PROP(BoxRichards, FluidSystem, FluidSystem2P<TypeTag>); + +/*! + * \brief The non-wetting phase used. + * + * By default we use gaseous nitrogen as non-wetting phase. Please be + * aware that you should be careful to use the Richards model in + * conjunction with liquid non-wetting phases. This is only meaningful + * if the viscosity of the liquid phase is _much_ lower than the + * viscosity of the wetting phase. + */ SET_PROP(BoxRichards, NonwettingPhase) { typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; typedef GasPhase<Scalar, N2<Scalar>> type; }; +//! The fluid state class SET_TYPE_PROP(BoxRichards, FluidState, RichardsFluidState<TypeTag>); // \} diff --git a/dumux/boxmodels/richards/richardsvolumevariables.hh b/dumux/boxmodels/richards/richardsvolumevariables.hh index 74094980c2daa9d4147e06f9952dc72543079a24..50c925feff7f9fb8344894dc4aa3246a844f459c 100644 --- a/dumux/boxmodels/richards/richardsvolumevariables.hh +++ b/dumux/boxmodels/richards/richardsvolumevariables.hh @@ -17,7 +17,7 @@ /*! * \file * - * \brief Quantities required by the richards box model defined on a vertex. + * \brief Volume averaged quantities required by the Richards model. */ #ifndef DUMUX_RICHARDS_VOLUME_VARIABLES_HH #define DUMUX_RICHARDS_VOLUME_VARIABLES_HH @@ -30,8 +30,10 @@ namespace Dumux /*! * \ingroup RichardsModel - * \brief Contains the quantities which are are constant within a - * finite volume in the Richards model. + * \brief Volume averaged quantities required by the Richards model. + * + * This contains the quantities which are are constant within a finite + * volume in the Richards model */ template <class TypeTag> class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag> @@ -65,6 +67,15 @@ class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag> public: /*! * \brief Update all quantities for a given control volume. + * + * \param priVars The primary variables as a vector for the finite + * volume. + * \param problem The physical problem which needs to be solved. + * \param element The DUNE Codim<0> enitity which intersects + * the control volume of the box method + * \param scvIdx The local index of the sub control volume inside the element + * \param isOldSol Specifies whether the solution is from + * the previous time step or from the current one */ void update(const PrimaryVariables &priVars, const Problem &problem, @@ -87,9 +98,10 @@ public: problem); // material law parameters + Scalar pnRef = problem.referencePressure(element, elemGeom, scvIdx); const MaterialLawParams &materialParams = problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx); - fluidState_.update(problem, materialParams, priVars, temperature_); + fluidState_.update(pnRef, materialParams, priVars, temperature_); mobility_[wPhaseIdx] = MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx)) / @@ -107,28 +119,52 @@ public: } /*! - * \brief Returns the effective saturation of a given phase within - * the control volume. + * \brief Returns the average porosity [] within the control volume. + * + * The porosity is defined as the ratio of the pore space to the + * total volume, i.e. \f[ \Phi := \frac{V_{pore}}{V_{pore} + V_{rock}} \f] + */ + Scalar porosity() const + { return porosity_; } + + /*! + * \brief Returns the average absolute saturation [] of a given + * fluid phase within the finite volume. + * + * The saturation of a fluid phase is defined as the fraction of + * the pore volume filled by it, i.e. + * \f[ S_\alpha := \frac{V_\alpha}{V_{pore}} = \phi \frac{V_\alpha}{V} \f] + * + * \param phaseIdx The index of the fluid phase */ Scalar saturation(int phaseIdx) const { return fluidState_.saturation(phaseIdx); } /*! - * \brief Returns the mass density of a given phase within the - * control volume. + * \brief Returns the average mass density [kg/m^3] of a given + * fluid phase within the control volume. + * + * \param phaseIdx The index of the fluid phase */ Scalar density(int phaseIdx) const { return fluidState_.density(phaseIdx); } /*! - * \brief Returns the effective pressure of a given phase within + * \brief Returns the effective pressure [Pa] of a given phase within * the control volume. + * + * For the non-wetting phase (i.e. the gas phase), we assume + * infinite mobility, which implies that the non-wetting phase + * pressure is equal to the finite volume's reference pressure + * defined by the problem. + * + * \param phaseIdx The index of the fluid phase */ Scalar pressure(int phaseIdx) const { return fluidState_.phasePressure(phaseIdx); } /*! - * \brief Returns temperature inside the sub-control volume. + * \brief Returns average temperature [K] inside the control volume. * * Note that we assume thermodynamic equilibrium, i.e. the * temperature of the rock matrix and of all fluid phases are @@ -138,24 +174,30 @@ public: { return fluidState_.temperature(); } /*! - * \brief Returns the effective mobility of a given phase within + * \brief Returns the effective mobility [1/(Pa s)] of a given phase within * the control volume. + * + * The mobility of a fluid phase is defined as the relative + * permeability of the phase (given by the chosen material law) + * divided by the dynamic viscosity of the fluid, i.e. + * \f[ \lambda_\alpha := \frac{k_{r\alpha}}{\mu_\alpha} \f] + * + * \param phaseIdx The index of the fluid phase */ Scalar mobility(int phaseIdx) const { return mobility_[phaseIdx]; } /*! - * \brief Returns the effective capillary pressure within the control volume. + * \brief Returns the effective capillary pressure [Pa] within the + * control volume. + * + * The capillary pressure is defined as the difference in + * pressures of the non-wetting and the wetting phase, i.e. + * \f[ p_c = p_n - p_w \f] */ Scalar capillaryPressure() const { return fluidState_.capillaryPressure(); } - /*! - * \brief Returns the average porosity within the control volume. - */ - Scalar porosity() const - { return porosity_; } - protected: void updateTemperature_(const PrimaryVariables &priVars, diff --git a/test/boxmodels/richards/richardslensproblem.hh b/test/boxmodels/richards/richardslensproblem.hh index cff1468081ca88e283df6d46f917e65d00021858..325ece4387213ac95bee75b16549cd0e77bc379e 100644 --- a/test/boxmodels/richards/richardslensproblem.hh +++ b/test/boxmodels/richards/richardslensproblem.hh @@ -1,6 +1,6 @@ /***************************************************************************** * Copyright (C) 2009 by Onur Dogan * - * Copyright (C) 2009 by Andreas Lauser * + * Copyright (C) 2009-2010 by Andreas Lauser * * Institute of Hydraulic Engineering * * University of Stuttgart, Germany * * email: <givenname>.<name>@iws.uni-stuttgart.de * @@ -13,6 +13,13 @@ * * * This program is distributed WITHOUT ANY WARRANTY. * *****************************************************************************/ +/*! + * \file + * + * \brief A water infiltration problem with a low-permeability lens + * embedded into a high- permeability domain which uses the + * Richards box model. + */ #ifndef DUMUX_RICHARDS_LENSPROBLEM_HH #define DUMUX_RICHARDS_LENSPROBLEM_HH @@ -39,8 +46,7 @@ namespace Properties { NEW_TYPE_TAG(RichardsLensProblem, INHERITS_FROM(BoxRichards)); -// Set the grid type -// Set the grid type +// Set the grid type. Use UG if available, else SGrid #if HAVE_UG SET_TYPE_PROP(RichardsLensProblem, Grid, Dune::UGGrid<2>); #else @@ -48,6 +54,8 @@ SET_PROP(RichardsLensProblem, Grid) { typedef Dune::SGrid<2, 2> type; }; //SET_TYPE_PROP(RichardsLensProblem, Grid, Dune::YaspGrid<2>); #endif +// set the local finite element space to be used to calculate the +// gradients in the flux calculation SET_PROP(RichardsLensProblem, LocalFEMSpace) { typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; @@ -59,11 +67,9 @@ public: // typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices }; -// Set the problem property +// Set the phsical problem to be solved SET_PROP(RichardsLensProblem, Problem) -{ - typedef Dumux::RichardsLensProblem<TTAG(RichardsLensProblem)> type; -}; +{ typedef Dumux::RichardsLensProblem<TypeTag> type; }; // Set the wetting phase SET_PROP(RichardsLensProblem, WettingPhase) @@ -76,11 +82,9 @@ public: // Set the spatial parameters SET_PROP(RichardsLensProblem, SpatialParameters) -{ - typedef Dumux::RichardsLensSpatialParameters<TypeTag> type; -}; +{ typedef Dumux::RichardsLensSpatialParameters<TypeTag> type; }; -// Enable gravity +// Enable gravity? SET_BOOL_PROP(RichardsLensProblem, EnableGravity, true); // Write the intermediate results of the newton method? @@ -88,17 +92,27 @@ SET_BOOL_PROP(RichardsLensProblem, NewtonWriteConvergence, false); } /*! - * \ingroup RichardsBoxProblems - * \brief Base class for defining an instance of a Richard`s problem, where water is infiltrating into an initially unsaturated zone. + * \ingroup RichardsModel + * + * \brief A water infiltration problem with a low-permeability lens + * embedded into a high- permeability domain which uses the + * Richards box model. * - * The domain is box shaped. Left and right boundaries are Dirichlet boundaries with fixed water pressure (fixed Saturation Sw = 0), - * bottom boundary is closed (Neumann 0 boundary), the top boundary (Neumann 0 boundary) is also closed except for infiltration section, - * where water is infiltrating into an initially unsaturated zone. Linear capillary pressure-saturation relationship is used. + * The domain is box shaped. Left and right boundaries are Dirichlet + * boundaries with fixed water pressure (fixed Saturation $S_w = 0$), + * bottom boundary is closed (Neumann 0 boundary), the top boundary + * (Neumann 0 boundary) is also closed except for infiltration + * section, where water is infiltrating into an initially unsaturated + * porous medium. This problem is very similar the the LensProblem + * which uses the TwoPBoxModel, with the main difference being that + * the domain is initally fully saturated by gas instead of water and + * water instead of a %DNAPL infiltrates from the top. * * To run the simulation execute the following line in shell: - * <tt>./lens_richards ./grids/lens_richards_2d.dgf 10000 1</tt> + * <tt>./test_richards grids/richardslens.dgf 10e6 100</tt> * - * where start simulation time = 1 second, end simulation time = 10000 seconds + * where the initial time step is 100 seconds, and the end of the + * simulation time is 10,000,000 seconds (115.7 days) */ template <class TypeTag = TTAG(RichardsLensProblem) > class RichardsLensProblem : public RichardsBoxProblem<TypeTag> @@ -130,6 +144,16 @@ class RichardsLensProblem : public RichardsBoxProblem<TypeTag> typedef Dune::FieldVector<Scalar, dim> GlobalPosition; public: + /*! + * \brief Constructor + * + * \param timeManager The Dumux TimeManager for simulation management. + * \param gridView The grid view on the spatial domain of the problem + * \param lensLowerLeft The lower left coordinate of the + * low-permeability lens + * \param lensUpperRight The upper right coordinate of the + * low-permeability lens + */ RichardsLensProblem(TimeManager &timeManager, const GridView &gridView, const GlobalPosition &lensLowerLeft, @@ -138,6 +162,7 @@ public: { lensLowerLeft_=lensLowerLeft; lensUpperRight_=lensUpperRight; + pnRef_ = 1e5; this->spatialParameters().setLensCoords(lensLowerLeft_, lensUpperRight_); } @@ -147,7 +172,7 @@ public: // \{ /*! - * \brief The problem name. + * \brief The problem name * * This is used as a prefix for files generated by the simulation. */ @@ -155,26 +180,37 @@ public: { return "richardslens"; } /*! - * \brief Returns the temperature within the domain. + * \brief Returns the temperature [K] within a finite volume * * This problem assumes a temperature of 10 degrees Celsius. + * + * \param element The DUNE Codim<0> entity which intersects with + * the finite volume in question + * \param fvElemGeom The finite volume geometry of the element + * \param scvIdx The sub control volume index inside the finite + * volume geometry */ Scalar temperature(const Element &element, const FVElementGeometry &fvElemGeom, int scvIdx) const - { - return 283.15; // -> 10°C - }; + { return 273.15 + 10; }; // -> 10°C /*! - * \brief Returns the reference pressure of the nonwetting phase. + * \brief Returns the reference pressure [Pa] of the non-wetting + * fluid phase within a finite volume + * + * This problem assumes a constant reference pressure of 1 bar. * - * This problem assumes a pressure of 1 bar. + * \param element The DUNE Codim<0> entity which intersects with + * the finite volume in question + * \param fvElemGeom The finite volume geometry of the element + * \param scvIdx The sub control volume index inside the finite + * volume geometry */ - Scalar pNreference() const - { - return 1e5; // reference non-wetting phase pressure [Pa] used for viscosity and density calculations - }; + Scalar referencePressure(const Element &element, + const FVElementGeometry &fvElemGeom, + int scvIdx) const + { return pnRef_; }; // \} @@ -225,6 +261,15 @@ public: * * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. + * + * \param element The DUNE Codim<0> entity which intersects with + * the finite volume in question + * \param fvElemGeom The finite volume geometry of the element + * \param is The DUNE boundary intersection of the boundary segment + * \param scvIdx The sub control volume index of the finite + * volume geometry + * \param boundaryFaceIdx The index of the boundary face of the + * finite volume geometry */ void neumann(PrimaryVariables &values, const Element &element, @@ -250,13 +295,20 @@ public: // \{ /*! - * \brief Evaluate the source term for all phases within a given - * sub-control-volume. + * \brief Specify the source term [kg/m^3] for the wetting phase + * within a given sub-control-volume. * * For this method, the \a values parameter stores the rate * wetting phase mass is generated or annihilate per volume * unit. Positive values mean that mass is created, negative ones * mean that it vanishes. + * + * \param values Storage for all values for the source terms + * \param element The DUNE Codim<0> entity which intersects with + * the finite volume in question + * \param fvElemGeom The finite volume geometry of the element + * \param scvIdx The sub control volume index of the finite + * volume geometry */ void source(PrimaryVariables &values, const Element &element, @@ -267,10 +319,17 @@ public: } /*! - * \brief Evaluate the initial value for a control volume. + * \brief Evaluate the initial values for a control volume. * * For this method, the \a values parameter stores primary * variables. + * + * \param values Storage for all primary variables of the initial condition + * \param element The DUNE Codim<0> entity which intersects with + * the finite volume in question + * \param fvElemGeom The finite volume geometry of the element + * \param scvIdx The sub control volume index of the finite + * volume geometry */ void initial(PrimaryVariables &values, const Element &element, @@ -291,7 +350,7 @@ private: Scalar pc = MaterialLaw::pC(this->spatialParameters().materialLawParams(pos), Sw); - values[pwIdx] = pNreference() - pc; + values[pwIdx] = pnRef_ - pc; } bool onLeftBoundary_(const GlobalPosition &globalPos) const @@ -322,6 +381,7 @@ private: } static const Scalar eps_ = 3e-6; + Scalar pnRef_; GlobalPosition lensLowerLeft_; GlobalPosition lensUpperRight_; }; diff --git a/test/boxmodels/richards/richardslensspatialparameters.hh b/test/boxmodels/richards/richardslensspatialparameters.hh index d8feb7c0141248e917016b02acc1246feb5890fd..7fa56e6350ce6294f448a3a91121fd6d8e4a9d68 100644 --- a/test/boxmodels/richards/richardslensspatialparameters.hh +++ b/test/boxmodels/richards/richardslensspatialparameters.hh @@ -3,7 +3,7 @@ * Copyright (C) 2010 by Markus Wolff * * Copyright (C) 2007-2008 by Klaus Mosthaf * * Copyright (C) 2007-2008 by Bernd Flemisch * - * Copyright (C) 2008-2009 by Andreas Lauser * + * Copyright (C) 2008-2010 by Andreas Lauser * * Institute of Hydraulic Engineering * * University of Stuttgart, Germany * * email: <givenname>.<name>@iws.uni-stuttgart.de * @@ -16,25 +16,23 @@ * * * This program is distributed WITHOUT ANY WARRANTY. * *****************************************************************************/ -#ifndef DUMUX_LENSSPATIALPARAMETERS_HH -#define DUMUX_LENSSPATIALPARAMETERS_HH +/*! + * \brief The spatial parameters for the RichardsLensProblem + */ +#ifndef DUMUX_RICHARDS_LENS_SPATIAL_PARAMETERS_HH +#define DUMUX_RICHARDS_LENS_SPATIAL_PARAMETERS_HH #include <dumux/material/spatialparameters/boxspatialparameters.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> -/** - * @file - * @brief Class for defining spatial parameters - * @author Bernd Flemisch, Klaus Mosthaf, Markus Wolff - */ - namespace Dumux { -/** \todo Please doc me! */ - +/*! + * \brief The spatial parameters for the RichardsLensProblem + */ template<class TypeTag> class RichardsLensSpatialParameters : public BoxSpatialParameters<TypeTag> { @@ -59,10 +57,21 @@ class RichardsLensSpatialParameters : public BoxSpatialParameters<TypeTag> typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; public: - // define the material law parameterized by absolute saturations + /*! + * \brief The material law to be used. + * + * This problem uses the RegularizedVanGenuchten material law. + */ typedef EffToAbsLaw<EffectiveLaw> MaterialLaw; + //! The parameters of the material law to be used typedef typename MaterialLaw::Params MaterialLawParams; + /*! + * \brief Constructor + * + * \param gridView The DUNE GridView representing the spatial + * domain of the problem. + */ RichardsLensSpatialParameters(const GridView& gridView) : ParentType(gridView) { @@ -91,13 +100,11 @@ public: } /*! - * \brief Apply the intrinsic permeability tensor to a pressure - * potential gradient. + * \brief Returns the intrinsic permeability tensor [m^2] at a given location * - * \param element The current finite element + * \param element An arbitrary DUNE Codim<0> entity of the grid view * \param fvElemGeom The current finite volume geometry of the element - * \param scvIdx The index sub-control volume face where the - * intrinsic velocity ought to be calculated. + * \param scvIdx The index of the sub-control volume */ Scalar intrinsicPermeability(const Element &element, const FVElementGeometry &fvElemGeom, @@ -108,13 +115,26 @@ public: return lensK_; return outerK_; } - + + /*! + * \brief Returns the porosity [] at a given location + * + * \param element An arbitrary DUNE Codim<0> entity of the grid view + * \param fvElemGeom The current finite volume geometry of the element + * \param scvIdx The index of the sub-control volume + */ Scalar porosity(const Element &element, const FVElementGeometry &fvElemGeom, int scvIdx) const { return 0.4; } - // return the brooks-corey context depending on the position + /*! + * \brief Returns the parameters for the material law at a given location + * + * \param element An arbitrary DUNE Codim<0> entity of the grid view + * \param fvElemGeom The current finite volume geometry of the element + * \param scvIdx The index of the sub-control volume + */ const MaterialLawParams& materialLawParams(const Element &element, const FVElementGeometry &fvElemGeom, int scvIdx) const @@ -122,7 +142,14 @@ public: return materialLawParams(fvElemGeom.subContVol[scvIdx].global); } - // return the brooks-corey context depending on the position + /*! + * \brief Returns the parameters for the material law at a given location + * + * This method is not actually required by the Richards model, but provided + * for the convenience of the RichardsLensProblem + * + * \param globalPos A global coordinate vector + */ const MaterialLawParams& materialLawParams(const GlobalPosition &globalPos) const { @@ -131,7 +158,14 @@ public: return outerMaterialParams_; } - //! Set the bounding box of the fine-sand lens + /*! + * \brief Set the bounding box of the low-permeability lens + * + * This method is not actually required by the Richards model, but provided + * for the convenience of the RichardsLensProblem + * + * \param globalPos A global coordinate vector + */ void setLensCoords(const GlobalPosition& lensLowerLeft, const GlobalPosition& lensUpperRight) { @@ -159,5 +193,6 @@ private: }; } // end namespace + #endif