diff --git a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh index 6f49c4ad03e72a2acb7e2ddd6a2d07b4de1eda34..c89d83adb7ab506f17cd527886767bb47af678d5 100644 --- a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh +++ b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh @@ -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 mole-/mass-fraction gradients, phase densities at + * This means pressure and mole-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 mole-/mass-fraction gradients, phase densities at + * This means pressure and mole-fraction gradients, phase densities at * the intergration point, etc. */ template <class TypeTag> @@ -58,7 +58,7 @@ class OnePTwoCFluxVariables typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { @@ -75,8 +75,8 @@ class OnePTwoCFluxVariables typedef typename GridView::ctype CoordScalar; typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; 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> DimMatrix; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename FVElementGeometry::SubControlVolumeFace SCVFace; @@ -87,33 +87,33 @@ public: * * \param problem The problem * \param element The finite element - * \param elemGeom The finite-volume geometry in the box scheme + * \param fvGeometry 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 elemVolVars 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, - const FVElementGeometry &elemGeom, - int faceIdx, - const ElementVolumeVariables &elemDat, - bool onBoundary = false) - : fvGeom_(elemGeom), faceIdx_(faceIdx), onBoundary_(onBoundary) + const FVElementGeometry &fvGeometry, + const int faceIdx, + const ElementVolumeVariables &elemVolVars, + const bool onBoundary = false) + : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary) { - viscosityAtIP_ = Scalar(0); - molarDensityAtIP_ = Scalar(0); - densityAtIP_ = Scalar(0); + viscosity_ = Scalar(0); + molarDensity_ = Scalar(0); + density_ = Scalar(0); potentialGrad_ = Scalar(0); concentrationGrad_ = Scalar(0); - moleFracGrad_ = Scalar(0); - massFracGrad_ = Scalar(0); - - calculateGradients_(problem, element, elemDat); - calculateK_(problem, element, elemDat); - calculateVelocities_(problem, element, elemDat); - calculateDiffCoeffPM_(problem, element, elemDat); - calculateDispersionTensor_(problem, element, elemDat); + moleFractionGrad_ = Scalar(0); + massFractionGrad_ = Scalar(0); + + calculateGradients_(problem, element, elemVolVars); + calculateK_(problem, element, elemVolVars); + calculateVelocities_(problem, element, elemVolVars); + calculatePorousDiffCoeff_(problem, element, elemVolVars); + calculateDispersionTensor_(problem, element, elemVolVars); }; public: @@ -136,7 +136,7 @@ public: * \brief Return the pressure potential multiplied with the * intrinsic permeability as vector (for velocity output). */ - Vector Kmvp() const + DimVector Kmvp() const { return Kmvp_; } /*! @@ -146,27 +146,27 @@ public: const SCVFace &face() const { if (onBoundary_) - return fvGeom_.boundaryFace[faceIdx_]; + return fvGeometry_.boundaryFace[faceIdx_]; else - return fvGeom_.subContVolFace[faceIdx_]; + return fvGeometry_.subContVolFace[faceIdx_]; } /*! * \brief Return the intrinsic permeability tensor \f$\mathrm{[m^2]}\f$. */ - const Tensor &intrinsicPermeability() const + const DimMatrix &intrinsicPermeability() const { return K_; } /*! * \brief Return the dispersion tensor \f$\mathrm{[m^2/s]}\f$. */ - const Tensor &dispersionTensor() const + const DimMatrix &dispersionTensor() const { return dispersionTensor_; } /*! * \brief Return the pressure potential gradient \f$\mathrm{[Pa/m]}\f$. */ - const Vector &potentialGrad() const + const DimVector &potentialGrad() const { return potentialGrad_; } /*! @@ -174,7 +174,8 @@ public: * * \param compIdx The index of the considered component */ - const Vector &concentrationGrad(int compIdx) const + DUMUX_DEPRECATED_MSG("use moleFractionGrad instead") + const DimVector &concentrationGrad(int compIdx) const { if (compIdx != 1) { DUNE_THROW(Dune::InvalidStateException, @@ -185,64 +186,105 @@ public: }; /*! - * \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 - */ - const Vector &moleFracGrad(int compIdx) const - { - if (compIdx != 1) - { DUNE_THROW(Dune::InvalidStateException, - "The 1p2c model is supposed to need " - "only the concentration gradient of " - "the second component!"); } - return moleFracGrad_; - }; - - /*! - * \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 - */ - const Vector &massFracGrad(int compIdx) const - { - if (compIdx != 1) - { DUNE_THROW(Dune::InvalidStateException, - "The 1p2c model is supposed to need " - "only the concentration gradient of " - "the second component!"); } - return massFracGrad_; - }; + * \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 + */ + const DimVector &moleFractionGrad(int compIdx) const + { + if (compIdx != 1) + { DUNE_THROW(Dune::InvalidStateException, + "The 1p2c model is supposed to need " + "only the concentration gradient of " + "the second component!"); } + return moleFractionGrad_; + }; + + /*! + * \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 + */ + DUMUX_DEPRECATED_MSG("use moleFractionGrad instead") + const DimVector &moleFracGrad(int compIdx) const + { + if (compIdx != 1) + { DUNE_THROW(Dune::InvalidStateException, + "The 1p2c model is supposed to need " + "only the concentration gradient of " + "the second component!"); } + return moleFractionGrad(compIdx); + }; + + /*! + * \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 + */ + DUMUX_DEPRECATED_MSG("use moleFractionGrad instead") + const DimVector &massFracGrad(int compIdx) const + { + if (compIdx != 1) + { DUNE_THROW(Dune::InvalidStateException, + "The 1p2c model is supposed to need " + "only the concentration gradient of " + "the second component!"); } + return massFractionGrad_; + }; /*! * \brief The binary diffusion coefficient for each fluid phase in the porous medium \f$\mathrm{[m^2/s]}\f$. */ Scalar porousDiffCoeff() const { - // TODO: tensorial diffusion coefficients - return diffCoeffPM_; + // TODO: tensorial porousDiffCoeff_usion coefficients + return porousDiffCoeff_; }; /*! * \brief Return viscosity \f$\mathrm{[Pa s]}\f$ of a phase at the integration * point. */ + Scalar viscosity() const + { return viscosity_;} + + /*! + * \brief Return viscosity \f$\mathrm{[Pa s]}\f$ of a phase at the integration + * point. + */ + DUMUX_DEPRECATED_MSG("use viscosity() instead") Scalar viscosityAtIP() const - { return viscosityAtIP_;} + { return viscosity();} + + /*! + * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration + * point. + */ + Scalar molarDensity() const + { return molarDensity_; } /*! * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration * point. */ + DUMUX_DEPRECATED_MSG("use molarDensity() instead") Scalar molarDensityAtIP() const - { return molarDensityAtIP_; } + { return molarDensity(); } + + /*! + * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration + * point. + */ + Scalar density() const + { return density_; } /*! * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration * point. */ + DUMUX_DEPRECATED_MSG("use density( instead") Scalar densityAtIP() const - { return densityAtIP_; } + { return density(); } /*! * \brief Given the intrinsic permeability times the pressure @@ -293,21 +335,21 @@ protected: const Element &element, const ElementVolumeVariables &elemVolVars) { - const VolumeVariables &vVars_i = elemVolVars[face().i]; - const VolumeVariables &vVars_j = elemVolVars[face().j]; + const VolumeVariables &volVarsI = elemVolVars[face().i]; + const VolumeVariables &volVarsJ = elemVolVars[face().j]; - Vector tmp; - //The decision of the if-statement depends on the function useTwoPointGradient(const Element &elem, + DimVector tmp; + //The decision of the if-statement depends on the function useTwoPointGradient(const Element &element, //int vertexI,int vertexJ) defined in test/tissue_tumor_spatialparameters.hh - if (!problem.spatialParameters().useTwoPointGradient(element, face().i, face().j)) { + if (!problem.spatialParams().useTwoPointGradient(element, face().i, face().j)) { // use finite-element gradients tmp = 0.0; for (int idx = 0; - idx < fvGeom_.numVertices; + idx < fvGeometry_.numVertices; idx++) // loop over adjacent vertices { // FE gradient at vertex idx - const Vector &feGrad = face().grad[idx]; + const DimVector &feGrad = face().grad[idx]; // the pressure gradient tmp = feGrad; @@ -322,40 +364,40 @@ protected: // the mole-fraction gradient tmp = feGrad; tmp *= elemVolVars[idx].moleFraction(comp1Idx); - moleFracGrad_ += tmp; + moleFractionGrad_ += tmp; // the mass-fraction gradient tmp = feGrad; tmp *= elemVolVars[idx].massFraction(comp1Idx); - massFracGrad_ += tmp; + massFractionGrad_ += tmp; // phase viscosity - viscosityAtIP_ += elemVolVars[idx].viscosity()*face().shapeValue[idx]; + viscosity_ += elemVolVars[idx].viscosity()*face().shapeValue[idx]; //phase molar density - molarDensityAtIP_ += elemVolVars[idx].molarDensity()*face().shapeValue[idx]; + molarDensity_ += elemVolVars[idx].molarDensity()*face().shapeValue[idx]; //phase density - densityAtIP_ += elemVolVars[idx].density()*face().shapeValue[idx]; + density_ += elemVolVars[idx].density()*face().shapeValue[idx]; } } else { // use two-point gradients - const GlobalPosition &posI = element.geometry().corner(face().i); - const GlobalPosition &posJ = element.geometry().corner(face().j); - tmp = posI; - tmp -= posJ; + const GlobalPosition &globalPosI = element.geometry().corner(face().i); + const GlobalPosition &globalPosJ = element.geometry().corner(face().j); + tmp = globalPosI; + tmp -= globalPosJ; Scalar dist = tmp.two_norm(); tmp = face().normal; tmp /= face().normal.two_norm()*dist; potentialGrad_ = tmp; - potentialGrad_ *= vVars_j.pressure() - vVars_i.pressure(); + potentialGrad_ *= volVarsJ.pressure() - volVarsI.pressure(); concentrationGrad_ = tmp; - concentrationGrad_ *= vVars_j.molarity(comp1Idx) - vVars_i.molarity(comp1Idx); - moleFracGrad_ = tmp; - moleFracGrad_ *= vVars_j.moleFraction(comp1Idx) - vVars_i.moleFraction(comp1Idx); + concentrationGrad_ *= volVarsJ.molarity(comp1Idx) - volVarsI.molarity(comp1Idx); + moleFractionGrad_ = tmp; + moleFractionGrad_ *= volVarsJ.moleFraction(comp1Idx) - volVarsI.moleFraction(comp1Idx); } /////////////// @@ -370,8 +412,8 @@ protected: // estimate the gravitational acceleration at a given SCV face // using the arithmetic mean - Vector f(problem.boxGravity(element, fvGeom_, face().i)); - f += problem.boxGravity(element, fvGeom_, face().j); + DimVector f(problem.boxGravity(element, fvGeometry_, face().i)); + f += problem.boxGravity(element, fvGeometry_, face().j); f /= 2; // make it a force @@ -389,19 +431,19 @@ protected: * * \param problem The considered problem file * \param element The considered element of the grid - * \param elemDat The parameters stored in the considered element + * \param elemVolVars The parameters stored in the considered element */ void calculateK_(const Problem &problem, const Element &element, - const ElementVolumeVariables &elemDat) + const ElementVolumeVariables &elemVolVars) { - const SpatialParameters &sp = problem.spatialParameters(); + const SpatialParams &sp = problem.spatialParams(); sp.meanK(K_, sp.intrinsicPermeability(element, - fvGeom_, + fvGeometry_, face().i), sp.intrinsicPermeability(element, - fvGeom_, + fvGeometry_, face().j)); } @@ -413,11 +455,11 @@ protected: * * \param problem The considered problem file * \param element The considered element of the grid - * \param elemDat The parameters stored in the considered element + * \param elemVolVars The parameters stored in the considered element */ void calculateVelocities_(const Problem &problem, const Element &element, - const ElementVolumeVariables &elemDat) + const ElementVolumeVariables &elemVolVars) { K_.mv(potentialGrad_, Kmvp_); KmvpNormal_ = -(Kmvp_*face().normal); @@ -437,21 +479,21 @@ protected: * * \param problem The considered problem file * \param element The considered element of the grid - * \param elemDat The parameters stored in the considered element + * \param elemVolVars The parameters stored in the considered element */ - void calculateDiffCoeffPM_(const Problem &problem, - const Element &element, - const ElementVolumeVariables &elemDat) + void calculatePorousDiffCoeff_(const Problem &problem, + const Element &element, + const ElementVolumeVariables &elemVolVars) { - const VolumeVariables &vDat_i = elemDat[face().i]; - const VolumeVariables &vDat_j = elemDat[face().j]; + const VolumeVariables &volVarsI = elemVolVars[face().i]; + const VolumeVariables &volVarsJ = elemVolVars[face().j]; // Diffusion coefficient in the porous medium - diffCoeffPM_ - = harmonicMean(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff(), - vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff()); -// = 1./2*(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff() + -// vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff()); + porousDiffCoeff_ + = harmonicMean(volVarsI.porosity() * volVarsI.tortuosity() * volVarsI.diffCoeff(), + volVarsJ.porosity() * volVarsJ.tortuosity() * volVarsJ.diffCoeff()); +// = 1./2*(volVarsI.porosity() * volVarsI.tortuosity() * volVarsI.diffCoeff() + +// volVarsJ.porosity() * volVarsJ.tortuosity() * volVarsJ.diffCoeff()); } /*! @@ -459,26 +501,26 @@ protected: * * \param problem The considered problem file * \param element The considered element of the grid - * \param elemDat The parameters stored in the considered element + * \param elemVolVars The parameters stored in the considered element */ void calculateDispersionTensor_(const Problem &problem, const Element &element, - const ElementVolumeVariables &elemDat) + const ElementVolumeVariables &elemVolVars) { - const VolumeVariables &vDat_i = elemDat[face().i]; - const VolumeVariables &vDat_j = elemDat[face().j]; + const VolumeVariables &volVarsI = elemVolVars[face().i]; + const VolumeVariables &volVarsJ = elemVolVars[face().j]; //calculate dispersivity at the interface: [0]: alphaL = longitudinal disp. [m], [1] alphaT = transverse disp. [m] Scalar dispersivity[2]; - dispersivity[0] = 0.5 * (vDat_i.dispersivity()[0] + vDat_j.dispersivity()[0]); - dispersivity[1] = 0.5 * (vDat_i.dispersivity()[1] + vDat_j.dispersivity()[1]); + dispersivity[0] = 0.5 * (volVarsI.dispersivity()[0] + volVarsJ.dispersivity()[0]); + dispersivity[1] = 0.5 * (volVarsI.dispersivity()[1] + volVarsJ.dispersivity()[1]); //calculate velocity at interface: v = -1/mu * vDarcy = -1/mu * K * grad(p) - Vector velocity; + DimVector velocity; Valgrind::CheckDefined(potentialGrad()); Valgrind::CheckDefined(K_); K_.mv(potentialGrad(), velocity); - velocity /= - 0.5 * (vDat_i.viscosity() + vDat_j.viscosity()); + velocity /= - 0.5 * (volVarsI.viscosity() + volVarsJ.viscosity()); //matrix multiplication of the velocity at the interface: vv^T dispersionTensor_ = 0; @@ -501,28 +543,28 @@ protected: dispersionTensor_[i][i] += vNorm*dispersivity[1]; } - const FVElementGeometry &fvGeom_; + const FVElementGeometry &fvGeometry_; const int faceIdx_; const bool onBoundary_; //! pressure potential gradient - Vector potentialGrad_; - //! concentratrion gradient - Vector concentrationGrad_; + DimVector potentialGrad_; + //! DEPRECATED concentration gradient + DimVector concentrationGrad_; //! mole-fraction gradient - Vector moleFracGrad_; - //! mass-fraction gradient - Vector massFracGrad_; + DimVector moleFractionGrad_; + //! DEPRECATED mass-fraction gradient + DimVector massFractionGrad_; //! the effective diffusion coefficent in the porous medium - Scalar diffCoeffPM_; + Scalar porousDiffCoeff_; //! the dispersion tensor in the porous medium - Tensor dispersionTensor_; + DimMatrix dispersionTensor_; //! the intrinsic permeability tensor - Tensor K_; + DimMatrix K_; // intrinsic permeability times pressure potential gradient - Vector Kmvp_; + DimVector Kmvp_; // projected on the face normal Scalar KmvpNormal_; @@ -532,10 +574,10 @@ protected: int downstreamIdx_; //! viscosity of the fluid at the integration point - Scalar viscosityAtIP_; + Scalar viscosity_; //! molar densities of the fluid at the integration point - Scalar molarDensityAtIP_, densityAtIP_; + Scalar molarDensity_, density_; }; } // end namepace diff --git a/dumux/boxmodels/1p2c/1p2cindices.hh b/dumux/boxmodels/1p2c/1p2cindices.hh index 2055fb3700ae854ca3ada77ea8e3ac1566659658..bbe796d812202ad85a386de36ce0dd2f24cad219 100644 --- a/dumux/boxmodels/1p2c/1p2cindices.hh +++ b/dumux/boxmodels/1p2c/1p2cindices.hh @@ -58,7 +58,7 @@ struct OnePTwoCIndices // primary variable indices static const int pressureIdx = PVOffset + 0; //!< pressure - static const int x1Idx = PVOffset + 1; //!< mole fraction of the second component + static const int massOrMoleFractionIdx = PVOffset + 1; //!< mole fraction of the second component }; // \} diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh index 0db5889660b93510c3f71cd67f0997b11a1bf891..9dbf104ec8f084edb3f4826717e6c0ce4d19a130 100644 --- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh +++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh @@ -109,11 +109,11 @@ public: * \brief Evaluate the amount of all conservation quantities * (e.g. phase mass) within a finite volume. * - * \param result The mass of the component within the sub-control volume + * \param storage The mass of the component within 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 + void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const { // if flag usePrevSol is set, the solution from the previous // time step is used, otherwise the current solution is @@ -123,23 +123,23 @@ public: const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); const VolumeVariables &volVars = elemVolVars[scvIdx]; - result = 0; + storage = 0; if(!useMoles) { // storage term of continuity equation - massfractions - result[contiEqIdx] += + storage[contiEqIdx] += volVars.fluidState().density(phaseIdx)*volVars.porosity(); //storage term of the transport equation - massfractions - result[transEqIdx] += + storage[transEqIdx] += volVars.fluidState().density(phaseIdx) * volVars.fluidState().massFraction(phaseIdx, comp1Idx) * volVars.porosity(); } else { // storage term of continuity equation- molefractions //careful: molarDensity changes with moleFrac! - result[contiEqIdx] += volVars.molarDensity()*volVars.porosity(); + storage[contiEqIdx] += volVars.molarDensity()*volVars.porosity(); // storage term of the transport equation - molefractions - result[transEqIdx] += + storage[transEqIdx] += volVars.fluidState().molarDensity(phaseIdx)*volVars.fluidState().moleFraction(phaseIdx, comp1Idx) * volVars.porosity(); } @@ -151,17 +151,17 @@ public: * 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 faceIdx The index of the considered face of the sub control volume * \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 + void computeFlux(PrimaryVariables &flux, const int faceIdx, const bool onBoundary=false) const { flux = 0; FluxVariables fluxVars(this->problem_(), - this->elem_(), - this->fvElemGeom_(), - faceId, + this->element_(), + this->fvGeometry_(), + faceIdx, this->curVolVars_(), onBoundary); @@ -241,22 +241,22 @@ public: if(!useMoles) { // diffusive flux of the second component - massfraction - tmp = -(fluxVars.moleFracGrad(comp1Idx)*fluxVars.face().normal); - tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensityAtIP(); + tmp = -(fluxVars.moleFractionGrad(comp1Idx)*fluxVars.face().normal); + tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensity(); // convert it to a mass flux and add it flux[transEqIdx] += tmp * FluidSystem::molarMass(comp1Idx); } else { // diffusive flux of the second component - molefraction - tmp = -(fluxVars.moleFracGrad(comp1Idx)*fluxVars.face().normal); - tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensityAtIP(); + tmp = -(fluxVars.moleFractionGrad(comp1Idx)*fluxVars.face().normal); + tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensity(); // dispersive flux of second component - molefraction // Vector normalDisp; // fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp); - // tmp -= fluxVars.molarDensityAtIP()* - // (normalDisp * fluxVars.moleFracGrad(comp1Idx)); + // tmp -= fluxVars.molarDensity()* + // (normalDisp * fluxVars.moleFractionGrad(comp1Idx)); flux[transEqIdx] += tmp; } @@ -264,16 +264,16 @@ public: /*! * \brief Calculate the source term of the equation - * \param q The source/sink in the SCV for each component - * \param localVertexIdx The index of the vertex of the sub control volume + * \param source The source/sink in the SCV for each component + * \param scvIdx The index of the vertex of the sub control volume * */ - void computeSource(PrimaryVariables &q, int localVertexIdx) + void computeSource(PrimaryVariables &source, const int scvIdx) { - this->problem_().boxSDSource(q, - this->elem_(), - this->fvElemGeom_(), - localVertexIdx, + this->problem_().boxSDSource(source, + this->element_(), + this->fvGeometry_(), + scvIdx, this->curVolVars_()); } @@ -286,8 +286,8 @@ public: * \param boundaryFaceIdx The index of the considered boundary face of the sub control volume */ void evalOutflowSegment(const IntersectionIterator &isIt, - int scvIdx, - int boundaryFaceIdx) + const int scvIdx, + const int boundaryFaceIdx) { const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx); // deal with outflow boundaries diff --git a/dumux/boxmodels/1p2c/1p2cmodel.hh b/dumux/boxmodels/1p2c/1p2cmodel.hh index 2a1e4ed5ec14adf1b15eff93032a46f4150db507..38e8ac6747efa81619ea6ff1f4cf0c1efbf4e0af 100644 --- a/dumux/boxmodels/1p2c/1p2cmodel.hh +++ b/dumux/boxmodels/1p2c/1p2cmodel.hh @@ -94,7 +94,7 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag> typedef typename GridView::template Codim<dim>::Iterator VertexIterator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: /*! @@ -124,18 +124,18 @@ public: unsigned numVertices = this->problem_().gridView().size(dim); ScalarField &pressure = *writer.allocateManagedBuffer(numVertices); ScalarField &delp = *writer.allocateManagedBuffer(numVertices); - ScalarField &moleFrac0 = *writer.allocateManagedBuffer(numVertices); - ScalarField &moleFrac1 = *writer.allocateManagedBuffer(numVertices); - ScalarField &massFrac0 = *writer.allocateManagedBuffer(numVertices); - ScalarField &massFrac1 = *writer.allocateManagedBuffer(numVertices); + ScalarField &moleFraction0 = *writer.allocateManagedBuffer(numVertices); + ScalarField &moleFraction1 = *writer.allocateManagedBuffer(numVertices); + ScalarField &massFraction0 = *writer.allocateManagedBuffer(numVertices); + ScalarField &massFraction1 = *writer.allocateManagedBuffer(numVertices); ScalarField &rho = *writer.allocateManagedBuffer(numVertices); ScalarField &mu = *writer.allocateManagedBuffer(numVertices); #ifdef VELOCITY_OUTPUT // check if velocity output is demanded ScalarField &velocityX = *writer.allocateManagedBuffer(numVertices); ScalarField &velocityY = *writer.allocateManagedBuffer(numVertices); ScalarField &velocityZ = *writer.allocateManagedBuffer(numVertices); - //use vertiacl faces for vx and horizontal faces for vy calculation - std::vector<GlobalPosition> boxSurface(numVertices); + //use vertical faces for vx and horizontal faces for vy calculation + std::vector<DimVector> boxSurface(numVertices); // initialize velocity fields for (int i = 0; i < numVertices; ++i) { @@ -156,7 +156,7 @@ public: ScalarField &rank = *writer.allocateManagedBuffer(numElements); - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeometry; VolumeVariables volVars; ElementBoundaryTypes elemBcTypes; @@ -167,8 +167,8 @@ public: int idx = this->problem_().model().elementMapper().map(*elemIt); rank[idx] = this->gridView_().comm().rank(); - fvElemGeom.update(this->gridView_(), *elemIt); - elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom); + fvGeometry.update(this->gridView_(), *elemIt); + elemBcTypes.update(this->problem_(), *elemIt, fvGeometry); int numVerts = elemIt->template count<dim> (); for (int i = 0; i < numVerts; ++i) @@ -177,16 +177,16 @@ public: volVars.update(sol[globalIdx], this->problem_(), *elemIt, - fvElemGeom, + fvGeometry, i, false); pressure[globalIdx] = volVars.pressure(); delp[globalIdx] = volVars.pressure() - 1e5; - moleFrac0[globalIdx] = volVars.moleFraction(0); - moleFrac1[globalIdx] = volVars.moleFraction(1); - massFrac0[globalIdx] = volVars.massFraction(0); - massFrac1[globalIdx] = volVars.massFraction(1); + moleFraction0[globalIdx] = volVars.moleFraction(0); + moleFraction1[globalIdx] = volVars.moleFraction(1); + massFraction0[globalIdx] = volVars.massFraction(0); + massFraction1[globalIdx] = volVars.massFraction(1); rho[globalIdx] = volVars.density(); mu[globalIdx] = volVars.viscosity(); }; @@ -194,29 +194,29 @@ public: #ifdef VELOCITY_OUTPUT // check if velocity output is demanded // In the box method, the velocity is evaluated on the FE-Grid. However, to get an // average apparent velocity at the vertex, all contributing velocities have to be interpolated. - GlobalPosition velocity; + DimVector velocity; ElementVolumeVariables elemVolVars; elemVolVars.update(this->problem_(), *elemIt, - fvElemGeom, + fvGeometry, false /* isOldSol? */); // loop over the phases - for (int faceIdx = 0; faceIdx < fvElemGeom.numEdges; faceIdx++) + for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++) { velocity = 0.0; //prepare the flux calculations (set up and prepare geometry, FE gradients) FluxVariables fluxVars(this->problem_(), *elemIt, - fvElemGeom, + fvGeometry, faceIdx, elemVolVars); - //use vertiacl faces for vx and horizontal faces for vy calculation - GlobalPosition xVector(0), yVector(0); + //use vertical faces for vx and horizontal faces for vy calculation + DimVector xVector(0), yVector(0); xVector[0] = 1; yVector[1] = 1; - Dune::SeqScalarProduct<GlobalPosition> sp; + Dune::SeqScalarProduct<DimVector> sp; Scalar xDir = std::abs(sp.dot(fluxVars.face().normal, xVector)); Scalar yDir = std::abs(sp.dot(fluxVars.face().normal, yVector)); @@ -286,7 +286,7 @@ public: { int i = this->problem_().vertexMapper().map(*vIt); - //use vertiacl faces for vx and horizontal faces for vy calculation + //use vertical faces for vx and horizontal faces for vy calculation velocityX[i] /= boxSurface[i][0]; if (dim >= 2) { @@ -306,18 +306,17 @@ public: if (dim > 2) writer.attachVertexData(velocityZ, "Vz"); #endif - char nameMoleFrac0[42], nameMoleFrac1[42]; - snprintf(nameMoleFrac0, 42, "x_%s", FluidSystem::componentName(0)); - snprintf(nameMoleFrac1, 42, "x_%s", FluidSystem::componentName(1)); - writer.attachVertexData(moleFrac0, nameMoleFrac0); - writer.attachVertexData(moleFrac1, nameMoleFrac1); - - char nameMassFrac0[42], nameMassFrac1[42]; - snprintf(nameMassFrac0, 42, "X_%s", FluidSystem::componentName(0)); - snprintf(nameMassFrac1, 42, "X_%s", FluidSystem::componentName(1)); - writer.attachVertexData(massFrac0, nameMassFrac0); - writer.attachVertexData(massFrac1, nameMassFrac1); -// writer.attachVertexData(delFrac, "delFrac_TRAIL"); + char nameMoleFraction0[42], nameMoleFraction1[42]; + snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0)); + snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1)); + writer.attachVertexData(moleFraction0, nameMoleFraction0); + writer.attachVertexData(moleFraction1, nameMoleFraction1); + + char nameMassFraction0[42], nameMassFraction1[42]; + snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0)); + snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1)); + writer.attachVertexData(massFraction0, nameMassFraction0); + writer.attachVertexData(massFraction1, nameMassFraction1); writer.attachVertexData(rho, "rho"); writer.attachVertexData(mu, "mu"); writer.attachCellData(rank, "process rank"); diff --git a/dumux/boxmodels/1p2c/1p2cproblem.hh b/dumux/boxmodels/1p2c/1p2cproblem.hh index 46be5c438a8112ba86eb944c45ced14de86290c9..1e7f73ecdee5b7b7c96f733a4eb83ee39d49e916 100644 --- a/dumux/boxmodels/1p2c/1p2cproblem.hh +++ b/dumux/boxmodels/1p2c/1p2cproblem.hh @@ -49,7 +49,7 @@ class OnePTwoCBoxProblem : public PorousMediaBoxProblem<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; public: /*! @@ -62,7 +62,7 @@ public: DUMUX_DEPRECATED_MSG("use PorousMediaBoxProblem instead") OnePTwoCBoxProblem(TimeManager &timeManager, const GridView &gridView, - bool verbose = true) + const bool verbose = true) : ParentType(timeManager, gridView) {} @@ -71,14 +71,14 @@ public: * * \param timeManager The time manager * \param gridView The grid view - * \param spatialParameters The spatial parameters object + * \param spatialParams The spatial parameters object * \param verbose Turn verbosity on or off */ DUMUX_DEPRECATED_MSG("use PorousMediaBoxProblem instead") OnePTwoCBoxProblem(TimeManager &timeManager, const GridView &gridView, - SpatialParameters &spatialParameters, - bool verbose = true) + SpatialParams &spatialParams, + const bool verbose = true) : ParentType(timeManager, gridView) {} diff --git a/dumux/boxmodels/1p2c/1p2cproperties.hh b/dumux/boxmodels/1p2c/1p2cproperties.hh index 80d682f631d468851bbb517f197958744779625f..3ce9b8f9807cad1123eda6cb3cc459ff5fea8dbd 100644 --- a/dumux/boxmodels/1p2c/1p2cproperties.hh +++ b/dumux/boxmodels/1p2c/1p2cproperties.hh @@ -58,7 +58,8 @@ NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system NEW_PROP_TAG(OnePTwoCIndices); //!< DEPRECATED Enumerations for the 1p2c models NEW_PROP_TAG(Indices); //!< Enumerations for the model -NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters +NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters +NEW_PROP_TAG(SpatialParameters); //!< DEPRECATED The type of the spatial parameters NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations NEW_PROP_TAG(UpwindWeight); //!< The default value of the upwind weight NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem diff --git a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh index 037c9a08b67702f8e2d6cabe9e449c430eaacccd..7cb25b070197e30204b3dedb23dfc09cc8691a2b 100644 --- a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh +++ b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh @@ -62,15 +62,15 @@ class OnePTwoCVolumeVariables : public BoxVolumeVariables<TypeTag> comp1Idx = Indices::comp1Idx, pressureIdx = Indices::pressureIdx, - x1Idx = Indices::x1Idx + massOrMoleFractionIdx = Indices::massOrMoleFractionIdx }; 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::FieldVector<Scalar,dim> DimVector; public: //! The type returned by the fluidState() method @@ -82,25 +82,25 @@ public: * \param priVars A vector containing the primary variables * \param problem The considered problem * \param element The considered element of the grid - * \param elemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param scvIdx The index of the considered sub-control volume * \param isOldSol Evaluate function with solution of current or previous time step */ void update(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int scvIdx, - bool isOldSol) + const FVElementGeometry &fvGeometry, + const int scvIdx, + const bool isOldSol) { - ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol); + ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol); //calculate all secondary variables from the primary variables and store results in fluidstate - completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_); + completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_); - porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx); - tortuosity_ = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx); - dispersivity_ = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx); + porosity_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx); + tortuosity_ = problem.spatialParams().tortuosity(element, fvGeometry, scvIdx); + dispersivity_ = problem.spatialParams().dispersivity(element, fvGeometry, scvIdx); // Second instance of a parameter cache. // Could be avoided if diffusion coefficients also @@ -120,7 +120,7 @@ public: Valgrind::CheckDefined(diffCoeff_); // energy related quantities not contained in the fluid state - asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol); + asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol); } /*! @@ -130,16 +130,16 @@ public: const Problem& problem, const Element& element, const FVElementGeometry& elementGeometry, - int scvIdx, + const int scvIdx, FluidState& fluidState) { - Scalar t = Implementation::temperature_(primaryVariables, problem, element, + Scalar T = Implementation::temperature_(primaryVariables, problem, element, elementGeometry, scvIdx); - fluidState.setTemperature(t); + fluidState.setTemperature(T); fluidState.setPressure(phaseIdx, primaryVariables[pressureIdx]); - Scalar x1 = primaryVariables[x1Idx]; //mole or mass fraction of component 1 + Scalar x1 = primaryVariables[massOrMoleFractionIdx]; //mole or mass fraction of component 1 if(!useMoles) //mass-fraction formulation { // convert mass to mole fractions @@ -228,7 +228,7 @@ public: /*! * \brief Returns the dispersivity of the fluid's streamlines. */ - const Vector &dispersivity() const + const DimVector &dispersivity() const { return dispersivity_; } /*! @@ -256,28 +256,28 @@ public: protected: static Scalar temperature_(const PrimaryVariables &priVars, - const Problem& problem, - const Element &element, - const FVElementGeometry &elemGeom, - int scvIdx) + const Problem& problem, + const Element &element, + const FVElementGeometry &fvGeometry, + const int scvIdx) { - return problem.boxTemperature(element, elemGeom, scvIdx); + return problem.boxTemperature(element, fvGeometry, scvIdx); } /*! * \brief Called by update() to compute the energy related quantities. */ - void updateEnergy_(const PrimaryVariables &sol, + void updateEnergy_(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int vertIdx, - bool isOldSol) + const FVElementGeometry &fvGeometry, + const int scvIdx, + const bool isOldSol) { } Scalar porosity_; //!< Effective porosity within the control volume Scalar tortuosity_; - Vector dispersivity_; + DimVector dispersivity_; Scalar diffCoeff_; FluidState fluidState_; diff --git a/test/boxmodels/1p2c/1p2coutflowproblem.hh b/test/boxmodels/1p2c/1p2coutflowproblem.hh index a75cd668e1e5fc3f7e817251ef44949a67a32533..600422eff68656abe9c5b023a6a8c7d5a1b49166 100644 --- a/test/boxmodels/1p2c/1p2coutflowproblem.hh +++ b/test/boxmodels/1p2c/1p2coutflowproblem.hh @@ -135,7 +135,7 @@ class OnePTwoCOutflowProblem : public PorousMediaBoxProblem<TypeTag> // indices of the primary variables pressureIdx = Indices::pressureIdx, - x1Idx = Indices::x1Idx, + massOrMoleFractionIdx = Indices::massOrMoleFractionIdx, // indices of the equations contiEqIdx = Indices::contiEqIdx, @@ -206,7 +206,7 @@ public: } /*! - * \brief Evaluate the boundary conditions for a dirichlet + * \brief Evaluate the boundary conditions for a Dirichlet * boundary segment. * * For this method, the \a values parameter stores primary variables. @@ -216,13 +216,13 @@ public: const GlobalPosition globalPos = vertex.geometry().center(); initial_(values, globalPos); - //condition for the trail molefraction at left boundary + //condition for the N2 molefraction at left boundary if(globalPos[0] < eps_) - values[x1Idx] = 2.0e-5; + values[massOrMoleFractionIdx] = 2.0e-5; } /*! - * \brief Evaluate the boundary conditions for a neumann + * \brief Evaluate the boundary conditions for a Neumann * boundary segment. * * For this method, the \a values parameter stores the mass flux @@ -231,10 +231,10 @@ public: */ void neumann(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, const Intersection &is, - int scvIdx, - int boundaryFaceIdx) const + const int scvIdx, + const int boundaryFaceIdx) const { //const GlobalPosition &globalPos = element.geometry().corner(scvIdx); values = 0; @@ -270,8 +270,8 @@ public: */ void initial(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + const FVElementGeometry &fvGeometry, + const int scvIdx) const { const GlobalPosition &globalPos = element.geometry().corner(scvIdx); @@ -287,7 +287,7 @@ private: const GlobalPosition &globalPos) const { values[pressureIdx] = 2e5 - 1e5*globalPos[0];//0.0; //initial condition for the pressure - values[x1Idx] = 0.0; //initial condition for the trail molefraction + values[massOrMoleFractionIdx] = 0.0; //initial condition for the N2 molefraction } const Scalar eps_; diff --git a/test/boxmodels/1p2c/1p2coutflowspatialparameters.hh b/test/boxmodels/1p2c/1p2coutflowspatialparameters.hh index 370296fac78e80a7dc00ff143838bddd65efc984..d193bace8198a22ae72934b6f86077c22bf11d52 100644 --- a/test/boxmodels/1p2c/1p2coutflowspatialparameters.hh +++ b/test/boxmodels/1p2c/1p2coutflowspatialparameters.hh @@ -68,8 +68,8 @@ class OnePTwoCOutflowSpatialParameters : public BoxSpatialParametersOneP<TypeTag //typedef LinearMaterial<Scalar> EffMaterialLaw; public: - OnePTwoCOutflowSpatialParameters(const GridView &gv) - : ParentType(gv) + OnePTwoCOutflowSpatialParameters(const GridView &gridView) + : ParentType(gridView) { permeability_ = 1e-10; porosity_ = 0.4; @@ -94,12 +94,12 @@ public: * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$. * * \param element The current finite element - * \param fvElemGeom The current finite volume geometry of the element + * \param fvGeometry The current finite volume geometry of the element * \param scvIdx The index of the sub-control volume */ const Scalar intrinsicPermeability(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + const FVElementGeometry &fvGeometry, + const int scvIdx) const { return permeability_; } @@ -108,12 +108,12 @@ public: * \brief Define the porosity \f$\mathrm{[-]}\f$. * * \param element The finite element - * \param fvElemGeom The finite volume geometry + * \param fvGeometry The finite volume geometry * \param scvIdx The local index of the sub-control volume where */ double porosity(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + const FVElementGeometry &fvGeometry, + const int scvIdx) const { return porosity_; } @@ -122,12 +122,12 @@ public: * \brief Define the tortuosity \f$\mathrm{[-]}\f$. * * \param element The finite element - * \param fvElemGeom The finite volume geometry + * \param fvGeometry The finite volume geometry * \param scvIdx The local index of the sub-control volume where */ double tortuosity(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + const FVElementGeometry &fvGeometry, + const int scvIdx) const { return tortuosity_; } @@ -136,19 +136,19 @@ public: * \brief Define the dispersivity. * * \param element The finite element - * \param fvElemGeom The finite volume geometry + * \param fvGeometry The finite volume geometry * \param scvIdx The local index of the sub-control volume where */ double dispersivity(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + const FVElementGeometry &fvGeometry, + const int scvIdx) const { return 0; } - bool useTwoPointGradient(const Element &elem, - int vertexI, - int vertexJ) const + bool useTwoPointGradient(const Element &element, + const int vertexI, + const int vertexJ) const { return false; }