diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 75f1a7ce4..7a92cd853 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -57,13 +57,14 @@ class BlackOilLocalResidualTPFA : public GetPropType; using RateVector = GetPropType; using FluidSystem = GetPropType; - + using GridView = GetPropType; enum { conti0EqIdx = Indices::conti0EqIdx }; enum { numEq = getPropValue() }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; - + + enum { dimWorld = GridView::dimensionworld }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; @@ -195,27 +196,117 @@ public: assert(timeIdx == 0); flux = 0.0; + // need for dary flux calculation + const auto& problem = elemCtx.problem(); + const auto& stencil = elemCtx.stencil(timeIdx); + const auto& scvf = stencil.interiorFace(scvfIdx); - const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + unsigned interiorDofIdx = scvf.interiorIndex(); + unsigned exteriorDofIdx = scvf.exteriorIndex(); + assert(interiorDofIdx != exteriorDofIdx); + + //unsigned I = stencil.globalSpaceIndex(interiorDofIdx); + //unsigned J = stencil.globalSpaceIndex(exteriorDofIdx); + Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0); + Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0); + const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx); + const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx); + Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); + Scalar faceArea = scvf.area(); + Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx); + + // estimate the gravity correction: for performance reasons we use a simplified + // approach for this flux module that assumes that gravity is constant and always + // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) + Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; + + const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx); + const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx); + + // this is quite hacky because the dune grid interface does not provide a + // cellCenterDepth() method (so we ask the problem to provide it). The "good" + // solution would be to take the Z coordinate of the element centroids, but since + // ECL seems to like to be inconsistent on that front, it needs to be done like + // here... + Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx); + Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx); + + // the distances from the DOF's depths. (i.e., the additional depth of the + // exterior DOF) + Scalar distZ = zIn - zEx; + + // + //const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); unsigned focusDofIdx = elemCtx.focusDofIndex(); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) continue; + // darcy flux calculation + short dnIdx; + // + short upIdx; + Evaluation pressureDifference; + ExtensiveQuantities::calculatePhasePressureDiff_(upIdx, + dnIdx, + pressureDifference, + intQuantsIn, + intQuantsEx, + scvfIdx,//input + timeIdx,//input + phaseIdx,//input + interiorDofIdx,//input + exteriorDofIdx,//intput + Vin, + Vex, + globalIndexIn, + globalIndexEx, + distZ*g, + thpres); + + + + const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx; + unsigned globalIndex; + if(upIdx == interiorDofIdx){ + //up = intQuantsIn; + globalIndex = globalIndexIn; + }else{ + //up = intQuantsEx; + globalIndex = globalIndexEx; + } + // TODO: should the rock compaction transmissibility multiplier be upstreamed + // or averaged? all fluids should see the same compaction?! + //const auto& globalIndex = stencil.globalSpaceIndex(upstreamIdx); + const Evaluation& transMult = + problem.template rockCompTransMultiplier(up, globalIndex); + Evaluation darcyFlux; + if(pressureDifference == 0){ + darcyFlux = 0.0; //NB maybe we could drop calculations + }else{ + if (upIdx == interiorDofIdx) + darcyFlux = + pressureDifference*up.mobility(phaseIdx)*transMult*(-trans/faceArea); + else + darcyFlux = + pressureDifference*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); + } + //const auto& darcyFlux = extQuants.volumeFlux(phaseIdx); + //unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); + + //const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + unsigned pvtRegionIdx = up.pvtRegionIndex(); + using FluidState = typename IntensiveQuantities::FluidState; + if (upIdx == focusDofIdx){ + const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); + const auto& surfaceVolumeFlux = invB*darcyFlux; + evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); + }else{ + const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); + const auto& surfaceVolumeFlux = invB*darcyFlux; + evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); + } + - const auto& darcyFlux = extQuants.volumeFlux(phaseIdx); - unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); - const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx); - unsigned pvtRegionIdx = up.pvtRegionIndex(); - using FluidState = typename IntensiveQuantities::FluidState; - if (upIdx == focusDofIdx){ - const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); - const auto& surfaceVolumeFlux = invB*darcyFlux; - evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); - }else{ - const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); - const auto& surfaceVolumeFlux = invB*darcyFlux; - evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); - } } // deal with solvents (if present) diff --git a/opm/models/discretization/common/fvbaseadlocallinearizertpfa.hh b/opm/models/discretization/common/fvbaseadlocallinearizertpfa.hh new file mode 100644 index 000000000..990b157c7 --- /dev/null +++ b/opm/models/discretization/common/fvbaseadlocallinearizertpfa.hh @@ -0,0 +1,308 @@ +// -*- 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::FvBaseAdLocalLinearizer + */ +#ifndef EWOMS_FV_BASE_AD_LOCAL_TPFA_LINEARIZER_HH +#define EWOMS_FV_BASE_AD_LOCAL_TPFA_LINEARIZER_HH + +#include "fvbaseproperties.hh" + +#include +#include +#include + +#include +#include + +#include +#include + +namespace Opm { +//forward declaration +template +class FvBaseAdLocalLinearizerTPFA; +} + +namespace Opm::Properties { + +//declare the property tags required for the finite differences local linearizer + +namespace TTag { +struct AutoDiffLocalLinearizerTPFA {}; +} // namespace TTag + +// set the properties to be spliced in +template +struct LocalLinearizer +{ using type = FvBaseAdLocalLinearizerTPFA; }; + +//! Set the function evaluation w.r.t. the primary variables +template +struct Evaluation +{ +private: + static const unsigned numEq = getPropValue(); + + using Scalar = GetPropType; + +public: + using type = DenseAd::Evaluation; +}; + +} // namespace Opm::Properties + +namespace Opm { + +/*! + * \ingroup FiniteVolumeDiscretizations + * + * \brief Calculates the local residual and its Jacobian for a single element of the grid. + * + * This class uses automatic differentiation to calculate the partial derivatives (the + * alternative is finite differences). + */ +template +class FvBaseAdLocalLinearizerTPFA +{ +private: + using Implementation = GetPropType; + using LocalResidual = GetPropType; + using Simulator = GetPropType; + using Problem = GetPropType; + using Model = GetPropType; + using PrimaryVariables = GetPropType; + using ElementContext = GetPropType; + using Scalar = GetPropType; + using GridView = GetPropType; + using Element = typename GridView::template Codim<0>::Entity; + + enum { numEq = getPropValue() }; + + using ScalarVectorBlock = Dune::FieldVector; + // extract local matrices from jacobian matrix for consistency + using ScalarMatrixBlock = typename GetPropType::MatrixBlock; + + using ScalarLocalBlockVector = Dune::BlockVector; + using ScalarLocalBlockMatrix = Dune::Matrix; + +public: + FvBaseAdLocalLinearizerTPFA() + : internalElemContext_(0) + { } + + // copying local linearizer objects around is a very bad idea, so we explicitly + // prevent it... + FvBaseAdLocalLinearizerTPFA(const FvBaseAdLocalLinearizerTPFA&) = delete; + + ~FvBaseAdLocalLinearizerTPFA() + { delete internalElemContext_; } + + /*! + * \brief Register all run-time parameters for the local jacobian. + */ + static void registerParameters() + { } + + /*! + * \brief Initialize the local Jacobian object. + * + * At this point we can assume that everything has been allocated, + * although some objects may not yet be completely initialized. + * + * \param simulator The simulator object of the simulation. + */ + void init(Simulator& simulator) + { + simulatorPtr_ = &simulator; + delete internalElemContext_; + internalElemContext_ = new ElementContext(simulator); + } + + /*! + * \brief Compute an element's local Jacobian matrix and evaluate its residual. + * + * The local Jacobian for a given context is defined as the derivatives of the + * residuals of all degrees of freedom featured by the stencil with regard to the + * primary variables of the stencil's "primary" degrees of freedom. Adding the local + * Jacobians for all elements in the grid will give the global Jacobian 'grad f(x)'. + * + * \param element The grid element for which the local residual and its local + * Jacobian should be calculated. + */ + void linearize(const Element& element) + { + linearize(*internalElemContext_, element); + } + + /*! + * \brief Compute an element's local Jacobian matrix and evaluate its residual. + * + * The local Jacobian for a given context is defined as the derivatives of the + * residuals of all degrees of freedom featured by the stencil with regard to the + * primary variables of the stencil's "primary" degrees of freedom. Adding the local + * Jacobians for all elements in the grid will give the global Jacobian 'grad f(x)'. + * + * After calling this method the ElementContext is in an undefined state, so do not + * use it anymore! + * + * \param elemCtx The element execution context for which the local residual and its + * local Jacobian should be calculated. + */ + void linearize(ElementContext& elemCtx, const Element& elem) + { + elemCtx.updateStencil(elem); + elemCtx.updateAllIntensiveQuantities(); + + // update the weights of the primary variables for the context + model_().updatePVWeights(elemCtx); + + // resize the internal arrays of the linearizer + resize_(elemCtx); + + // compute the local residual and its Jacobian + unsigned numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); + for (unsigned focusDofIdx = 0; focusDofIdx < numPrimaryDof; focusDofIdx++) { + elemCtx.setFocusDofIndex(focusDofIdx); + //elemCtx.updateAllExtensiveQuantities();//NB should not be need anymore + + // calculate the local residual + localResidual_.eval(elemCtx); + + // convert the local Jacobian matrix and the right hand side from the data + // structures used by the automatic differentiation code to the conventional + // ones used by the linear solver. + updateLocalLinearization_(elemCtx, focusDofIdx); + } + } + + /*! + * \brief Return reference to the local residual. + */ + LocalResidual& localResidual() + { return localResidual_; } + + /*! + * \brief Return reference to the local residual. + */ + const LocalResidual& localResidual() const + { return localResidual_; } + + /*! + * \brief Returns the local Jacobian matrix of the residual of a sub-control volume. + * + * \param domainScvIdx The local index of the sub control volume to which the primary + * variables are associated with + * \param rangeScvIdx The local index of the sub control volume which contains the + * local residual + */ + const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const + { return jacobian_[domainScvIdx][rangeScvIdx]; } + + /*! + * \brief Returns the local residual of a sub-control volume. + * + * \param dofIdx The local index of the sub control volume + */ + const ScalarVectorBlock& residual(unsigned dofIdx) const + { return residual_[dofIdx]; } + +protected: + Implementation& asImp_() + { return *static_cast(this); } + const Implementation& asImp_() const + { return *static_cast(this); } + + const Simulator& simulator_() const + { return *simulatorPtr_; } + const Problem& problem_() const + { return simulatorPtr_->problem(); } + const Model& model_() const + { return simulatorPtr_->model(); } + + /*! + * \brief Resize all internal attributes to the size of the + * element. + */ + void resize_(const ElementContext& elemCtx) + { + size_t numDof = elemCtx.numDof(/*timeIdx=*/0); + size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0); + + residual_.resize(numDof); + if (jacobian_.N() != numDof || jacobian_.M() != numPrimaryDof) + jacobian_.setSize(numDof, numPrimaryDof); + } + + /*! + * \brief Reset the all relevant internal attributes to 0 + */ + void reset_(const ElementContext&) + { + residual_ = 0.0; + jacobian_ = 0.0; + } + + /*! + * \brief Updates the current local Jacobian matrix with the partial derivatives of + * all equations for the degree of freedom associated with 'focusDofIdx'. + */ + void updateLocalLinearization_(const ElementContext& elemCtx, + unsigned focusDofIdx) + { + const auto& resid = localResidual_.residual(); + + for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) + residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value(); + + size_t numDof = elemCtx.numDof(/*timeIdx=*/0); + for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) { + for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) { + for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) { + // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of + // the residual function 'eqIdx' for the degree of freedom 'dofIdx' + // with regard to the focus variable 'pvIdx' of the degree of freedom + // 'focusDofIdx' + jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx); + Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]); + } + } + } + } + + Simulator *simulatorPtr_; + Model *modelPtr_; + + ElementContext *internalElemContext_; + + LocalResidual localResidual_; + + ScalarLocalBlockVector residual_; + ScalarLocalBlockMatrix jacobian_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/discretization/common/fvbaselocalresidualtpfa.hh b/opm/models/discretization/common/fvbaselocalresidualtpfa.hh new file mode 100644 index 000000000..d8abbf5ca --- /dev/null +++ b/opm/models/discretization/common/fvbaselocalresidualtpfa.hh @@ -0,0 +1,648 @@ +// -*- 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(); + // 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