From f21ca0c63e6176615be33ceb28de1836977a9528 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Jun 2022 14:05:02 +0200 Subject: [PATCH] Remove unneeded class. --- .../common/fvbaselocalresidualtpfa.hh | 649 ------------------ 1 file changed, 649 deletions(-) delete mode 100644 opm/models/discretization/common/fvbaselocalresidualtpfa.hh diff --git a/opm/models/discretization/common/fvbaselocalresidualtpfa.hh b/opm/models/discretization/common/fvbaselocalresidualtpfa.hh deleted file mode 100644 index e1a363846..000000000 --- a/opm/models/discretization/common/fvbaselocalresidualtpfa.hh +++ /dev/null @@ -1,649 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/* - This file is part of the Open Porous Media project (OPM). - - OPM 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. - - OPM 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 OPM. If not, see . - - Consult the COPYING file in the top-level source directory of this - module for the precise wording of the license and the list of - copyright holders. -*/ -/*! - * \file - * - * \copydoc Opm::FvBaseLocalResidual - */ -#ifndef EWOMS_FV_BASE_LOCAL_TPFA_RESIDUAL_HH -#define EWOMS_FV_BASE_LOCAL_TPFA_RESIDUAL_HH - -#include "fvbaseproperties.hh" - -#include -#include - -#include -#include - -#include -#include - -#include - -#include - -#include - -namespace Opm { -/*! - * \ingroup FiniteVolumeDiscretizations - * - * \brief Element-wise caculation of the residual matrix for models based on a finite - * volume spatial discretization. - * - * \copydetails Doxygen::typeTagTParam - */ -template -class FvBaseLocalResidualTPFA -{ -private: - using Implementation = GetPropType; - - using GridView = GetPropType; - using Element = typename GridView::template Codim<0>::Entity; - - using Problem = GetPropType; - using Scalar = GetPropType; - using Evaluation = GetPropType; - using BoundaryRateVector = GetPropType; - using RateVector = GetPropType; - using EqVector = GetPropType; - using PrimaryVariables = GetPropType; - using ElementContext = GetPropType; - using BoundaryContext = GetPropType; - - static constexpr bool useVolumetricResidual = getPropValue(); - - enum { numEq = getPropValue() }; - enum { extensiveStorageTerm = getPropValue() }; - - using Toolbox = MathToolbox; - using EvalVector = Dune::FieldVector; - - // copying the local residual class is not a good idea - FvBaseLocalResidualTPFA(const FvBaseLocalResidualTPFA& ) - {} - -public: - using LocalEvalBlockVector = Dune::BlockVector >; - - FvBaseLocalResidualTPFA() - { } - - ~FvBaseLocalResidualTPFA() - { } - - /*! - * \brief Register all run-time parameters for the local residual. - */ - static void registerParameters() - { } - - /*! - * \brief Return the result of the eval() call using internal - * storage. - */ - const LocalEvalBlockVector& residual() const - { return internalResidual_; } - - /*! - * \brief Return the result of the eval() call using internal - * storage. - * - * \copydetails Doxygen::ecfvScvIdxParam - */ - const EvalVector& residual(unsigned dofIdx) const - { return internalResidual_[dofIdx]; } - - /*! - * \brief Compute the local residual, i.e. the deviation of the - * conservation equations from zero and store the results - * internally. - * - * The results can be requested afterwards using the residual() method. - * - * \copydetails Doxygen::problemParam - * \copydetails Doxygen::elementParam - */ - void eval(const Problem& problem, const Element& element) - { - ElementContext elemCtx(problem); - elemCtx.updateAll(element); - eval(elemCtx); - } - - /*! - * \brief Compute the local residual, i.e. the deviation of the - * conservation equations from zero and store the results - * internally. - * - * The results can be requested afterwards using the residual() method. - * - * \copydetails Doxygen::ecfvElemCtxParam - */ - void eval(ElementContext& elemCtx) - { - size_t numDof = elemCtx.numDof(/*timeIdx=*/0); - internalResidual_.resize(numDof); - asImp_().eval(internalResidual_, elemCtx); - } - - /*! - * \brief Compute the local residual, i.e. the deviation of the - * conservation equations from zero. - * - * \copydetails Doxygen::residualParam - * \copydetails Doxygen::ecfvElemCtxParam - */ - void eval(LocalEvalBlockVector& residual, - ElementContext& elemCtx) const - { - assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0)); - - residual = 0.0; - - // evaluate the flux terms - asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0); - - // evaluate the storage and the source terms - asImp_().evalVolumeTerms_(residual, elemCtx); - - // evaluate the boundary conditions - //asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0); - - if (useVolumetricResidual) { - // make the residual volume specific (i.e., make it incorrect mass per cubic - // meter instead of total mass) - size_t numDof = elemCtx.numDof(/*timeIdx=*/0); - for (unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) { - if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0.0) { - // interior DOF - Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0); - - assert(std::isfinite(dofVolume)); - Valgrind::CheckDefined(dofVolume); - - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - residual[dofIdx][eqIdx] /= dofVolume; - } - } - } - } - - /*! - * \brief Calculate the amount of all conservation quantities stored in all element's - * sub-control volumes for a given history index. - * - * This is used to figure out how much of each conservation quantity is inside the - * element. - * - * \copydetails Doxygen::storageParam - * \copydetails Doxygen::ecfvElemCtxParam - * \copydetails Doxygen::timeIdxParam - */ - void evalStorage(LocalEvalBlockVector& storage, - const ElementContext& elemCtx, - unsigned timeIdx) const - { - // the derivative of the storage term depends on the current primary variables; - // for time indices != 0, the storage term is constant (because these solutions - // are not changed by the Newton method!) - if (timeIdx == 0) { - // calculate the amount of conservation each quantity inside - // all primary sub control volumes - size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); - for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) { - storage[dofIdx] = 0.0; - - // the volume of the associated DOF - Scalar alpha = - elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume(); - //* elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor(); - - // If the degree of freedom which we currently look at is the one at the - // center of attention, we need to consider the derivatives for the - // storage term, else the storage term is constant w.r.t. the primary - // variables of the focused DOF. - if (dofIdx == elemCtx.focusDofIndex()) { - asImp_().computeStorage(storage[dofIdx], - elemCtx, - dofIdx, - timeIdx); - - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - storage[dofIdx][eqIdx] *= alpha; - } - else { - Dune::FieldVector tmp; - asImp_().computeStorage(tmp, - elemCtx, - dofIdx, - timeIdx); - - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - storage[dofIdx][eqIdx] = tmp[eqIdx]*alpha; - } - } - } - else { - // for all previous solutions, the storage term does _not_ depend on the - // current primary variables, so we use scalars to store it. - if (elemCtx.enableStorageCache()) { - size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx); - for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) { - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); - const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx); - for (unsigned eqIdx=0; eqIdx < numEq; eqIdx++) - storage[dofIdx][eqIdx] = cachedStorage[eqIdx]; - } - } - else { - // calculate the amount of conservation each quantity inside - // all primary sub control volumes - Dune::FieldVector tmp; - size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); - for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) { - tmp = 0.0; - asImp_().computeStorage(tmp, - elemCtx, - dofIdx, - timeIdx); - tmp *= - elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume() - * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor(); - - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - storage[dofIdx][eqIdx] = tmp[eqIdx]; - } - } - } - -#ifndef NDEBUG - size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); - for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) { - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - Valgrind::CheckDefined(storage[dofIdx][eqIdx]); - assert(isfinite(storage[dofIdx][eqIdx])); - } - } -#endif - } - - /*! - * \brief Add the flux term to a local residual. - * - * \copydetails Doxygen::residualParam - * \copydetails Doxygen::ecfvElemCtxParam - * \copydetails Doxygen::timeIdxParam - */ - void evalFluxes(LocalEvalBlockVector& residual, - const ElementContext& elemCtx, - unsigned timeIdx) const - { - RateVector flux; - - const auto& stencil = elemCtx.stencil(timeIdx); - // calculate the mass flux over the sub-control volume faces - size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx); - for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) { - const auto& face = stencil.interiorFace(scvfIdx); - unsigned i = face.interiorIndex(); - unsigned j = face.exteriorIndex(); - - Valgrind::SetUndefined(flux); - asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx); - Valgrind::CheckDefined(flux); -#ifndef NDEBUG - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - assert(isfinite(flux[eqIdx])); -#endif - - // Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor(); - //Scalar alpha = face.area(); - Scalar alpha = 1.0;//NB tpfa has - // Valgrind::CheckDefined(alpha); - // assert(alpha > 0.0); - // assert(isfinite(alpha)); - - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - flux[eqIdx] *= alpha; - - // The balance equation for a finite volume is given by - // - // dStorage/dt + Flux = Source - // - // where the 'Flux' and the 'Source' terms represent the - // mass per second which leaves the finite - // volume. Re-arranging this, we get - // - // dStorage/dt + Flux - Source = 0 - // - // Since the mass flux as calculated by computeFlux() goes out of sub-control - // volume i and into sub-control volume j, we need to add the flux to finite - // volume i and subtract it from finite volume j - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - assert(isfinite(flux[eqIdx])); - residual[i][eqIdx] += flux[eqIdx]; - residual[j][eqIdx] -= flux[eqIdx]; - } - } - -#if !defined NDEBUG - // in debug mode, ensure that the residual is well-defined - size_t numDof = elemCtx.numDof(timeIdx); - for (unsigned i=0; i < numDof; i++) { - for (unsigned j = 0; j < numEq; ++ j) { - assert(isfinite(residual[i][j])); - Valgrind::CheckDefined(residual[i][j]); - } - } -#endif - - } - - ///////////////////////////// - // The following methods _must_ be overloaded by the actual flow - // models! - ///////////////////////////// - - /*! - * \brief Evaluate the amount all conservation quantities - * (e.g. phase mass) within a finite sub-control volume. - * - * \copydetails Doxygen::storageParam - * \copydetails Doxygen::ecfvScvCtxParams - */ - void computeStorage(EqVector&, - const ElementContext&, - unsigned, - unsigned) const - { - throw std::logic_error("Not implemented: The local residual "+Dune::className() - +" does not implement the required method 'computeStorage()'"); - } - - /*! - * \brief Evaluates the total mass flux of all conservation - * quantities over a face of a sub-control volume. - * - * \copydetails Doxygen::areaFluxParam - * \copydetails Doxygen::ecfvScvfCtxParams - */ - void computeFlux(RateVector&, - const ElementContext&, - unsigned, - unsigned) const - { - throw std::logic_error("Not implemented: The local residual "+Dune::className() - +" does not implement the required method 'computeFlux()'"); - } - - /*! - * \brief Calculate the source term of the equation - * - * \copydoc Doxygen::sourceParam - * \copydoc Doxygen::ecfvScvCtxParams - */ - void computeSource(RateVector&, - const ElementContext&, - unsigned, - unsigned) const - { - throw std::logic_error("Not implemented: The local residual "+Dune::className() - +" does not implement the required method 'computeSource()'"); - } - -protected: - /*! - * \brief Evaluate the boundary conditions of an element. - */ - void evalBoundary_(LocalEvalBlockVector& residual, - const ElementContext& elemCtx, - unsigned timeIdx) const - { - - if (!elemCtx.onBoundary()) - return; - throw std::logic_error("Not implemented: Boundary??? "+Dune::className() - +" does not implement the required method 'computeSource()'"); - - BoundaryContext boundaryCtx(elemCtx); - // move the iterator to the first boundary - if(boundaryCtx.intersection(0).neighbor()) - boundaryCtx.increment(); - - // evaluate the boundary for all boundary faces of the current context - size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); - for (unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx, boundaryCtx.increment()) { - // add the residual of all vertices of the boundary - // segment - evalBoundarySegment_(residual, - boundaryCtx, - faceIdx, - timeIdx); - } - -#if !defined NDEBUG - // in debug mode, ensure that the residual and the storage terms are well-defined - size_t numDof = elemCtx.numDof(/*timeIdx=*/0); - for (unsigned i=0; i < numDof; i++) { - for (unsigned j = 0; j < numEq; ++ j) { - assert(isfinite(residual[i][j])); - Valgrind::CheckDefined(residual[i][j]); - } - } -#endif - - } - - /*! - * \brief Evaluate all boundary conditions for a single - * sub-control volume face to the local residual. - */ - void evalBoundarySegment_(LocalEvalBlockVector& residual, - const BoundaryContext& boundaryCtx, - unsigned boundaryFaceIdx, - unsigned timeIdx) const - { - throw std::logic_error("Not implemented: Boundary??? "+Dune::className() - +" does not implement the required method 'computeSource()'"); - - BoundaryRateVector values; - - Valgrind::SetUndefined(values); - boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx); - Valgrind::CheckDefined(values); - - const auto& stencil = boundaryCtx.stencil(timeIdx); - unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex(); - const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx); - for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx) { - values[eqIdx] *= - stencil.boundaryFace(boundaryFaceIdx).area() - * insideIntQuants.extrusionFactor(); - - Valgrind::CheckDefined(values[eqIdx]); - assert(isfinite(values[eqIdx])); - } - - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - residual[dofIdx][eqIdx] += values[eqIdx]; - } - - /*! - * \brief Add the change in the storage terms and the source term - * to the local residual of all sub-control volumes of the - * current element. - */ - void evalVolumeTerms_(LocalEvalBlockVector& residual, - ElementContext& elemCtx) const - { - EvalVector tmp; - EqVector tmp2; - RateVector sourceRate; - - tmp = 0.0; - tmp2 = 0.0; - - // evaluate the volumetric terms (storage + source terms) - size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); - for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) { - // Scalar extrusionFactor = - // elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor(); - // Valgrind::CheckDefined(extrusionFactor); - // assert(isfinite(extrusionFactor)); - // assert(extrusionFactor > 0.0); - Scalar scvVolume = - elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume();// * extrusionFactor; - Valgrind::CheckDefined(scvVolume); - assert(isfinite(scvVolume)); - assert(scvVolume > 0.0); - - // if the model uses extensive quantities in its storage term, and we use - // automatic differention and current DOF is also not the one we currently - // focus on, the storage term does not need any derivatives! - // if (!extensiveStorageTerm && - // !std::is_same::value && - // dofIdx != elemCtx.focusDofIndex()) - // { - // asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0); - // for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - // tmp[eqIdx] = tmp2[eqIdx]; - // } - // else - asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0); - -#ifndef NDEBUG - Valgrind::CheckDefined(tmp); - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - assert(isfinite(tmp[eqIdx])); -#endif - - if (elemCtx.enableStorageCache()) { - const auto& model = elemCtx.model(); - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - if (model.newtonMethod().numIterations() == 0 && - !elemCtx.haveStashedIntensiveQuantities()) - { - if (!elemCtx.problem().recycleFirstIterationStorage()) { - // we re-calculate the storage term for the solution of the - // previous time step from scratch instead of using the one of - // the first iteration of the current time step. - tmp2 = 0.0; - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1); - asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1); - } - else { - // if the storage term is cached and we're in the first iteration - // of the time step, use the storage term of the first iteration - // as the one as the solution of the last time step (this assumes - // that the initial guess for the solution at the end of the time - // step is the same as the solution at the beginning of the time - // step. This is usually true, but some fancy preprocessing - // scheme might invalidate that assumption.) - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); - } - - Valgrind::CheckDefined(tmp2); - - model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); - } - else { - // if the mass storage at the beginning of the time step is not cached, - // if the storage term is cached and we're not looking at the first - // iteration of the time step, we take the cached data. - tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1); - Valgrind::CheckDefined(tmp2); - } - } - else { - // if the mass storage at the beginning of the time step is not cached, - // we re-calculate it from scratch. - tmp2 = 0.0; - asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1); - Valgrind::CheckDefined(tmp2); - } - - // Use the implicit Euler time discretization - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - double dt = elemCtx.simulator().timeStepSize(); - assert(dt > 0); - tmp[eqIdx] -= tmp2[eqIdx]; - tmp[eqIdx] *= scvVolume / dt; - - residual[dofIdx][eqIdx] += tmp[eqIdx]; - } - - Valgrind::CheckDefined(residual[dofIdx]); - - // deal with the source term - asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0); - - // if the model uses extensive quantities in its storage term, and we use - // automatic differention and current DOF is also not the one we currently - // focus on, the storage term does not need any derivatives! - if (!extensiveStorageTerm && - !std::is_same::value && - dofIdx != elemCtx.focusDofIndex()) - { - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - residual[dofIdx][eqIdx] -= scalarValue(sourceRate[eqIdx])*scvVolume; - } - else { - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - sourceRate[eqIdx] *= scvVolume; - residual[dofIdx][eqIdx] -= sourceRate[eqIdx]; - } - } - - Valgrind::CheckDefined(residual[dofIdx]); - } - -#if !defined NDEBUG - // in debug mode, ensure that the residual is well-defined - size_t numDof = elemCtx.numDof(/*timeIdx=*/0); - for (unsigned i=0; i < numDof; i++) { - for (unsigned j = 0; j < numEq; ++ j) { - assert(isfinite(residual[i][j])); - Valgrind::CheckDefined(residual[i][j]); - } - } -#endif - } - - -private: - Implementation& asImp_() - { return *static_cast(this); } - - const Implementation& asImp_() const - { return *static_cast(this); } - - LocalEvalBlockVector internalResidual_; -}; - -} // namespace Opm - -#endif