diff --git a/dumux/decoupled/1p/1pproperties.hh b/dumux/decoupled/1p/1pproperties.hh index 447f95ebc8c50673ea40435a77386d48bbc1d539..5b95ee51ce1dab8a7dc08ec0199e06cf3eee5bbf 100644 --- a/dumux/decoupled/1p/1pproperties.hh +++ b/dumux/decoupled/1p/1pproperties.hh @@ -63,6 +63,7 @@ NEW_TYPE_TAG(DecoupledOneP, INHERITS_FROM(DecoupledModel)); ////////////////////////////////////////////////////////////////// NEW_PROP_TAG( SpatialParameters ); //!< The type of the spatial parameters object +NEW_PROP_TAG( SpatialParams ); //!< The type of the spatial parameters object NEW_PROP_TAG( EnableGravity); //!< Returns whether gravity is considered in the problem NEW_PROP_TAG( Fluid ); //!< The fluid for one-phase models NEW_PROP_TAG( Indices ); //!< Set of indices for the one-phase model @@ -100,6 +101,9 @@ SET_TYPE_PROP(DecoupledOneP, Variables, VariableClass<TypeTag>); //! Set standart CellData of immiscible one-phase models as default SET_TYPE_PROP(DecoupledOneP, CellData, CellData1P<TypeTag>); + +//! DEPRECATED SpatialParameters property +SET_TYPE_PROP(DecoupledOneP, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams)); } } #endif diff --git a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh index 4fa32eb771267b34b7f79929f19de8652f487a0b..bb5a970c5ea7332e815a7e29bfde781bf0ba2c65 100644 --- a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh +++ b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh @@ -61,7 +61,7 @@ class DiffusionProblem1P: public OneModelProblem<TypeTag> // material properties typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GridView::Traits::template Codim<0>::Entity Element; @@ -82,7 +82,7 @@ public: DiffusionProblem1P(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView), gravity_(0) { - spatialParameters_ = new SpatialParameters(gridView); + spatialParams_ = new SpatialParams(gridView); newSpatialParams_ = true; gravity_ = 0; if (GET_PARAM(TypeTag, bool, EnableGravity)) @@ -93,10 +93,10 @@ public: * * \param timeManager the time manager * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ - DiffusionProblem1P(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters) - : ParentType(timeManager, gridView), gravity_(0), spatialParameters_(&spatialParameters) + DiffusionProblem1P(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams) + : ParentType(timeManager, gridView), gravity_(0), spatialParams_(&spatialParams) { newSpatialParams_ = false; gravity_ = 0; @@ -111,7 +111,7 @@ public: DiffusionProblem1P(const GridView &gridView) : ParentType(gridView, false), gravity_(0) { - spatialParameters_ = new SpatialParameters(gridView); + spatialParams_ = new SpatialParams(gridView); newSpatialParams_ = true; gravity_ = 0; if (GET_PARAM(TypeTag, bool, EnableGravity)) @@ -121,10 +121,10 @@ public: * \brief The constructor * * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ - DiffusionProblem1P(const GridView &gridView, SpatialParameters &spatialParameters) - : ParentType(gridView, false), gravity_(0), spatialParameters_(&spatialParameters) + DiffusionProblem1P(const GridView &gridView, SpatialParams &spatialParams) + : ParentType(gridView, false), gravity_(0), spatialParams_(&spatialParams) { newSpatialParams_ = false; gravity_ = 0; @@ -136,7 +136,7 @@ public: { if (newSpatialParams_) { - delete spatialParameters_; + delete spatialParams_; } } @@ -229,19 +229,27 @@ public: /*! * \brief Returns the spatial parameters object. */ - SpatialParameters &spatialParameters() + SpatialParams &spatialParams() { - return *spatialParameters_; + return *spatialParams_; } + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + SpatialParams &spatialParameters() + { return *spatialParams_; } + /*! * \brief Returns the spatial parameters object. */ - const SpatialParameters &spatialParameters() const + const SpatialParams &spatialParams() const { - return *spatialParameters_; + return *spatialParams_; } + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + const SpatialParams &spatialParameters() const + { return *spatialParams_; } + // \} private: @@ -256,7 +264,7 @@ private: GlobalPosition gravity_; // fluids and material properties - SpatialParameters* spatialParameters_; + SpatialParams* spatialParams_; bool newSpatialParams_; }; diff --git a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh index 196beab4b6cd7bcecabfc1aed056881f40aaad26..7abaebbc51f480ca184ebfe43994ecebb11a046e 100644 --- a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh +++ b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh @@ -65,7 +65,7 @@ template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid; @@ -98,7 +98,7 @@ template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag> typedef typename GridView::Intersection Intersection; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; public: @@ -271,10 +271,10 @@ void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entry, const I Scalar dist = distVec.two_norm(); // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI), - problem_.spatialParameters().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI), + problem_.spatialParams().intrinsicPermeability(*elementJ)); Dune::FieldVector<Scalar, dim> permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -331,10 +331,10 @@ const Intersection& intersection, const CellData& cellData, const bool first) //permeability vector at boundary // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element)); Dune::FieldVector<Scalar, dim> permeability(0); meanPermeability.mv(unitOuterNormal, permeability); diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh index 44c7b0b2d1b1733e2c7b0ebf7b597f9031eee109..928f04308f3630bf2d21088e7326153b618a3b8b 100644 --- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh +++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh @@ -55,7 +55,7 @@ class FVVelocity1P typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; @@ -80,7 +80,7 @@ typedef typename GridView::IntersectionIterator IntersectionIterator; }; typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix; public: //! Constructs a FVVelocity1P object @@ -150,8 +150,8 @@ public: 0); // get the transposed Jacobian of the element mapping - const FieldMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos); - FieldMatrix jacobianT(jacobianInv); + const DimMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos); + DimMatrix jacobianT(jacobianInv); jacobianT.invert(); // calculate the element velocity by the Piola transformation @@ -208,10 +208,10 @@ void FVVelocity1P<TypeTag>::calculateVelocity(const Intersection& intersection, Scalar dist = distVec.two_norm(); // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI), - problem_.spatialParameters().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI), + problem_.spatialParams().intrinsicPermeability(*elementJ)); Dune::FieldVector < Scalar, dim > permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -286,9 +286,9 @@ void FVVelocity1P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte //permeability vector at boundary // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element)); //multiply with normal vector at the boundary Dune::FieldVector < Scalar, dim > permeability(0); diff --git a/dumux/decoupled/1p/fluxData1p.hh b/dumux/decoupled/1p/fluxData1p.hh index 60a820bf0c371eb82789ec3ae3a8078f550127c6..ca75c55b48f00a66887578d2a79c26b2746dff4d 100644 --- a/dumux/decoupled/1p/fluxData1p.hh +++ b/dumux/decoupled/1p/fluxData1p.hh @@ -54,8 +54,8 @@ private: dim = GridView::dimension, dimWorld = GridView::dimensionworld }; - typedef Dune::FieldVector<Scalar, dim> FieldVector; - typedef Dune::FieldVector<FieldVector, 2 * dim> VelocityVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; + typedef Dune::FieldVector<DimVector, 2 * dim> VelocityVector; VelocityVector velocity_; Scalar potential_[2 * dim]; @@ -68,7 +68,7 @@ public: { for (int face = 0; face < 2*dim; face++) { - velocity_[face] = FieldVector(0.0); + velocity_[face] = DimVector(0.0); potential_[face] = 0.0; velocityMarker_[face] = false; } @@ -82,7 +82,7 @@ public: * * \param indexInInside Index of the cell-cell interface in this cell */ - const FieldVector& velocity(int indexInInside) + const DimVector& velocity(int indexInInside) { return velocity_[indexInInside]; } @@ -90,7 +90,7 @@ public: * * \param indexInInside Index of the cell-cell interface in this cell */ - const FieldVector& velocity(int indexInInside) const + const DimVector& velocity(int indexInInside) const { return velocity_[indexInInside]; } @@ -99,7 +99,7 @@ public: * \param indexInInside Index of the cell-cell interface in this cell * \param velocity Velocity vector which is stored */ - void setVelocity(int indexInInside, FieldVector& velocity) + void setVelocity(int indexInInside, DimVector& velocity) { velocity_[indexInInside] = velocity; } diff --git a/dumux/decoupled/2p/2pproperties.hh b/dumux/decoupled/2p/2pproperties.hh index 1847b21915f47e63d3098a965a0839f251d41c52..67cf09a6bdf5c583694a50b7e798826f553e94c2 100644 --- a/dumux/decoupled/2p/2pproperties.hh +++ b/dumux/decoupled/2p/2pproperties.hh @@ -132,7 +132,7 @@ typedef DecoupledTwoPIndices<GET_PROP_VALUE(TypeTag, Formulation), 0> type; }; //! \cond \private -// keep only for compatibility with box models +//! DEPRECATED TwoPIndices property SET_TYPE_PROP(DecoupledTwoP, TwoPIndices, typename GET_PROP_TYPE(TypeTag, Indices)); //! \endcond @@ -173,12 +173,7 @@ public: typedef IsothermalImmiscibleFluidState<Scalar, FluidSystem> type; }; //! DEPRECATED SpatialParameters property -#warning Please use SpatialParams instead of SpatialParameters -//TODO: next line enables old Models using SpatialParameters to work -// with base class impesproblem2p. If models are adapted, -// l180 should be replaced by l181 to deprecate old problems using SpatialParameters. -SET_TYPE_PROP(DecoupledTwoP, SpatialParams, typename GET_PROP_TYPE(TypeTag, SpatialParameters)); -//SET_TYPE_PROP(DecoupledTwoP, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams)); +SET_TYPE_PROP(DecoupledTwoP, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams)); /*! * \brief Set the property for the material parameters by extracting diff --git a/dumux/decoupled/2p/diffusion/diffusionproblem2p.hh b/dumux/decoupled/2p/diffusion/diffusionproblem2p.hh index 636272a4b1df7dbf3c2b0be7d0b8370da09f6618..316ee8548c1bc6f2328ec69c6e1f2266be751e27 100644 --- a/dumux/decoupled/2p/diffusion/diffusionproblem2p.hh +++ b/dumux/decoupled/2p/diffusion/diffusionproblem2p.hh @@ -57,7 +57,7 @@ class DiffusionProblem2P: public OneModelProblem<TypeTag> // material properties typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GridView::Traits::template Codim<0>::Entity Element; @@ -82,7 +82,7 @@ public: DiffusionProblem2P(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView), gravity_(0) { - spatialParameters_ = new SpatialParameters(gridView); + spatialParams_ = new SpatialParams(gridView); newSpatialParams_ = true; gravity_ = 0; if (GET_PARAM(TypeTag, bool, EnableGravity)) @@ -93,10 +93,10 @@ public: * * \param timeManager the time manager * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ - DiffusionProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters) - : ParentType(timeManager, gridView), gravity_(0), spatialParameters_(&spatialParameters) + DiffusionProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams) + : ParentType(timeManager, gridView), gravity_(0), spatialParams_(&spatialParams) { newSpatialParams_ = false; gravity_ = 0; @@ -112,7 +112,7 @@ public: DiffusionProblem2P(const GridView &gridView) : ParentType(gridView, false), gravity_(0) { - spatialParameters_ = new SpatialParameters(gridView); + spatialParams_ = new SpatialParams(gridView); newSpatialParams_ = true; gravity_ = 0; if (GET_PARAM(TypeTag, bool, EnableGravity)) @@ -122,10 +122,10 @@ public: * \brief Constructs a DiffusionProblem2P object * * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ - DiffusionProblem2P(const GridView &gridView, SpatialParameters &spatialParameters) - : ParentType(gridView, false), gravity_(0), spatialParameters_(&spatialParameters) + DiffusionProblem2P(const GridView &gridView, SpatialParams &spatialParams) + : ParentType(gridView, false), gravity_(0), spatialParams_(&spatialParams) { newSpatialParams_ = false; gravity_ = 0; @@ -138,7 +138,7 @@ public: { if (newSpatialParams_) { - delete spatialParameters_; + delete spatialParams_; } } @@ -223,19 +223,29 @@ public: /*! * \brief Returns the spatial parameters object. */ - SpatialParameters &spatialParameters() + SpatialParams &spatialParams() { - return *spatialParameters_; + return *spatialParams_; } + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + SpatialParams &spatialParameters() + { return *spatialParams_; } + /*! * \brief Returns the spatial parameters object. */ - const SpatialParameters &spatialParameters() const + const SpatialParams &spatialParams() const { - return *spatialParameters_; + return *spatialParams_; } + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + const SpatialParams &spatialParameters() const + { return *spatialParams_; } + + + // \} private: @@ -250,7 +260,7 @@ private: GlobalPosition gravity_; // fluids and material properties - SpatialParameters* spatialParameters_; + SpatialParams* spatialParams_; bool newSpatialParams_; }; diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh index 2b75c58351129b3a6a4e41335fef899590e32830..4cd710d4353e6c400cd31a231d6a7a4d69c28ed7 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh @@ -95,8 +95,8 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -137,7 +137,7 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag> typedef typename GridView::Intersection Intersection; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; protected: //! \cond \private @@ -494,7 +494,7 @@ void FVPressure2P<TypeTag>::getStorage(EntryType& entry, const Element& element // cell volume, assume linear map here Scalar volume = element.geometry().volume(); - Scalar porosity = problem_.spatialParameters().porosity(element); + Scalar porosity = problem_.spatialParams().porosity(element); switch (saturationType_) { @@ -592,10 +592,10 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters Scalar dist = distVec.two_norm(); // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI), - problem_.spatialParameters().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI), + problem_.spatialParams().intrinsicPermeability(*elementJ)); Dune::FieldVector<Scalar, dim> permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -743,10 +743,10 @@ const Intersection& intersection, const CellData& cellData, const bool first) //permeability vector at boundary // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element)); Dune::FieldVector<Scalar, dim> permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -783,7 +783,7 @@ const Intersection& intersection, const CellData& cellData, const bool first) Scalar pressBound = boundValues[pressureIdx]; //calculate consitutive relations depending on the kind of saturation used - Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*element), satW); + Scalar pcBound = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*element), satW); //determine phase pressures from primary pressure variable Scalar pressW = 0; @@ -829,9 +829,9 @@ const Intersection& intersection, const CellData& cellData, const bool first) rhoMeanNW = 0.5 * (cellData.density(nPhaseIdx) + densityNWBound); } - Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satW) + Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satW) / viscosityWBound; - Scalar lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satW) + Scalar lambdaNWBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satW) / viscosityNWBound; Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNWBound); @@ -981,7 +981,7 @@ void FVPressure2P<TypeTag>::updateMaterialLaws() Scalar satW = cellData.saturation(wPhaseIdx); - Scalar pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW); + Scalar pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt), satW); //determine phase pressures from primary pressure variable Scalar pressW = 0; @@ -1036,8 +1036,8 @@ void FVPressure2P<TypeTag>::updateMaterialLaws() } // initialize mobilities - Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; - Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; + Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; + Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; if (compressibility_) { diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh index 7c92cc2d277ed130c62ce6109324feb986f1c9b9..5335b9c0a0f8c47365524b8156fb414e164da963 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh @@ -290,12 +290,12 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection FieldMatrix permeabilityJ(0); FieldMatrix permeabilityK(0); - problem_.spatialParameters().meanK(permeabilityI, - problem_.spatialParameters().intrinsicPermeability(*elementI)); - problem_.spatialParameters().meanK(permeabilityJ, - problem_.spatialParameters().intrinsicPermeability(*elementJ)); - problem_.spatialParameters().meanK(permeabilityK, - problem_.spatialParameters().intrinsicPermeability(*elementK)); + problem_.spatialParams().meanK(permeabilityI, + problem_.spatialParams().intrinsicPermeability(*elementI)); + problem_.spatialParams().meanK(permeabilityJ, + problem_.spatialParams().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(permeabilityK, + problem_.spatialParams().intrinsicPermeability(*elementK)); // Calculate permeablity component normal to interface Scalar kI, kJ, kK, kMean, ng; diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh index a9982413e4f81dcae7bbc5dda3cb4fb295ae4ce9..716332d2f3f1d91e0f1581145a288b87d1fce583 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh @@ -61,8 +61,8 @@ class FVVelocity2P typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -108,7 +108,7 @@ class FVVelocity2P }; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; public: /*! \brief Constructs a FVVelocity2P object @@ -208,8 +208,8 @@ public: ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0); // get the transposed Jacobian of the element mapping - const FieldMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos); - FieldMatrix jacobianT(jacobianInv); + const DimMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos); + DimMatrix jacobianT(jacobianInv); jacobianT.invert(); // calculate the element velocity by the Piola transformation @@ -320,10 +320,10 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, Scalar dist = distVec.two_norm(); // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI), - problem_.spatialParameters().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI), + problem_.spatialParams().intrinsicPermeability(*elementJ)); Dune::FieldVector < Scalar, dim > permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -493,9 +493,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte //permeability vector at boundary // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element)); Dune::FieldVector < Scalar, dim > permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -528,7 +528,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte } Scalar pressBound = boundValues[pressureIdx]; - Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*element), satW); + Scalar pcBound = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*element), satW); //determine phase pressures from primary pressure variable Scalar pressWBound = 0; @@ -572,9 +572,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte viscosityNWBound = FluidSystem::viscosity(fluidState, nPhaseIdx) / densityNWBound; } - Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satW) + Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satW) / viscosityWBound; - Scalar lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satW) + Scalar lambdaNWBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satW) / viscosityNWBound; Scalar potentialW = 0; diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh index c30874a454125b5be8309309604e7371bf82322f..93312c189c9c74e9246c8801f4fb66b1adbc9cbf 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh @@ -245,12 +245,12 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters FieldMatrix permeabilityJ(0); FieldMatrix permeabilityK(0); - problem_.spatialParameters().meanK(permeabilityI, - problem_.spatialParameters().intrinsicPermeability(*elementI)); - problem_.spatialParameters().meanK(permeabilityJ, - problem_.spatialParameters().intrinsicPermeability(*elementJ)); - problem_.spatialParameters().meanK(permeabilityK, - problem_.spatialParameters().intrinsicPermeability(*elementK)); + problem_.spatialParams().meanK(permeabilityI, + problem_.spatialParams().intrinsicPermeability(*elementI)); + problem_.spatialParams().meanK(permeabilityJ, + problem_.spatialParams().intrinsicPermeability(*elementJ)); + problem_.spatialParams().meanK(permeabilityK, + problem_.spatialParams().intrinsicPermeability(*elementK)); // Calculate permeablity component normal to interface Scalar kI, kJ, kK, ng, kMean; //, gI, gJ, gK; diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh index 31828434e88c8343e362f889212ade1c5449d84a..72678fc84269180ed300f5c1bed0080fd6193e14 100644 --- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh +++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh @@ -70,8 +70,8 @@ class FVMPFAOPressure2P: public FVPressure<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -114,11 +114,11 @@ class FVMPFAOPressure2P: public FVPressure<TypeTag> typedef typename GridView::IntersectionIterator IntersectionIterator; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix; typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) Vector; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; //initializes the matrix to store the system of equations friend class FVPressure<TypeTag>; @@ -506,7 +506,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() this->f_=0; // introduce matrix R for vector rotation and R is initialized as zero matrix - FieldMatrix R(0); + DimMatrix R(0); // evaluate matrix R if (dim==2) @@ -543,7 +543,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() this->f_[globalIdx1] = volume1*(source[wPhaseIdx]/density_[wPhaseIdx] + source[nPhaseIdx]/density_[nPhaseIdx]); // get absolute permeability of cell 1 - FieldMatrix K1(problem_.spatialParameters().intrinsicPermeability(*eIt)); + DimMatrix K1(problem_.spatialParams().intrinsicPermeability(*eIt)); //compute total mobility of cell 1 Scalar lambda1 = 0; @@ -693,7 +693,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() globalPos2 = outside->geometry().center(); // get absolute permeability of neighbor cell 2 - FieldMatrix K2(problem_.spatialParameters().intrinsicPermeability(*outside)); + DimMatrix K2(problem_.spatialParams().intrinsicPermeability(*outside)); // get total mobility of neighbor cell 2 Scalar lambda2 = 0; @@ -718,7 +718,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = 0; @@ -726,7 +726,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() // neighbor cell 4 GlobalPosition globalPos4(0); - FieldMatrix K4(0); + DimMatrix K4(0); Scalar lambda4 = 0; int globalIdx4 = 0; @@ -754,7 +754,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() globalPos4 = innerisItoutside->geometry().center(); // get absolute permeability of neighbor cell 4 - K4 += problem_.spatialParameters().intrinsicPermeability(*innerisItoutside); + K4 += problem_.spatialParams().intrinsicPermeability(*innerisItoutside); // get total mobility of neighbor cell 4 lambda4 = cellData4.mobility(wPhaseIdx) + cellData4.mobility(nPhaseIdx); @@ -846,63 +846,63 @@ void FVMPFAOPressure2P<TypeTag>::assemble() *= face34vol/2.0; // compute normal vectors nu11,nu21; nu12, nu22; nu13, nu23; nu14, nu24; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24-globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12-globalPos2, nu22); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3-globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3-globalPosFace34, nu23); - FieldVector nu14(0); + DimVector nu14(0); R.umv(globalPos4-globalPosFace24, nu14); - FieldVector nu24(0); + DimVector nu24(0); R.umv(globalPosFace34-globalPos4, nu24); // compute dF1, dF2, dF3, dF4 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); - FieldVector Rnu24(0); + DimVector Rnu24(0); R.umv(nu24, Rnu24); Scalar dF4 = fabs(nu14 * Rnu24); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); - FieldVector K4nu14(0); + DimVector K4nu14(0); K4.umv(nu14, K4nu14); - FieldVector K4nu24(0); + DimVector K4nu24(0); K4.umv(nu24, K4nu24); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1042,35 +1042,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar J4 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24-globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12-globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1151,10 +1151,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda2 = lambdaWBound + lambdaNWBound; } @@ -1164,35 +1164,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() } // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24-globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12-globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1202,8 +1202,8 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g122 = alambda2 * (integrationOuterNormaln1 * K2nu22)/dF2; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g111 + g112; @@ -1222,7 +1222,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // assemble the global matrix this->A_ and right hand side f @@ -1267,10 +1267,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1287,35 +1287,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar J4 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24-globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12-globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1327,8 +1327,8 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g222 = lambda2 * (integrationOuterNormaln4 * K2nu22)/dF2; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g111 + g112; @@ -1347,7 +1347,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // assemble the global matrix this->A_ and right hand side this->f_ @@ -1391,10 +1391,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda2 = lambdaWBound + lambdaNWBound; } @@ -1404,35 +1404,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() } // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24-globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12-globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1442,8 +1442,8 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g122 = alambda2 * (integrationOuterNormaln1 * K2nu22)/dF2; // compute the matrix T & vector r - FieldMatrix T(0); - FieldVector r(0); + DimMatrix T(0); + DimVector r(0); Scalar coe = g111 + g112; @@ -1522,10 +1522,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1538,21 +1538,21 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g3 = boundValues[pressureIdx]; // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1584,7 +1584,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = 0; @@ -1639,35 +1639,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar J2 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3-globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3-globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1764,10 +1764,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda3 = lambdaWBound + lambdaNWBound; } @@ -1777,35 +1777,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() } // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3-globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3-globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -1815,7 +1815,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g223 = alambda3 * (integrationOuterNormaln3 * K3nu23)/dF3; // compute transmissibility matrix T = CA^{-1}B+F - FieldMatrix C(0), A(0), F(0), B(0); + DimMatrix C(0), A(0), F(0), B(0); // evaluate matrix C, F, A, B C[0][0] = -g111; @@ -1836,19 +1836,19 @@ void FVMPFAOPressure2P<TypeTag>::assemble() // compute T A.invert(); - FieldMatrix CAinv(C.rightmultiply(A)); + DimMatrix CAinv(C.rightmultiply(A)); F += B.leftmultiply(CAinv); - FieldMatrix T(F); + DimMatrix T(F); // compute vector r // evaluate r1, r2 - FieldVector r1(0), r2(0); + DimVector r1(0), r2(0); r1[1] = -g213 * g2; r2[0] = -J1 * face12vol/2.0; r2[1] = g213 * g2; // compute r = CA^{-1}r1 - FieldVector r(0); + DimVector r(0); CAinv.umv(r2, r); r += r1; @@ -1895,10 +1895,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1948,10 +1948,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1961,21 +1961,21 @@ void FVMPFAOPressure2P<TypeTag>::assemble() } // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); // compute dF1 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; @@ -2002,21 +2002,21 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar J3 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); // compute dF1 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -2051,7 +2051,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = 0; @@ -2132,10 +2132,10 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda3 = lambdaWBound + lambdaNWBound; } @@ -2145,35 +2145,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() } // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3-globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3-globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -2183,8 +2183,8 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g223 = alambda3 * (integrationOuterNormaln3 * K3nu23)/dF3; // compute the matrix T & vector r - FieldMatrix T(0); - FieldVector r(0); + DimMatrix T(0); + DimVector r(0); Scalar coe = g221 + g223; @@ -2212,35 +2212,35 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar J2 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13-globalPos1 ,nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1-globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3-globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3-globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11)/dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21)/dF1; @@ -2252,8 +2252,8 @@ void FVMPFAOPressure2P<TypeTag>::assemble() Scalar g223 = lambda3 * (integrationOuterNormaln3 * K3nu23)/dF3; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g113; @@ -2272,7 +2272,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // assemble the global matrix this->A_ and right hand side this->f_ @@ -2435,7 +2435,7 @@ void FVMPFAOPressure2P<TypeTag>::updateMaterialLaws() Scalar satW = cellData.saturation(wPhaseIdx); Scalar satNW = cellData.saturation(nPhaseIdx); - Scalar pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW); + Scalar pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt), satW); cellData.setSaturation(wPhaseIdx, satW); cellData.setSaturation(nPhaseIdx, satNW); @@ -2443,8 +2443,8 @@ void FVMPFAOPressure2P<TypeTag>::updateMaterialLaws() // initialize mobilities - Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; - Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; + Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; + Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; // initialize mobilities cellData.setMobility(wPhaseIdx, mobilityW); diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh index 92e6a7ce11894b71c40fb7e7ef056e68f4939815..fce17a7732c482acbc3b61a9830126139500da88 100644 --- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh +++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh @@ -59,8 +59,8 @@ template<class TypeTag> class FVMPFAOVelocity2P:public FVMPFAOPressure2P<TypeTag typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -106,8 +106,8 @@ template<class TypeTag> class FVMPFAOVelocity2P:public FVMPFAOPressure2P<TypeTag }; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: //! Constructs a FVMPFAOVelocity2P object @@ -187,7 +187,7 @@ template<class TypeTag> void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() { // introduce matrix R for vector rotation and R is initialized as zero matrix - FieldMatrix R(0); + DimMatrix R(0); // evaluate matrix R if (dim == 2) @@ -228,7 +228,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar q1 = source[wPhaseIdx] + source[nPhaseIdx]; // get absolute permeability of cell 1 - FieldMatrix K1(problem_.spatialParameters().intrinsicPermeability(*eIt)); + DimMatrix K1(problem_.spatialParams().intrinsicPermeability(*eIt)); // compute total mobility of cell 1 Scalar lambda1 = cellData1.mobility(wPhaseIdx) @@ -388,7 +388,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() GlobalPosition globalPos2 = outside->geometry().center(); // get absolute permeability of neighbor cell 2 - FieldMatrix K2(problem_.spatialParameters().intrinsicPermeability(*outside)); + DimMatrix K2(problem_.spatialParams().intrinsicPermeability(*outside)); // get total mobility of neighbor cell 2 Scalar lambda2 = cellData2.mobility(wPhaseIdx) @@ -415,7 +415,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() GlobalPosition globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = cellData3.mobility(wPhaseIdx) @@ -423,7 +423,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() // neighbor cell 4 GlobalPosition globalPos4(0); - FieldMatrix K4(0); + DimMatrix K4(0); Scalar lambda4 = 0; int globalIdx4 = 0; Scalar press4 = 0; @@ -455,7 +455,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() globalPos4 = innerisItoutside->geometry().center(); // get absolute permeability of neighbor cell 4 - K4 += problem_.spatialParameters().intrinsicPermeability(*innerisItoutside); + K4 += problem_.spatialParams().intrinsicPermeability(*innerisItoutside); lambda4 = cellData4.mobility(wPhaseIdx) + cellData4.mobility(nPhaseIdx); @@ -546,63 +546,63 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() integrationOuterNormaln2 *= face34vol / 2.0; // compute normal vectors nu11,nu21; nu12, nu22; nu13, nu23; nu14, nu24; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24 - globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12 - globalPos2, nu22); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3 - globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3 - globalPosFace34, nu23); - FieldVector nu14(0); + DimVector nu14(0); R.umv(globalPos4 - globalPosFace24, nu14); - FieldVector nu24(0); + DimVector nu24(0); R.umv(globalPosFace34 - globalPos4, nu24); // compute dF1, dF2, dF3, dF4 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); - FieldVector Rnu24(0); + DimVector Rnu24(0); R.umv(nu24, Rnu24); Scalar dF4 = fabs(nu14 * Rnu24); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); - FieldVector K4nu14(0); + DimVector K4nu14(0); K4.umv(nu14, K4nu14); - FieldVector K4nu24(0); + DimVector K4nu24(0); K4.umv(nu24, K4nu24); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -677,13 +677,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() T.umv(u, Tu); // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= Tu[0] / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= Tu[2] / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -754,35 +754,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J4 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24 - globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12 - globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -827,7 +827,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() + g121 * T[1][1]) * press2 - (g111 * r[0] + g121 * r[1]); // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); @@ -867,10 +867,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda2 = lambdaWBound + lambdaNWBound; } @@ -880,35 +880,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() } // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24 - globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12 - globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -918,8 +918,8 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g122 = alambda2 * (integrationOuterNormaln1 * K2nu22) / dF2; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g111 + g112; @@ -938,7 +938,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // use the pressure values to compute the fluxes @@ -946,7 +946,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() + g121 * T[1][1]) * press2 - (g111 * r[0] + g121 * r[1]); // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); @@ -987,10 +987,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1007,35 +1007,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J4 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24 - globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12 - globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1047,8 +1047,8 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g222 = lambda2 * (integrationOuterNormaln4 * K2nu22) / dF2; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g111 + g112; @@ -1067,7 +1067,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // use the pressure values to compute the fluxes @@ -1077,13 +1077,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() - g211 * r[0]; // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1123,10 +1123,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda2 = lambdaWBound + lambdaNWBound; } @@ -1136,35 +1136,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() } // compute normal vectors nu11,nu21; nu12, nu22; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu12(0); + DimVector nu12(0); R.umv(globalPosFace24 - globalPos2, nu12); - FieldVector nu22(0); + DimVector nu22(0); R.umv(globalPosFace12 - globalPos2, nu22); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu22(0); + DimVector Rnu22(0); R.umv(nu22, Rnu22); Scalar dF2 = fabs(nu12 * Rnu22); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K2nu12(0); + DimVector K2nu12(0); K2.umv(nu12, K2nu12); - FieldVector K2nu22(0); + DimVector K2nu22(0); K2.umv(nu22, K2nu22); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1174,8 +1174,8 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g122 = alambda2 * (integrationOuterNormaln1 * K2nu22) / dF2; // compute the matrix T & vector r - FieldMatrix T(0); - FieldVector r(0); + DimMatrix T(0); + DimVector r(0); Scalar coe = g111 + g112; @@ -1194,13 +1194,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f3 = T[1][0] * press1 + T[1][1] * press2 + r[1]; // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1226,7 +1226,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J1 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= -J1; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); @@ -1269,10 +1269,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1285,21 +1285,21 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g3 = boundValues[pressureIdx]; // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1311,7 +1311,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() * (-J1) * face12vol) / (2.0 * g111); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1338,7 +1338,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() GlobalPosition globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = cellData3.mobility(wPhaseIdx) @@ -1391,35 +1391,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J2 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3 - globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3 - globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1479,7 +1479,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f3 = T[2][0] * press1 + T[2][1] * press3 + r[2]; // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1519,10 +1519,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda3 = lambdaWBound + lambdaNWBound; } @@ -1532,35 +1532,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() } // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3 - globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3 - globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1570,7 +1570,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g223 = alambda3 * (integrationOuterNormaln3 * K3nu23) / dF3; // compute transmissibility matrix T = CA^{-1}B+F - FieldMatrix C(0), A(0), F(0), B(0); + DimMatrix C(0), A(0), F(0), B(0); // evaluate matrix C, F, A, B C[0][0] = -g111; @@ -1591,19 +1591,19 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() // compute T A.invert(); - FieldMatrix CAinv(C.rightmultiply(A)); + DimMatrix CAinv(C.rightmultiply(A)); F += B.leftmultiply(CAinv); - FieldMatrix T(F); + DimMatrix T(F); // compute vector r // evaluate r1, r2 - FieldVector r1(0), r2(0); + DimVector r1(0), r2(0); r1[1] = -g213 * g2; r2[0] = -J1 * face12vol / 2.0; r2[1] = g213 * g2; // compute r = CA^{-1}r1 - FieldVector r(0); + DimVector r(0); CAinv.umv(r2, r); r += r1; @@ -1611,7 +1611,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f3 = T[1][0] * press1 + T[1][1] * press3 + r[1]; // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1653,10 +1653,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1706,10 +1706,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda1 = lambdaWBound + lambdaNWBound; } @@ -1719,21 +1719,21 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() } // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); // compute dF1 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1751,13 +1751,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f3 = T3 * press1 - r3; // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1771,21 +1771,21 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J3 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); // compute dF1 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1800,7 +1800,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f1 = T * press1 + r; // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); @@ -1827,7 +1827,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() GlobalPosition globalPos3 = nextisItoutside->geometry().center(); // get absolute permeability of neighbor cell 3 - FieldMatrix K3(problem_.spatialParameters().intrinsicPermeability(*nextisItoutside)); + DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside)); // get total mobility of neighbor cell 3 Scalar lambda3 = cellData3.mobility(wPhaseIdx) @@ -1906,10 +1906,10 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar lambdaNWBound = 0; lambdaWBound = MaterialLaw::krw( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; lambdaNWBound = MaterialLaw::krn( - problem_.spatialParameters().materialLawParams(*eIt), satW) + problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; alambda3 = lambdaWBound + lambdaNWBound; } @@ -1919,35 +1919,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() } // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3 - globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3 - globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -1957,8 +1957,8 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g223 = alambda3 * (integrationOuterNormaln3 * K3nu23) / dF3; // compute the matrix T & vector r - FieldMatrix T(0); - FieldVector r(0); + DimMatrix T(0); + DimVector r(0); Scalar coe = g221 + g223; @@ -1977,13 +1977,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar f3 = T[1][0] * press1 + T[1][1] * press3 + r[1]; // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); @@ -1997,35 +1997,35 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar J2 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]); // compute normal vectors nu11,nu21; nu13, nu23; - FieldVector nu11(0); + DimVector nu11(0); R.umv(globalPosFace13 - globalPos1, nu11); - FieldVector nu21(0); + DimVector nu21(0); R.umv(globalPos1 - globalPosFace12, nu21); - FieldVector nu13(0); + DimVector nu13(0); R.umv(globalPos3 - globalPosFace13, nu13); - FieldVector nu23(0); + DimVector nu23(0); R.umv(globalPos3 - globalPosFace34, nu23); // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu' - FieldVector Rnu21(0); + DimVector Rnu21(0); R.umv(nu21, Rnu21); Scalar dF1 = fabs(nu11 * Rnu21); - FieldVector Rnu23(0); + DimVector Rnu23(0); R.umv(nu23, Rnu23); Scalar dF3 = fabs(nu13 * Rnu23); // compute components needed for flux calculation, denoted as 'g' - FieldVector K1nu11(0); + DimVector K1nu11(0); K1.umv(nu11, K1nu11); - FieldVector K1nu21(0); + DimVector K1nu21(0); K1.umv(nu21, K1nu21); - FieldVector K3nu13(0); + DimVector K3nu13(0); K3.umv(nu13, K3nu13); - FieldVector K3nu23(0); + DimVector K3nu23(0); K3.umv(nu23, K3nu23); Scalar g111 = alambda1 * (integrationOuterNormaln1 * K1nu11) / dF1; Scalar g121 = alambda1 * (integrationOuterNormaln1 * K1nu21) / dF1; @@ -2037,8 +2037,8 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() Scalar g223 = lambda3 * (integrationOuterNormaln3 * K3nu23) / dF3; // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r - FieldMatrix A(0), B(0); - FieldVector r1(0), r(0); + DimMatrix A(0), B(0); + DimVector r1(0), r(0); // evaluate matrix A, B A[0][0] = g113; @@ -2057,7 +2057,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() // compute T and r A.invert(); B.leftmultiply(A); - FieldMatrix T(B); + DimMatrix T(B); A.umv(r1, r); // use the pressure values to compute the fluxes @@ -2067,13 +2067,13 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity() + g221 * r[1]); // evaluate velocity of facet 'isIt' - FieldVector vector1 = unitOuterNormaln1; + DimVector vector1 = unitOuterNormaln1; vector1 *= f1 / face12vol; vector1 += cellData1.fluxData().velocityTotal(indexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1); // evaluate velocity of facet 'nextisIt' - FieldVector vector3 = unitOuterNormaln3; + DimVector vector3 = unitOuterNormaln3; vector3 *= f3 / face13vol; vector3 += cellData1.fluxData().velocityTotal(nextindexInInside); cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3); diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh index cf80f525f2b71d53b6cb8ae83075f92fd814bcd0..41ef58e47cd2d1921f9a4116ceb478d9a5c87ff7 100644 --- a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh +++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh @@ -203,7 +203,7 @@ public: // eval diffusion tensor, ASSUMING to be constant over each cell Dune::FieldMatrix<Scalar,dim,dim> K(0); - K = problem_.spatialParameters().intrinsicPermeability(element); + K = problem_.spatialParams().intrinsicPermeability(element); K *= (cellData.mobility(wPhaseIdx) + cellData.mobility(nPhaseIdx)); diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh index 4671aa2f4f593bb235158fd453a6ee2da10f4807..ea45381800be372291138a31f612d40a4fab7783 100644 --- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh +++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh @@ -72,8 +72,8 @@ template<class TypeTag> class MimeticPressure2P typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, Variables) Variables; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -107,7 +107,6 @@ template<class TypeTag> class MimeticPressure2P typedef typename GridView::IntersectionIterator IntersectionIterator; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; typedef typename GET_PROP_TYPE(TypeTag, LocalStiffness) LocalStiffness; typedef Dune::BlockVector< Dune::FieldVector<Scalar, 1> > TraceType; @@ -353,13 +352,13 @@ void MimeticPressure2P<TypeTag>::updateMaterialLaws() //determine phase saturations from primary saturation variable Scalar satW = cellData.saturation(wPhaseIdx); - Scalar pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW); + Scalar pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt), satW); cellData.setCapillaryPressure(pc); // initialize mobilities - Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; - Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; + Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; + Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; // initialize mobilities cellData.setMobility(wPhaseIdx, mobilityW); diff --git a/dumux/decoupled/2p/fluxData2p.hh b/dumux/decoupled/2p/fluxData2p.hh index 27ff6c5bfd43d46c45ec67c995647acec45a194b..58c3857c8aa635ed3169969f9309f9f4fd63fc62 100644 --- a/dumux/decoupled/2p/fluxData2p.hh +++ b/dumux/decoupled/2p/fluxData2p.hh @@ -66,8 +66,8 @@ private: numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; - typedef Dune::FieldVector<Scalar, dim> FieldVector; - typedef Dune::FieldVector<FieldVector, 2 * dim> VelocityVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; + typedef Dune::FieldVector<DimVector, 2 * dim> VelocityVector; VelocityVector velocity_[numPhases]; Scalar potential_[2 * dim][numPhases]; @@ -82,7 +82,7 @@ public: { for (int phase = 0; phase < numPhases; phase++) { - velocity_[phase][face] = FieldVector(0.0); + velocity_[phase][face] = DimVector(0.0); potential_[face][phase] = 0.0; } @@ -99,7 +99,7 @@ public: * \param phaseIdx Index of a fluid phase * \param indexInInside Index of the cell-cell interface in this cell */ - const FieldVector& velocity(int phaseIdx, int indexInInside) + const DimVector& velocity(int phaseIdx, int indexInInside) { return velocity_[phaseIdx][indexInInside]; } @@ -109,7 +109,7 @@ public: * \param phaseIdx Index of a fluid phase * \param indexInInside Index of the cell-cell interface in this cell */ - const FieldVector& velocity(int phaseIdx, int indexInInside) const + const DimVector& velocity(int phaseIdx, int indexInInside) const { return velocity_[phaseIdx][indexInInside]; } @@ -120,7 +120,7 @@ public: * \param indexInInside Index of the cell-cell interface in this cell * \param velocity Phase velocity vector which is stored */ - void setVelocity(int phaseIdx, int indexInInside, const FieldVector& velocity) + void setVelocity(int phaseIdx, int indexInInside, const DimVector& velocity) { velocity_[phaseIdx][indexInInside] = velocity; } @@ -131,7 +131,7 @@ public: * \param indexInInside Index of the cell-cell interface in this cell * \param velocity Phase velocity vector which is added */ - void addVelocity(int phaseIdx, int indexInInside, const FieldVector& velocity) + void addVelocity(int phaseIdx, int indexInInside, const DimVector& velocity) { velocity_[phaseIdx][indexInInside] += velocity; } @@ -154,7 +154,7 @@ public: * * \param indexInInside Index of the cell-cell interface in this cell */ - FieldVector velocityTotal(int indexInInside) + DimVector velocityTotal(int indexInInside) { return velocity_[wPhaseIdx][indexInInside] + velocity_[nPhaseIdx][indexInInside]; @@ -164,7 +164,7 @@ public: * * \param indexInInside Index of the cell-cell interface in this cell */ - FieldVector velocityTotal(int indexInInside) const + DimVector velocityTotal(int indexInInside) const { return velocity_[wPhaseIdx][indexInInside] + velocity_[nPhaseIdx][indexInInside]; diff --git a/dumux/decoupled/2p/impes/gridadaptionindicator2p.hh b/dumux/decoupled/2p/impes/gridadaptionindicator2p.hh index 2c3a99e8fb0f354e5d9e8df6345fc732f8117a51..d67a0d5f31f984d39a2250f8f405099adb894874 100644 --- a/dumux/decoupled/2p/impes/gridadaptionindicator2p.hh +++ b/dumux/decoupled/2p/impes/gridadaptionindicator2p.hh @@ -176,6 +176,12 @@ public: return (indicatorVector_[problem_.elementMapper().map(element)] < coarsenBound_); } + /*! \brief Initializes the adaption indicator class*/ + void init() + { + refineBound_ = 0.; + coarsenBound_ = 0.; + }; /*! @brief Constructs a GridAdaptionIndicator instance * @@ -189,8 +195,6 @@ public: { refinetol_ = GET_PARAM(TypeTag, Scalar, RefineTolerance); coarsentol_ = GET_PARAM(TypeTag, Scalar, CoarsenTolerance); - refineBound_ = 0.; - coarsenBound_ = 0.; } private: diff --git a/dumux/decoupled/2p/impes/impesproblem2p.hh b/dumux/decoupled/2p/impes/impesproblem2p.hh index e3c7f304a4c1f018a5efbc486f3db0cf50d7381c..378e0f3cb2916cff9f5c554e8ed32fce8875bc8b 100644 --- a/dumux/decoupled/2p/impes/impesproblem2p.hh +++ b/dumux/decoupled/2p/impes/impesproblem2p.hh @@ -97,7 +97,7 @@ public: * * \param timeManager The time manager * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ IMPESProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams) : ParentType(timeManager, gridView), diff --git a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh index bd7d5461132be2f2ce9f7fa7be4c626e34790a30..f2d863f5da8b64f87a9de92a64c1149aea315041 100644 --- a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh +++ b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh @@ -55,8 +55,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState; @@ -76,9 +76,9 @@ private: typedef typename GridView::Traits::template Codim<0>::Entity Element; typedef typename GridView::template Codim<0>::EntityPointer ElementPointer; typedef typename GridView::Intersection Intersection; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix; public: //! Returns capillary diffusion term @@ -89,7 +89,7 @@ public: * \param satJ saturation of neighbor element * \param pcGradient gradient of capillary pressure between element I and J */ - void getFlux (FieldVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const + void getFlux (DimVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const DimVector& pcGradient) const { ElementPointer element = intersection.inside(); // get global coordinate of cell center @@ -120,13 +120,13 @@ public: fluidState.setPressure(wPhaseIdx, referencePressure); fluidState.setPressure(nPhaseIdx, referencePressure); fluidState.setTemperature(temperature); - mobilityWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satI); + mobilityWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satI); mobilityWI /= FluidSystem::viscosity(fluidState, wPhaseIdx); - mobilityNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satI); + mobilityNWI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satI); mobilityNWI /= FluidSystem::viscosity(fluidState, nPhaseIdx); } - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); if (intersection.neighbor()) { @@ -140,18 +140,18 @@ public: const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center(); // distance vector between barycenters - FieldVector distVec = globalPosNeighbor - globalPos; + DimVector distVec = globalPosNeighbor - globalPos; // compute distance between cell centers Scalar dist = distVec.two_norm(); - FieldVector unitDistVec(distVec); + DimVector unitDistVec(distVec); unitDistVec /= dist; // get permeability - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element), - problem_.spatialParameters().intrinsicPermeability(*neighborPointer)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element), + problem_.spatialParams().intrinsicPermeability(*neighborPointer)); Scalar mobilityWJ = 0; @@ -169,9 +169,9 @@ public: fluidState.setPressure(nPhaseIdx, referencePressure); fluidState.setTemperature(temperature); - mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ); + mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighborPointer), satJ); mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx); - mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ); + mobilityNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighborPointer), satJ); mobilityNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx); } Scalar mobilityWMean = 0.5*(mobilityWI + mobilityWJ); @@ -181,8 +181,8 @@ public: else { // get permeability - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element)); Scalar mobilityWJ = 0; Scalar mobilityNWJ = 0; @@ -192,9 +192,9 @@ public: fluidState.setPressure(wPhaseIdx, referencePressure); fluidState.setPressure(nPhaseIdx, referencePressure); fluidState.setTemperature(temperature); - mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satJ); + mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satJ); mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx); - mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satJ); + mobilityNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satJ); mobilityNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx); Scalar mobWMean = 0.5 * (mobilityWI + mobilityWJ); diff --git a/dumux/decoupled/2p/transport/fv/convectivepart.hh b/dumux/decoupled/2p/transport/fv/convectivepart.hh index 59ec527c8a27504e112f44b0b1238775dad1fa60..2572f07ff8250c407195b3da60cabe23ff614ffb 100644 --- a/dumux/decoupled/2p/transport/fv/convectivepart.hh +++ b/dumux/decoupled/2p/transport/fv/convectivepart.hh @@ -49,7 +49,7 @@ private: enum{dim = GridView::dimension, dimWorld = GridView::dimensionworld}; typedef typename GridView::Traits::template Codim<0>::Entity Element; typedef typename GridView::Intersection Intersection; - typedef Dune::FieldVector<Scalar, dimWorld> FieldVector; + typedef Dune::FieldVector<Scalar, dimWorld> DimVector; public: /*! \brief Returns convective term for current element face @@ -69,7 +69,7 @@ public: * \param satI Saturation of current element * \param satJ Saturation of neighbor element */ - void getFlux(FieldVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const + void getFlux(DimVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const {} //! Constructs a ConvectivePart object diff --git a/dumux/decoupled/2p/transport/fv/diffusivepart.hh b/dumux/decoupled/2p/transport/fv/diffusivepart.hh index da262e560922151788886eb36305f5c14a112227..7dee5ba53c9c7a4b947ac2cff6d0aef2b843ac7b 100644 --- a/dumux/decoupled/2p/transport/fv/diffusivepart.hh +++ b/dumux/decoupled/2p/transport/fv/diffusivepart.hh @@ -48,7 +48,7 @@ private: enum{dim = GridView::dimension}; typedef typename GridView::Traits::template Codim<0>::Entity Element; typedef typename GridView::Intersection Intersection; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: /*! \brief Returns diffusive term for current element face @@ -59,7 +59,7 @@ public: * \param satJ saturation of neighbor element * \param pcGradient gradient of capillary pressure between element I and J */ - void getFlux(FieldVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const + void getFlux(DimVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const DimVector& pcGradient) const {} /*! \brief Returns diffusive term for current element face @@ -70,8 +70,8 @@ public: * \param satGradient gradient of saturation between element I and J * \param time time */ - void getFlux(FieldVector& flux, const Intersection& intersection, - const Scalar satIntersection, const FieldVector& satGradient, const Scalar time) const + void getFlux(DimVector& flux, const Intersection& intersection, + const Scalar satIntersection, const DimVector& satGradient, const Scalar time) const {} /*! \brief Returns diffusive term for current element face @@ -84,8 +84,8 @@ public: * \param satI saturation of current element * \param satJ saturation of neighbor element */ - void getFlux(FieldVector& flux, const Intersection& intersection, - const Scalar satIntersection, const FieldVector& satGradient, const Scalar time, + void getFlux(DimVector& flux, const Intersection& intersection, + const Scalar satIntersection, const DimVector& satGradient, const Scalar time, Scalar satI, Scalar satJ) const {} diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh index b6973da1bf413a4ae44d45e8a8b26cd4df440903..eede6bf47b2f39eaed84ed7130f83f18708ed354 100644 --- a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh +++ b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh @@ -45,8 +45,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -85,7 +85,7 @@ private: typedef typename GridView::Intersection Intersection; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; public: @@ -96,10 +96,10 @@ public: void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux, const Element& element, int phaseIdx = -1) { - if (hasHangingNode_) - { - return; - } +// if (hasHangingNode_) +// { +// return; +// } ParentType::addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, element, phaseIdx); } @@ -111,16 +111,16 @@ public: void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux, const Intersection& intersection, int phaseIdx = -1) { - if (hasHangingNode_) - { - return; - } - Scalar parentMob = 0.5; Scalar parentVis = 1; // ParentType::addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, intersection, phaseIdx); ParentType::addFlux(parentMob, parentMob, parentVis, parentVis, flux, intersection, phaseIdx); + if (hasHangingNode_) + { + return; + } + ElementPointer element = intersection.inside(); //coordinates of cell center @@ -137,7 +137,7 @@ public: Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar dPcdSI = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(*element), satI); + Scalar dPcdSI = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satI); const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal(); @@ -170,14 +170,14 @@ public: Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx); Scalar lambdaNWJ = cellDataI.mobility(nPhaseIdx); - Scalar dPcdSJ = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(*neighbor), satJ); + Scalar dPcdSJ = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*neighbor), satJ); // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element), - problem_.spatialParameters().intrinsicPermeability(*neighbor)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element), + problem_.spatialParams().intrinsicPermeability(*neighbor)); Dune::FieldVector < Scalar, dim > permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -204,8 +204,8 @@ public: dS += eps_; } - Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighbor), std::abs(satPlus)) / viscosityW; - dLambdaWdS -= MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighbor), std::abs(satMinus)) / viscosityW; + Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satPlus)) / viscosityW; + dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satMinus)) / viscosityW; dLambdaWdS /= (dS); if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside)) @@ -227,8 +227,8 @@ public: dS += eps_; } - Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighbor), satPlus) / viscosityNW; - dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighbor), satMinus) / viscosityNW; + Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satPlus) / viscosityNW; + dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satMinus) / viscosityNW; dLambdaNWdS /= (dS); Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ); @@ -288,10 +288,10 @@ public: //permeability vector at boundary // compute vectorized permeabilities - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element)); Dune::FieldVector<Scalar, dim> permeability(0); meanPermeability.mv(unitOuterNormal, permeability); @@ -322,7 +322,7 @@ public: } - Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(*element), satWBound); + Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satWBound); Scalar lambdaWBound = 0; Scalar lambdaNWBound = 0; @@ -337,8 +337,8 @@ public: Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx); Scalar viscosityNWBound = FluidSystem::viscosity(fluidState, nPhaseIdx); - lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satWBound) / viscosityWBound; - lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satWBound) / viscosityNWBound; + lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityWBound; + lambdaNWBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityNWBound; Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist; @@ -362,8 +362,8 @@ public: dS += eps_; } - Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satPlus) / viscosityW; - dLambdaWdS -= MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satMinus) / viscosityW; + Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityW; + dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityW; dLambdaWdS /= (dS); if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside)) @@ -385,8 +385,8 @@ public: dS += eps_; } - Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satPlus) / viscosityNW; - dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satMinus) / viscosityNW; + Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityNW; + dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityNW; dLambdaNWdS /= (dS); Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound); @@ -458,7 +458,9 @@ public: return (returnValue == cflFluxDefault) ? 0.95 / returnValue : 1 / returnValue; } else - return 1e100; + { + return cflFluxDefault; + } } /*! \brief Returns the CFL time-step @@ -467,7 +469,7 @@ public: */ Scalar getDt(const Element& element) { - return getCFLFluxFunction(element) * problem_.spatialParameters().porosity(element) * element.geometry().volume(); + return getCFLFluxFunction(element) * problem_.spatialParams().porosity(element) * element.geometry().volume(); } /*! \brief Constructs an EvalCflFluxDefault object diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh index fbde8410096b45905cb16f5511576c63ec7b0fdb..de5fe8edda1463a2a389d925f43c861fa2eea906 100644 --- a/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh +++ b/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh @@ -108,7 +108,7 @@ public: */ Scalar getDt(const Element& element) { - return (getCFLFluxFunction(element) * problem_.spatialParameters().porosity(element) * element.geometry().volume()); + return (getCFLFluxFunction(element) * problem_.spatialParams().porosity(element) * element.geometry().volume()); } /*! \brief Constructs an EvalCflFluxDefault object @@ -235,8 +235,8 @@ private: template<class TypeTag> typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCFLFluxFunction(const Element& element) { - Scalar residualSatW = problem_.spatialParameters().materialLawParams(element).Swr(); - Scalar residualSatNW = problem_.spatialParameters().materialLawParams(element).Snr(); + Scalar residualSatW = problem_.spatialParams().materialLawParams(element).Swr(); + Scalar residualSatNW = problem_.spatialParams().materialLawParams(element).Snr(); // compute dt restriction Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNW; diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh index e5ce3576152326e59967b6552af6bd2f3a856ac3..9a824d7eac1a8f7d09d6f535d342b9ea9304b892 100644 --- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh +++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh @@ -82,8 +82,8 @@ class FVSaturation2P: public FVTransport<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, CapillaryFlux) CapillaryFlux; typedef typename GET_PROP_TYPE(TypeTag, GravityFlux) GravityFlux; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; @@ -129,7 +129,7 @@ class FVSaturation2P: public FVTransport<TypeTag> typedef typename GridView::Intersection Intersection; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; Velocity& velocity() { @@ -435,7 +435,7 @@ void FVSaturation2P<TypeTag>::getFlux(Scalar& update, const Intersection& inters // cell volume, assume linear map here Scalar volume = elementI->geometry().volume(); - Scalar porosity = problem_.spatialParameters().porosity(*elementI); + Scalar porosity = problem_.spatialParams().porosity(*elementI); if (compressibility_) { @@ -531,7 +531,7 @@ void FVSaturation2P<TypeTag>::getFlux(Scalar& update, const Intersection& inters pcGradient *= (pcI - pcJ) / dist; // get the diffusive part - FieldVector flux(0.); + DimVector flux(0.); capillaryFlux().getFlux(flux, intersection, satWI, satWJ, pcGradient); Scalar capillaryFlux = (flux * unitOuterNormal * faceArea); @@ -620,7 +620,7 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti // cell volume, assume linear map here Scalar volume = elementI->geometry().volume(); - Scalar porosity = problem_.spatialParameters().porosity(*elementI); + Scalar porosity = problem_.spatialParams().porosity(*elementI); if (compressibility_) { @@ -679,7 +679,7 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti } } - Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*elementI), satWBound); + Scalar pcBound = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*elementI), satWBound); Scalar lambdaW = 0; Scalar lambdaNW = 0; @@ -697,12 +697,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti { if (compressibility_) { - lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*elementI), satWBound) + lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*elementI), satWBound) / FluidSystem::viscosity(cellDataI.fluidState(), wPhaseIdx); } else { - lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*elementI), satWBound) + lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*elementI), satWBound) / viscosity_[wPhaseIdx]; } } @@ -719,12 +719,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti { if (compressibility_) { - lambdaNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*elementI), satWBound) + lambdaNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*elementI), satWBound) / FluidSystem::viscosity(cellDataI.fluidState(), nPhaseIdx); } else { - lambdaNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*elementI), satWBound) + lambdaNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*elementI), satWBound) / viscosity_[nPhaseIdx]; } } @@ -742,7 +742,7 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti pcGradient *= (pcI - pcBound) / dist; // get the diffusive part -> give 1-sat because sat = S_n and lambda = lambda(S_w) and pc = pc(S_w) - FieldVector flux(0.); + DimVector flux(0.); capillaryFlux().getFlux(flux, intersection, satWI, satWBound, pcGradient); Scalar capillaryFlux = flux * unitOuterNormal * faceArea; @@ -908,7 +908,7 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element, // cell volume, assume linear map here Scalar volume = element.geometry().volume(); - Scalar porosity = problem_.spatialParameters().porosity(element); + Scalar porosity = problem_.spatialParams().porosity(element); if (compressibility_) { @@ -1037,7 +1037,7 @@ void FVSaturation2P<TypeTag>::updateMaterialLaws() Scalar satW = cellData.saturation(wPhaseIdx); Scalar satNW = cellData.saturation(nPhaseIdx); - Scalar pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW); + Scalar pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt), satW); cellData.setSaturation(wPhaseIdx, satW); cellData.setSaturation(nPhaseIdx, satNW); @@ -1045,8 +1045,8 @@ void FVSaturation2P<TypeTag>::updateMaterialLaws() cellData.setCapillaryPressure(pc); // initialize mobilities - Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; - Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; + Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx]; + Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx]; // initialize mobilities cellData.setMobility(wPhaseIdx, mobilityW); diff --git a/dumux/decoupled/2p/transport/fv/gravitypart.hh b/dumux/decoupled/2p/transport/fv/gravitypart.hh index 1f875909092ebb2601ea2561f43e796f8ea57789..8de8f7011fad20529bede3203d9c9271e2e39550 100644 --- a/dumux/decoupled/2p/transport/fv/gravitypart.hh +++ b/dumux/decoupled/2p/transport/fv/gravitypart.hh @@ -56,8 +56,8 @@ private: typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; - typedef typename SpatialParameters::MaterialLaw MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; + typedef typename SpatialParams::MaterialLaw MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState; @@ -77,9 +77,9 @@ private: typedef typename GridView::template Codim<0>::EntityPointer ElementPointer; typedef typename GridView::template Codim<0>::Iterator ElementIterator; typedef typename GridView::Intersection Intersection; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; - typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; + typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix; public: /*! \brief Returns convective term for current element face @@ -89,7 +89,7 @@ public: * \param satI saturation of current element * \param satJ saturation of neighbor element */ - void getFlux(FieldVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const + void getFlux(DimVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const { ElementPointer element = intersection.inside(); @@ -111,13 +111,13 @@ public: } else { - lambdaWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satI); + lambdaWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satI); lambdaWI /= viscosity_[wPhaseIdx]; - lambdaNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satI); + lambdaNWI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satI); lambdaNWI /= viscosity_[nPhaseIdx]; } - FieldMatrix meanPermeability(0); + DimMatrix meanPermeability(0); if (intersection.neighbor()) { @@ -128,9 +128,9 @@ public: CellData& cellDataJ = problem_.variables().cellData(globalIdxJ); // get permeability - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element), - problem_.spatialParameters().intrinsicPermeability(*neighborPointer)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element), + problem_.spatialParams().intrinsicPermeability(*neighborPointer)); //get lambda_bar = lambda_n*f_w if (preComput_) @@ -140,22 +140,22 @@ public: } else { - lambdaWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ); + lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighborPointer), satJ); lambdaWJ /= viscosity_[wPhaseIdx]; - lambdaNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ); + lambdaNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighborPointer), satJ); lambdaNWJ /= viscosity_[nPhaseIdx]; } } else { // get permeability - problem_.spatialParameters().meanK(meanPermeability, - problem_.spatialParameters().intrinsicPermeability(*element)); + problem_.spatialParams().meanK(meanPermeability, + problem_.spatialParams().intrinsicPermeability(*element)); //calculate lambda_n*f_w at the boundary - lambdaWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satJ); + lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satJ); lambdaWJ /= viscosity_[wPhaseIdx]; - lambdaNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satJ); + lambdaNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satJ); lambdaNWJ /= viscosity_[nPhaseIdx]; } diff --git a/dumux/decoupled/2p/transport/transportproblem2p.hh b/dumux/decoupled/2p/transport/transportproblem2p.hh index 24cae45d44a2080bb17cf2e969a99bd0cd220cec..6ff25839ca2f8e97e75d4be8d72372ac428d6c19 100644 --- a/dumux/decoupled/2p/transport/transportproblem2p.hh +++ b/dumux/decoupled/2p/transport/transportproblem2p.hh @@ -65,7 +65,7 @@ class TransportProblem2P : public OneModelProblem<TypeTag> // material properties typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters; + typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams; typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes; typedef typename SolutionTypes::ScalarSolution Solution; @@ -104,7 +104,7 @@ public: cFLFactor_ = GET_PARAM(TypeTag, Scalar, CFLFactor); newSpatialParams_ = true; - spatialParameters_ = new SpatialParameters(gridView); + spatialParams_ = new SpatialParams(gridView); gravity_ = 0; if (GET_PARAM(TypeTag, bool, EnableGravity)) @@ -116,11 +116,11 @@ public: * * \param timeManager The time manager * \param gridView The grid view - * \param spatialParameters SpatialParameters instantiation + * \param spatialParams SpatialParams instantiation */ - TransportProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters) + TransportProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams) : ParentType(timeManager, gridView), - gravity_(0),spatialParameters_(spatialParameters) + gravity_(0),spatialParams_(spatialParams) { cFLFactor_ = GET_PARAM(TypeTag, Scalar, CFLFactor); @@ -135,7 +135,7 @@ public: { if (newSpatialParams_) { - delete spatialParameters_; + delete spatialParams_; } } @@ -206,14 +206,22 @@ public: /*! * \brief Returns the spatial parameters object. */ - SpatialParameters &spatialParameters() - { return *spatialParameters_; } + SpatialParams &spatialParams() + { return *spatialParams_; } + + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + SpatialParams &spatialParameters() + { return *spatialParams_; } /*! * \brief Returns the spatial parameters object. */ - const SpatialParameters &spatialParameters() const - { return *spatialParameters_; } + const SpatialParams &spatialParams() const + { return *spatialParams_; } + + DUMUX_DEPRECATED_MSG("use spatialParams() method instead") + const SpatialParams &spatialParameters() const + { return *spatialParams_; } void timeIntegration() { @@ -248,7 +256,7 @@ private: GlobalPosition gravity_; // material properties - SpatialParameters* spatialParameters_; + SpatialParams* spatialParams_; bool newSpatialParams_; Scalar cFLFactor_;