diff --git a/dumux/boxmodels/richards/richardsfluxvariables.hh b/dumux/boxmodels/richards/richardsfluxvariables.hh index 403ebd91faf18e2c4290286826688fc269bfdf59..cc4ba20247b2c8e2517aa105025c5cf08ab05233 100644 --- a/dumux/boxmodels/richards/richardsfluxvariables.hh +++ b/dumux/boxmodels/richards/richardsfluxvariables.hh @@ -46,7 +46,7 @@ class RichardsFluxVariables { typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParams; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { @@ -55,11 +55,11 @@ class RichardsFluxVariables typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GridView::template Codim<0>::Entity Element; - enum { dimWorld = GridView::dimensionworld}; + enum { dim = GridView::dimension}; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef Dune::FieldVector<Scalar, dimWorld> Vector; - typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor; + typedef Dune::FieldVector<Scalar, dim> DimVector; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimTensor; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename FVElementGeometry::SubControlVolumeFace SCVFace; @@ -80,13 +80,13 @@ public: */ RichardsFluxVariables(const Problem &problem, const Element &element, - const FVElementGeometry &fvElemGeom, - int scvfIdx, + const FVElementGeometry &fvGeometry, + const int faceIdx, const ElementVolumeVariables &elemVolVars, const bool onBoundary = false) - : fvElemGeom_(fvElemGeom), onBoundary_(onBoundary) + : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary) { - scvfIdx_ = scvfIdx; + calculateGradients_(problem, element, elemVolVars); calculateK_(problem, element, elemVolVars); @@ -95,13 +95,13 @@ public: /*! * \brief Return the intrinsic permeability \f$\mathrm{[m^2]}\f$. */ - const Tensor &intrinsicPermeability() const + const DimTensor &intrinsicPermeability() const { return K_; } /*! * \brief Return the pressure potential gradient \f$\mathrm{[Pa/m]}\f$ */ - const Vector &potentialGradW() const + const DimVector &potentialGradW() const { return potentialGrad_; } /*! @@ -128,18 +128,22 @@ public: 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 + * \brief The face of the current sub-control volume. This may be either + * an inner sub-control-volume face or a face on the boundary. */ const SCVFace &face() const { if (onBoundary_) - return fvElemGeom_.boundaryFace[scvfIdx_]; + return fvGeometry_.boundaryFace[faceIdx_]; else - return fvElemGeom_.subContVolFace[scvfIdx_]; + return fvGeometry_.subContVolFace[faceIdx_]; } + + protected: void calculateGradients_(const Problem &problem, const Element &element, @@ -148,15 +152,16 @@ protected: potentialGrad_ = 0.0; // calculate gradients for (int idx = 0; - idx < fvElemGeom_.numVertices; + idx < fvGeometry_.numVertices; idx++) // loop over adjacent vertices { // FE gradient at vertex index - const Vector &feGrad = face().grad[idx]; + const DimVector &feGrad = face().grad[idx]; // the pressure gradient - Vector tmp(feGrad); + DimVector tmp(feGrad); tmp *= elemVolVars[idx].pressure(wPhaseIdx); + potentialGrad_ += tmp; } @@ -180,8 +185,8 @@ protected: // estimate the gravitational acceleration at a given SCV face // using the arithmetic mean - Vector f(problem.boxGravity(element, fvElemGeom_, face().i)); - f += problem.boxGravity(element, fvElemGeom_, face().j); + DimVector f(problem.boxGravity(element, fvGeometry_, face().i)); + f += problem.boxGravity(element, fvGeometry_, face().j); f /= 2; // make it a force @@ -196,26 +201,26 @@ protected: const Element &element, const ElementVolumeVariables &elemDat) { - const SpatialParameters &sp = problem.spatialParameters(); + const SpatialParams &spatialParams = problem.spatialParams(); // calculate the intrinsic permeability - sp.meanK(K_, - sp.intrinsicPermeability(element, - fvElemGeom_, + spatialParams.meanK(K_, + spatialParams.intrinsicPermeability(element, + fvGeometry_, face().i), - sp.intrinsicPermeability(element, - fvElemGeom_, + spatialParams.intrinsicPermeability(element, + fvGeometry_, face().j)); } - const FVElementGeometry &fvElemGeom_; - int scvfIdx_; + const FVElementGeometry &fvGeometry_; + const int faceIdx_; const bool onBoundary_; // gradients - Vector potentialGrad_; + DimVector potentialGrad_; // intrinsic permeability - Tensor K_; + DimTensor K_; }; } // end namepace diff --git a/dumux/boxmodels/richards/richardslocalresidual.hh b/dumux/boxmodels/richards/richardslocalresidual.hh index 4ee30cb185092285d0e1fd8a16135d9f6fa7befa..30304a2716dfc45ef35fc7146fa864bae83ec1b5 100644 --- a/dumux/boxmodels/richards/richardslocalresidual.hh +++ b/dumux/boxmodels/richards/richardslocalresidual.hh @@ -55,10 +55,10 @@ class RichardsLocalResidual : public BoxLocalResidual<TypeTag> }; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - enum { dimWorld = GridView::dimensionworld}; + enum { dim = GridView::dimension}; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef Dune::FieldVector<Scalar, dimWorld> Vector; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: /*! @@ -80,12 +80,12 @@ public: * * This function should not include the source and sink terms. * - * \param result Stores the average mass per unit volume for each phase \f$\mathrm{[kg/m^3]}\f$ + * \param storage Stores the average mass per unit volume for each phase \f$\mathrm{[kg/m^3]}\f$ * \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 + void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const { // if flag usePrevSol is set, the solution from the previous // time step is used, otherwise the current solution is @@ -98,7 +98,7 @@ public: this->curVolVars_(scvIdx); // partial time derivative of the wetting phase mass - result[contiEqIdx] = + storage[contiEqIdx] = volVars.density(wPhaseIdx) * volVars.saturation(wPhaseIdx) * volVars.porosity(); @@ -117,18 +117,18 @@ public: * \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 scvfIdx, const bool onBoundary=false) const + void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const { FluxVariables fluxVars(this->problem_(), - this->elem_(), - this->fvElemGeom_(), - scvfIdx, + this->element_(), + this->fvGeometry_(), + faceIdx, this->curVolVars_(), onBoundary); // calculate the flux in the normal direction of the // current sub control volume face - Vector tmpVec; + DimVector tmpVec; fluxVars.intrinsicPermeability().mv(fluxVars.potentialGradW(), tmpVec); Scalar normalFlux = -(tmpVec*fluxVars.face().normal); @@ -149,16 +149,16 @@ public: /*! * \brief Calculate the source term of the equation * - * \param q Stores the average source term of all phases inside a + * \param source Stores the average source term of all phases inside a * sub-control volume of the current element \f$\mathrm{[kg/(m^3 * s)]}\f$ * \param scvIdx The sub control volume index inside the current * element */ - void computeSource(PrimaryVariables &q, int scvIdx) + void computeSource(PrimaryVariables &source, int scvIdx) { - this->problem_().boxSDSource(q, - this->elem_(), - this->fvElemGeom_(), + this->problem_().boxSDSource(source, + this->element_(), + this->fvGeometry_(), scvIdx, this->curVolVars_()); } diff --git a/dumux/boxmodels/richards/richardsnewtoncontroller.hh b/dumux/boxmodels/richards/richardsnewtoncontroller.hh index 241f1872d154b626874acd9c704c80ac47043efa..fbe1451472a7f8017b42ab9849d65f901ee40bb2 100644 --- a/dumux/boxmodels/richards/richardsnewtoncontroller.hh +++ b/dumux/boxmodels/richards/richardsnewtoncontroller.hh @@ -46,7 +46,7 @@ class RichardsNewtonController : public NewtonController<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParams; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; @@ -90,21 +90,21 @@ public: return; // clamp saturation change to at most 20% per iteration - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeomtry; const GridView &gv = this->problem_().gridView(); ElementIterator eIt = gv.template begin<0>(); const ElementIterator eEndIt = gv.template end<0>(); for (; eIt != eEndIt; ++eIt) { - fvElemGeom.update(gv, *eIt); - for (int i = 0; i < fvElemGeom.numVertices; ++i) { + fvGeomtry.update(gv, *eIt); + for (int i = 0; i < fvGeomtry.numVertices; ++i) { int globI = this->problem_().vertexMapper().map(*eIt, i, dim); // calculate the old wetting phase saturation - const SpatialParameters &sp = this->problem_().spatialParameters(); - const MaterialLawParams &mp = sp.materialLawParams(*eIt, fvElemGeom, i); + const SpatialParams &spatialParams = this->problem_().spatialParams(); + const MaterialLawParams &mp = spatialParams.materialLawParams(*eIt, fvGeomtry, i); Scalar pcMin = MaterialLaw::pC(mp, 1.0); Scalar pW = uLastIter[globI][pwIdx]; - Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvElemGeom, i), + Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvGeomtry, i), pW + pcMin); Scalar pcOld = pN - pW; Scalar SwOld = std::max<Scalar>(0.0, MaterialLaw::Sw(mp, pcOld)); diff --git a/dumux/boxmodels/richards/richardsproblem.hh b/dumux/boxmodels/richards/richardsproblem.hh index e3e3c9c1200c99710ed6b329d685f1ca64a95c00..000bbf9466bc58ab888554a89b11adc698b6bae5 100644 --- a/dumux/boxmodels/richards/richardsproblem.hh +++ b/dumux/boxmodels/richards/richardsproblem.hh @@ -84,7 +84,7 @@ public: * \param scvIdx The local index of the sub control volume inside the element */ Scalar referencePressure(const Element &element, - const FVElementGeometry fvGeom, + const FVElementGeometry fvGeometry, int scvIdx) const { DUNE_THROW(Dune::NotImplemented, "referencePressure() method not implemented by the actual problem"); }; // \} diff --git a/dumux/boxmodels/richards/richardsvolumevariables.hh b/dumux/boxmodels/richards/richardsvolumevariables.hh index 27f587fb9a117c9073061638a54836b5040b4e78..52e98f4ed2a0521e43199a1c627166482477d333 100644 --- a/dumux/boxmodels/richards/richardsvolumevariables.hh +++ b/dumux/boxmodels/richards/richardsvolumevariables.hh @@ -85,37 +85,37 @@ public: * \param isOldSol Specifies whether the solution is from * the previous time step or from the current one */ - void update(const PrimaryVariables &priVars, + void update(const PrimaryVariables &primaryVariables, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, + const FVElementGeometry &fvGeometry, int scvIdx, bool isOldSol) { assert(!FluidSystem::isLiquid(nPhaseIdx)); - ParentType::update(priVars, + ParentType::update(primaryVariables, problem, element, - elemGeom, + fvGeometry, scvIdx, isOldSol); - completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_); + completeFluidState(primaryVariables, problem, element, fvGeometry, scvIdx, fluidState_); ////////// // specify the other parameters ////////// const MaterialLawParams &matParams = - problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx); + problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx); relativePermeabilityWetting_ = MaterialLaw::krw(matParams, fluidState_.saturation(wPhaseIdx)); - porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx); + porosity_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx); // energy related quantities not contained in the fluid state - asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol); + asImp_().updateEnergy_(primaryVariables, problem, element, fvGeometry, scvIdx, isOldSol); } /*! @@ -124,19 +124,19 @@ public: static void completeFluidState(const PrimaryVariables& primaryVariables, const Problem& problem, const Element& element, - const FVElementGeometry& elementGeometry, + const FVElementGeometry& fvGeometry, int scvIdx, FluidState& fluidState) { // temperature Scalar t = Implementation::temperature_(primaryVariables, problem, element, - elementGeometry, scvIdx); + fvGeometry, scvIdx); fluidState.setTemperature(t); // pressures - Scalar pnRef = problem.referencePressure(element, elementGeometry, scvIdx); + Scalar pnRef = problem.referencePressure(element, fvGeometry, scvIdx); const MaterialLawParams &matParams = - problem.spatialParameters().materialLawParams(element, elementGeometry, scvIdx); + problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx); Scalar minPc = MaterialLaw::pC(matParams, 1.0); fluidState.setPressure(wPhaseIdx, primaryVariables[pwIdx]); fluidState.setPressure(nPhaseIdx, std::max(pnRef, primaryVariables[pwIdx] + minPc)); @@ -257,13 +257,13 @@ public: { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); } protected: - static Scalar temperature_(const PrimaryVariables &priVars, + static Scalar temperature_(const PrimaryVariables &primaryVariables, const Problem& problem, const Element &element, - const FVElementGeometry &elemGeom, + const FVElementGeometry &fvGeometry, int scvIdx) { - return problem.boxTemperature(element, elemGeom, scvIdx); + return problem.boxTemperature(element, fvGeometry, scvIdx); } /*! @@ -272,8 +272,8 @@ protected: void updateEnergy_(const PrimaryVariables &sol, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int vertIdx, + const FVElementGeometry &fvGeometry, + int scvIdx, bool isOldSol) { }