From 6c8c7ce5aea9cd3cdb3df8042a4ecae3fecd7ed3 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 16:47:04 +0200 Subject: [PATCH 01/46] [implicit][localjacobian] adjust assembly map setup The new type of the stencils requires a different way of finding the flux vars required to compute the derivative w.r.t. a neighbor. --- dumux/implicit/cellcentered/localjacobian.hh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index efb9210922..8eca5f7b9c 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -128,8 +128,10 @@ public: auto fluxVarsIdx = scvFaceJ.index(); // if globalI is in flux var stencil, add to list - if (this->model_().fluxVars(fluxVarsIdx).stencil().count(globalI)) - fluxVarIndices.push_back(fluxVarsIdx); + const auto& fluxStencil = this->model_().fluxVars(fluxVarsIdx).stencil(); + for (auto it = fluxStencil.begin(); it != fluxStencil.end(); ++it) + if (*it == globalI) + fluxVarIndices.push_back(fluxVarsIdx); } assemblyMap_[globalI].emplace_back(std::move(fluxVarIndices)); -- GitLab From bcce2a4db7823fb18f8b7cf38a7ed28f3159ecd8 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 16:55:42 +0200 Subject: [PATCH 02/46] [implicit][localjacobian] replace loop by range-based loop --- dumux/implicit/cellcentered/localjacobian.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 8eca5f7b9c..898cc7f4d7 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -129,8 +129,8 @@ public: // if globalI is in flux var stencil, add to list const auto& fluxStencil = this->model_().fluxVars(fluxVarsIdx).stencil(); - for (auto it = fluxStencil.begin(); it != fluxStencil.end(); ++it) - if (*it == globalI) + for (auto globalIdx : fluxStencil) + if (globalIdx == globalI) fluxVarIndices.push_back(fluxVarsIdx); } -- GitLab From 7438918f9b213607e771884fd56c92209faeefdb Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 16:57:05 +0200 Subject: [PATCH 03/46] [implicit][localjacobian] move temporary volvars into vector --- dumux/implicit/cellcentered/localjacobian.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 898cc7f4d7..f362b3d8ec 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -296,7 +296,7 @@ public: neighborDeriv /= delta; // restore the original state of the scv's volume variables - this->model_().curVolVars_(scv) = origVolVars; + this->model_().curVolVars_(scv) = std::move(origVolVars); #if HAVE_VALGRIND for (unsigned i = 0; i < partialDeriv.size(); ++i) -- GitLab From 0f85262bc8ac31b10df694e360e028f65845d219 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:00:21 +0200 Subject: [PATCH 04/46] [implicit] change stencil type The stencil type is now std::vector. The flux stencil is now stored in the fluxvariables class and no longer in darcy's law etc. --- dumux/implicit/cellcentered/stencils.hh | 18 +++++++++++++++--- dumux/implicit/fluxvariables.hh | 6 ++++-- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index ad6f15b3e3..c70bfdbdf1 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -40,17 +40,29 @@ class CCElementStencils using IndexType = typename GridView::IndexSet::IndexType; using Element = typename GridView::template Codim<0>::Entity; // TODO a separate stencil class all stencils can derive from? - using Stencil = std::set; + using Stencil = std::vector; public: void update(const Problem& problem, const Element& element) { + elementStencil_.clear(); for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) { auto&& fluxStencil = problem.model().fluxVars(scvf).stencil(); - elementStencil_.insert(fluxStencil.begin(), fluxStencil.end()); + elementStencil_.insert(elementStencil_.end(), fluxStencil.begin(), fluxStencil.end()); } + // make values in elementstencil unique + std::sort(elementStencil_.begin(), elementStencil_.end()); + elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); + neighborStencil_ = elementStencil_; - neighborStencil_.erase(problem.elementMapper().index(element)); + for (auto it = neighborStencil_.begin(); it != neighborStencil_.end(); ++it) + { + if (*it == problem.elementMapper().index(element)) + { + neighborStencil_.erase(it); + break; + } + } } //! The full element stencil (all element this element is interacting with) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 3e1d681471..5c144517fb 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -51,7 +51,8 @@ class FluxVariables using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = std::set; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Stencil = std::vector; using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); @@ -103,7 +104,8 @@ class FluxVariables using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; - using Stencil = std::set; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Stencil = std::vector; using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); -- GitLab From f0d5599f7dcfbd22ba7b5601309b606dc06b99bb Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:03:22 +0200 Subject: [PATCH 05/46] [implicit][tpfa][darcyslaw] change to static member functions The darcys law class has no private variables any longer. This substantially reduces the memory used by the program. Member functions are now static and transmissibilities are computed fresh at every call. --- .../implicit/cellcentered/tpfa/darcyslaw.hh | 239 ++++++------------ 1 file changed, 73 insertions(+), 166 deletions(-) diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh index 41f8016f2c..6408ce9fab 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh @@ -72,234 +72,141 @@ class CCTpfaDarcysLaw public: - void update(const Problem &problem, - const Element& element, - const SubControlVolumeFace &scvFace) + static Scalar flux(const Problem& problem, + const SubControlVolumeFace& scvFace, + const IndexType phaseIdx, + VolumeVariables* boundaryVolVars, + const bool boundaryVolVarsUpdated) { - problemPtr_ = &problem; - scvFacePtr_ = &scvFace; - enableGravity_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity); + Scalar tij = calculateTransmissibility_(problem, scvFace); - updateTransmissibilities_(); - updateStencil_(); - } - - void update(const Problem& problem, - const Element& element, - const SubControlVolumeFace &scvFace, - VolumeVariables* boundaryVolVars) - { - boundaryVolVars_ = boundaryVolVars; - update(problem, element, scvFace); - } - - void updateTransmissibilities(const Problem &problem, const SubControlVolumeFace &scvFace) - { - updateTransmissibilities_(); - } - - void beginFluxComputation(bool boundaryVolVarsUpdated = false) - { // Get the inside volume variables - const auto insideScvIdx = scvFace_().insideScvIdx(); - const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto* insideVolVars = &problem_().model().curVolVars(insideScv); + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto* insideVolVars = &problem.model().curVolVars(insideScv); // and the outside volume variables const VolumeVariables* outsideVolVars; - if (!scvFace_().boundary()) - outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx()); + if (!scvFace.boundary()) + outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); else { - outsideVolVars = boundaryVolVars_; if (!boundaryVolVarsUpdated) { // update the boudary volvars for Dirichlet boundaries - const auto element = problem_().model().fvGeometries().element(insideScv); - const auto dirichletPriVars = problem_().dirichlet(element, scvFace_()); - boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv); + const auto element = problem.model().fvGeometries().element(insideScv); + const auto dirichletPriVars = problem.dirichlet(element, scvFace); + boundaryVolVars->update(dirichletPriVars, problem, element, insideScv); } + outsideVolVars = boundaryVolVars; } - // loop over all phases to compute the volume flux - for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) + auto hInside = insideVolVars->pressure(phaseIdx); + auto hOutside = outsideVolVars->pressure(phaseIdx); + + if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity)) { - auto hInside = insideVolVars->pressure(phaseIdx); - auto hOutside = outsideVolVars->pressure(phaseIdx); + // do averaging for the density + const auto rhoInside = insideVolVars->density(phaseIdx); + const auto rhoOutide = outsideVolVars->density(phaseIdx); + const auto rho = (rhoInside + rhoOutide)*0.5; - if (enableGravity_) - { - // do averaging for the density - const auto rhoInside = insideVolVars->density(phaseIdx); - const auto rhoOutide = outsideVolVars->density(phaseIdx); - const auto rho = (rhoInside + rhoOutide)*0.5; - - - // ask for the gravitational acceleration in the inside neighbor - const auto xInside = insideScv.center(); - const auto gInside = problem_().gravityAtPos(xInside); - - hInside -= rho*(gInside*xInside); - - // and the outside neighbor - if (scvFace_().boundary()) - { - const auto xOutside = scvFace_().center(); - const auto gOutside = problem_().gravityAtPos(xOutside); - hOutside -= rho*(gOutside*xOutside); - } - else - { - const auto outsideScvIdx = scvFace_().outsideScvIdx(); - const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx); - const auto xOutside = outsideScv.center(); - const auto gOutside = problem_().gravityAtPos(xOutside); - hOutside -= rho*(gOutside*xOutside); - } - } - kGradPNormal_[phaseIdx] = tij_*(hInside - hOutside); + // ask for the gravitational acceleration in the inside neighbor + const auto xInside = insideScv.center(); + const auto gInside = problem.gravityAtPos(xInside); + + hInside -= rho*(gInside*xInside); - if (std::signbit(kGradPNormal_[phaseIdx])) + // and the outside neighbor + if (scvFace.boundary()) { - if (scvFace_().boundary()) - upWindIndices_[phaseIdx] = std::make_pair(-1, scvFace_().insideScvIdx()); - else - upWindIndices_[phaseIdx] = std::make_pair(scvFace_().outsideScvIdx(), scvFace_().insideScvIdx()); + const auto xOutside = scvFace.center(); + const auto gOutside = problem.gravityAtPos(xOutside); + hOutside -= rho*(gOutside*xOutside); } else { - if (scvFace_().boundary()) - upWindIndices_[phaseIdx] = std::make_pair(scvFace_().insideScvIdx(), -1); - else - upWindIndices_[phaseIdx] = std::make_pair(scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()); + const auto outsideScvIdx = scvFace.outsideScvIdx(); + const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx); + const auto xOutside = outsideScv.center(); + const auto gOutside = problem.gravityAtPos(xOutside); + hOutside -= rho*(gOutside*xOutside); } } - } - - /*! - * \brief A function to calculate the mass flux over a sub control volume face - * - * \param phaseIdx The index of the phase of which the flux is to be calculated - * \param upwindFunction A function which does the upwinding - */ - template - Scalar flux(IndexType phaseIdx, FunctionType upwindFunction) const - { - return kGradPNormal_[phaseIdx]*upwindFunction(upVolVars(phaseIdx), dnVolVars(phaseIdx)); - } - const std::set& stencil() const - { - return stencil_; + return tij*(hInside - hOutside); } - const VolumeVariables& upVolVars(IndexType phaseIdx) const + static std::vector stencil(const SubControlVolumeFace& scvFace) { - if(upWindIndices_[phaseIdx].first != -1) - return problem_().model().curVolVars(upWindIndices_[phaseIdx].first); + std::vector stencil; + stencil.clear(); + if (!scvFace.boundary()) + { + stencil.push_back(scvFace.insideScvIdx()); + stencil.push_back(scvFace.outsideScvIdx()); + } else - return *boundaryVolVars_; - } + stencil.push_back(scvFace.insideScvIdx()); - const VolumeVariables& dnVolVars(IndexType phaseIdx) const - { - if(upWindIndices_[phaseIdx].second != -1) - return problem_().model().curVolVars(upWindIndices_[phaseIdx].second); - else - return *boundaryVolVars_; + return stencil; } private: - void updateTransmissibilities_() + static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace) { - if (!scvFace_().boundary()) - { - const auto insideScvIdx = scvFace_().insideScvIdx(); - const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto insideK = problem_().spatialParams().intrinsicPermeability(insideScv); - Scalar ti = calculateOmega_(insideK, insideScv); + Scalar tij; + + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv); + Scalar ti = calculateOmega_(problem, scvFace, insideK, insideScv); - const auto outsideScvIdx = scvFace_().outsideScvIdx(); - const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx); - const auto outsideK = problem_().spatialParams().intrinsicPermeability(outsideScv); - Scalar tj = -1.0*calculateOmega_(outsideK, outsideScv); + if (!scvFace.boundary()) + { + const auto outsideScvIdx = scvFace.outsideScvIdx(); + const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx); + const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv); + Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideScv); - tij_ = scvFace_().area()*(ti * tj)/(ti + tj); + tij = scvFace.area()*(ti * tj)/(ti + tj); } else { - const auto insideScvIdx = scvFace_().insideScvIdx(); - const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto insideK = problem_().spatialParams().intrinsicPermeability(insideScv); - Scalar ti = calculateOmega_(insideK, insideScv); - - tij_ = scvFace_().area()*ti; + tij = scvFace.area()*ti; } - } - void updateStencil_() - { - // fill the stencil - if (!scvFace_().boundary()) - stencil_.insert({scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()}); - else - // fill the stencil - stencil_.insert(scvFace_().insideScvIdx()); + return tij; } - Scalar calculateOmega_(const DimWorldMatrix &K, const SubControlVolume &scv) const + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &K, const SubControlVolume &scv) { GlobalPosition Knormal; - K.mv(scvFace_().unitOuterNormal(), Knormal); + K.mv(scvFace.unitOuterNormal(), Knormal); - auto distanceVector = scvFace_().center(); + auto distanceVector = scvFace.center(); distanceVector -= scv.center(); distanceVector /= distanceVector.two_norm2(); Scalar omega = Knormal * distanceVector; - omega *= problem_().model().curVolVars(scv).extrusionFactor(); + omega *= problem.model().curVolVars(scv).extrusionFactor(); return omega; } - Scalar calculateOmega_(Scalar K, const SubControlVolume &scv) const + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar K, const SubControlVolume &scv) { - auto distanceVector = scvFace_().center(); + auto distanceVector = scvFace.center(); distanceVector -= scv.center(); distanceVector /= distanceVector.two_norm2(); - Scalar omega = K * (distanceVector * scvFace_().unitOuterNormal()); - omega *= problem_().model().curVolVars(scv).extrusionFactor(); + Scalar omega = K * (distanceVector * scvFace.unitOuterNormal()); + omega *= problem.model().curVolVars(scv).extrusionFactor(); return omega; } - - const Problem &problem_() const - { - return *problemPtr_; - } - - const SubControlVolumeFace& scvFace_() const - { - return *scvFacePtr_; - } - - const Problem *problemPtr_; //! Pointer to the problem - const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created - bool enableGravity_; //! If we have a problem considering gravitational effects - std::set stencil_; //! Indices of the cells of which the pressure is needed for the flux calculation - - //! Boundary volume variables (they only get updated on Dirichlet boundaries) - VolumeVariables* boundaryVolVars_; - - //! The upstream (first) and downstream (second) volume variable indices - std::array, numPhases> upWindIndices_; - - //! Precomputed values - Scalar tij_; //! transmissibility for the flux calculation tij(ui - uj) - std::array kGradPNormal_; //! K(grad(p) - rho*g)*n }; } // end namespace -- GitLab From 768c765efddf56119dc55671105c26858617d2a6 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:06:38 +0200 Subject: [PATCH 06/46] [implicit][tpfa][fickslaw] remove private variables Fick's law has no more private variables, which reduces the memory consumption of the program substantially. The member functions are now static. --- .../implicit/cellcentered/tpfa/fickslaw.hh | 155 ++++++------------ 1 file changed, 51 insertions(+), 104 deletions(-) diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh index 0636512688..841f018282 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh @@ -71,56 +71,35 @@ class CCTpfaFicksLaw public: - void update(const Problem& problem, - const Element& element, - const SubControlVolumeFace &scvFace, - int phaseIdx, int compIdx) - { - problemPtr_ = &problem; - scvFacePtr_ = &scvFace; - - phaseIdx_ = phaseIdx; - compIdx_ = compIdx; - - updateStencil_(); - - // TODO for non solution dependent diffusion tensors... - } - - void update(const Problem &problem, - const Element &element, - const SubControlVolumeFace &scvFace, - int phaseIdx, int compIdx, - VolumeVariables* boundaryVolVars) - { - boundaryVolVars_ = boundaryVolVars; - update(problem, element, scvFace, phaseIdx, compIdx); - } - - void beginFluxComputation(bool boundaryVolVarsUpdated = false) + static Scalar flux(const Problem& problem, + const SubControlVolumeFace& scvFace, + const int phaseIdx_, + const int compIdx_, + VolumeVariables* boundaryVolVars, + bool boundaryVolVarsUpdated) { // diffusion tensors are always solution dependent - updateTransmissibilities_(); + Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx_, compIdx_); // Get the inside volume variables - const auto insideScvIdx = scvFace_().insideScvIdx(); - const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto* insideVolVars = &problem_().model().curVolVars(insideScv); + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto* insideVolVars = &problem.model().curVolVars(insideScv); // and the outside volume variables const VolumeVariables* outsideVolVars; - if (!scvFace_().boundary()) - outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx()); + if (!scvFace.boundary()) + outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); else { - outsideVolVars = boundaryVolVars_; if (!boundaryVolVarsUpdated) { // update the boudary volvars for Dirichlet boundaries - const auto element = problem_().model().fvGeometries().element(insideScv); - const auto dirichletPriVars = problem_().dirichlet(element, scvFace_()); - boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv); + const auto element = problem.model().fvGeometries().element(insideScv); + const auto dirichletPriVars = problem.dirichlet(element, scvFace); + boundaryVolVars->update(dirichletPriVars, problem, element, insideScv); } + outsideVolVars = boundaryVolVars; } // compute the diffusive flux @@ -128,117 +107,85 @@ public: const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_); const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_)); - rhoDGradXNormal_ = rho*tij_*(xInside - xOutside); + return rho*tij*(xInside - xOutside); } - - - /*! - * \brief A function to calculate the mass flux over a sub control volume face - * - * \param phaseIdx The index of the phase of which the flux is to be calculated - * \param compIdx The index of the transported component - */ - Scalar flux() const + static std::vector stencil(const SubControlVolumeFace& scvFace) { - return rhoDGradXNormal_; - } + std::vector stencil; + stencil.clear(); + if (!scvFace.boundary()) + { + stencil.push_back(scvFace.insideScvIdx()); + stencil.push_back(scvFace.outsideScvIdx()); + } + else + stencil.push_back(scvFace.insideScvIdx()); - std::set stencil() const - { - return stencil_; + return stencil; } protected: - void updateTransmissibilities_() + static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx_, const int compIdx_) { - const auto insideScvIdx = scvFace_().insideScvIdx(); - const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto& insideVolVars = problem_().model().curVolVars(insideScvIdx); + Scalar tij; + + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto& insideVolVars = problem.model().curVolVars(insideScvIdx); auto insideD = insideVolVars.diffusionCoefficient(phaseIdx_, compIdx_); insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx_), insideD); - Scalar ti = calculateOmega_(insideD, insideScv); + Scalar ti = calculateOmega_(problem, scvFace, insideD, insideScv); - if (!scvFace_().boundary()) + if (!scvFace.boundary()) { - const auto outsideScvIdx = scvFace_().outsideScvIdx(); - const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx); - const auto& outsideVolVars = problem_().model().curVolVars(outsideScvIdx); + const auto outsideScvIdx = scvFace.outsideScvIdx(); + const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx); + const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx); auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx_, compIdx_); outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx_), outsideD); - Scalar tj = -1.0*calculateOmega_(outsideD, outsideScv); + Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideD, outsideScv); - tij_ = scvFace_().area()*(ti * tj)/(ti + tj); + tij = scvFace.area()*(ti * tj)/(ti + tj); } else { - tij_ = scvFace_().area()*ti; + tij = scvFace.area()*ti; } - } - void updateStencil_() - { - // fill the stencil - if (!scvFace_().boundary()) - stencil_= {scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()}; - else - // fill the stencil - stencil_ = {scvFace_().insideScvIdx()}; + return tij; } - Scalar calculateOmega_(const DimWorldMatrix &D, const SubControlVolume &scv) const + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv) { GlobalPosition Dnormal; - D.mv(scvFace_().unitOuterNormal(), Dnormal); + D.mv(scvFace.unitOuterNormal(), Dnormal); - auto distanceVector = scvFace_().center(); + auto distanceVector = scvFace.center(); distanceVector -= scv.center(); distanceVector /= distanceVector.two_norm2(); Scalar omega = Dnormal * distanceVector; - omega *= problem_().model().curVolVars(scv).extrusionFactor(); + omega *= problem.model().curVolVars(scv).extrusionFactor(); return omega; } - Scalar calculateOmega_(Scalar D, const SubControlVolume &scv) const + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv) { - auto distanceVector = scvFace_().center(); + auto distanceVector = scvFace.center(); distanceVector -= scv.center(); distanceVector /= distanceVector.two_norm2(); - Scalar omega = D * (distanceVector * scvFace_().unitOuterNormal()); - omega *= problem_().model().curVolVars(scv).extrusionFactor(); + Scalar omega = D * (distanceVector * scvFace.unitOuterNormal()); + omega *= problem.model().curVolVars(scv).extrusionFactor(); return omega; } - - const Problem &problem_() const - { - return *problemPtr_; - } - - const SubControlVolumeFace& scvFace_() const - { - return *scvFacePtr_; - } - - const Problem *problemPtr_; - const SubControlVolumeFace *scvFacePtr_; //!< Pointer to the sub control volume face for which the flux variables are created - std::set stencil_; //!< Indices of the cells of which the pressure is needed for the flux calculation - - //! Boundary volume variables (they only get updated on Dirichlet boundaries) - VolumeVariables* boundaryVolVars_; - - IndexType phaseIdx_; - IndexType compIdx_; - //! Precomputed values - Scalar tij_; //!< transmissibility for the flux calculation - Scalar rhoDGradXNormal_; //! rho*D(grad(x))*n }; } // end namespace -- GitLab From 820584e1d4e8cf6e26d721e60f9f8649038c08e5 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:08:30 +0200 Subject: [PATCH 07/46] [3p3c][localresidual] adapt to new flux var layout --- .../porousmediumflow/3p3c/implicit/localresidual.hh | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh index ffe92cbabc..a7c8ebfa7c 100644 --- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh +++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh @@ -141,17 +141,17 @@ public: // get equation index auto eqIdx = conti0EqIdx + compIdx; - flux[eqIdx] += fluxVars.advection().flux(phaseIdx, upwindRule); + flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindRule); } } // diffusive fluxes - Scalar jGW = fluxVars.molecularDiffusion(wPhaseIdx, gCompIdx).flux(); - Scalar jNW = fluxVars.molecularDiffusion(wPhaseIdx, nCompIdx).flux(); + Scalar jGW = fluxVars.molecularDiffusionFlux(wPhaseIdx, gCompIdx); + Scalar jNW = fluxVars.molecularDiffusionFlux(wPhaseIdx, nCompIdx); Scalar jWW = -(jGW+jNW); - Scalar jWG = fluxVars.molecularDiffusion(gPhaseIdx, wCompIdx).flux(); - Scalar jNG = fluxVars.molecularDiffusion(gPhaseIdx, nCompIdx).flux(); + Scalar jWG = fluxVars.molecularDiffusionFlux(gPhaseIdx, wCompIdx); + Scalar jNG = fluxVars.molecularDiffusionFlux(gPhaseIdx, nCompIdx); Scalar jGG = -(jWG+jNG); // Scalar jWN = fluxVars.molecularDiffusion().flux(nPhaseIdx, wCompIdx); @@ -167,6 +167,8 @@ public: flux[conti1EqIdx] += jNW+jNG+jNN; flux[conti2EqIdx] += jGW+jGG+jGN; + fluxVars.endFluxComputation(); + return flux; } -- GitLab From a70f20e872ebc548aaeacec3e62b56be946b32f0 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:09:14 +0200 Subject: [PATCH 08/46] [immiscible][localresidual] adapt to new fluxvar layout --- dumux/porousmediumflow/immiscible/localresidual.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh index 2fe3347df6..eaa9862f80 100644 --- a/dumux/porousmediumflow/immiscible/localresidual.hh +++ b/dumux/porousmediumflow/immiscible/localresidual.hh @@ -112,8 +112,9 @@ public: + (dn.density(phaseIdx)*dn.mobility(phaseIdx))*(1-w); }; auto eqIdx = conti0EqIdx + phaseIdx; - flux[eqIdx] = fluxVars.advection().flux(phaseIdx, upwindRule); + flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule); } + fluxVars.endFluxComputation(); return flux; } -- GitLab From 80151b6c3ff3c05ea4ae5aaf3131bfe4f26a974d Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 12 May 2016 17:10:35 +0200 Subject: [PATCH 09/46] [implicit][fluxvars] adapt to static constitutive laws --- dumux/implicit/fluxvariables.hh | 154 ++++++++++++++++++-------------- 1 file changed, 88 insertions(+), 66 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 5c144517fb..f9caa53b9f 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -60,39 +60,66 @@ class FluxVariables public: void update(const Problem& problem, const Element& element, - const SubControlVolumeFace &scvf) + const SubControlVolumeFace &scvFace) { - if (scvf.boundary()) + problemPtr_ = &problem; + scvFacePtr_ = &scvFace; + + if (scvFace.boundary()) { if(!boundaryVolVars_) boundaryVolVars_ = Dune::Std::make_unique(); - advection_.update(problem, element, scvf, boundaryVolVars_.get()); - } - else - { - advection_.update(problem, element, scvf); } + + stencil_ = AdvectionType::stencil(*scvFacePtr_); } void beginFluxComputation() { - advection_.beginFluxComputation(); + // TODO if constant BC, do not set to false + boundaryVolVarsUpdated_ = false; } - const AdvectionType& advection() const + void endFluxComputation() {} + + template + Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - return advection_; + Scalar flux = AdvectionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + + const auto* insideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->insideScvIdx()); + const VolumeVariables* outsideVolVars; + + if (scvFacePtr_->boundary()) + outsideVolVars = boundaryVolVars_.get(); + else + outsideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->outsideScvIdx()); + + if (std::signbit(flux)) + return flux*upwindFunction(*outsideVolVars, *insideVolVars); + else + return flux*upwindFunction(*insideVolVars, *outsideVolVars); + + // we are sure the boundary vol vars have been updated at this point (if existing) + boundaryVolVarsUpdated_ = true; + + return flux; } const Stencil& stencil() const { - return advection_.stencil(); + return stencil_; } private: - AdvectionType advection_; + const Problem *problemPtr_; //! Pointer to the problem + const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created + // boundary volume variables in case of Dirichlet boundaries std::unique_ptr boundaryVolVars_; + bool boundaryVolVarsUpdated_; + // the flux stencil + Stencil stencil_; }; @@ -118,93 +145,88 @@ class FluxVariables }; public: + + FluxVariables() + {} + void update(const Problem& problem, const Element& element, - const SubControlVolumeFace &scvf) + const SubControlVolumeFace &scvFace) { - if (scvf.boundary()) + problemPtr_ = &problem; + scvFacePtr_ = &scvFace; + + if (scvFace.boundary()) { if(!boundaryVolVars_) boundaryVolVars_ = Dune::Std::make_unique(); - - advection_.update(problem, element, scvf, boundaryVolVars_.get()); - for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - for (int compIdx = 0; compIdx < numComponents; ++compIdx) - { - if (phaseIdx != compIdx) - molecularDiffusion(phaseIdx, compIdx).update(problem, element, scvf, phaseIdx, compIdx, boundaryVolVars_.get()); - } - } - else - { - advection_.update(problem, element, scvf); - for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - for (int compIdx = 0; compIdx < numComponents; ++compIdx) - { - if (phaseIdx != compIdx) - molecularDiffusion(phaseIdx, compIdx).update(problem, element, scvf, phaseIdx, compIdx); - } } + + boundaryVolVarsUpdated_ = false; + + // TODO compute stencil as unity of advective and diffusive flux stencils + stencil_ = AdvectionType::stencil(*scvFacePtr_); } void beginFluxComputation() { - advection_.beginFluxComputation(); + // TODO if constant BC, do not set to false + boundaryVolVarsUpdated_ = false; + + // calculate advective fluxes and store them + advectiveVolFluxes_ = new std::array(); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - for (int compIdx = 0; compIdx < numComponents; ++compIdx) - { - if (phaseIdx != compIdx) - molecularDiffusion(phaseIdx, compIdx).beginFluxComputation(true); - } + (*advectiveVolFluxes_)[phaseIdx] = AdvectionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + boundaryVolVarsUpdated_ = true; } } - // TODO update transmissibilities? - - const AdvectionType& advection() const + void endFluxComputation() { - return advection_; + delete advectiveVolFluxes_; } - AdvectionType& advection() + template + Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - return advection_; - } + const auto* insideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->insideScvIdx()); + const VolumeVariables* outsideVolVars; - const MolecularDiffusionType& molecularDiffusion(const int phaseIdx, const int compIdx) const - { - if (compIdx < phaseIdx) - return molecularDiffusion_[phaseIdx][compIdx]; - else if (compIdx > phaseIdx) - return molecularDiffusion_[phaseIdx][compIdx-1]; + if (scvFacePtr_->boundary()) + outsideVolVars = boundaryVolVars_.get(); + else + outsideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->outsideScvIdx()); + + if (std::signbit((*advectiveVolFluxes_)[phaseIdx])) + return (*advectiveVolFluxes_)[phaseIdx]*upwindFunction(*outsideVolVars, *insideVolVars); else - DUNE_THROW(Dune::InvalidStateException, "Diffusion flux called for phaseIdx = compIdx"); + return (*advectiveVolFluxes_)[phaseIdx]*upwindFunction(*insideVolVars, *outsideVolVars); } - MolecularDiffusionType& molecularDiffusion(const int phaseIdx, const int compIdx) + Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - if (compIdx < phaseIdx) - return molecularDiffusion_[phaseIdx][compIdx]; - else if (compIdx > phaseIdx) - return molecularDiffusion_[phaseIdx][compIdx-1]; - else - DUNE_THROW(Dune::InvalidStateException, "Diffusion flux called for phaseIdx = compIdx"); + Scalar flux = MolecularDiffusionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, compIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + boundaryVolVarsUpdated_ = true; + return flux; } - Stencil stencil() const + const Stencil& stencil() const { - // std::set stencil(advection().stencil().begin(), advection().stencil().end()); - // TODO: insert all molecularDiffusion stencils of all components in all phases - // stencil.insert(); - return advection().stencil(); + return stencil_; } private: - AdvectionType advection_; - std::array< std::array, numPhases> molecularDiffusion_; + const Problem *problemPtr_; //! Pointer to the problem + const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created + + // storage for calculated advective fluxes to not having to calculate them again + std::array* advectiveVolFluxes_; // boundary volume variables in case of Dirichlet boundaries std::unique_ptr boundaryVolVars_; + bool boundaryVolVarsUpdated_; + // the flux stencil + Stencil stencil_; }; -- GitLab From a8445f51d6c5f220df883a94e81023c665af55f2 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 13 May 2016 09:22:19 +0200 Subject: [PATCH 10/46] [implicit][cc][stencils] use std::remove to eliminate own element index from stencil --- dumux/implicit/cellcentered/stencils.hh | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index c70bfdbdf1..44389eb041 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -54,15 +54,9 @@ public: std::sort(elementStencil_.begin(), elementStencil_.end()); elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end()); + auto globalI = problem.elementMapper().index(element); neighborStencil_ = elementStencil_; - for (auto it = neighborStencil_.begin(); it != neighborStencil_.end(); ++it) - { - if (*it == problem.elementMapper().index(element)) - { - neighborStencil_.erase(it); - break; - } - } + neighborStencil_.erase(std::remove(neighborStencil_.begin(), neighborStencil_.end(), globalI), neighborStencil_.end()); } //! The full element stencil (all element this element is interacting with) -- GitLab From 0a9b485f1e8b686783a660e30221e7391b11caa9 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 13 May 2016 09:30:06 +0200 Subject: [PATCH 11/46] [implicit][fluxvars] use return function instead of dereferencing pointers --- dumux/implicit/fluxvariables.hh | 42 ++++++++++++++++++++++++--------- 1 file changed, 31 insertions(+), 11 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index f9caa53b9f..8834a5cb28 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -71,7 +71,7 @@ public: boundaryVolVars_ = Dune::Std::make_unique(); } - stencil_ = AdvectionType::stencil(*scvFacePtr_); + stencil_ = AdvectionType::stencil(scvFace); } void beginFluxComputation() @@ -85,15 +85,15 @@ public: template Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - Scalar flux = AdvectionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + Scalar flux = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); - const auto* insideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->insideScvIdx()); + const auto* insideVolVars = &problem().model().curVolVars(scvFace().insideScvIdx()); const VolumeVariables* outsideVolVars; - if (scvFacePtr_->boundary()) + if (scvFace().boundary()) outsideVolVars = boundaryVolVars_.get(); else - outsideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->outsideScvIdx()); + outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); if (std::signbit(flux)) return flux*upwindFunction(*outsideVolVars, *insideVolVars); @@ -111,6 +111,16 @@ public: return stencil_; } + const Problem& problem() const + { + return *problemPtr_; + } + + const SubControlVolumeFace& scvFace() const + { + return *scvFacePtr_; + } + private: const Problem *problemPtr_; //! Pointer to the problem const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created @@ -165,7 +175,7 @@ public: boundaryVolVarsUpdated_ = false; // TODO compute stencil as unity of advective and diffusive flux stencils - stencil_ = AdvectionType::stencil(*scvFacePtr_); + stencil_ = AdvectionType::stencil(scvFace); } void beginFluxComputation() @@ -177,7 +187,7 @@ public: advectiveVolFluxes_ = new std::array(); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - (*advectiveVolFluxes_)[phaseIdx] = AdvectionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + (*advectiveVolFluxes_)[phaseIdx] = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); boundaryVolVarsUpdated_ = true; } } @@ -190,13 +200,13 @@ public: template Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - const auto* insideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->insideScvIdx()); + const auto* insideVolVars = &problem().model().curVolVars(scvFace().insideScvIdx()); const VolumeVariables* outsideVolVars; - if (scvFacePtr_->boundary()) + if (scvFace().boundary()) outsideVolVars = boundaryVolVars_.get(); else - outsideVolVars = &problemPtr_->model().curVolVars(scvFacePtr_->outsideScvIdx()); + outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); if (std::signbit((*advectiveVolFluxes_)[phaseIdx])) return (*advectiveVolFluxes_)[phaseIdx]*upwindFunction(*outsideVolVars, *insideVolVars); @@ -206,7 +216,7 @@ public: Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - Scalar flux = MolecularDiffusionType::flux(*problemPtr_, *scvFacePtr_, phaseIdx, compIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + Scalar flux = MolecularDiffusionType::flux(problem(), scvFace(), phaseIdx, compIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); boundaryVolVarsUpdated_ = true; return flux; } @@ -216,6 +226,16 @@ public: return stencil_; } + const Problem& problem() const + { + return *problemPtr_; + } + + const SubControlVolumeFace& scvFace() const + { + return *scvFacePtr_; + } + private: const Problem *problemPtr_; //! Pointer to the problem const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created -- GitLab From f31f4f133dda8b92b8433d6658cf3af173adeb60 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 13 May 2016 09:55:00 +0200 Subject: [PATCH 12/46] [implicit][tpfa][fickslaw] make protected member functions private --- dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh index 841f018282..f1740637ef 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh @@ -125,7 +125,7 @@ public: return stencil; } -protected: +private: static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx_, const int compIdx_) -- GitLab From edc401d14a2cc87691010eacbee576a29305005d Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 13 May 2016 10:09:27 +0200 Subject: [PATCH 13/46] [implicit][scvface] fix comment on outsidescvindex() --- dumux/implicit/subcontrolvolumeface.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/implicit/subcontrolvolumeface.hh b/dumux/implicit/subcontrolvolumeface.hh index ed2f9292eb..d9261e819a 100644 --- a/dumux/implicit/subcontrolvolumeface.hh +++ b/dumux/implicit/subcontrolvolumeface.hh @@ -109,7 +109,7 @@ public: } //! index of the outside sub control volume for spatial param evaluation - // This results in undefined behaviour if boundary is false + // This results in undefined behaviour if boundary is true IndexType outsideScvIdx() const { return scvIndices_[1]; -- GitLab From 44db3bd3ee189e82db5eb78a8e160ac9ef86bb01 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 13:44:24 +0200 Subject: [PATCH 14/46] [implicit] introduce flux var cache The flux variables do not exist as a global vector in the model anymore. Instead, a cache class can be saved per scv face in order to store data that should be recycled. This can be switched on via the property EnableFluxVariablesCache and leads to a better performance on small grids. On large grids memory usage becomes an issue and it should be switched off. If the property ConstantBoundaryConditions is switched on, the volume variables at the boundary, used by cc models, will be stored in the cache and recomputation does not occur. --- dumux/implicit/cellcentered/localjacobian.hh | 4 +- dumux/implicit/cellcentered/stencils.hh | 6 +- dumux/implicit/fluxvariablescache.hh | 121 ++++++++++++++++++ dumux/implicit/fluxvariablescachevector.hh | 76 +++++++++++ dumux/implicit/model.hh | 30 ++--- dumux/implicit/properties.hh | 7 +- dumux/implicit/propertydefaults.hh | 22 +++- .../immiscible/localresidual.hh | 5 +- 8 files changed, 245 insertions(+), 26 deletions(-) create mode 100644 dumux/implicit/fluxvariablescache.hh create mode 100644 dumux/implicit/fluxvariablescachevector.hh diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index f362b3d8ec..73903e36a2 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -80,6 +80,7 @@ class CCLocalJacobian : public ImplicitLocalJacobian typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; @@ -128,7 +129,8 @@ public: auto fluxVarsIdx = scvFaceJ.index(); // if globalI is in flux var stencil, add to list - const auto& fluxStencil = this->model_().fluxVars(fluxVarsIdx).stencil(); + FluxVariables fluxVars; + const auto& fluxStencil = fluxVars.stencil(problem, scvFaceJ); for (auto globalIdx : fluxStencil) if (globalIdx == globalI) fluxVarIndices.push_back(fluxVarsIdx); diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 44389eb041..27f6d05965 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -37,6 +37,7 @@ class CCElementStencils { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); using IndexType = typename GridView::IndexSet::IndexType; using Element = typename GridView::template Codim<0>::Entity; // TODO a separate stencil class all stencils can derive from? @@ -47,8 +48,9 @@ public: elementStencil_.clear(); for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) { - auto&& fluxStencil = problem.model().fluxVars(scvf).stencil(); - elementStencil_.insert(elementStencil_.end(), fluxStencil.begin(), fluxStencil.end()); + FluxVariables fluxVars; + auto&& stencil = fluxVars.stencil(problem, scvf); + elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); } // make values in elementstencil unique std::sort(elementStencil_.begin(), elementStencil_.end()); diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh new file mode 100644 index 0000000000..156b576a57 --- /dev/null +++ b/dumux/implicit/fluxvariablescache.hh @@ -0,0 +1,121 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Base class for the flux variables + */ +#ifndef DUMUX_IMPLICIT_FLUXVARIABLESCACHE_HH +#define DUMUX_IMPLICIT_FLUXVARIABLESCACHE_HH + +#include + +namespace Dumux +{ + +namespace Properties +{ +NEW_PROP_TAG(NumPhases); +NEW_PROP_TAG(NumComponents); +} + +/*! + * \ingroup ImplicitModel + * \brief The flux variables cache class + * stores the transmissibilities and stencils for the darcy flux calculation on a scv face + */ +template +class FluxVariablesCache +{}; + +// specialization for the Box Method +template +class FluxVariablesCache +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Stencil = std::vector; + using TransmissibilityVector = std::vector; + +public: + void update(const Problem& problem, + const Element& element, + const SubControlVolumeFace &scvFace) + { + darcyStencil_ = AdvectionType::stencil(problem, scvFace); + tij_ = AdvectionType::calculateTransmissibilities(problem, scvFace); + } + + const Stencil& stencil() const + { return darcyStencil_; } + + const Stencil& darcyStencil() const + { return darcyStencil_; } + + const TransmissibilityVector& tij() const + { return tij_; } + +private: + Stencil darcyStencil_; + TransmissibilityVector tij_; +}; + +// specialization for cell centered methods +template +class FluxVariablesCache : public FluxVariablesCache +{ + using ParentType = FluxVariablesCache; + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Stencil = std::vector; + +public: + void update(const Problem& problem, + const Element& element, + const SubControlVolumeFace &scvFace) + { + ParentType::update(problem, element, scvFace); + FluxVariables fluxVars; + stencil_ = fluxVars.computeFluxStencil(problem, scvFace); + boundaryVolVars_ = std::make_unique(std::move(fluxVars.getBoundaryVolumeVariables(problem, element, scvFace))); + } + + const Stencil& stencil() const + { return stencil_; } + + const VolumeVariables& boundaryVolumeVariables() const + { return *boundaryVolVars_.get(); } + +private: + Stencil stencil_; + std::unique_ptr boundaryVolVars_; +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/fluxvariablescachevector.hh b/dumux/implicit/fluxvariablescachevector.hh new file mode 100644 index 0000000000..159768f1a8 --- /dev/null +++ b/dumux/implicit/fluxvariablescachevector.hh @@ -0,0 +1,76 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Base class for the volume variables vector + */ +#ifndef DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH +#define DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH + +#include + +namespace Dumux +{ + +/*! + * \ingroup ImplicitModel + * \brief Base class for the flux variables cache vector, we store one cache per face + */ +template +class FluxVariablesCacheVector +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using IndexType = typename GridView::IndexSet::IndexType; + using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + +public: + template + typename std::enable_if::type update(const Problem& problem) + { + fluxVarsCache_.resize(problem.model().fvGeometries().numScvf()); + for (const auto& element : elements(problem.gridView())) + { + for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) + { + (*this)[scvf.index()].update(problem, element, scvf); + } + } + } + + template + typename std::enable_if::type update(const Problem& problem) + {} + + template + const typename std::enable_if::type& operator [](IndexType scvfIdx) const + { return fluxVarsCache_[scvfIdx]; } + + template + typename std::enable_if::type& operator [](IndexType scvfIdx) + { return fluxVarsCache_[scvfIdx]; } + +private: + std::vector fluxVarsCache_; +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 1979e3eaa8..5b209f2b41 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -57,8 +57,8 @@ class ImplicitModel typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; - typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; - typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesVector) FluxVariablesVector; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCacheVector) FluxVariablesCacheVector; typedef typename GET_PROP_TYPE(TypeTag, StencilsVector) StencilsVector; enum { @@ -119,8 +119,8 @@ public: curVolVarsVector_.resize(fvGeometries().numScv()); curVolVarsVector_.update(problem_(), curSol()); - // update the flux vars (precompute transmissibilities) - fluxVarsVector_.update(problem_()); + // update the flux variables caches + fluxVarsCacheVector_.update(problem_()); // update stencils stencilsVector_.update(problem_()); @@ -740,7 +740,7 @@ public: */ bool onBoundary(const SubControlVolume &scv) const { - return asImp_().onBoundary(onBoundary(scv.dofIndex())); + return asImp_().onBoundary(scv.dofIndex()); } /*! @@ -769,11 +769,11 @@ public: // fvGeometry, scvIdx, fluidState); // } - const FluxVariables& fluxVars(const SubControlVolumeFace& scvf) const - { return fluxVarsVector_[scvf.index()]; } + const FluxVariablesCache& fluxVarsCache(const SubControlVolumeFace& scvf) const + { return fluxVarsCacheVector_[scvf.index()]; } - const FluxVariables& fluxVars(unsigned int scvfIdx) const - { return fluxVarsVector_[scvfIdx]; } + const FluxVariablesCache& fluxVarsCache(unsigned int scvfIdx) const + { return fluxVarsCacheVector_[scvfIdx]; } const VolumeVariables& curVolVars(const SubControlVolume& scv) const { return curVolVarsVector_[scv.index()]; } @@ -814,11 +814,11 @@ protected: VolumeVariables& prevVolVars_(unsigned int scvIdx) { return prevVolVarsVector_[scvIdx]; } - FluxVariables& fluxVars_(const SubControlVolumeFace& scvf) - { return fluxVarsVector_[scvf.index()]; } + FluxVariablesCache& fluxVarsCache_(const SubControlVolumeFace& scvf) + { return fluxVarsCacheVector_[scvf.index()]; } - FluxVariables& fluxVars_(unsigned int scvfIdx) - { return fluxVarsVector_[scvfIdx]; } + FluxVariablesCache& fluxVarsCache_(unsigned int scvfIdx) + { return fluxVarsCacheVector_[scvfIdx]; } /*! * \brief A reference to the problem on which the model is applied. @@ -971,8 +971,8 @@ protected: VolumeVariablesVector curVolVarsVector_; VolumeVariablesVector prevVolVarsVector_; - // the flux variables vector - FluxVariablesVector fluxVarsVector_; + // the flux variables cache vector vector + FluxVariablesCacheVector fluxVarsCacheVector_; // the stencils vector StencilsVector stencilsVector_; diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index 41e22f039f..2666ca42e5 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -74,11 +74,14 @@ NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a NEW_PROP_TAG(VolumeVariables); //!< The secondary variables within a sub-control volume NEW_PROP_TAG(VolumeVariablesVector); //!< The type for a container of volume variables NEW_PROP_TAG(FluxVariables); //!< Container storing the different types of flux variables -NEW_PROP_TAG(FluxVariablesVector); //!< The global vector of flux variable containers -NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces (outflow) +NEW_PROP_TAG(FluxVariablesCache); //!< Stores data associated with flux vars (if enabled) +NEW_PROP_TAG(FluxVariablesCacheVector); //!< The global vector of flux variable containers +NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces in cc models(outflow) +NEW_PROP_TAG(ConstantBoundaryConditions); //!< boundary data is stored in case the BC are constant // Specify the forms of fluxes that should be considered in the model // also, specify their corresponding flux variables +NEW_PROP_TAG(EnableFluxVariablesCache); //! specifies if data on flux vars should be saved (faster, but more memory consuming) NEW_PROP_TAG(EnableAdvection); //! specifies if advection is considered in the model NEW_PROP_TAG(AdvectionType); //! The type for the calculation the advective fluxes NEW_PROP_TAG(EnableMolecularDiffusion); //! specifies if molecular diffusive fluxes are considered in the model diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 4a413661a1..8d9af7d5b1 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -42,7 +42,8 @@ #include "volumevariables.hh" #include "volumevariablesvector.hh" #include "fluxvariables.hh" -#include "fluxvariablesvector.hh" +#include "fluxvariablescache.hh" +#include "fluxvariablescachevector.hh" #include "fvelementgeometry.hh" namespace Dumux { @@ -98,8 +99,17 @@ SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables); -//! The global volume variables vector class -SET_TYPE_PROP(ImplicitBase, FluxVariablesVector, Dumux::FluxVariablesVector); +//! The flux variables cache class +SET_PROP(ImplicitBase, FluxVariablesCache) +{ +private: + enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + public: + typedef typename Dumux::FluxVariablesCache type; +}; + +//! The global flux variables cache vector class +SET_TYPE_PROP(ImplicitBase, FluxVariablesCacheVector, Dumux::FluxVariablesCacheVector); // //! The class that contains the different flux variables (i.e. darcy, diffusion, energy) SET_PROP(ImplicitBase, FluxVariables) @@ -142,6 +152,12 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false); //! disable partial reassembling by default SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false); +//! disable flux variables data caching by default +SET_BOOL_PROP(ImplicitBase, EnableFluxVariablesCache, false); + +//! boundary conditions are not stationary by default +SET_BOOL_PROP(ImplicitBase, ConstantBoundaryConditions, false); + //! Set the type of a global jacobian matrix from the solution types SET_PROP(ImplicitBase, JacobianMatrix) { diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh index eaa9862f80..efe919230b 100644 --- a/dumux/porousmediumflow/immiscible/localresidual.hh +++ b/dumux/porousmediumflow/immiscible/localresidual.hh @@ -97,8 +97,8 @@ public: */ PrimaryVariables computeFlux(const SubControlVolumeFace& scvf) { - auto& fluxVars = this->model_().fluxVars_(scvf); - fluxVars.beginFluxComputation(); + FluxVariables fluxVars; + fluxVars.update(asImp_()->problem_(), this->element_(), scvf); // copy weight to local scope for use in lambda expression auto w = upwindWeight_; @@ -114,7 +114,6 @@ public: auto eqIdx = conti0EqIdx + phaseIdx; flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule); } - fluxVars.endFluxComputation(); return flux; } -- GitLab From f71c51ecc593955fb0bb8cf1a6057f7300611f12 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 13:49:32 +0200 Subject: [PATCH 15/46] [implicit][cc][stencils] remove obsolete comment characters --- dumux/implicit/cellcentered/stencils.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 27f6d05965..bc311716bd 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -67,7 +67,7 @@ public: return elementStencil_; } - //! //! The full element stencil without this element + //! The full element stencil without this element const Stencil& neighborStencil() const { return neighborStencil_; -- GitLab From ceba93af76bff1047bb84fefa1555e91049fd2a9 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 13:50:50 +0200 Subject: [PATCH 16/46] [implicit][cc][stencils] introduce typedef for stencil type --- dumux/implicit/cellcentered/stencils.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index bc311716bd..84d3915586 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -87,6 +87,7 @@ class CCStencilsVector { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using StencilType = CCElementStencils; public: void update(const Problem& problem) @@ -102,7 +103,7 @@ public: //! overload for elements template - typename std::enable_if&>::type + typename std::enable_if::type get(const Entity& entity) const { return elementStencils_[problemPtr_->elementMapper().index(entity)]; -- GitLab From 689bd2de8bf7526026b04469011d684cc89fdf31 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 13:52:46 +0200 Subject: [PATCH 17/46] [implicit][fluxvars] introduce cache handling --- dumux/implicit/fluxvariables.hh | 118 ++++++++++++++++++++++++++------ 1 file changed, 96 insertions(+), 22 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 8834a5cb28..7dee525e52 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -57,7 +57,17 @@ class FluxVariables using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); + enum + { + isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox), + enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache), + constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) + }; + public: + FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr), boundaryVolVars_(nullptr) + {} + void update(const Problem& problem, const Element& element, const SubControlVolumeFace &scvFace) @@ -65,33 +75,31 @@ public: problemPtr_ = &problem; scvFacePtr_ = &scvFace; - if (scvFace.boundary()) - { - if(!boundaryVolVars_) - boundaryVolVars_ = Dune::Std::make_unique(); - } + // boundary vol vars only need to be handled at boundaries for cc methods + if (scvFace.boundary() && !isBox) + setBoundaryVolumeVariables_(problem, element, scvFace); - stencil_ = AdvectionType::stencil(scvFace); + // update the stencil if needed + if (!enableFluxVarsCache && stencil_.empty()) + stencil_ = stencil(problem, scvFace); } - void beginFluxComputation() - { - // TODO if constant BC, do not set to false - boundaryVolVarsUpdated_ = false; - } - void endFluxComputation() {} template Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - Scalar flux = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); + Scalar flux = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_); const auto* insideVolVars = &problem().model().curVolVars(scvFace().insideScvIdx()); const VolumeVariables* outsideVolVars; if (scvFace().boundary()) - outsideVolVars = boundaryVolVars_.get(); + { + if (boundaryVolVars_ == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); + outsideVolVars = boundaryVolVars_; + } else outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); @@ -100,34 +108,100 @@ public: else return flux*upwindFunction(*insideVolVars, *outsideVolVars); - // we are sure the boundary vol vars have been updated at this point (if existing) - boundaryVolVarsUpdated_ = true; - return flux; } - const Stencil& stencil() const + // interface allowing for stencil information without having to update the flux vars. + // this becomes useful when the element is not at hand for a call to update(...), e.g. during the assembly - localjacobian.hh + template + const typename std::enable_if::type& stencil(const Problem& problem, const SubControlVolumeFace& scvFace) const + { return problem.model().fluxVarsCache(scvFace).stencil(); } + + // provide interface in case caching is disabled + template + const typename std::enable_if::type& stencil(const Problem& problem, const SubControlVolumeFace& scvFace) { + if (stencil_.empty()) + stencil_ = computeFluxStencil(problem, scvFace); return stencil_; } + // returns the boundary vol vars belonging to this flux variables object + VolumeVariables getBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + { + VolumeVariables boundaryVolVars; + + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto dirichletPriVars = problem.dirichlet(element, scvFace); + boundaryVolVars.update(dirichletPriVars, problem, element, insideScv); + + return boundaryVolVars; + } + const Problem& problem() const { + if (problemPtr_ == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Problem pointer not valid. Call update before using the flux variables class."); return *problemPtr_; } const SubControlVolumeFace& scvFace() const { + if (scvFacePtr_ == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Scv face pointer not valid. Call update before using the flux variables class."); return *scvFacePtr_; } + // when caching is enabled, get the stencil from the cache class + template + const typename std::enable_if::type& stencil() const + { + if (problemPtr_ == nullptr || scvFacePtr_ == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + return problem().model().fluxVarsCache(scvFace()).stencil(); + } + + // when caching is disabled, return the private stencil variable. The update(...) routine has to be called beforehand. + template + const typename std::enable_if::type& stencil() + { + if (stencil_.empty()) + DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + return stencil_; + } + + Stencil computeFluxStencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { return AdvectionType::stencil(problem, scvFace); } + private: + + // if flux variables caching is enabled, we use the boundary volume variables stored in the cache class + template + typename std::enable_if::type + setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + { + boundaryVolVars_ = &problem.model().fluxVarsCache(scvFace).boundaryVolumeVariables(); + } + + // if flux variables caching is disabled, we update the boundary volume variables + template + typename std::enable_if::type + setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + { + VolumeVariables tmp; + + const auto dirichletPriVars = problem.dirichlet(element, scvFace); + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + tmp.update(dirichletPriVars, problem, element, insideScv); + + boundaryVolVars_ = new VolumeVariables(std::move(tmp)); + } + const Problem *problemPtr_; //! Pointer to the problem const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created - - // boundary volume variables in case of Dirichlet boundaries - std::unique_ptr boundaryVolVars_; - bool boundaryVolVarsUpdated_; + const VolumeVariables* boundaryVolVars_; //! boundary volume variables in case of Dirichlet boundaries // the flux stencil Stencil stencil_; }; -- GitLab From 1ed84dd42727a2609f2fc1fc06646de744c514d9 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 13:54:01 +0200 Subject: [PATCH 18/46] [implicit][tpfa][darcyslaw] introduce cache handling --- .../implicit/cellcentered/tpfa/darcyslaw.hh | 43 +++++++++++-------- 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh index 6408ce9fab..11c2fb878d 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh @@ -53,16 +53,19 @@ NEW_PROP_TAG(ProblemEnableGravity); template class CCTpfaDarcysLaw { - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GridView::IndexSet::IndexType IndexType; - using Element = typename GridView::template Codim<0>::Entity; + typedef std::vector TransmissivityVector; + typedef std::vector Stencil; + using Element = typename GridView::template Codim<0>::Entity; enum { dim = GridView::dimension} ; enum { dimWorld = GridView::dimensionworld} ; enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ; @@ -75,10 +78,9 @@ public: static Scalar flux(const Problem& problem, const SubControlVolumeFace& scvFace, const IndexType phaseIdx, - VolumeVariables* boundaryVolVars, - const bool boundaryVolVarsUpdated) + const VolumeVariables* boundaryVolVars) { - Scalar tij = calculateTransmissibility_(problem, scvFace); + const auto& tij = getTransmissibilities(problem, scvFace); // Get the inside volume variables const auto insideScvIdx = scvFace.insideScvIdx(); @@ -91,13 +93,8 @@ public: outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); else { - if (!boundaryVolVarsUpdated) - { - // update the boudary volvars for Dirichlet boundaries - const auto element = problem.model().fvGeometries().element(insideScv); - const auto dirichletPriVars = problem.dirichlet(element, scvFace); - boundaryVolVars->update(dirichletPriVars, problem, element, insideScv); - } + if (boundaryVolVars == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); outsideVolVars = boundaryVolVars; } @@ -135,10 +132,10 @@ public: } } - return tij*(hInside - hOutside); + return tij[0]*(hInside - hOutside); } - static std::vector stencil(const SubControlVolumeFace& scvFace) + static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace) { std::vector stencil; stencil.clear(); @@ -153,9 +150,19 @@ public: return stencil; } -private: + template + static const typename std::enable_if::type& getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) + { + return problem.model().fluxVarsCache(scvFace).tij(); + } - static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace) + template + static const typename std::enable_if::type getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) + { + return calculateTransmissibilities(problem, scvFace); + } + + static TransmissivityVector calculateTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) { Scalar tij; @@ -178,9 +185,11 @@ private: tij = scvFace.area()*ti; } - return tij; + return TransmissivityVector(1, tij); } +private: + static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &K, const SubControlVolume &scv) { GlobalPosition Knormal; -- GitLab From 5b2d464da666700926514bc5aec7dfda3c1fc0a2 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 14:47:42 +0200 Subject: [PATCH 19/46] [implicit][fluxvariables] add destructor In case flux variable caching is enabled we need a delete statement for the previously allocated memory in the setBoundaryVolumeVariables_() routine. --- dumux/implicit/fluxvariables.hh | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 7dee525e52..6c3817e926 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -68,6 +68,13 @@ public: FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr), boundaryVolVars_(nullptr) {} + // if flux variables caching is disabled, we need to delete an eventually allocated boundaryVolVars_ pointer + ~FluxVariables() + { + if (boundaryVolVars_ != nullptr) + delete boundaryVolVars_; + } + void update(const Problem& problem, const Element& element, const SubControlVolumeFace &scvFace) -- GitLab From 787118f11d3f8a78c839ff6e8d0f9d4c6abf689e Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 15:03:05 +0200 Subject: [PATCH 20/46] [implicit][fluxvars] fix if statement in destructor --- dumux/implicit/fluxvariables.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 6c3817e926..a394b8abf3 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -68,10 +68,10 @@ public: FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr), boundaryVolVars_(nullptr) {} - // if flux variables caching is disabled, we need to delete an eventually allocated boundaryVolVars_ pointer + // if boundary volume variables have been allocated previously, we need to delete an allocated memory ~FluxVariables() { - if (boundaryVolVars_ != nullptr) + if (boundaryVolVars_ != nullptr && (!constantBC || !enableFluxVarsCache)) delete boundaryVolVars_; } -- GitLab From 80e537446d7fb66c25115eb897df1ef0575e917f Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 15:04:04 +0200 Subject: [PATCH 21/46] [implicit][fluxvarcache] only store boundary volvars when property on --- dumux/implicit/fluxvariablescache.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh index 156b576a57..efb0514d85 100644 --- a/dumux/implicit/fluxvariablescache.hh +++ b/dumux/implicit/fluxvariablescache.hh @@ -102,7 +102,8 @@ public: ParentType::update(problem, element, scvFace); FluxVariables fluxVars; stencil_ = fluxVars.computeFluxStencil(problem, scvFace); - boundaryVolVars_ = std::make_unique(std::move(fluxVars.getBoundaryVolumeVariables(problem, element, scvFace))); + if (GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions)) + boundaryVolVars_ = std::make_unique(std::move(fluxVars.getBoundaryVolumeVariables(problem, element, scvFace))); } const Stencil& stencil() const -- GitLab From e26736cf956b617d4d7a82841f50057d0d8b66c5 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 15:11:46 +0200 Subject: [PATCH 22/46] [implicit][fluxvars] make setboundaryvolvars() public we need this function to be public because the different specializations of the flux vars will now inherit from each other, so access needs to be guaranteed. --- dumux/implicit/fluxvariables.hh | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index a394b8abf3..a572170a75 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -84,7 +84,7 @@ public: // boundary vol vars only need to be handled at boundaries for cc methods if (scvFace.boundary() && !isBox) - setBoundaryVolumeVariables_(problem, element, scvFace); + setBoundaryVolumeVariables(problem, element, scvFace); // update the stencil if needed if (!enableFluxVarsCache && stencil_.empty()) @@ -181,12 +181,11 @@ public: Stencil computeFluxStencil(const Problem& problem, const SubControlVolumeFace& scvFace) { return AdvectionType::stencil(problem, scvFace); } -private: // if flux variables caching is enabled, we use the boundary volume variables stored in the cache class template typename std::enable_if::type - setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + setBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) { boundaryVolVars_ = &problem.model().fluxVarsCache(scvFace).boundaryVolumeVariables(); } @@ -194,7 +193,7 @@ private: // if flux variables caching is disabled, we update the boundary volume variables template typename std::enable_if::type - setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + setBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) { VolumeVariables tmp; -- GitLab From d00f6fc444932690b622dd3851694b0ff6b293a2 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 17:04:46 +0200 Subject: [PATCH 23/46] [implicit][fluxvars] make setBoundaryVolumeVariables_ function protected Previous commit making it public is not necessary since protected members can be acesses from the child class --- dumux/implicit/fluxvariables.hh | 46 ++++++++++----------------------- 1 file changed, 13 insertions(+), 33 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index a572170a75..fbee6fe9fb 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -84,7 +84,7 @@ public: // boundary vol vars only need to be handled at boundaries for cc methods if (scvFace.boundary() && !isBox) - setBoundaryVolumeVariables(problem, element, scvFace); + setBoundaryVolumeVariables_(problem, element, scvFace); // update the stencil if needed if (!enableFluxVarsCache && stencil_.empty()) @@ -160,63 +160,43 @@ public: return *scvFacePtr_; } - // when caching is enabled, get the stencil from the cache class - template - const typename std::enable_if::type& stencil() const - { - if (problemPtr_ == nullptr || scvFacePtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); - return problem().model().fluxVarsCache(scvFace()).stencil(); - } - - // when caching is disabled, return the private stencil variable. The update(...) routine has to be called beforehand. - template - const typename std::enable_if::type& stencil() - { - if (stencil_.empty()) - DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); - return stencil_; - } - Stencil computeFluxStencil(const Problem& problem, const SubControlVolumeFace& scvFace) { return AdvectionType::stencil(problem, scvFace); } +protected: // if flux variables caching is enabled, we use the boundary volume variables stored in the cache class template typename std::enable_if::type - setBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) { - boundaryVolVars_ = &problem.model().fluxVarsCache(scvFace).boundaryVolumeVariables(); + boundaryVolVars_ = std::shared_ptr(problem.model().fluxVarsCache(scvFace).boundaryVolumeVariables()); } // if flux variables caching is disabled, we update the boundary volume variables template typename std::enable_if::type - setBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) { - VolumeVariables tmp; + boundaryVolVars_ = std::make_shared(std::move(computeBoundaryVolumeVariables(problem, element, scvFace))); + } - const auto dirichletPriVars = problem.dirichlet(element, scvFace); - const auto insideScvIdx = scvFace.insideScvIdx(); - const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); - tmp.update(dirichletPriVars, problem, element, insideScv); + std::shared_ptr boundaryVolVars_; //! boundary volume variables in case of Dirichlet boundaries + // the flux stencil + Stencil stencil_; - boundaryVolVars_ = new VolumeVariables(std::move(tmp)); - } +private: const Problem *problemPtr_; //! Pointer to the problem const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created - const VolumeVariables* boundaryVolVars_; //! boundary volume variables in case of Dirichlet boundaries - // the flux stencil - Stencil stencil_; }; // specialization for isothermal advection molecularDiffusion equations template -class FluxVariables +class FluxVariables : public FluxVariables { + using ParentType = FluxVariables; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; -- GitLab From 62ad25332141a56c140c27799d0b07221411f120 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 17:07:44 +0200 Subject: [PATCH 24/46] [implicit][fluxvars] make boundaryvolvars a shared pointer --- dumux/implicit/fluxvariables.hh | 120 +++++++++++++------------------- 1 file changed, 48 insertions(+), 72 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index fbee6fe9fb..a35c9483f7 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -65,16 +65,9 @@ class FluxVariables }; public: - FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr), boundaryVolVars_(nullptr) + FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr) {} - // if boundary volume variables have been allocated previously, we need to delete an allocated memory - ~FluxVariables() - { - if (boundaryVolVars_ != nullptr && (!constantBC || !enableFluxVarsCache)) - delete boundaryVolVars_; - } - void update(const Problem& problem, const Element& element, const SubControlVolumeFace &scvFace) @@ -103,9 +96,9 @@ public: if (scvFace().boundary()) { - if (boundaryVolVars_ == nullptr) + if (!boundaryVolVars_) DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = boundaryVolVars_; + outsideVolVars = boundaryVolVars_.get(); } else outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); @@ -133,8 +126,26 @@ public: return stencil_; } + // when caching is enabled, get the stencil from the cache class + template + const typename std::enable_if::type& stencil() const + { + if (problemPtr_ == nullptr || scvFacePtr_ == nullptr) + DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + return problem().model().fluxVarsCache(scvFace()).stencil(); + } + + // when caching is disabled, return the private stencil variable. The update(...) routine has to be called beforehand. + template + const typename std::enable_if::type& stencil() + { + if (stencil_.empty()) + DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + return stencil_; + } + // returns the boundary vol vars belonging to this flux variables object - VolumeVariables getBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + VolumeVariables computeBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) { VolumeVariables boundaryVolVars; @@ -223,90 +234,55 @@ public: const Element& element, const SubControlVolumeFace &scvFace) { - problemPtr_ = &problem; - scvFacePtr_ = &scvFace; + ParentType::update(problem, element, scvFace); - if (scvFace.boundary()) - { - if(!boundaryVolVars_) - boundaryVolVars_ = Dune::Std::make_unique(); - } - - boundaryVolVarsUpdated_ = false; - - // TODO compute stencil as unity of advective and diffusive flux stencils - stencil_ = AdvectionType::stencil(scvFace); - } + // Stencil calculated in the base class is only the stencil for the advective flux. + // get stencil of the diffusive flux calculation + Stencil diffStencil = diffusionStencil(problem, scvFace); + // this has to be united with the diffusion stencil now + Stencil& stencil = this->stencil_; + stencil.insert(stencil.end(), diffStencil.begin(), diffStencil.end()); + std::sort(stencil.begin(), stencil.end()); + stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); - void beginFluxComputation() - { - // TODO if constant BC, do not set to false - boundaryVolVarsUpdated_ = false; - - // calculate advective fluxes and store them - advectiveVolFluxes_ = new std::array(); + // precompute advective volume fluxes for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - { - (*advectiveVolFluxes_)[phaseIdx] = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); - boundaryVolVarsUpdated_ = true; - } + advectiveVolumeFluxes_[phaseIdx] = AdvectionType::flux(problem, scvFace, phaseIdx, this->boundaryVolVars_); } - void endFluxComputation() - { - delete advectiveVolFluxes_; - } + Stencil diffusionStencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { return MolecularDiffusionType::stencil(problem, scvFace); } template Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - const auto* insideVolVars = &problem().model().curVolVars(scvFace().insideScvIdx()); + const auto* insideVolVars = &this->problem().model().curVolVars(this->scvFace().insideScvIdx()); const VolumeVariables* outsideVolVars; - if (scvFace().boundary()) - outsideVolVars = boundaryVolVars_.get(); + if (this->scvFace().boundary()) + { + if (!this->boundaryVolVars_) + DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); + outsideVolVars = this->boundaryVolVars_.get(); + } else - outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); + outsideVolVars = &this->problem().model().curVolVars(this->scvFace().outsideScvIdx()); - if (std::signbit((*advectiveVolFluxes_)[phaseIdx])) - return (*advectiveVolFluxes_)[phaseIdx]*upwindFunction(*outsideVolVars, *insideVolVars); + if (std::signbit(advectiveVolumeFluxes_[phaseIdx])) + return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(*outsideVolVars, *insideVolVars); else - return (*advectiveVolFluxes_)[phaseIdx]*upwindFunction(*insideVolVars, *outsideVolVars); + return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(*insideVolVars, *outsideVolVars); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - Scalar flux = MolecularDiffusionType::flux(problem(), scvFace(), phaseIdx, compIdx, boundaryVolVars_.get(), boundaryVolVarsUpdated_); - boundaryVolVarsUpdated_ = true; + Scalar flux = MolecularDiffusionType::flux(this->problem(), this->scvFace(), phaseIdx, compIdx, this->boundaryVolVars_); return flux; } - const Stencil& stencil() const - { - return stencil_; - } - - const Problem& problem() const - { - return *problemPtr_; - } - - const SubControlVolumeFace& scvFace() const - { - return *scvFacePtr_; - } - private: - const Problem *problemPtr_; //! Pointer to the problem - const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created - // storage for calculated advective fluxes to not having to calculate them again - std::array* advectiveVolFluxes_; - // boundary volume variables in case of Dirichlet boundaries - std::unique_ptr boundaryVolVars_; - bool boundaryVolVarsUpdated_; - // the flux stencil - Stencil stencil_; + std::array advectiveVolumeFluxes_; }; -- GitLab From b55f1033cd02199850a8e103aa791b3b7e848a0d Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 17:42:17 +0200 Subject: [PATCH 25/46] [implicit][fluxvarcache] make boundary volvars a shared ptr --- dumux/implicit/fluxvariablescache.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh index efb0514d85..e2bcc34ceb 100644 --- a/dumux/implicit/fluxvariablescache.hh +++ b/dumux/implicit/fluxvariablescache.hh @@ -103,18 +103,18 @@ public: FluxVariables fluxVars; stencil_ = fluxVars.computeFluxStencil(problem, scvFace); if (GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions)) - boundaryVolVars_ = std::make_unique(std::move(fluxVars.getBoundaryVolumeVariables(problem, element, scvFace))); + boundaryVolVars_ = std::make_shared(std::move(fluxVars.computeBoundaryVolumeVariables(problem, element, scvFace))); } const Stencil& stencil() const { return stencil_; } - const VolumeVariables& boundaryVolumeVariables() const - { return *boundaryVolVars_.get(); } + std::shared_ptr boundaryVolumeVariables() const + { return boundaryVolVars_; } private: Stencil stencil_; - std::unique_ptr boundaryVolVars_; + std::shared_ptr boundaryVolVars_; }; } // end namespace -- GitLab From ae7aedbec69cb4a63fc548b4f84280ed70184760 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 17:44:46 +0200 Subject: [PATCH 26/46] [3p3c][localresidual] adapt to new flux var structure --- dumux/porousmediumflow/3p3c/implicit/localresidual.hh | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh index a7c8ebfa7c..acfc1c28df 100644 --- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh +++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh @@ -47,6 +47,7 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { @@ -117,8 +118,8 @@ public: */ PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace) { - auto& fluxVars = this->model_().fluxVars_(scvFace); - fluxVars.beginFluxComputation(); + FluxVariables fluxVars; + fluxVars.update(asImp_()->problem_(), this->element_(), scvFace); // get upwind weights into local scope auto massWeight = massWeight_; @@ -167,8 +168,6 @@ public: flux[conti1EqIdx] += jNW+jNG+jNN; flux[conti2EqIdx] += jGW+jGG+jGN; - fluxVars.endFluxComputation(); - return flux; } -- GitLab From 34566cd33b4732c55f29f18aa025f4f3c8173bc9 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Wed, 18 May 2016 17:48:01 +0200 Subject: [PATCH 27/46] [tpfa][darcyslaw/fickslaw] accept shared pointer to vol vars --- .../implicit/cellcentered/tpfa/darcyslaw.hh | 6 +++--- .../implicit/cellcentered/tpfa/fickslaw.hh | 18 +++++++----------- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh index 11c2fb878d..d1ecb461df 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh @@ -78,7 +78,7 @@ public: static Scalar flux(const Problem& problem, const SubControlVolumeFace& scvFace, const IndexType phaseIdx, - const VolumeVariables* boundaryVolVars) + std::shared_ptr boundaryVolVars) { const auto& tij = getTransmissibilities(problem, scvFace); @@ -93,9 +93,9 @@ public: outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); else { - if (boundaryVolVars == nullptr) + if (!boundaryVolVars) DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = boundaryVolVars; + outsideVolVars = boundaryVolVars.get(); } auto hInside = insideVolVars->pressure(phaseIdx); diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh index f1740637ef..0caf2ad56d 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh @@ -60,6 +60,8 @@ class CCTpfaFicksLaw typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GridView::IndexSet::IndexType IndexType; + typedef typename std::vector Stencil; + using Element = typename GridView::template Codim<0>::Entity; enum { dim = GridView::dimension} ; @@ -75,8 +77,7 @@ public: const SubControlVolumeFace& scvFace, const int phaseIdx_, const int compIdx_, - VolumeVariables* boundaryVolVars, - bool boundaryVolVarsUpdated) + std::shared_ptr boundaryVolVars) { // diffusion tensors are always solution dependent Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx_, compIdx_); @@ -92,14 +93,9 @@ public: outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); else { - if (!boundaryVolVarsUpdated) - { - // update the boudary volvars for Dirichlet boundaries - const auto element = problem.model().fvGeometries().element(insideScv); - const auto dirichletPriVars = problem.dirichlet(element, scvFace); - boundaryVolVars->update(dirichletPriVars, problem, element, insideScv); - } - outsideVolVars = boundaryVolVars; + if (!boundaryVolVars) + DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); + outsideVolVars = boundaryVolVars.get(); } // compute the diffusive flux @@ -110,7 +106,7 @@ public: return rho*tij*(xInside - xOutside); } - static std::vector stencil(const SubControlVolumeFace& scvFace) + static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace) { std::vector stencil; stencil.clear(); -- GitLab From 807834416036d845dd72f7a7d6afc3f9f9eae49c Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 15:52:04 +0200 Subject: [PATCH 28/46] [implicit] fluxvars have to be updated before using stencil method --- dumux/implicit/cellcentered/localjacobian.hh | 4 +++- dumux/implicit/cellcentered/stencils.hh | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 73903e36a2..0200b16258 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -127,10 +127,12 @@ public: for (auto&& scvFaceJ : this->model_().fvGeometries(globalJ).scvfs()) { auto fluxVarsIdx = scvFaceJ.index(); + const auto& elementJ = this->model_().fvGeometries().element(globalJ); // if globalI is in flux var stencil, add to list FluxVariables fluxVars; - const auto& fluxStencil = fluxVars.stencil(problem, scvFaceJ); + fluxVars.update(problem, elementJ, scvFaceJ, false); + const auto& fluxStencil = fluxVars.stencil(); for (auto globalIdx : fluxStencil) if (globalIdx == globalI) fluxVarIndices.push_back(fluxVarsIdx); diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 84d3915586..e9ffd79756 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -49,7 +49,8 @@ public: for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) { FluxVariables fluxVars; - auto&& stencil = fluxVars.stencil(problem, scvf); + fluxVars.update(problem, element, scvf, false); + auto&& stencil = fluxVars.stencil(); elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); } // make values in elementstencil unique -- GitLab From 69155bf8c4fd366496dd4a2bdf5ca1c5eb79383f Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 15:54:21 +0200 Subject: [PATCH 29/46] [implicit][cc] introduce boundary volume variables --- dumux/implicit/cellcentered/localresidual.hh | 14 +++++ .../tpfa/fvelementgeometryvector.hh | 12 ++++- dumux/implicit/volumevariablesvector.hh | 53 +++++++++++++++++++ .../implicit/cellcentered/tpfa/darcyslaw.hh | 25 +++------ .../implicit/cellcentered/tpfa/fickslaw.hh | 37 +++++-------- 5 files changed, 97 insertions(+), 44 deletions(-) diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index a543c8932c..dee2e65c52 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -51,6 +51,8 @@ class CCLocalResidual : public ImplicitLocalResidual numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) }; + typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; @@ -73,7 +75,19 @@ protected: return this->asImp_().computeFlux(scvf); } else + { + if (!constantBC) + { + // update corresponding boundary volume variables before flux calculation + const auto insideScvIdx = scvf.insideScvIdx(); + const auto& insideScv = this->problem_().model().fvGeometries().subControlVolume(insideScvIdx); + const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf); + + this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv); + } + return this->asImp_().evalBoundary_(scvf); + } } PrimaryVariables evalBoundary_(const SubControlVolumeFace &scvf) diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh index 83a77b0bc5..88332de8cd 100644 --- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh +++ b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh @@ -82,12 +82,18 @@ public: return scvs_.size(); } - //! The total number of sun control volume faces + //! The total number of sub control volume faces std::size_t numScvf() const { return scvfs_.size(); } + //! The total number of boundary sub control volume faces + std::size_t numBoundaryScvf() const + { + return numBoundaryScvf_; + } + // Get an element from a sub control volume contained in it Element element(const SubControlVolume& scv) const { return elementMap_.element(scv.elementIndex()); } @@ -106,6 +112,7 @@ public: // Build the SCV and SCV faces IndexType scvfIdx = 0; + numBoundaryScvf_ = 0; elementMap_.resize(gridView_.size(0)); scvs_.resize(gridView_.size(0)); for (const auto& element : elements(gridView_)) @@ -137,7 +144,7 @@ public: intersection.geometry().center(), intersection.centerUnitOuterNormal(), scvfIdx, - std::vector({eIdx}), + std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}), true)); scvfsIndexSet.push_back(scvfIdx++); } @@ -154,6 +161,7 @@ private: std::vector> scvs_; std::vector> scvfs_; std::vector fvGeometries_; + IndexType numBoundaryScvf_; }; } // end namespace diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh index 6a3db2f03e..fb9555d2f9 100644 --- a/dumux/implicit/volumevariablesvector.hh +++ b/dumux/implicit/volumevariablesvector.hh @@ -36,13 +36,26 @@ template class VolumeVariablesVector : public std::vector { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using IndexType = typename GridView::IndexSet::IndexType; + + enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; public: void update(const Problem& problem, const SolutionVector& sol) { + numScvs_ = problem.model().fvGeometries().numScv(); + if (!isBox) + numBoundaryScvs_ = problem.model().fvGeometries().numBoundaryScvf(); + else + numBoundaryScvs_ = 0; + + volumeVariables_.resize(numScvs_); + boundaryVolumeVariables_.resize(numBoundaryScvs_); for (const auto& element : elements(problem.gridView())) { for (auto&& scv : problem.model().fvGeometries(element).scvs()) @@ -53,8 +66,48 @@ public: scv); } + if (!isBox) + { + for (auto&& scvFace : problem.model().fvGeometries(element).scvfs()) + { + if (!scvFace.boundary()) + continue; + + const auto insideScvIdx = scvFace.insideScvIdx(); + const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); + const auto dirichletPriVars = problem.dirichlet(element, scvFace); + + (*this)[scvFace.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv); + } + } } } + + const VolumeVariables& operator [](IndexType scvIdx) const + { + assert(scvIdx < numScvs_ + numBoundaryScvs_); + + if (scvIdx < numScvs_) + return volumeVariables_[scvIdx]; + else + return boundaryVolumeVariables_[scvIdx-numScvs_]; + } + + VolumeVariables& operator [](IndexType scvIdx) + { + assert(scvIdx < numScvs_ + numBoundaryScvs_); + + if (scvIdx < numScvs_) + return volumeVariables_[scvIdx]; + else + return boundaryVolumeVariables_[scvIdx-numScvs_]; + } + +private: + IndexType numScvs_; + IndexType numBoundaryScvs_; + std::vector volumeVariables_; + std::vector boundaryVolumeVariables_; }; } // end namespace diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh index d1ecb461df..6ed9bef547 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh @@ -56,8 +56,6 @@ class CCTpfaDarcysLaw typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; - typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; @@ -77,35 +75,26 @@ public: static Scalar flux(const Problem& problem, const SubControlVolumeFace& scvFace, - const IndexType phaseIdx, - std::shared_ptr boundaryVolVars) + const IndexType phaseIdx) { const auto& tij = getTransmissibilities(problem, scvFace); // Get the inside volume variables const auto insideScvIdx = scvFace.insideScvIdx(); const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); - const auto* insideVolVars = &problem.model().curVolVars(insideScv); + const auto& insideVolVars = problem.model().curVolVars(insideScv); // and the outside volume variables - const VolumeVariables* outsideVolVars; - if (!scvFace.boundary()) - outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); - else - { - if (!boundaryVolVars) - DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = boundaryVolVars.get(); - } + const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx()); - auto hInside = insideVolVars->pressure(phaseIdx); - auto hOutside = outsideVolVars->pressure(phaseIdx); + auto hInside = insideVolVars.pressure(phaseIdx); + auto hOutside = outsideVolVars.pressure(phaseIdx); if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity)) { // do averaging for the density - const auto rhoInside = insideVolVars->density(phaseIdx); - const auto rhoOutide = outsideVolVars->density(phaseIdx); + const auto rhoInside = insideVolVars.density(phaseIdx); + const auto rhoOutide = outsideVolVars.density(phaseIdx); const auto rho = (rhoInside + rhoOutide)*0.5; diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh index 0caf2ad56d..309e438e81 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh @@ -56,8 +56,6 @@ class CCTpfaFicksLaw typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; - typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GridView::IndexSet::IndexType IndexType; typedef typename std::vector Stencil; @@ -75,33 +73,24 @@ public: static Scalar flux(const Problem& problem, const SubControlVolumeFace& scvFace, - const int phaseIdx_, - const int compIdx_, - std::shared_ptr boundaryVolVars) + const int phaseIdx, + const int compIdx) { // diffusion tensors are always solution dependent - Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx_, compIdx_); + Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx); // Get the inside volume variables const auto insideScvIdx = scvFace.insideScvIdx(); const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); - const auto* insideVolVars = &problem.model().curVolVars(insideScv); + const auto& insideVolVars = problem.model().curVolVars(insideScv); // and the outside volume variables - const VolumeVariables* outsideVolVars; - if (!scvFace.boundary()) - outsideVolVars = &problem.model().curVolVars(scvFace.outsideScvIdx()); - else - { - if (!boundaryVolVars) - DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = boundaryVolVars.get(); - } + const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx()); // compute the diffusive flux - const auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_); - const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_); - const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_)); + const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx); + const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx); + const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx)); return rho*tij*(xInside - xOutside); } @@ -124,7 +113,7 @@ public: private: - static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx_, const int compIdx_) + static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx) { Scalar tij; @@ -132,8 +121,8 @@ private: const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); const auto& insideVolVars = problem.model().curVolVars(insideScvIdx); - auto insideD = insideVolVars.diffusionCoefficient(phaseIdx_, compIdx_); - insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx_), insideD); + auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx); + insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD); Scalar ti = calculateOmega_(problem, scvFace, insideD, insideScv); if (!scvFace.boundary()) @@ -142,8 +131,8 @@ private: const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx); const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx); - auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx_, compIdx_); - outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx_), outsideD); + auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx); + outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD); Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideD, outsideScv); tij = scvFace.area()*(ti * tj)/(ti + tj); -- GitLab From f5c4b90880521cf7e0ca6b7d4b9430f029dd7f7b Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 15:55:18 +0200 Subject: [PATCH 30/46] [implicit][fluxvars] reimplement flux variables --- dumux/implicit/fluxvariables.hh | 238 +++++++++++++++----------------- 1 file changed, 108 insertions(+), 130 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index a35c9483f7..a242b0c679 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -37,15 +37,10 @@ NEW_PROP_TAG(NumComponents); /*! * \ingroup ImplicitModel * \brief Base class for the flux variables - * specializations are provided for combinations of physical processes + * Actual flux variables inherit from this class */ -template -class FluxVariables {}; - - -// specialization for pure advective flow (e.g. 1p/2p/3p immiscible darcy flow) -template -class FluxVariables +template +class FluxVariablesBase { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); @@ -53,19 +48,12 @@ class FluxVariables using IndexType = typename GridView::IndexSet::IndexType; using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Stencil = std::vector; - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); - enum - { - isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox), - enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache), - constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) - }; + enum{ enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache) }; public: - FluxVariables() : problemPtr_(nullptr), scvFacePtr_(nullptr) + FluxVariablesBase() : problemPtr_(nullptr), scvFacePtr_(nullptr) {} void update(const Problem& problem, @@ -75,55 +63,9 @@ public: problemPtr_ = &problem; scvFacePtr_ = &scvFace; - // boundary vol vars only need to be handled at boundaries for cc methods - if (scvFace.boundary() && !isBox) - setBoundaryVolumeVariables_(problem, element, scvFace); - // update the stencil if needed - if (!enableFluxVarsCache && stencil_.empty()) - stencil_ = stencil(problem, scvFace); - } - - - - template - Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) - { - Scalar flux = AdvectionType::flux(problem(), scvFace(), phaseIdx, boundaryVolVars_); - - const auto* insideVolVars = &problem().model().curVolVars(scvFace().insideScvIdx()); - const VolumeVariables* outsideVolVars; - - if (scvFace().boundary()) - { - if (!boundaryVolVars_) - DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = boundaryVolVars_.get(); - } - else - outsideVolVars = &problem().model().curVolVars(scvFace().outsideScvIdx()); - - if (std::signbit(flux)) - return flux*upwindFunction(*outsideVolVars, *insideVolVars); - else - return flux*upwindFunction(*insideVolVars, *outsideVolVars); - - return flux; - } - - // interface allowing for stencil information without having to update the flux vars. - // this becomes useful when the element is not at hand for a call to update(...), e.g. during the assembly - localjacobian.hh - template - const typename std::enable_if::type& stencil(const Problem& problem, const SubControlVolumeFace& scvFace) const - { return problem.model().fluxVarsCache(scvFace).stencil(); } - - // provide interface in case caching is disabled - template - const typename std::enable_if::type& stencil(const Problem& problem, const SubControlVolumeFace& scvFace) - { - if (stencil_.empty()) - stencil_ = computeFluxStencil(problem, scvFace); - return stencil_; + if (!enableFluxVarsCache) + stencil_ = asImp_().computeStencil(problem, scvFace); } // when caching is enabled, get the stencil from the cache class @@ -131,7 +73,7 @@ public: const typename std::enable_if::type& stencil() const { if (problemPtr_ == nullptr || scvFacePtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); return problem().model().fluxVarsCache(scvFace()).stencil(); } @@ -140,74 +82,114 @@ public: const typename std::enable_if::type& stencil() { if (stencil_.empty()) - DUNE_THROW(Dune::InvalidStateException, "Calling the stencil() method before update of the FluxVariables object."); + DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); return stencil_; } - // returns the boundary vol vars belonging to this flux variables object - VolumeVariables computeBoundaryVolumeVariables(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) - { - VolumeVariables boundaryVolVars; - - const auto insideScvIdx = scvFace.insideScvIdx(); - const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); - const auto dirichletPriVars = problem.dirichlet(element, scvFace); - boundaryVolVars.update(dirichletPriVars, problem, element, insideScv); - - return boundaryVolVars; - } - const Problem& problem() const { if (problemPtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "Problem pointer not valid. Call update before using the flux variables class."); + DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); return *problemPtr_; } const SubControlVolumeFace& scvFace() const { if (scvFacePtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "Scv face pointer not valid. Call update before using the flux variables class."); + DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); return *scvFacePtr_; } - Stencil computeFluxStencil(const Problem& problem, const SubControlVolumeFace& scvFace) - { return AdvectionType::stencil(problem, scvFace); } + Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); } -protected: +private: - // if flux variables caching is enabled, we use the boundary volume variables stored in the cache class - template - typename std::enable_if::type - setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + Implementation &asImp_() { - boundaryVolVars_ = std::shared_ptr(problem.model().fluxVarsCache(scvFace).boundaryVolumeVariables()); + assert(static_cast(this) != 0); + return *static_cast(this); } - // if flux variables caching is disabled, we update the boundary volume variables - template - typename std::enable_if::type - setBoundaryVolumeVariables_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace) + const Implementation &asImp_() const { - boundaryVolVars_ = std::make_shared(std::move(computeBoundaryVolumeVariables(problem, element, scvFace))); + assert(static_cast(this) != 0); + return *static_cast(this); } - std::shared_ptr boundaryVolVars_; //! boundary volume variables in case of Dirichlet boundaries - // the flux stencil - Stencil stencil_; - -private: - const Problem *problemPtr_; //! Pointer to the problem const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created + Stencil stencil_; //! The flux stencil +}; + + + +/*! + * \ingroup ImplicitModel + * \brief the flux variables class + * specializations are provided for combinations of physical processes + */ +template +class FluxVariables {}; + + +// specialization for pure advective flow (e.g. 1p/2p/3p immiscible darcy flow) +template +class FluxVariables : public FluxVariablesBase> +{ + using ParentType = FluxVariablesBase>; + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Stencil = std::vector; + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); + + enum + { + isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox), + enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache), + constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) + }; + +public: + FluxVariables() + {} + + void update(const Problem& problem, + const Element& element, + const SubControlVolumeFace &scvFace, + const bool preComputation = true) + { + ParentType::update(problem, element, scvFace); + } + + template + Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) + { + Scalar flux = AdvectionType::flux(this->problem(), this->scvFace(), phaseIdx); + + const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx()); + const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx()); + + if (std::signbit(flux)) + return flux*upwindFunction(outsideVolVars, insideVolVars); + else + return flux*upwindFunction(insideVolVars, outsideVolVars); + } + + Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { return AdvectionType::stencil(problem, scvFace); } }; // specialization for isothermal advection molecularDiffusion equations template -class FluxVariables : public FluxVariables +class FluxVariables : public FluxVariablesBase> { - using ParentType = FluxVariables; + using ParentType = FluxVariablesBase>; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; @@ -232,51 +214,47 @@ public: void update(const Problem& problem, const Element& element, - const SubControlVolumeFace &scvFace) + const SubControlVolumeFace &scvFace, + const bool preComputation = true) { ParentType::update(problem, element, scvFace); - // Stencil calculated in the base class is only the stencil for the advective flux. - // get stencil of the diffusive flux calculation - Stencil diffStencil = diffusionStencil(problem, scvFace); - // this has to be united with the diffusion stencil now - Stencil& stencil = this->stencil_; - stencil.insert(stencil.end(), diffStencil.begin(), diffStencil.end()); + // precompute advective volume fluxes + if (preComputation) + { + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + advectiveVolumeFluxes_[phaseIdx] = AdvectionType::flux(problem, scvFace, phaseIdx); + } + } + + Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace) + { + // unifiy advective and diffusive stencil + Stencil stencil = AdvectionType::stencil(problem, scvFace); + Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, scvFace); + + stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end()); std::sort(stencil.begin(), stencil.end()); stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); - // precompute advective volume fluxes - for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - advectiveVolumeFluxes_[phaseIdx] = AdvectionType::flux(problem, scvFace, phaseIdx, this->boundaryVolVars_); + return stencil; } - Stencil diffusionStencil(const Problem& problem, const SubControlVolumeFace& scvFace) - { return MolecularDiffusionType::stencil(problem, scvFace); } - template Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction) { - const auto* insideVolVars = &this->problem().model().curVolVars(this->scvFace().insideScvIdx()); - const VolumeVariables* outsideVolVars; - - if (this->scvFace().boundary()) - { - if (!this->boundaryVolVars_) - DUNE_THROW(Dune::InvalidStateException, "Trying to access invalid boundary volume variables."); - outsideVolVars = this->boundaryVolVars_.get(); - } - else - outsideVolVars = &this->problem().model().curVolVars(this->scvFace().outsideScvIdx()); + const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx()); + const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx()); if (std::signbit(advectiveVolumeFluxes_[phaseIdx])) - return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(*outsideVolVars, *insideVolVars); + return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); else - return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(*insideVolVars, *outsideVolVars); + return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) { - Scalar flux = MolecularDiffusionType::flux(this->problem(), this->scvFace(), phaseIdx, compIdx, this->boundaryVolVars_); + Scalar flux = MolecularDiffusionType::flux(this->problem(), this->scvFace(), phaseIdx, compIdx); return flux; } -- GitLab From 095b94c9fa1968570a9ecb4e0ae7ec05066850b6 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 15:56:01 +0200 Subject: [PATCH 31/46] [implicit][fluxvarscache] remove cached boundary vol vars --- dumux/implicit/fluxvariablescache.hh | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh index e2bcc34ceb..5f2eb9c6cc 100644 --- a/dumux/implicit/fluxvariablescache.hh +++ b/dumux/implicit/fluxvariablescache.hh @@ -101,20 +101,14 @@ public: { ParentType::update(problem, element, scvFace); FluxVariables fluxVars; - stencil_ = fluxVars.computeFluxStencil(problem, scvFace); - if (GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions)) - boundaryVolVars_ = std::make_shared(std::move(fluxVars.computeBoundaryVolumeVariables(problem, element, scvFace))); + stencil_ = fluxVars.computeStencil(problem, scvFace); } const Stencil& stencil() const { return stencil_; } - std::shared_ptr boundaryVolumeVariables() const - { return boundaryVolVars_; } - private: Stencil stencil_; - std::shared_ptr boundaryVolVars_; }; } // end namespace -- GitLab From a20f0f92c0b7ec84ef9ceacda2ae6765d1de7ff6 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 15:56:48 +0200 Subject: [PATCH 32/46] [implicit][model] make localresiduals friend classes The local residuals need to be able to update the boundary volume variables, thus they need to be friend classes in order to have access to the protected non constant return function for the volvarsvector. --- dumux/implicit/model.hh | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 5b209f2b41..22cb8403b7 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -34,6 +34,10 @@ namespace Dumux { +// forward declaration of the friend localresidual classes +template class CCLocalResidual; +template class BoxLocalResidual; + /*! * \ingroup ImplicitModel * \brief The base class for the vertex centered finite volume @@ -43,7 +47,8 @@ template class ImplicitModel { friend typename GET_PROP_TYPE(TypeTag, LocalJacobian); - friend typename GET_PROP_TYPE(TypeTag, LocalResidual); + friend CCLocalResidual; + friend BoxLocalResidual; typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; @@ -116,7 +121,6 @@ public: asImp_().applyInitialSolution_(); // resize and update the volVars with the initial solution - curVolVarsVector_.resize(fvGeometries().numScv()); curVolVarsVector_.update(problem_(), curSol()); // update the flux variables caches -- GitLab From da71f73f63834ca3b43d253cb093a2ece80af867 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 16:34:00 +0200 Subject: [PATCH 33/46] [implicit] only update boundary volvars on pure Dirichlet boundaries When complex boundary handling is inactive, we only use the boundary volume variables on pure Dirichlet boundaries. Thus, updating becomes obsolete. --- dumux/implicit/cellcentered/localresidual.hh | 23 ++++++++------------ dumux/implicit/volumevariablesvector.hh | 6 +++++ 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index dee2e65c52..829f6c135b 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -71,23 +71,9 @@ protected: PrimaryVariables computeFlux_(const SubControlVolumeFace &scvf) { if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/) - { return this->asImp_().computeFlux(scvf); - } else - { - if (!constantBC) - { - // update corresponding boundary volume variables before flux calculation - const auto insideScvIdx = scvf.insideScvIdx(); - const auto& insideScv = this->problem_().model().fvGeometries().subControlVolume(insideScvIdx); - const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf); - - this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv); - } - return this->asImp_().evalBoundary_(scvf); - } } PrimaryVariables evalBoundary_(const SubControlVolumeFace &scvf) @@ -229,6 +215,15 @@ protected: // temporary vector to store the Dirichlet boundary fluxes PrimaryVariables flux(0); + if (!constantBC) + { + // update corresponding boundary volume variables before flux calculation + const auto insideScvIdx = scvf.insideScvIdx(); + const auto& insideScv = this->problem_().model().fvGeometries().subControlVolume(insideScvIdx); + const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf); + this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv); + } + auto dirichletFlux = this->asImp_().computeFlux(scvf); // add fluxes to the residual diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh index fb9555d2f9..ef0bea9ee2 100644 --- a/dumux/implicit/volumevariablesvector.hh +++ b/dumux/implicit/volumevariablesvector.hh @@ -70,9 +70,15 @@ public: { for (auto&& scvFace : problem.model().fvGeometries(element).scvfs()) { + // if we are not on a boundary, skip the rest if (!scvFace.boundary()) continue; + // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries + auto bcTypes = problem.boundaryTypes(element, scvFace); + if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann())) + continue; + const auto insideScvIdx = scvFace.insideScvIdx(); const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx); const auto dirichletPriVars = problem.dirichlet(element, scvFace); -- GitLab From 58f80a6d1bf36baa182c5ad74554b8129bda3b34 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 19:10:25 +0200 Subject: [PATCH 34/46] [implicit][fluxVars] rename update function The flux variables base class now has the init() routine which takes care of setting the pointers and eventually the stencil. The actual implementation of the flux vars calls this method in its method initAndComputeFluxes(). This eventually precomputes the advective fluxes. Wherever only the stencil is desired (localjacobian, stencilvector), only the init() method from the base class is called in order to save computational time. --- dumux/implicit/cellcentered/localjacobian.hh | 2 +- dumux/implicit/cellcentered/stencils.hh | 2 +- dumux/implicit/fluxvariables.hh | 37 +++++++------------ .../3p3c/implicit/localresidual.hh | 2 +- .../immiscible/localresidual.hh | 2 +- 5 files changed, 17 insertions(+), 28 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 0200b16258..40b2494272 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -131,7 +131,7 @@ public: // if globalI is in flux var stencil, add to list FluxVariables fluxVars; - fluxVars.update(problem, elementJ, scvFaceJ, false); + fluxVars.init(problem, elementJ, scvFaceJ); const auto& fluxStencil = fluxVars.stencil(); for (auto globalIdx : fluxStencil) if (globalIdx == globalI) diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index e9ffd79756..24a2deab7d 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -49,7 +49,7 @@ public: for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) { FluxVariables fluxVars; - fluxVars.update(problem, element, scvf, false); + fluxVars.init(problem, element, scvf); auto&& stencil = fluxVars.stencil(); elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); } diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index a242b0c679..5e074a82a9 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -56,9 +56,9 @@ public: FluxVariablesBase() : problemPtr_(nullptr), scvFacePtr_(nullptr) {} - void update(const Problem& problem, - const Element& element, - const SubControlVolumeFace &scvFace) + void init(const Problem& problem, + const Element& element, + const SubControlVolumeFace &scvFace) { problemPtr_ = &problem; scvFacePtr_ = &scvFace; @@ -155,15 +155,12 @@ class FluxVariables : public FluxVariablesBase @@ -209,22 +206,14 @@ class FluxVariables : public FluxVariablesBaseproblem_(), this->element_(), scvFace); + fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvFace); // get upwind weights into local scope auto massWeight = massWeight_; diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh index efe919230b..ee999d5aae 100644 --- a/dumux/porousmediumflow/immiscible/localresidual.hh +++ b/dumux/porousmediumflow/immiscible/localresidual.hh @@ -98,7 +98,7 @@ public: PrimaryVariables computeFlux(const SubControlVolumeFace& scvf) { FluxVariables fluxVars; - fluxVars.update(asImp_()->problem_(), this->element_(), scvf); + fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvf); // copy weight to local scope for use in lambda expression auto w = upwindWeight_; -- GitLab From e5b34429db47683b15ad1e9b66a5f754331bf7a4 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 19:20:39 +0200 Subject: [PATCH 35/46] [implicit][fluxvarscache] method-specific flux cache classes The flux var cache classes are specific to the different methods now in order to customize the stored variables and types. For tpfa models the transmissivity can thus be stored as a scalar instead of a vector. This reduces the memory usage and increases the computational speed. --- .../cellcentered/tpfa/propertydefaults.hh | 4 ++ dumux/implicit/fluxvariablescache.hh | 38 ++++++++++--------- dumux/implicit/propertydefaults.hh | 9 ----- .../implicit/cellcentered/tpfa/darcyslaw.hh | 12 +++--- 4 files changed, 29 insertions(+), 34 deletions(-) diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index ed73e4960b..3259a4bf45 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -29,6 +29,7 @@ #include #include + #include #include #include #include @@ -64,6 +65,9 @@ public: typedef Dumux::SubControlVolumeFace type; }; +//! The flux variables cache class +SET_TYPE_PROP(CCTpfaModel, FluxVariablesCache, Dumux::CCTpfaFluxVariablesCache); + //! The darcy flux variables SET_TYPE_PROP(CCTpfaModel, AdvectionType, Dumux::CCTpfaDarcysLaw); diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh index 5f2eb9c6cc..153b249912 100644 --- a/dumux/implicit/fluxvariablescache.hh +++ b/dumux/implicit/fluxvariablescache.hh @@ -36,16 +36,12 @@ NEW_PROP_TAG(NumComponents); /*! * \ingroup ImplicitModel - * \brief The flux variables cache class - * stores the transmissibilities and stencils for the darcy flux calculation on a scv face + * \brief The flux variables cache classes + * stores the transmissibilities and stencils */ -template -class FluxVariablesCache -{}; - // specialization for the Box Method template -class FluxVariablesCache +class BoxFluxVariablesCache { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); @@ -55,41 +51,43 @@ class FluxVariablesCache using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; using Stencil = std::vector; - using TransmissibilityVector = std::vector; + using TransmissibilityVector = std::vector; public: void update(const Problem& problem, const Element& element, const SubControlVolumeFace &scvFace) { - darcyStencil_ = AdvectionType::stencil(problem, scvFace); + stencil_ = AdvectionType::stencil(problem, scvFace); + volVarsStencil_ = AdvectionType::volVarsStencil(problem, element, scvFace); tij_ = AdvectionType::calculateTransmissibilities(problem, scvFace); } const Stencil& stencil() const - { return darcyStencil_; } + { return stencil_; } const Stencil& darcyStencil() const - { return darcyStencil_; } + { return volVarsStencil_; } const TransmissibilityVector& tij() const { return tij_; } private: - Stencil darcyStencil_; + Stencil volVarsStencil_; + Stencil stencil_; TransmissibilityVector tij_; }; -// specialization for cell centered methods +// specialization for the cell centered tpfa method template -class FluxVariablesCache : public FluxVariablesCache +class CCTpfaFluxVariablesCache { - using ParentType = FluxVariablesCache; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType); using Element = typename GridView::template Codim<0>::Entity; using IndexType = typename GridView::IndexSet::IndexType; using Stencil = std::vector; @@ -99,16 +97,20 @@ public: const Element& element, const SubControlVolumeFace &scvFace) { - ParentType::update(problem, element, scvFace); FluxVariables fluxVars; stencil_ = fluxVars.computeStencil(problem, scvFace); + tij_ = AdvectionType::calculateTransmissibilities(problem, scvFace); } const Stencil& stencil() const { return stencil_; } + const Scalar& tij() const + { return tij_; } + private: Stencil stencil_; + Scalar tij_; }; } // end namespace diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 8d9af7d5b1..6a500fc2b1 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -99,15 +99,6 @@ SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables); -//! The flux variables cache class -SET_PROP(ImplicitBase, FluxVariablesCache) -{ -private: - enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - public: - typedef typename Dumux::FluxVariablesCache type; -}; - //! The global flux variables cache vector class SET_TYPE_PROP(ImplicitBase, FluxVariablesCacheVector, Dumux::FluxVariablesCacheVector); diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh index 6ed9bef547..49a35b14d0 100644 --- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh +++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh @@ -60,7 +60,6 @@ class CCTpfaDarcysLaw typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GridView::IndexSet::IndexType IndexType; - typedef std::vector TransmissivityVector; typedef std::vector Stencil; using Element = typename GridView::template Codim<0>::Entity; @@ -121,13 +120,12 @@ public: } } - return tij[0]*(hInside - hOutside); + return tij*(hInside - hOutside); } static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace) { std::vector stencil; - stencil.clear(); if (!scvFace.boundary()) { stencil.push_back(scvFace.insideScvIdx()); @@ -140,18 +138,18 @@ public: } template - static const typename std::enable_if::type& getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) + static const typename std::enable_if::type& getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) { return problem.model().fluxVarsCache(scvFace).tij(); } template - static const typename std::enable_if::type getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) + static const typename std::enable_if::type getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) { return calculateTransmissibilities(problem, scvFace); } - static TransmissivityVector calculateTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) + static Scalar calculateTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace) { Scalar tij; @@ -174,7 +172,7 @@ public: tij = scvFace.area()*ti; } - return TransmissivityVector(1, tij); + return tij; } private: -- GitLab From 0b7ece2498f997196c37c615039bbf7175aa2a5c Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 19:23:43 +0200 Subject: [PATCH 36/46] [implicit][fluxvars] remove assertion of pointers In order to reduce the computational cost the pointers are no longer checked for their validity. The user has to make sure that the flux variables objects are properly initialized, otherwise a segmentation fault will follow. --- dumux/implicit/fluxvariables.hh | 20 +++----------------- 1 file changed, 3 insertions(+), 17 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index 5e074a82a9..ac940158b1 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -71,32 +71,18 @@ public: // when caching is enabled, get the stencil from the cache class template const typename std::enable_if::type& stencil() const - { - if (problemPtr_ == nullptr || scvFacePtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); - return problem().model().fluxVarsCache(scvFace()).stencil(); - } + { return problem().model().fluxVarsCache(scvFace()).stencil(); } // when caching is disabled, return the private stencil variable. The update(...) routine has to be called beforehand. template const typename std::enable_if::type& stencil() - { - if (stencil_.empty()) - DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); - return stencil_; - } + { return stencil_; } const Problem& problem() const - { - if (problemPtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); - return *problemPtr_; - } + { return *problemPtr_; } const SubControlVolumeFace& scvFace() const { - if (scvFacePtr_ == nullptr) - DUNE_THROW(Dune::InvalidStateException, "FluxVariables object has not been properly initialized. Call the update(...) method before using it."); return *scvFacePtr_; } -- GitLab From 3e460d3fa88686ad0a5751975174a8dddeadc09c Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Thu, 19 May 2016 19:27:03 +0200 Subject: [PATCH 37/46] [implicit][fluxvars] rename advective volume fluxes the advective volume fluxes are now simply called advective fluxes as they are no volume fluxes regarding their unit. --- dumux/implicit/fluxvariables.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh index ac940158b1..e7a4513161 100644 --- a/dumux/implicit/fluxvariables.hh +++ b/dumux/implicit/fluxvariables.hh @@ -221,10 +221,10 @@ public: const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx()); const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx()); - if (std::signbit(advectiveVolumeFluxes_[phaseIdx])) - return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); + if (std::signbit(advectiveFluxes_[phaseIdx])) + return advectiveFluxes_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars); else - return advectiveVolumeFluxes_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); + return advectiveFluxes_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars); } Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx) @@ -235,7 +235,7 @@ public: private: // storage for calculated advective fluxes to not having to calculate them again - std::array advectiveVolumeFluxes_; + std::array advectiveFluxes_; }; -- GitLab From b16c900ec39818c36180d2c9c98cbd186bc5759f Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 May 2016 14:47:00 +0200 Subject: [PATCH 38/46] [implicit] make vol vars cache optional The volume variables caching can now be disabled in order for the simulation to use less memory. The FVGeometry can be bound to an element which then creates the geometries necessary for operations on the element or the complete stencil. Several classes now need to be friend with the model, otherwise the FVGeometry can not be bound to the element. Maybe the bind function could return a geometry instead of a reference, then access rights can be limited better. It would decrease the efficiency though. --- dumux/implicit/cellcentered/localjacobian.hh | 3 + dumux/implicit/cellcentered/localresidual.hh | 4 +- dumux/implicit/fluxvariablescachevector.hh | 11 +- dumux/implicit/localresidual.hh | 12 + dumux/implicit/model.hh | 33 +- dumux/implicit/properties.hh | 4 +- dumux/implicit/propertydefaults.hh | 10 +- dumux/implicit/volumevariablesvector.hh | 284 +++++++++++++++++- dumux/porousmediumflow/1p/implicit/model.hh | 5 +- dumux/porousmediumflow/2p/implicit/model.hh | 5 +- dumux/porousmediumflow/3p/implicit/model.hh | 5 +- dumux/porousmediumflow/3p3c/implicit/model.hh | 10 +- .../compositional/primaryvariableswitch.hh | 12 +- 13 files changed, 365 insertions(+), 33 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 40b2494272..53230f0b91 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -151,6 +151,9 @@ public: */ void assemble(const Element& element, JacobianMatrix& matrix) { + this->model_().curVolVars_().bind(element); + this->model_().prevVolVars_().bindElement(element); + // set the current grid element and update the element's // finite volume geometry globalI_ = this->problem_().elementMapper().index(element); diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index 829f6c135b..d599045821 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -53,6 +53,8 @@ class CCLocalResidual : public ImplicitLocalResidual enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) }; + enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache) }; + typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; @@ -215,7 +217,7 @@ protected: // temporary vector to store the Dirichlet boundary fluxes PrimaryVariables flux(0); - if (!constantBC) + if (!constantBC || !enableVolVarsCache) { // update corresponding boundary volume variables before flux calculation const auto insideScvIdx = scvf.insideScvIdx(); diff --git a/dumux/implicit/fluxvariablescachevector.hh b/dumux/implicit/fluxvariablescachevector.hh index 159768f1a8..e426ac6dd2 100644 --- a/dumux/implicit/fluxvariablescachevector.hh +++ b/dumux/implicit/fluxvariablescachevector.hh @@ -43,12 +43,17 @@ class FluxVariablesCacheVector public: template - typename std::enable_if::type update(const Problem& problem) + typename std::enable_if::type update(Problem& problem) { fluxVarsCache_.resize(problem.model().fvGeometries().numScvf()); for (const auto& element : elements(problem.gridView())) { - for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) + // bind the geometries and volume variables to the element (all the elements in stencil) + problem.model().fvGeometries_().bind(element); + problem.model().curVolVars_().bind(element); + + const auto& fvGeometry = problem.model().fvGeometries(element); + for (const auto& scvf : fvGeometry.scvfs()) { (*this)[scvf.index()].update(problem, element, scvf); } @@ -56,7 +61,7 @@ public: } template - typename std::enable_if::type update(const Problem& problem) + typename std::enable_if::type update(Problem& problem) {} template diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index 1ab9f085ae..a5a4c79319 100644 --- a/dumux/implicit/localresidual.hh +++ b/dumux/implicit/localresidual.hh @@ -91,6 +91,10 @@ public: */ void eval(const Element &element) { + // make sure FVElementGeometry and volume variables are bound to the element + model_().curVolVars_().bind(element); + model_().prevVolVars_().bindElement(element); + ElementBoundaryTypes bcTypes; bcTypes.update(problem_(), element, fvGeometry_()); @@ -110,6 +114,10 @@ public: { elemPtr_ = &element; + // make sure FVElementGeometry and volume variables are bound to the element + model_().curVolVars_().bindElement(element); + model_().prevVolVars_().bindElement(element); + ElementBoundaryTypes bcTypes; bcTypes.update(problem_(), element, fvGeometry_()); bcTypesPtr_ = &bcTypes; @@ -129,6 +137,10 @@ public: { elemPtr_ = &element; + // make sure FVElementGeometry and volume variables are bound to the element + model_().curVolVars_().bind(element); + model_().prevVolVars_().bindElement(element); + ElementBoundaryTypes bcTypes; bcTypes.update(problem_(), element, fvGeometry_()); diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 22cb8403b7..8d09e4de0d 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -27,9 +27,11 @@ #include #include "properties.hh" +#include "localresidual.hh" #include #include #include + #include namespace Dumux { @@ -46,9 +48,17 @@ template class BoxLocalResidual; template class ImplicitModel { - friend typename GET_PROP_TYPE(TypeTag, LocalJacobian); friend CCLocalResidual; friend BoxLocalResidual; + friend ImplicitLocalResidual; + friend PrimaryVariableSwitch; + friend typename GET_PROP_TYPE(TypeTag, Problem); + friend typename GET_PROP_TYPE(TypeTag, LocalJacobian); + friend typename GET_PROP_TYPE(TypeTag, StencilsVector); + friend typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector); + friend typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector); + friend typename GET_PROP_TYPE(TypeTag, FluxVariablesCacheVector); + typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; @@ -58,7 +68,8 @@ class ImplicitModel typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler; - typedef typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector) VolumeVariablesVector; + typedef typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector) CurrentVolumeVariablesVector; + typedef typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector) PreviousVolumeVariablesVector; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; @@ -123,12 +134,12 @@ public: // resize and update the volVars with the initial solution curVolVarsVector_.update(problem_(), curSol()); - // update the flux variables caches - fluxVarsCacheVector_.update(problem_()); - // update stencils stencilsVector_.update(problem_()); + // update the flux variables caches + fluxVarsCacheVector_.update(problem_()); + // initialize assembler and create matrix localJacobian_.init(problem_()); jacAsm_ = std::make_shared(); @@ -801,17 +812,19 @@ public: { return fvGeometries_->fvGeometry(eIdx); } protected: - // only to be called by friends and deriving classes - VolumeVariablesVector& curVolVars_() + + CurrentVolumeVariablesVector& curVolVars_() { return curVolVarsVector_; } - // only to be called by friends and deriving classes VolumeVariables& curVolVars_(const SubControlVolume& scv) { return curVolVarsVector_[scv.index()]; } VolumeVariables& curVolVars_(unsigned int scvIdx) { return curVolVarsVector_[scvIdx]; } + PreviousVolumeVariablesVector& prevVolVars_() + { return prevVolVarsVector_; } + VolumeVariables& prevVolVars_(const SubControlVolume& scv) { return prevVolVarsVector_[scv.index()]; } @@ -972,8 +985,8 @@ protected: SolutionVector uPrev_; // the current and previous variables (primary and secondary variables) - VolumeVariablesVector curVolVarsVector_; - VolumeVariablesVector prevVolVarsVector_; + CurrentVolumeVariablesVector curVolVarsVector_; + PreviousVolumeVariablesVector prevVolVarsVector_; // the flux variables cache vector vector FluxVariablesCacheVector fluxVarsCacheVector_; diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index 2666ca42e5..f9079f30f3 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -71,8 +71,10 @@ NEW_PROP_TAG(PrimaryVariables); //!< A vector of primary variables within a sub- NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector of the grid NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume +NEW_PROP_TAG(EnableVolumeVariablesCache); //!< If disabled, the volume variables are not stored (reduces memory, but is slower) NEW_PROP_TAG(VolumeVariables); //!< The secondary variables within a sub-control volume -NEW_PROP_TAG(VolumeVariablesVector); //!< The type for a container of volume variables +NEW_PROP_TAG(CurrentVolumeVariablesVector); //!< The type for a container of current volume variables +NEW_PROP_TAG(PreviousVolumeVariablesVector); //!< The type for a container of previous volume variables NEW_PROP_TAG(FluxVariables); //!< Container storing the different types of flux variables NEW_PROP_TAG(FluxVariablesCache); //!< Stores data associated with flux vars (if enabled) NEW_PROP_TAG(FluxVariablesCacheVector); //!< The global vector of flux variable containers diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 6a500fc2b1..98ac0c8906 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -96,8 +96,11 @@ SET_TYPE_PROP(ImplicitBase, FVElementGeometry, Dumux::FVElementGeometry //! The volume variable class, to be overloaded by the model SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables); -//! The global volume variables vector class -SET_TYPE_PROP(ImplicitBase, VolumeVariablesVector, Dumux::VolumeVariablesVector); +//! The global current volume variables vector class +SET_TYPE_PROP(ImplicitBase, CurrentVolumeVariablesVector, Dumux::VolumeVariablesVector); + +//! The global previous volume variables vector class +SET_TYPE_PROP(ImplicitBase, PreviousVolumeVariablesVector, Dumux::VolumeVariablesVector); //! The global flux variables cache vector class SET_TYPE_PROP(ImplicitBase, FluxVariablesCacheVector, Dumux::FluxVariablesCacheVector); @@ -143,6 +146,9 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false); //! disable partial reassembling by default SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false); +//! We do not store the volume variables by default +SET_BOOL_PROP(ImplicitBase, EnableVolumeVariablesCache, false); + //! disable flux variables data caching by default SET_BOOL_PROP(ImplicitBase, EnableFluxVariablesCache, false); diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh index ef0bea9ee2..99af3b4132 100644 --- a/dumux/implicit/volumevariablesvector.hh +++ b/dumux/implicit/volumevariablesvector.hh @@ -32,21 +32,42 @@ namespace Dumux * \ingroup ImplicitModel * \brief Base class for the volume variables vector */ -template -class VolumeVariablesVector : public std::vector +template +class VolumeVariablesVector +{}; + +// specialization in case of storing the volume variables +template +class VolumeVariablesVector : public std::vector { + friend VolumeVariablesVector; + friend ImplicitModel; using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); - using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + VolumeVariablesVector& operator= (const VolumeVariablesVector& other) = default; + + VolumeVariablesVector& operator= (const VolumeVariablesVector& other) + { + // do the copy + numScvs_ = other.numScvs_; + numBoundaryScvs_ = other.numBoundaryScvs_; + volumeVariables_ = other.volumeVariables_; + boundaryVolumeVariables_ = other.boundaryVolumeVariables_; + + // return the existing object + return *this; + }; + public: - void update(const Problem& problem, const SolutionVector& sol) + void update(Problem& problem, const SolutionVector& sol) { numScvs_ = problem.model().fvGeometries().numScv(); if (!isBox) @@ -109,6 +130,12 @@ public: return boundaryVolumeVariables_[scvIdx-numScvs_]; } + // For compatibility reasons with the case of not storing the vol vars. + // function to be called before assembling an element, preparing the vol vars within the stencil + void bind(const Element& element) const {} + // function to prepare the vol vars within the element + void bindElement(const Element& element) const {} + private: IndexType numScvs_; IndexType numBoundaryScvs_; @@ -116,6 +143,255 @@ private: std::vector boundaryVolumeVariables_; }; + +// Specialization when the current volume variables are not stored +template +class VolumeVariablesVector +{ + // prev vol vars have to be a friend class in order for the assignment operator to work + friend VolumeVariablesVector; + friend ImplicitModel; + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; + + VolumeVariablesVector& operator= (const VolumeVariablesVector& other) = default; + + void operator= (const VolumeVariablesVector& other) + { + eIdxBound_ = -1; + } + + VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {} + + +public: + + void update(Problem& problem, const SolutionVector& sol) + { + problemPtr_ = &problem; + } + + + // Binding of an element, prepares the volume variables within the element stencil + // called by the local jacobian to prepare element assembly + // specialization for cc models + template + typename std::enable_if::type + bind(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + if (eIdx == eIdxBound_ && volVarIndices_.size() > 1) + return; + + eIdxBound_ = eIdx; + + const auto& fvGeometry = problem_().model().fvGeometries(eIdx); + + // get stencil information + const auto& elementStencil = problem_().model().stencils(element).elementStencil(); + const auto& neighborStencil = problem_().model().stencils(element).neighborStencil(); + const auto numDofs = elementStencil.size(); + + // resize volume variables to the required size + volumeVariables_.resize(numDofs); + volVarIndices_.resize(numDofs); + int localIdx = 0; + + // update the volume variables of the element at hand + const auto& scvI = problem_().model().fvGeometries().subControlVolume(eIdx); + const auto& solI = problem_().model().curSol()[eIdx]; + volumeVariables_[localIdx].update(solI, problem_(), element, scvI); + volVarIndices_[localIdx] = scvI.index(); + localIdx++; + + // Update the volume variables of the neighboring elements and the boundary + for (auto globalJ : neighborStencil) + { + const auto& elementJ = problem_().model().fvGeometries().element(globalJ); + const auto& scvJ = problem_().model().fvGeometries().subControlVolume(globalJ); + const auto& solJ = problem_().model().curSol()[globalJ]; + volumeVariables_[localIdx].update(solJ, problem_(), elementJ, scvJ); + volVarIndices_[localIdx] = scvJ.index(); + localIdx++; + } + + // Check for boundary volume variables + for (const auto& scvFace : fvGeometry.scvfs()) + { + // if we are not on a boundary, skip the rest + if (!scvFace.boundary()) + continue; + + // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries + const auto bcTypes = problem_().boundaryTypes(element, scvFace); + if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann())) + continue; + + volumeVariables_.resize(localIdx+1); + volVarIndices_.resize(localIdx+1); + const auto dirichletPriVars = problem_().dirichlet(element, scvFace); + volumeVariables_[localIdx].update(dirichletPriVars, problem_(), element, scvI); + volVarIndices_[localIdx] = scvFace.outsideScvIdx(); + localIdx++; + } + } + + // Binding of an element, prepares only the volume variables of the element + // specialization for cc models + template + typename std::enable_if::type + bindElement(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + if (eIdx == eIdxBound_ || std::find(volVarIndices_.begin(), volVarIndices_.end(), eIdx) != volVarIndices_.end()) + return; + + volumeVariables_.resize(1); + volVarIndices_.resize(1); + eIdxBound_ = eIdx; + + // update the volume variables of the element at hand + const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx); + const auto& sol = problem_().model().curSol()[eIdx]; + volumeVariables_[0].update(sol, problem_(), element, scv); + volVarIndices_[0] = scv.index(); + } + + const VolumeVariables& operator [](IndexType scvIdx) const + { + return volumeVariables_[getLocalIdx_(scvIdx)]; + } + + VolumeVariables& operator [](IndexType scvIdx) + { + return volumeVariables_[getLocalIdx_(scvIdx)]; + } + +private: + + void release_() + { + volumeVariables_.clear(); + volVarIndices_.clear(); + } + + const int getLocalIdx_(const int volVarIdx) const + { + auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx); + + if (it != volVarIndices_.end()) + return std::distance(volVarIndices_.begin(), it); + else + DUNE_THROW(Dune::InvalidStateException, "Could not find the current volume variables for volVarIdx = " << volVarIdx << + ", make sure to properly bind the volume variables to the element before using them"); + } + + Problem& problem_() const + { return *problemPtr_;} + + Problem* problemPtr_; + IndexType eIdxBound_; + std::vector volVarIndices_; + std::vector volumeVariables_; +}; + +// Specialization when the previous volume variables are not stored +template +class VolumeVariablesVector +{ + // current vol vars have to be a friend class in order for the assignment operator to work + friend VolumeVariablesVector; + friend ImplicitModel; + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; + + VolumeVariablesVector& operator= (const VolumeVariablesVector& other) + { + release_(); + problemPtr_ = other.problemPtr_; + eIdxBound_ = -1; + return *this; + } + + VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {} + + +public: + + void update(const Problem& problem, const SolutionVector& sol) + { + problemPtr_ = &problem; + } + + // Binding of an element, prepares the volume variables of only the element + // specialization for cc models + template + typename std::enable_if::type + bindElement(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + if (eIdx == eIdxBound_) + return; + + volumeVariables_.resize(1); + volVarIndices_.resize(1); + eIdxBound_ = eIdx; + + // update the volume variables of the element at hand + const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx); + const auto& sol = problem_().model().prevSol()[eIdx]; + volumeVariables_[0].update(sol, problem_(), element, scv); + volVarIndices_[0] = scv.index(); + } + + const VolumeVariables& operator [](IndexType scvIdx) const + { + return volumeVariables_[getLocalIdx_(scvIdx)]; + } + + VolumeVariables& operator [](IndexType scvIdx) + { + return volumeVariables_[getLocalIdx_(scvIdx)]; + } + +private: + + void release_() + { + volumeVariables_.clear(); + volVarIndices_.clear(); + } + + const int getLocalIdx_(const int volVarIdx) const + { + auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx); + + if (it != volVarIndices_.end()) + return std::distance(volVarIndices_.begin(), it); + else + DUNE_THROW(Dune::InvalidStateException, "Could not find the previous volume variables for volVarIdx = " << volVarIdx << + ", make sure to properly bind the volume variables to the element before using them"); + } + + Problem& problem_() const + { return *problemPtr_;} + + Problem* problemPtr_; + IndexType eIdxBound_; + std::vector volVarIndices_; + std::vector volumeVariables_; +}; + } // end namespace #endif diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh index 50c585d679..f90607ada6 100644 --- a/dumux/porousmediumflow/1p/implicit/model.hh +++ b/dumux/porousmediumflow/1p/implicit/model.hh @@ -107,7 +107,10 @@ public: int eIdx = this->elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); - for (auto&& scv : this->fvGeometries(element).scvs()) + this->curVolVars_().bindElement(element); + + const auto& fvGeometry = this->fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { const auto& spatialParams = this->problem_().spatialParams(); auto dofIdxGlobal = scv.dofIndex(); diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh index d0ad7eebdd..43670abb45 100644 --- a/dumux/porousmediumflow/2p/implicit/model.hh +++ b/dumux/porousmediumflow/2p/implicit/model.hh @@ -144,7 +144,10 @@ public: int eIdx = this->elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); - for (auto&& scv : this->fvGeometries(element).scvs()) + this->curVolVars_().bindElement(element); + + const auto& fvGeometry = this->fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { auto dofIdxGlobal = scv.dofIndex(); const auto& volVars = this->curVolVars(scv); diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh index 0048e07f59..bb7423673f 100644 --- a/dumux/porousmediumflow/3p/implicit/model.hh +++ b/dumux/porousmediumflow/3p/implicit/model.hh @@ -139,9 +139,10 @@ public: { if(element.partitionType() == Dune::InteriorEntity) { - const auto& fvGeometry = this->fvGeometries(element); + this->curVolVars_().bindElement(element); - for (auto&& scv : fvGeometry.scvs()) + const auto& fvGeometry = this->fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { auto eIdx = scv.elementIndex(); auto dofIdxGlobal = scv.dofIndex(); diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh index 16fbf7bfec..08f59ae950 100644 --- a/dumux/porousmediumflow/3p3c/implicit/model.hh +++ b/dumux/porousmediumflow/3p3c/implicit/model.hh @@ -159,7 +159,10 @@ public: { for (const auto& element : elements(this->problem_().gridView())) { - for (auto&& scv : this->fvGeometries(element).scvs()) + this->curVolVars_().bindElement(element); + + const auto& fvGeometry = this->fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { auto dofIdxGlobal = scv.dofIndex(); if (priVarSwitch_().wasSwitched(dofIdxGlobal)) @@ -274,9 +277,10 @@ public: int eIdx = this->problem_().elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); - const auto& fvGeometry = this->fvGeometries(element); + this->curVolVars_().bindElement(element); - for (auto&& scv : fvGeometry.scvs()) + const auto& fvGeometry = this->fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { const auto& volVars = this->curVolVars(scv); int dofIdxGlobal = scv.dofIndex(); diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh index 799a2a5ecf..4770e083df 100644 --- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh +++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh @@ -45,11 +45,11 @@ class PrimaryVariableSwitch using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using VolumeVariablesVector = typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector); + using CurrentVolumeVariablesVector = typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector); public: - void init(const Problem& problem) + void init(Problem& problem) { auto numDofs = problem.model().numDofs(); phasePresence_.resize(numDofs); @@ -113,15 +113,17 @@ public: * \param problem The problem * \param curSol The current solution to be updated / modified */ - bool update(const Problem& problem, SolutionVector& curSol, - VolumeVariablesVector& volVarsVector) + bool update(Problem& problem, SolutionVector& curSol, + CurrentVolumeVariablesVector& volVarsVector) { bool switched = false; visited_.assign(phasePresence_.size(), false); for (const auto& element : elements(problem.gridView())) { + volVarsVector.bindElement(element); + const auto& fvGeometry = problem.model().fvGeometries(element); - for (auto&& scv : fvGeometry.scvs()) + for (const auto& scv : fvGeometry.scvs()) { auto dofIdxGlobal = scv.dofIndex(); if (!visited_[dofIdxGlobal]) -- GitLab From 66266bc805e6a6503b3cafe51ef006681f3abe09 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 May 2016 15:02:52 +0200 Subject: [PATCH 39/46] [implicit] make FvGeometry caching optional the storing of the FVElementGeometry can be disabled which requires their on-the-fly calculation. This is done by binding the FVGeometry to an element by calling the bindElement() function. The bind() function precomputes all the geometries within the stencil. --- dumux/implicit/cellcentered/localjacobian.hh | 10 +- dumux/implicit/cellcentered/stencils.hh | 5 +- .../tpfa/fvelementgeometryvector.hh | 276 +++++++++++++++++- .../cellcentered/tpfa/propertydefaults.hh | 2 +- dumux/implicit/localresidual.hh | 6 +- dumux/implicit/model.hh | 23 +- dumux/implicit/properties.hh | 1 + dumux/implicit/propertydefaults.hh | 3 + dumux/implicit/volumevariablesvector.hh | 16 +- dumux/porousmediumflow/1p/implicit/model.hh | 2 + dumux/porousmediumflow/2p/implicit/model.hh | 2 + dumux/porousmediumflow/3p/implicit/model.hh | 2 + dumux/porousmediumflow/3p3c/implicit/model.hh | 4 + .../compositional/primaryvariableswitch.hh | 5 + .../3p3c/implicit/infiltrationproblem.hh | 7 +- 15 files changed, 331 insertions(+), 33 deletions(-) diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 53230f0b91..ceb6e1fba8 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -115,6 +115,9 @@ public: assemblyMap_.resize(problem.gridView().size(0)); for (const auto& element : elements(problem.gridView())) { + // bind FVGeometry to the element to prepare all the geometries of the stencil + problem.model().fvGeometries_().bind(element); + auto globalI = problem.elementMapper().index(element); const auto& neighborStencil = this->model_().stencils(element).neighborStencil(); @@ -151,6 +154,8 @@ public: */ void assemble(const Element& element, JacobianMatrix& matrix) { + // prepare the volvars/fvGeometries in case caching is disabled + this->model_().fvGeometries_().bind(element); this->model_().curVolVars_().bind(element); this->model_().prevVolVars_().bindElement(element); @@ -158,13 +163,14 @@ public: // finite volume geometry globalI_ = this->problem_().elementMapper().index(element); - bcTypes_.update(this->problem_(), element, fvElemGeom_()); + const auto& fvGeometry = fvElemGeom_(); + bcTypes_.update(this->problem_(), element, fvGeometry); // calculate the local residual this->localResidual().eval(element, bcTypes_); this->residual_ = this->localResidual().residual(); - this->model_().updatePVWeights(fvElemGeom_()); + this->model_().updatePVWeights(fvGeometry); // calculate derivatives of all dofs in stencil with respect to the dofs in the element evalPartialDerivatives_(element, matrix); diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 24a2deab7d..4ef579ee69 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -91,12 +91,15 @@ class CCStencilsVector using StencilType = CCElementStencils; public: - void update(const Problem& problem) + void update(Problem& problem) { problemPtr_ = &problem; elementStencils_.resize(problem.gridView().size(0)); for (const auto& element : elements(problem.gridView())) { + // bind the FvGeometry to the element before using it + problem.model().fvGeometries_().bind(element); + auto eIdx = problem.elementMapper().index(element); elementStencils_[eIdx].update(problem, element); } diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh index 88332de8cd..4d9825551a 100644 --- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh +++ b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh @@ -39,8 +39,13 @@ namespace Dumux * This builds up the sub control volumes and sub control volume faces * for each element. */ -template +template class CCTpfaFVElementGeometryVector +{}; + +// specialization in case the FVElementGeometries are stored +template +class CCTpfaFVElementGeometryVector { using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); @@ -73,7 +78,7 @@ public: //! Get a sub control volume face with a global scvf index const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const { - return *scvfs_[scvfIdx]; + return scvfs_[scvfIdx]; } //! The total number of sub control volumes @@ -130,22 +135,22 @@ public: if (!intersection.boundary()) { auto nIdx = problem.elementMapper().index(intersection.outside()); - scvfs_.push_back(std::make_shared(intersection.geometry(), - intersection.geometry().center(), - intersection.centerUnitOuterNormal(), - scvfIdx, - std::vector({eIdx, nIdx}), - false)); + scvfs_.push_back(SubControlVolumeFace(intersection.geometry(), + intersection.geometry().center(), + intersection.centerUnitOuterNormal(), + scvfIdx, + std::vector({eIdx, nIdx}), + false)); scvfsIndexSet.push_back(scvfIdx++); } else { - scvfs_.push_back(std::make_shared(intersection.geometry(), - intersection.geometry().center(), - intersection.centerUnitOuterNormal(), - scvfIdx, - std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}), - true)); + scvfs_.push_back(SubControlVolumeFace(intersection.geometry(), + intersection.geometry().center(), + intersection.centerUnitOuterNormal(), + scvfIdx, + std::vector({eIdx, gridView_.size(0) + numBoundaryScvf_++}), + true)); scvfsIndexSet.push_back(scvfIdx++); } } @@ -155,15 +160,256 @@ public: } } + // Binding of an element, called by the local jacobian to prepare element assembly + // For compatibility reasons with the FVGeometry cache disabled + void bind(const Element& element) {} + void bindElement(const Element& element) {} + private: GridView gridView_; Dumux::ElementMap elementMap_; std::vector> scvs_; - std::vector> scvfs_; + std::vector scvfs_; std::vector fvGeometries_; IndexType numBoundaryScvf_; }; +// specialization in case the FVElementGeometries are not stored +template +class CCTpfaFVElementGeometryVector +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using IndexType = typename GridView::IndexSet::IndexType; + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using Element = typename GridView::template Codim<0>::Entity; + +public: + //! Constructor + CCTpfaFVElementGeometryVector(const GridView gridView) + : gridView_(gridView), elementMap_(gridView), eIdxBound_(-1) {} + + /* \brief Get the finite volume geometry of an element + * \note The finite volume geometry offers iterators over the sub control volumes + * and the sub control volume faces of an element. + */ + const FVElementGeometry& fvGeometry(IndexType eIdx) const + { + return fvGeometries_[getStencilScvIdx_(eIdx)]; + } + + const FVElementGeometry& fvGeometry(const Element& element) const + { + return fvGeometries_[getStencilScvIdx_(problem_().elementMapper().index(element))]; + } + + //! Get a sub control volume with a global scv index + const SubControlVolume& subControlVolume(IndexType scvIdx) const + { + return stencilScvs_[getStencilScvIdx_(scvIdx)]; + } + + //! Get a sub control volume face with a global scvf index + const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const + { + return stencilScvfs_[getStencilScvFaceIdx_(scvfIdx)]; + } + + //! The total number of sub control volumes + std::size_t numScv() const + { + return numScvs_; + } + + //! The total number of sub control volume faces + std::size_t numScvf() const + { + return numScvf_; + } + + //! The total number of boundary sub control volume faces + std::size_t numBoundaryScvf() const + { + return numBoundaryScvf_; + } + + // Get an element from a sub control volume contained in it + Element element(const SubControlVolume& scv) const + { return elementMap_.element(scv.elementIndex()); } + + // Get an element from a global element index + Element element(IndexType eIdx) const + { return elementMap_.element(eIdx); } + + //! update all fvElementGeometries (do this again after grid adaption) + void update(const Problem& problem) + { + problemPtr_ = &problem; + elementMap_.clear(); + scvFaceToElementMap_.clear(); + eIdxBound_ = -1; + + // Build the SCV and SCV faces + numScvf_ = 0; + numBoundaryScvf_ = 0; + numScvs_ = gridView_.size(0); + elementMap_.resize(numScvs_); + scvFaceIndices_.resize(numScvs_); + neighborVolVarIndices_.resize(numScvs_); + for (const auto& element : elements(gridView_)) + { + auto eIdx = problem.elementMapper().index(element); + + // fill the element map with seeds + elementMap_[eIdx] = element.seed(); + + // the element-wise index sets for finite volume geometry + std::vector scvfsIndexSet; + std::vector neighborVolVarIndexSet; + for (const auto& intersection : intersections(gridView_, element)) + { + scvFaceToElementMap_.push_back(eIdx); + scvfsIndexSet.push_back(numScvf_++); + if (!intersection.boundary()) + { + auto nIdx = problem.elementMapper().index(intersection.outside()); + neighborVolVarIndexSet.push_back(nIdx); + } + else + neighborVolVarIndexSet.push_back(numScvs_ + numBoundaryScvf_++); + } + + // store the sets of indices in the data container + scvFaceIndices_[eIdx] = scvfsIndexSet; + neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet; + } + } + + // Binding of an element preparing the geometries of the whole stencil + // called by the local jacobian to prepare element assembly + void bind(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + if (eIdxBound_ == eIdx && fvGeometries_.size() > 1) + return; + + release_(); + eIdxBound_ = eIdx; + + makeElementGeometries_(element); + for (const auto& intersection : intersections(gridView_, element)) + { + if (!intersection.boundary()) + makeElementGeometries_(intersection.outside()); + } + } + + // Binding of an element preparing the geometries only inside the element + void bindElement(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + if (eIdx == eIdxBound_ || std::find(stencilScvIndices_.begin(), stencilScvIndices_.end(), eIdx) != stencilScvIndices_.end() ) + return; + + release_(); + eIdxBound_ = eIdx; + + makeElementGeometries_(element); + } + + +private: + + void release_() + { + stencilScvs_.clear(); + stencilScvIndices_.clear(); + stencilScvfs_.clear(); + stencilScvFaceIndices_.clear(); + fvGeometries_.clear(); + } + + void makeElementGeometries_(const Element& element) + { + const auto eIdx = problem_().elementMapper().index(element); + + std::vector scv({SubControlVolume(element.geometry(), eIdx)}); + stencilScvs_.push_back(scv[0]); + stencilScvIndices_.push_back(eIdx); + + const auto& scvFaceIndices = scvFaceIndices_[eIdx]; + const auto& neighborVolVarIndices = neighborVolVarIndices_[eIdx]; + + int scvfCounter = 0; + for (const auto& intersection : intersections(gridView_, element)) + { + stencilScvfs_.push_back(SubControlVolumeFace(intersection.geometry(), + intersection.geometry().center(), + intersection.centerUnitOuterNormal(), + scvFaceIndices[scvfCounter], + std::vector({eIdx, neighborVolVarIndices[scvfCounter]}), + intersection.boundary())); + + stencilScvFaceIndices_.push_back(scvFaceIndices[scvfCounter]); + scvfCounter++; + } + + // Compute the finite volume element geometries + fvGeometries_.push_back(FVElementGeometry(*this, {eIdx}, scvFaceIndices)); + } + + const int getStencilScvIdx_(const int scvIdx) const + { + auto it = std::find(stencilScvIndices_.begin(), stencilScvIndices_.end(), scvIdx); + + if (it != stencilScvIndices_.end()) + return std::distance(stencilScvIndices_.begin(), it); + else + DUNE_THROW(Dune::InvalidStateException, + "Could not find the sub control volume for scvIdx = " << scvIdx << + ", make sure to properly bind the FVGeometries to the element before using them."); + } + + const int getStencilScvFaceIdx_(const int scvfIdx) const + { + auto it = std::find(stencilScvFaceIndices_.begin(), stencilScvFaceIndices_.end(), scvfIdx); + + if (it != stencilScvFaceIndices_.end()) + return std::distance(stencilScvFaceIndices_.begin(), it); + else + DUNE_THROW(Dune::InvalidStateException, + "Could not find the sub control volume face for scvfIdx = " << scvfIdx << + ", make sure to properly bind the FVGeometries to the to the element before using them."); + } + + const Problem& problem_() const + { return *problemPtr_; } + + const Problem* problemPtr_; + GridView gridView_; + + // Information on the global number of geometries + IndexType numScvs_; + IndexType numScvf_; + IndexType numBoundaryScvf_; + + // vectors that stores global data + Dumux::ElementMap elementMap_; + std::vector> scvFaceIndices_; + std::vector> neighborVolVarIndices_; + std::vector scvFaceToElementMap_; + + // vectors to store the geometries temporarily after binding an element + IndexType eIdxBound_; + std::vector stencilScvs_; + std::vector stencilScvIndices_; + std::vector stencilScvfs_; + std::vector stencilScvFaceIndices_; + std::vector fvGeometries_; +}; + } // end namespace #endif diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index 3259a4bf45..348ad210a3 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -42,7 +42,7 @@ template class CCElementBoundaryTypes; namespace Properties { //! Set the default for the FVElementGeometry vector -SET_TYPE_PROP(CCTpfaModel, FVElementGeometryVector, CCTpfaFVElementGeometryVector); +SET_TYPE_PROP(CCTpfaModel, FVElementGeometryVector, CCTpfaFVElementGeometryVector); //! The sub control volume SET_PROP(CCTpfaModel, SubControlVolume) diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index a5a4c79319..7073f0220b 100644 --- a/dumux/implicit/localresidual.hh +++ b/dumux/implicit/localresidual.hh @@ -92,6 +92,7 @@ public: void eval(const Element &element) { // make sure FVElementGeometry and volume variables are bound to the element + model_().fvGeometries_().bind(element); model_().curVolVars_().bind(element); model_().prevVolVars_().bindElement(element); @@ -115,6 +116,7 @@ public: elemPtr_ = &element; // make sure FVElementGeometry and volume variables are bound to the element + model_().fvGeometries_().bindElement(element); model_().curVolVars_().bindElement(element); model_().prevVolVars_().bindElement(element); @@ -138,6 +140,7 @@ public: elemPtr_ = &element; // make sure FVElementGeometry and volume variables are bound to the element + model_().fvGeometries_().bind(element); model_().curVolVars_().bind(element); model_().prevVolVars_().bindElement(element); @@ -170,6 +173,7 @@ public: void eval(const Element &element, const ElementBoundaryTypes &bcTypes) { + // At this point the fv geometry has to be bound already elemPtr_ = &element; bcTypesPtr_ = &bcTypes; @@ -451,7 +455,7 @@ protected: * \brief Returns a reference to the current element's finite * volume geometry. */ - const FVElementGeometry &fvGeometry_() const + const FVElementGeometry& fvGeometry_() const { return model_().fvGeometries(element_()); } diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 8d09e4de0d..d666958efe 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -119,8 +119,8 @@ public: updateBoundaryIndices_(); - fvGeometries_ = std::make_shared(gridView_()); - fvGeometries_->update(problem_()); + fvGeometryVector_ = std::make_shared(gridView_()); + fvGeometryVector_->update(problem_()); int numDofs = asImp_().numDofs(); uCur_.resize(numDofs); @@ -803,13 +803,13 @@ public: { return prevVolVarsVector_[scvIdx]; } const FVElementGeometryVector& fvGeometries() const - { return *fvGeometries_; } + { return *fvGeometryVector_; } const FVElementGeometry& fvGeometries(const Element& element) const - { return fvGeometries_->fvGeometry(elementMapper().index(element)); } + { return fvGeometryVector_->fvGeometry(elementMapper().index(element)); } const FVElementGeometry& fvGeometries(unsigned int eIdx) const - { return fvGeometries_->fvGeometry(eIdx); } + { return fvGeometryVector_->fvGeometry(eIdx); } protected: @@ -837,6 +837,9 @@ protected: FluxVariablesCache& fluxVarsCache_(unsigned int scvfIdx) { return fluxVarsCacheVector_[scvfIdx]; } + FVElementGeometryVector& fvGeometries_() + { return *fvGeometryVector_; } + /*! * \brief A reference to the problem on which the model is applied. */ @@ -876,8 +879,12 @@ protected: // condition at the center of each sub control volume for (const auto& element : elements(gridView_())) { - // deal with the current element - for(auto&& scv : fvGeometries(element).scvs()) + // deal with the current element, bind FVGeometry to it + fvGeometries_().bindElement(element); + const auto& fvGeometry = fvGeometries(element); + + // loop over sub control volumes + for(const auto& scv : fvGeometry.scvs()) { // let the problem do the dirty work of nailing down // the initial solution. @@ -995,7 +1002,7 @@ protected: StencilsVector stencilsVector_; // the finite volume element geometries - std::shared_ptr fvGeometries_; + std::shared_ptr fvGeometryVector_; Dune::BlockVector > boxVolume_; diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index f9079f30f3..6218686aa6 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -61,6 +61,7 @@ NEW_PROP_TAG(SubControlVolume);//!< The type of the sub control volume NEW_PROP_TAG(SubControlVolumeFace); //!< The type of the sub control volume face NEW_PROP_TAG(FVElementGeometry); //!< The type of the finite volume geometry (iterators over scvs, scvfs) NEW_PROP_TAG(FVElementGeometryVector); //!< The type of the finite volume geometry vector +NEW_PROP_TAG(EnableFVElementGeometryCache); //! specifies if geometric data is be saved (faster, but more memory consuming) NEW_PROP_TAG(JacobianAssembler); //!< Assembles the global jacobian matrix NEW_PROP_TAG(JacobianMatrix); //!< Type of the global jacobian matrix diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 98ac0c8906..375b0c46ef 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -146,6 +146,9 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false); //! disable partial reassembling by default SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false); +//! We do not store the FVGeometry by default +SET_BOOL_PROP(ImplicitBase, EnableFVElementGeometryCache, false); + //! We do not store the volume variables by default SET_BOOL_PROP(ImplicitBase, EnableVolumeVariablesCache, false); diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh index 99af3b4132..567e560d54 100644 --- a/dumux/implicit/volumevariablesvector.hh +++ b/dumux/implicit/volumevariablesvector.hh @@ -79,7 +79,9 @@ public: boundaryVolumeVariables_.resize(numBoundaryScvs_); for (const auto& element : elements(problem.gridView())) { - for (auto&& scv : problem.model().fvGeometries(element).scvs()) + problem.model().fvGeometries_().bind(element); + const auto& fvGeometry = problem.model().fvGeometries(element); + for (const auto& scv : fvGeometry.scvs()) { (*this)[scv.index()].update(sol[scv.dofIndex()], problem, @@ -89,14 +91,14 @@ public: if (!isBox) { - for (auto&& scvFace : problem.model().fvGeometries(element).scvfs()) + for (const auto& scvFace : fvGeometry.scvfs()) { // if we are not on a boundary, skip the rest if (!scvFace.boundary()) continue; // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries - auto bcTypes = problem.boundaryTypes(element, scvFace); + const auto bcTypes = problem.boundaryTypes(element, scvFace); if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann())) continue; @@ -190,6 +192,8 @@ public: eIdxBound_ = eIdx; + // make sure the FVElementGeometry is bound to the element + problem_().model().fvGeometries_().bind(element); const auto& fvGeometry = problem_().model().fvGeometries(eIdx); // get stencil information @@ -255,6 +259,9 @@ public: volVarIndices_.resize(1); eIdxBound_ = eIdx; + // make sure the FVElementGeometry is bound to the element + problem_().model().fvGeometries_().bindElement(element); + // update the volume variables of the element at hand const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx); const auto& sol = problem_().model().curSol()[eIdx]; @@ -347,6 +354,9 @@ public: volVarIndices_.resize(1); eIdxBound_ = eIdx; + // make sure FVElementGeometry is bound to the element + problem_().model().fvGeometries_().bindElement(element); + // update the volume variables of the element at hand const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx); const auto& sol = problem_().model().prevSol()[eIdx]; diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh index f90607ada6..6b6684a34f 100644 --- a/dumux/porousmediumflow/1p/implicit/model.hh +++ b/dumux/porousmediumflow/1p/implicit/model.hh @@ -107,6 +107,8 @@ public: int eIdx = this->elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); + // make sure FVElementGeometry and the volume variables are bound to the element + this->fvGeometries_().bindElement(element); this->curVolVars_().bindElement(element); const auto& fvGeometry = this->fvGeometries(element); diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh index 43670abb45..4014bece2f 100644 --- a/dumux/porousmediumflow/2p/implicit/model.hh +++ b/dumux/porousmediumflow/2p/implicit/model.hh @@ -144,6 +144,8 @@ public: int eIdx = this->elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); + // make sure FVElementGeometry & vol vars are bound to the element + this->fvGeometries_().bindElement(element); this->curVolVars_().bindElement(element); const auto& fvGeometry = this->fvGeometries(element); diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh index bb7423673f..697e796f03 100644 --- a/dumux/porousmediumflow/3p/implicit/model.hh +++ b/dumux/porousmediumflow/3p/implicit/model.hh @@ -139,6 +139,8 @@ public: { if(element.partitionType() == Dune::InteriorEntity) { + // make sure FVElementGeometry & vol vars are bound to the element + this->fvGeometries_().bindElement(element); this->curVolVars_().bindElement(element); const auto& fvGeometry = this->fvGeometries(element); diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh index 08f59ae950..5fdaea88fc 100644 --- a/dumux/porousmediumflow/3p3c/implicit/model.hh +++ b/dumux/porousmediumflow/3p3c/implicit/model.hh @@ -159,6 +159,8 @@ public: { for (const auto& element : elements(this->problem_().gridView())) { + // make sure FVElementGeometry & vol vars are bound to the element + this->fvGeometries_().bindElement(element); this->curVolVars_().bindElement(element); const auto& fvGeometry = this->fvGeometries(element); @@ -277,6 +279,8 @@ public: int eIdx = this->problem_().elementMapper().index(element); (*rank)[eIdx] = this->gridView_().comm().rank(); + // make sure FVElementGeometry & vol vars are bound to the element + this->fvGeometries_().bindElement(element); this->curVolVars_().bindElement(element); const auto& fvGeometry = this->fvGeometries(element); diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh index 4770e083df..1eae984eeb 100644 --- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh +++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh @@ -57,6 +57,9 @@ public: for (const auto& element : elements(problem.gridView())) { + // make sure FVElementGeometry is bound to the element + problem.model().fvGeometries_().bindElement(element); + const auto& fvGeometry = problem.model().fvGeometries(element); for (auto&& scv : fvGeometry.scvs()) { @@ -120,6 +123,8 @@ public: visited_.assign(phasePresence_.size(), false); for (const auto& element : elements(problem.gridView())) { + // make sure FVElementGeometry is bound to the element + problem.model().fvGeometries_().bindElement(element); volVarsVector.bindElement(element); const auto& fvGeometry = problem.model().fvGeometries(element); diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh index d639edadbe..786220b7d3 100644 --- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh +++ b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh @@ -345,8 +345,11 @@ public: for (const auto& element : elements(this->gridView())) { - auto fvGeometry = this->model().fvGeometries(element); - for (auto&& scv : fvGeometry.scvs()) + // make sure the FVElementGeometry is bound to the element + this->model().fvGeometries_().bindElement(element); + const auto& fvGeometry = this->model().fvGeometries(element); + + for (const auto& scv : fvGeometry.scvs()) { auto dofIdxGlobal = scv.dofIndex(); (*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(scv); -- GitLab From 6d2405e6b96992677bd0a2fba9a6ebce03daaba8 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 May 2016 15:06:27 +0200 Subject: [PATCH 40/46] [implicit] change type in range-based loops over geometries --- .../implicit/cellcentered/elementboundarytypes.hh | 4 ++-- dumux/implicit/cellcentered/localjacobian.hh | 7 ++++--- dumux/implicit/cellcentered/stencils.hh | 9 ++++++--- dumux/implicit/localresidual.hh | 15 ++++++++++----- 4 files changed, 22 insertions(+), 13 deletions(-) diff --git a/dumux/implicit/cellcentered/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh index de787f440a..db59ac6170 100644 --- a/dumux/implicit/cellcentered/elementboundarytypes.hh +++ b/dumux/implicit/cellcentered/elementboundarytypes.hh @@ -87,12 +87,12 @@ public: (*this)[0].reset(); - for (auto&& scv : fvGeometry.scvs()) + for (const auto& scv : fvGeometry.scvs()) { if (!problem.model().onBoundary(scv)) return; - for (auto&& scvFace : fvGeometry.scvfs()) + for (const auto& scvFace : fvGeometry.scvfs()) { if (!scvFace.boundary()) continue; diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index ceb6e1fba8..9aa326132f 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -122,12 +122,13 @@ public: const auto& neighborStencil = this->model_().stencils(element).neighborStencil(); assemblyMap_[globalI].reserve(neighborStencil.size()); - - for (auto&& globalJ : neighborStencil) + for (auto globalJ : neighborStencil) { // find the flux vars needed for the calculation of the flux into element std::vector fluxVarIndices; - for (auto&& scvFaceJ : this->model_().fvGeometries(globalJ).scvfs()) + + const auto& fvGeometry = this->model_().fvGeometries(globalJ); + for (const auto& scvFaceJ : fvGeometry.scvfs()) { auto fluxVarsIdx = scvFaceJ.index(); const auto& elementJ = this->model_().fvGeometries().element(globalJ); diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 4ef579ee69..65e506dab1 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -46,11 +46,14 @@ public: void update(const Problem& problem, const Element& element) { elementStencil_.clear(); - for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) + + const auto& fvGeometry = problem.model().fvGeometries(element); + // loop over sub control faces + for (const auto& scvf : fvGeometry.scvfs()) { FluxVariables fluxVars; - fluxVars.init(problem, element, scvf); - auto&& stencil = fluxVars.stencil(); + const auto& stencil = fluxVars.computeStencil(problem, scvf); + elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end()); } // make values in elementstencil unique diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index 7073f0220b..10467405d3 100644 --- a/dumux/implicit/localresidual.hh +++ b/dumux/implicit/localresidual.hh @@ -266,7 +266,8 @@ protected: void evalFluxes_() { // calculate the mass flux over the scv faces and subtract - for (auto&& scvf : fvGeometry_().scvfs()) + const auto& fvGeometry = fvGeometry_(); + for (const auto& scvf : fvGeometry.scvfs()) { auto flux = asImp_().computeFlux_(scvf); @@ -322,7 +323,8 @@ protected: // calculate the amount of conservation each quantity inside // all sub control volumes - for (auto&& scv : fvGeometry_().scvs()) + const auto& fvGeometry = fvGeometry_(); + for (const auto& scv : fvGeometry.scvs()) { auto scvIdx = scv.indexInElement(); @@ -334,7 +336,8 @@ protected: PrimaryVariables evalSource_() { PrimaryVariables source(0); - for (auto&& scv : fvGeometry_().scvs()) + const auto& fvGeometry = fvGeometry_(); + for (const auto& scv : fvGeometry) { source += this->problem_().source(element_(), scv); @@ -355,7 +358,8 @@ protected: evalVolumeTerms_() { // evaluate the volume terms (storage + source terms) - for (auto&& scv : fvGeometry_().scvs()) + const auto& fvGeometry = fvGeometry_(); + for (const auto& scv : fvGeometry.scvs()) { auto scvIdx = scv.indexInElement(); auto curExtrusionFactor = model_().curVolVars(scv).extrusionFactor(); @@ -378,7 +382,8 @@ protected: evalVolumeTerms_() { // evaluate the volume terms (storage + source terms) - for (auto&& scv : fvGeometry_().scvs()) + const auto& fvGeometry = fvGeometry_(); + for (const auto& scv : fvGeometry.scvs()) { auto scvIdx = scv.indexInElement(); auto prevExtrusionFactor = model_().prevVolVars(scv).extrusionFactor(); -- GitLab From 3c1de044fda4761eb38d90e99e6440607254ab03 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 May 2016 15:07:40 +0200 Subject: [PATCH 41/46] [implicit][cc][localresidual] use of convenience function --- dumux/implicit/cellcentered/localresidual.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index d599045821..c9e0b970b9 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -221,7 +221,7 @@ protected: { // update corresponding boundary volume variables before flux calculation const auto insideScvIdx = scvf.insideScvIdx(); - const auto& insideScv = this->problem_().model().fvGeometries().subControlVolume(insideScvIdx); + const auto& insideScv = this->model_().fvGeometries().subControlVolume(insideScvIdx); const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf); this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv); } @@ -248,7 +248,7 @@ protected: PrimaryVariables dirichletValues = this->problem_().dirichlet(this->element_(), scvf); // get the primary variables - const auto& priVars = this->problem_().model().curVolVars(scvf.insideScvIdx()).priVars(); + const auto& priVars = this->model_().curVolVars(scvf.insideScvIdx()).priVars(); // set Dirichlet conditions in a strong sense for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) -- GitLab From 10b5f7de3c32ae2598df8322546949dbb7142055 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Mon, 23 May 2016 15:08:20 +0200 Subject: [PATCH 42/46] [implicit][properties] rearrange property declaration --- dumux/implicit/properties.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index 6218686aa6..cb9e36510f 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -79,12 +79,12 @@ NEW_PROP_TAG(PreviousVolumeVariablesVector); //!< The type for a container of p NEW_PROP_TAG(FluxVariables); //!< Container storing the different types of flux variables NEW_PROP_TAG(FluxVariablesCache); //!< Stores data associated with flux vars (if enabled) NEW_PROP_TAG(FluxVariablesCacheVector); //!< The global vector of flux variable containers +NEW_PROP_TAG(EnableFluxVariablesCache); //! specifies if data on flux vars should be saved (faster, but more memory consuming) NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces in cc models(outflow) NEW_PROP_TAG(ConstantBoundaryConditions); //!< boundary data is stored in case the BC are constant // Specify the forms of fluxes that should be considered in the model // also, specify their corresponding flux variables -NEW_PROP_TAG(EnableFluxVariablesCache); //! specifies if data on flux vars should be saved (faster, but more memory consuming) NEW_PROP_TAG(EnableAdvection); //! specifies if advection is considered in the model NEW_PROP_TAG(AdvectionType); //! The type for the calculation the advective fluxes NEW_PROP_TAG(EnableMolecularDiffusion); //! specifies if molecular diffusive fluxes are considered in the model -- GitLab From 64fd3e40f9ac194f6bc8613b323fb208cd68d19b Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 27 May 2016 13:23:20 +0200 Subject: [PATCH 43/46] [cc] Fix cc assembly for parallel runs --- dumux/implicit/assembler.hh | 28 +------- dumux/implicit/cellcentered/localjacobian.hh | 64 +++++++++++-------- dumux/implicit/cellcentered/stencils.hh | 12 +++- .../tpfa/fvelementgeometryvector.hh | 20 ++++-- 4 files changed, 65 insertions(+), 59 deletions(-) diff --git a/dumux/implicit/assembler.hh b/dumux/implicit/assembler.hh index 2193aaaf96..dd5ffdd40a 100644 --- a/dumux/implicit/assembler.hh +++ b/dumux/implicit/assembler.hh @@ -152,31 +152,9 @@ protected: { asImp_().resetSystem_(); - // reassemble the elements... - for (const auto& element : elements(gridView_())) { - if (element.partitionType() == Dune::GhostEntity) - { - asImp_().assembleGhostElement_(element); - } - else - { - asImp_().assembleElement_(element); - } - } - } - // "assemble" a ghost element - void assembleGhostElement_(const Element &element) - { - int globalI = elementMapper_().index(element); - - // update the right hand side - residual_[globalI] = 0.0; - - // update the diagonal entry - typedef typename JacobianMatrix::block_type BlockType; - BlockType &J = (*matrix_)[globalI][globalI]; - for (int j = 0; j < BlockType::rows; ++j) - J[j][j] = 1.0; + // assemble the elements... + for (const auto& element : elements(gridView_())) + asImp_().assembleElement_(element); } protected: diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh index 9aa326132f..0cf4858490 100644 --- a/dumux/implicit/cellcentered/localjacobian.hh +++ b/dumux/implicit/cellcentered/localjacobian.hh @@ -124,14 +124,14 @@ public: assemblyMap_[globalI].reserve(neighborStencil.size()); for (auto globalJ : neighborStencil) { + const auto& fvGeometry = this->model_().fvGeometries(globalJ); + const auto& elementJ = this->model_().fvGeometries().element(globalJ); + // find the flux vars needed for the calculation of the flux into element std::vector fluxVarIndices; - - const auto& fvGeometry = this->model_().fvGeometries(globalJ); for (const auto& scvFaceJ : fvGeometry.scvfs()) { auto fluxVarsIdx = scvFaceJ.index(); - const auto& elementJ = this->model_().fvGeometries().element(globalJ); // if globalI is in flux var stencil, add to list FluxVariables fluxVars; @@ -155,6 +155,8 @@ public: */ void assemble(const Element& element, JacobianMatrix& matrix) { + const bool isGhost = (element.partitionType() == Dune::GhostEntity); + // prepare the volvars/fvGeometries in case caching is disabled this->model_().fvGeometries_().bind(element); this->model_().curVolVars_().bind(element); @@ -168,13 +170,20 @@ public: bcTypes_.update(this->problem_(), element, fvGeometry); // calculate the local residual - this->localResidual().eval(element, bcTypes_); - this->residual_ = this->localResidual().residual(); + if (isGhost) + { + this->residual_ = 0.0; + } + else + { + this->localResidual().eval(element, bcTypes_); + this->residual_ = this->localResidual().residual(); + } this->model_().updatePVWeights(fvGeometry); // calculate derivatives of all dofs in stencil with respect to the dofs in the element - evalPartialDerivatives_(element, matrix); + evalPartialDerivatives_(element, matrix, isGhost); // TODO: calculate derivatives in the case of an extended source stencil // const auto& extendedSourceStencil = model_().stencils(element).extendedSourceStencil(); @@ -191,14 +200,14 @@ public: // } } - void evalPartialDerivatives_(const Element& element, JacobianMatrix& matrix) + void evalPartialDerivatives_(const Element& element, JacobianMatrix& matrix, const bool isGhost) { // get stencil informations const auto& neighborStencil = this->model_().stencils(element).neighborStencil(); const auto numNeighbors = neighborStencil.size(); // derivatives in the element and the neighbors - PrimaryVariables partialDeriv(0.0); + PrimaryVariables partialDeriv(isGhost ? 1.0 : 0.0); Dune::BlockVector neighborDeriv(numNeighbors); // the localresidual class used for the flux calculations @@ -248,11 +257,14 @@ public: // update the volume variables this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv); - // calculate the residual with the deflected primary variables - this->localResidual().eval(element, bcTypes_); + if (!isGhost) + { + // calculate the residual with the deflected primary variables + this->localResidual().eval(element, bcTypes_); - // store the residual and the storage term - partialDeriv = this->localResidual().residual(0); + // store the residual and the storage term + partialDeriv = this->localResidual().residual(0); + } // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) @@ -266,7 +278,8 @@ public: // we are using backward differences, i.e. we don't need // to calculate f(x + \epsilon) and we can recycle the // (already calculated) residual f(x) - partialDeriv = this->residual_[0]; + if (!isGhost) + partialDeriv = this->residual_[0]; neighborDeriv = origFlux; } @@ -282,11 +295,14 @@ public: // update the volume variables this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv); - // calculate the residual with the deflected primary variables - this->localResidual().eval(element, bcTypes_); + if (!isGhost) + { + // calculate the residual with the deflected primary variables + this->localResidual().eval(element, bcTypes_); - // subtract the residual from the derivative storage - partialDeriv -= this->localResidual().residual(0); + // subtract the residual from the derivative storage + partialDeriv -= this->localResidual().residual(0); + } // calculate the fluxes into element with the deflected primary variables for (std::size_t k = 0; k < numNeighbors; ++k) @@ -300,26 +316,20 @@ public: // we are using forward differences, i.e. we don't need to // calculate f(x - \epsilon) and we can recycle the // (already calculated) residual f(x) - partialDeriv -= this->residual_[0]; + if (!isGhost) + partialDeriv -= this->residual_[0]; neighborDeriv -= origFlux; } // divide difference in residuals by the magnitude of the // deflections between the two function evaluation - partialDeriv /= delta; + if (!isGhost) + partialDeriv /= delta; neighborDeriv /= delta; // restore the original state of the scv's volume variables this->model_().curVolVars_(scv) = std::move(origVolVars); -#if HAVE_VALGRIND - for (unsigned i = 0; i < partialDeriv.size(); ++i) - Valgrind::CheckDefined(partialDeriv[i]); - for (int i = 0; i < neighborDeriv.size(); ++idx) - for(unsigned k = 0; k < neighborDeriv[i].size(); ++k) - Valgrind::CheckDefined(neighborDeriv[i][k]); -#endif - // update the global jacobian matrix with the current partial derivatives this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv); diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh index 65e506dab1..c0dcfea856 100644 --- a/dumux/implicit/cellcentered/stencils.hh +++ b/dumux/implicit/cellcentered/stencils.hh @@ -62,7 +62,17 @@ public: auto globalI = problem.elementMapper().index(element); neighborStencil_ = elementStencil_; - neighborStencil_.erase(std::remove(neighborStencil_.begin(), neighborStencil_.end(), globalI), neighborStencil_.end()); + + //remove the element itself and possible ghost neighbors from the stencil + auto pred = [&problem, globalI](const int i) -> bool + { + if (i == globalI) + return true; + if (problem.model().fvGeometries().element(i).partitionType() == Dune::GhostEntity) + return true; + return false; + }; + neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), pred), neighborStencil_.end()); } //! The full element stencil (all element this element is interacting with) diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh index 4d9825551a..083f088488 100644 --- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh +++ b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh @@ -132,7 +132,8 @@ public: std::vector scvfsIndexSet; for (const auto& intersection : intersections(gridView_, element)) { - if (!intersection.boundary()) + // inner sub control volume faces + if (intersection.neighbor()) { auto nIdx = problem.elementMapper().index(intersection.outside()); scvfs_.push_back(SubControlVolumeFace(intersection.geometry(), @@ -143,7 +144,8 @@ public: false)); scvfsIndexSet.push_back(scvfIdx++); } - else + // boundary sub control volume faces + else if (intersection.boundary()) { scvfs_.push_back(SubControlVolumeFace(intersection.geometry(), intersection.geometry().center(), @@ -270,15 +272,21 @@ public: std::vector neighborVolVarIndexSet; for (const auto& intersection : intersections(gridView_, element)) { - scvFaceToElementMap_.push_back(eIdx); - scvfsIndexSet.push_back(numScvf_++); - if (!intersection.boundary()) + // inner sub control volume faces + if (intersection.neighbor()) { + scvFaceToElementMap_.push_back(eIdx); + scvfsIndexSet.push_back(numScvf_++); auto nIdx = problem.elementMapper().index(intersection.outside()); neighborVolVarIndexSet.push_back(nIdx); } - else + // boundary sub control volume faces + else if (intersection.boundary()) + { + scvFaceToElementMap_.push_back(eIdx); + scvfsIndexSet.push_back(numScvf_++); neighborVolVarIndexSet.push_back(numScvs_ + numBoundaryScvf_++); + } } // store the sets of indices in the data container -- GitLab From bf85a7046d21f3a82495fe8742ba9fef0cf8ef74 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 27 May 2016 13:23:46 +0200 Subject: [PATCH 44/46] Remove some unnecessary debug statements --- dumux/implicit/localresidual.hh | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index 10467405d3..4f9502b9be 100644 --- a/dumux/implicit/localresidual.hh +++ b/dumux/implicit/localresidual.hh @@ -186,19 +186,7 @@ public: storageTerm_ = 0.0; asImp_().evalFluxes_(); - -#if !defined NDEBUG && HAVE_VALGRIND - for (int i=0; i < fvGeometry_().numScv; i++) - Valgrind::CheckDefined(residual_[i]); -#endif // HAVE_VALGRIND - - asImp_().evalVolumeTerms_(); - -#if !defined NDEBUG && HAVE_VALGRIND - for (int i=0; i < fvGeometry_().numScv; i++) { - Valgrind::CheckDefined(residual_[i]); - } -#endif // HAVE_VALGRIND + asImp_().evalVolumeTerms_(); } /*! -- GitLab From ff5c916b538319cbb5d33172ffcde03991a850c1 Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 27 May 2016 15:27:06 +0200 Subject: [PATCH 45/46] [implicit][tpfa] fix parallel for geometry caching off --- .../tpfa/fvelementgeometryvector.hh | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh index 083f088488..96081fd91c 100644 --- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh +++ b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh @@ -309,7 +309,7 @@ public: makeElementGeometries_(element); for (const auto& intersection : intersections(gridView_, element)) { - if (!intersection.boundary()) + if (intersection.neighbor()) makeElementGeometries_(intersection.outside()); } } @@ -353,15 +353,18 @@ private: int scvfCounter = 0; for (const auto& intersection : intersections(gridView_, element)) { - stencilScvfs_.push_back(SubControlVolumeFace(intersection.geometry(), - intersection.geometry().center(), - intersection.centerUnitOuterNormal(), - scvFaceIndices[scvfCounter], - std::vector({eIdx, neighborVolVarIndices[scvfCounter]}), - intersection.boundary())); - - stencilScvFaceIndices_.push_back(scvFaceIndices[scvfCounter]); - scvfCounter++; + if (intersection.neighbor() || intersection.boundary()) + { + stencilScvfs_.push_back(SubControlVolumeFace(intersection.geometry(), + intersection.geometry().center(), + intersection.centerUnitOuterNormal(), + scvFaceIndices[scvfCounter], + std::vector({eIdx, neighborVolVarIndices[scvfCounter]}), + intersection.boundary())); + + stencilScvFaceIndices_.push_back(scvFaceIndices[scvfCounter]); + scvfCounter++; + } } // Compute the finite volume element geometries -- GitLab From e86615c3e3e517ea3474446c768e450c12fda06e Mon Sep 17 00:00:00 2001 From: DennisGlaeser Date: Fri, 27 May 2016 15:28:28 +0200 Subject: [PATCH 46/46] [implicit] minimize fvgeometry binding in update routine --- dumux/implicit/volumevariablesvector.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh index 567e560d54..45d99cd44a 100644 --- a/dumux/implicit/volumevariablesvector.hh +++ b/dumux/implicit/volumevariablesvector.hh @@ -79,7 +79,7 @@ public: boundaryVolumeVariables_.resize(numBoundaryScvs_); for (const auto& element : elements(problem.gridView())) { - problem.model().fvGeometries_().bind(element); + problem.model().fvGeometries_().bindElement(element); const auto& fvGeometry = problem.model().fvGeometries(element); for (const auto& scv : fvGeometry.scvs()) { -- GitLab