| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Assembly | ||
| 10 | * \ingroup CCDiscretization | ||
| 11 | * \ingroup MultiDomain | ||
| 12 | * \brief A multidomain local assembler for Jacobian and residual contribution per element (cell-centered methods) | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_MULTIDOMAIN_CC_LOCAL_ASSEMBLER_HH | ||
| 15 | #define DUMUX_MULTIDOMAIN_CC_LOCAL_ASSEMBLER_HH | ||
| 16 | |||
| 17 | #include <dune/common/indices.hh> | ||
| 18 | #include <dune/common/reservedvector.hh> | ||
| 19 | #include <dune/common/hybridutilities.hh> | ||
| 20 | #include <dune/grid/common/gridenums.hh> // for GhostEntity | ||
| 21 | #include <dune/istl/matrixindexset.hh> | ||
| 22 | |||
| 23 | #include <dumux/common/reservedblockvector.hh> | ||
| 24 | #include <dumux/common/properties.hh> | ||
| 25 | #include <dumux/common/parameters.hh> | ||
| 26 | #include <dumux/common/numericdifferentiation.hh> | ||
| 27 | #include <dumux/common/numeqvector.hh> | ||
| 28 | #include <dumux/discretization/elementsolution.hh> | ||
| 29 | #include <dumux/discretization/extrusion.hh> | ||
| 30 | #include <dumux/assembly/numericepsilon.hh> | ||
| 31 | #include <dumux/assembly/diffmethod.hh> | ||
| 32 | #include <dumux/assembly/fvlocalassemblerbase.hh> | ||
| 33 | |||
| 34 | namespace Dumux { | ||
| 35 | |||
| 36 | /*! | ||
| 37 | * \ingroup Assembly | ||
| 38 | * \ingroup CCDiscretization | ||
| 39 | * \ingroup MultiDomain | ||
| 40 | * \brief A base class for all multidomain local assemblers | ||
| 41 | * \tparam id the id of the sub domain | ||
| 42 | * \tparam TypeTag the TypeTag | ||
| 43 | * \tparam Assembler the assembler type | ||
| 44 | * \tparam Implementation the actual assembler implementation | ||
| 45 | * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit) | ||
| 46 | */ | ||
| 47 | template<std::size_t id, class TypeTag, class Assembler, class Implementation, bool implicit> | ||
| 48 | 53291655 | class SubDomainCCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler,Implementation, implicit> | |
| 49 | { | ||
| 50 | using ParentType = FVLocalAssemblerBase<TypeTag, Assembler,Implementation, implicit>; | ||
| 51 | |||
| 52 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 53 | using LocalResidualValues = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 54 | using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>; | ||
| 55 | using SolutionVector = typename Assembler::SolutionVector; | ||
| 56 | using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>; | ||
| 57 | |||
| 58 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 59 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 60 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 61 | using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView; | ||
| 62 | using Scalar = typename GridVariables::Scalar; | ||
| 63 | |||
| 64 | using GridGeometry = typename GridVariables::GridGeometry; | ||
| 65 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 66 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 67 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 68 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 69 | using GridView = typename GridGeometry::GridView; | ||
| 70 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 71 | |||
| 72 | using CouplingManager = typename Assembler::CouplingManager; | ||
| 73 | using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector; | ||
| 74 | |||
| 75 | public: | ||
| 76 | //! export the domain id of this sub-domain | ||
| 77 | static constexpr auto domainId = typename Dune::index_constant<id>(); | ||
| 78 | //! the local residual type of this domain | ||
| 79 | using LocalResidual = typename ParentType::LocalResidual; | ||
| 80 | //! pull up constructor of parent class | ||
| 81 | using ParentType::ParentType; | ||
| 82 | |||
| 83 | //! the constructor | ||
| 84 | 56919862 | explicit SubDomainCCLocalAssemblerBase(const Assembler& assembler, | |
| 85 | const Element& element, | ||
| 86 | const SolutionVector& curSol, | ||
| 87 | CouplingManager& couplingManager) | ||
| 88 | : ParentType(assembler, | ||
| 89 | element, | ||
| 90 | curSol, | ||
| 91 | localView(assembler.gridGeometry(domainId)), | ||
| 92 |
1/2✓ Branch 1 taken 30056867 times.
✗ Branch 2 not taken.
|
56919862 | localView(assembler.gridVariables(domainId).curGridVolVars()), |
| 93 |
1/2✓ Branch 1 taken 30056867 times.
✗ Branch 2 not taken.
|
56919862 | localView(assembler.gridVariables(domainId).prevGridVolVars()), |
| 94 |
1/2✓ Branch 1 taken 30056867 times.
✗ Branch 2 not taken.
|
56919862 | localView(assembler.gridVariables(domainId).gridFluxVarsCache()), |
| 95 |
1/2✓ Branch 1 taken 30056867 times.
✗ Branch 2 not taken.
|
56919862 | assembler.localResidual(domainId), |
| 96 |
1/2✓ Branch 1 taken 3242051 times.
✗ Branch 2 not taken.
|
56919862 | (element.partitionType() == Dune::GhostEntity)) |
| 97 |
3/5✗ Branch 0 not taken.
✓ Branch 1 taken 3719133 times.
✓ Branch 3 taken 477082 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 26337734 times.
|
62078476 | , couplingManager_(couplingManager) |
| 98 | 56919862 | {} | |
| 99 | |||
| 100 | /*! | ||
| 101 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 102 | * to the global matrix. The element residual is written into the right hand side. | ||
| 103 | */ | ||
| 104 | template<class JacobianMatrixRow, class SubResidualVector, class GridVariablesTuple> | ||
| 105 | 32166307 | void assembleJacobianAndResidual(JacobianMatrixRow& jacRow, SubResidualVector& res, GridVariablesTuple& gridVariables) | |
| 106 | { | ||
| 107 | 32166307 | this->asImp_().bindLocalViews(); | |
| 108 | |||
| 109 | // for the diagonal jacobian block | ||
| 110 | 32166307 | const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element()); | |
| 111 |
1/2✓ Branch 2 taken 5781 times.
✗ Branch 3 not taken.
|
32166307 | res[globalI] = this->asImp_().assembleJacobianAndResidualImplInverse(jacRow[domainId], *std::get<domainId>(gridVariables)); |
| 112 | |||
| 113 | // for the coupling blocks | ||
| 114 | using namespace Dune::Hybrid; | ||
| 115 |
1/2✓ Branch 1 taken 5781 times.
✗ Branch 2 not taken.
|
63267983 | forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i) |
| 116 | { | ||
| 117 | if (i != id) | ||
| 118 | 20508553 | this->assembleJacobianCoupling(i, jacRow, res[globalI], gridVariables); | |
| 119 | }); | ||
| 120 | 32166307 | } | |
| 121 | |||
| 122 | /*! | ||
| 123 | * \brief Assemble the entries in a coupling block of the jacobian. | ||
| 124 | * There is no coupling block between a domain and itself. | ||
| 125 | */ | ||
| 126 | template<std::size_t otherId, class JacRow, class GridVariables, | ||
| 127 | typename std::enable_if_t<(otherId == id), int> = 0> | ||
| 128 | void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow, | ||
| 129 | const LocalResidualValues& res, GridVariables& gridVariables) | ||
| 130 | {} | ||
| 131 | |||
| 132 | /*! | ||
| 133 | * \brief Assemble the entries in a coupling block of the jacobian. | ||
| 134 | */ | ||
| 135 | template<std::size_t otherId, class JacRow, class GridVariables, | ||
| 136 | typename std::enable_if_t<(otherId != id), int> = 0> | ||
| 137 | 16452553 | void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow, | |
| 138 | const LocalResidualValues& res, GridVariables& gridVariables) | ||
| 139 | { | ||
| 140 | 16452553 | this->asImp_().assembleJacobianCoupling(domainJ, jacRow[domainJ], res, *std::get<domainId>(gridVariables)); | |
| 141 | } | ||
| 142 | |||
| 143 | /*! | ||
| 144 | * \brief Assemble the residual only | ||
| 145 | */ | ||
| 146 | template<class SubResidualVector> | ||
| 147 | 24753555 | void assembleResidual(SubResidualVector& res) | |
| 148 | { | ||
| 149 | 24753555 | this->asImp_().bindLocalViews(); | |
| 150 | 24753555 | const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element()); | |
| 151 | 24753555 | res[globalI] = this->evalLocalResidual()[0]; // forward to the internal implementation | |
| 152 | |||
| 153 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 154 | { | ||
| 155 | 705110 | const auto applyDirichlet = [&] (const auto& scvI, | |
| 156 | const auto& dirichletValues, | ||
| 157 | const auto eqIdx, | ||
| 158 | const auto pvIdx) | ||
| 159 | { | ||
| 160 | 110 | res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx]; | |
| 161 | }; | ||
| 162 | |||
| 163 | 705000 | this->asImp_().enforceInternalDirichletConstraints(applyDirichlet); | |
| 164 | } | ||
| 165 | 24753553 | } | |
| 166 | |||
| 167 | /*! | ||
| 168 | * \brief Evaluates the local source term for an element and given element volume variables | ||
| 169 | */ | ||
| 170 | 2859586 | ElementResidualVector evalLocalSourceResidual(const Element& element, const ElementVolumeVariables& elemVolVars) const | |
| 171 | { | ||
| 172 | // initialize the residual vector for all scvs in this element | ||
| 173 | 2859586 | ElementResidualVector residual(this->fvGeometry().numScv()); | |
| 174 | |||
| 175 | // evaluate the volume terms (storage + source terms) | ||
| 176 | // forward to the local residual specialized for the discretization methods | ||
| 177 |
2/2✓ Branch 0 taken 1430633 times.
✓ Branch 1 taken 1430633 times.
|
5719172 | for (auto&& scv : scvs(this->fvGeometry())) |
| 178 | { | ||
| 179 | 2859586 | const auto& curVolVars = elemVolVars[scv]; | |
| 180 | 5575896 | auto source = this->localResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv); | |
| 181 | 2859586 | source *= -Extrusion::volume(this->fvGeometry(), scv)*curVolVars.extrusionFactor(); | |
| 182 | 2859586 | residual[scv.indexInElement()] = std::move(source); | |
| 183 | } | ||
| 184 | |||
| 185 | 2859586 | return residual; | |
| 186 | } | ||
| 187 | |||
| 188 | /*! | ||
| 189 | * \brief Evaluates the local source term depending on time discretization scheme | ||
| 190 | */ | ||
| 191 | 1430633 | ElementResidualVector evalLocalSourceResidual(const Element& neighbor) const | |
| 192 | 1430633 | { return this->evalLocalSourceResidual(neighbor, implicit ? this->curElemVolVars() : this->prevElemVolVars()); } | |
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Evaluates the storage terms within the element | ||
| 196 | */ | ||
| 197 | 4056000 | LocalResidualValues evalLocalStorageResidual() const | |
| 198 | { | ||
| 199 | 8112000 | return this->localResidual().evalStorage(this->element(), this->fvGeometry(), this->prevElemVolVars(), this->curElemVolVars())[0]; | |
| 200 | } | ||
| 201 | |||
| 202 | /*! | ||
| 203 | * \brief Evaluates the fluxes depending on the chose time discretization scheme | ||
| 204 | */ | ||
| 205 | 169984356 | LocalResidualValues evalFluxResidual(const Element& neighbor, | |
| 206 | const SubControlVolumeFace& scvf) const | ||
| 207 | { | ||
| 208 | 169984356 | const auto& elemVolVars = implicit ? this->curElemVolVars() : this->prevElemVolVars(); | |
| 209 | 169984356 | return this->localResidual().evalFlux(problem(), neighbor, this->fvGeometry(), elemVolVars, this->elemFluxVarsCache(), scvf); | |
| 210 | } | ||
| 211 | |||
| 212 | /*! | ||
| 213 | * \brief Prepares all local views necessary for local assembly. | ||
| 214 | */ | ||
| 215 | 56919862 | void bindLocalViews() | |
| 216 | { | ||
| 217 | // get some references for convenience | ||
| 218 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | const auto& element = this->element(); |
| 219 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | const auto& curSol = this->curSol()[domainId]; |
| 220 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | auto&& fvGeometry = this->fvGeometry(); |
| 221 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | auto&& curElemVolVars = this->curElemVolVars(); |
| 222 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | auto&& elemFluxVarsCache = this->elemFluxVarsCache(); |
| 223 | |||
| 224 | // bind the caches | ||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 173500 times.
|
56919862 | couplingManager_.bindCouplingContext(domainId, element, this->assembler()); |
| 226 | 56919862 | fvGeometry.bind(element); | |
| 227 | |||
| 228 | if (implicit) | ||
| 229 | { | ||
| 230 | 52863862 | curElemVolVars.bind(element, fvGeometry, curSol); | |
| 231 | 2674414 | elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars); | |
| 232 |
2/2✓ Branch 0 taken 721180 times.
✓ Branch 1 taken 828455 times.
|
52863862 | if (!this->assembler().isStationaryProblem()) |
| 233 | 1258140 | this->prevElemVolVars().bindElement(element, fvGeometry, this->assembler().prevSol()[domainId]); | |
| 234 | } | ||
| 235 | else | ||
| 236 | { | ||
| 237 | 4056000 | auto& prevElemVolVars = this->prevElemVolVars(); | |
| 238 | 4056000 | const auto& prevSol = this->assembler().prevSol()[domainId]; | |
| 239 | |||
| 240 | 4056000 | curElemVolVars.bindElement(element, fvGeometry, curSol); | |
| 241 | 4056000 | prevElemVolVars.bind(element, fvGeometry, prevSol); | |
| 242 | 4056000 | elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars); | |
| 243 | } | ||
| 244 | 56919862 | } | |
| 245 | |||
| 246 | //! return reference to the underlying problem | ||
| 247 | 197594759 | const Problem& problem() const | |
| 248 |
10/16✓ Branch 0 taken 706680 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 75299 times.
✓ Branch 3 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✓ Branch 11 taken 12 times.
✓ Branch 12 taken 1 times.
✓ Branch 14 taken 12 times.
✗ Branch 15 not taken.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✓ Branch 10 taken 7 times.
✓ Branch 13 taken 7 times.
|
172120141 | { return this->assembler().problem(domainId); } |
| 249 | |||
| 250 | //! return reference to the coupling manager | ||
| 251 | 103221183 | CouplingManager& couplingManager() | |
| 252 | 132133900 | { return couplingManager_; } | |
| 253 | |||
| 254 | private: | ||
| 255 | CouplingManager& couplingManager_; //!< the coupling manager | ||
| 256 | }; | ||
| 257 | |||
| 258 | /*! | ||
| 259 | * \ingroup Assembly | ||
| 260 | * \ingroup CCDiscretization | ||
| 261 | * \ingroup MultiDomain | ||
| 262 | * \brief The cell-centered scheme multidomain local assembler | ||
| 263 | * \tparam id the id of the sub domain | ||
| 264 | * \tparam TypeTag the TypeTag | ||
| 265 | * \tparam Assembler the assembler type | ||
| 266 | * \tparam DM the numeric differentiation method | ||
| 267 | * \tparam implicit whether the assembler is explicit or implicit in time | ||
| 268 | */ | ||
| 269 | template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true> | ||
| 270 | class SubDomainCCLocalAssembler; | ||
| 271 | |||
| 272 | /*! | ||
| 273 | * \ingroup Assembly | ||
| 274 | * \ingroup CCDiscretization | ||
| 275 | * \ingroup MultiDomain | ||
| 276 | * \brief Cell-centered scheme multidomain local assembler using numeric differentiation and implicit time discretization | ||
| 277 | */ | ||
| 278 | template<std::size_t id, class TypeTag, class Assembler> | ||
| 279 | 51263654 | class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true> | |
| 280 | : public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, | ||
| 281 | SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true>, true> | ||
| 282 | { | ||
| 283 | using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>; | ||
| 284 | using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/true>; | ||
| 285 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 286 | |||
| 287 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 288 | using LocalResidualValues = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 289 | |||
| 290 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 291 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 292 | using GridView = typename GridGeometry::GridView; | ||
| 293 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 294 | |||
| 295 | enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() }; | ||
| 296 | enum { dim = GridView::dimension }; | ||
| 297 | |||
| 298 | static constexpr bool enableGridFluxVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridFluxVariablesCache::cachingEnabled; | ||
| 299 | static constexpr bool enableGridVolVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridVolumeVariables::cachingEnabled; | ||
| 300 | static constexpr int maxElementStencilSize = GridGeometry::maxElementStencilSize; | ||
| 301 | static constexpr auto domainI = Dune::index_constant<id>(); | ||
| 302 | |||
| 303 | public: | ||
| 304 | 28028867 | using ParentType::ParentType; | |
| 305 | |||
| 306 | /*! | ||
| 307 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 308 | * to the global matrix. | ||
| 309 | * | ||
| 310 | * \return The element residual at the current solution. | ||
| 311 | */ | ||
| 312 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 313 | 28110307 | LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | |
| 314 | { | ||
| 315 | ////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 316 | // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the // | ||
| 317 | // neighboring elements we do so by computing the derivatives of the fluxes which depend on the // | ||
| 318 | // actual element. In the actual element we evaluate the derivative of the entire residual. // | ||
| 319 | ////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 320 | |||
| 321 | // get some aliases for convenience | ||
| 322 | 28110307 | const auto& element = this->element(); | |
| 323 | 28110307 | const auto& fvGeometry = this->fvGeometry(); | |
| 324 | 28110307 | const auto& gridGeometry = fvGeometry.gridGeometry(); | |
| 325 | 28110307 | auto&& curElemVolVars = this->curElemVolVars(); | |
| 326 | 28110307 | auto&& elemFluxVarsCache = this->elemFluxVarsCache(); | |
| 327 | |||
| 328 | // get stencil information | ||
| 329 | 28110307 | const auto globalI = gridGeometry.elementMapper().index(element); | |
| 330 | 28110307 | const auto& connectivityMap = gridGeometry.connectivityMap(); | |
| 331 | 28110307 | const auto numNeighbors = connectivityMap[globalI].size(); | |
| 332 | |||
| 333 | // container to store the neighboring elements | ||
| 334 | 28110307 | Dune::ReservedVector<Element, maxElementStencilSize> neighborElements; | |
| 335 | 28110307 | neighborElements.resize(numNeighbors); | |
| 336 | |||
| 337 | // assemble the undeflected residual | ||
| 338 | using Residuals = ReservedBlockVector<LocalResidualValues, maxElementStencilSize>; | ||
| 339 | 28110467 | Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0; | |
| 340 |
1/2✓ Branch 1 taken 572897 times.
✗ Branch 2 not taken.
|
28110307 | origResiduals[0] = this->evalLocalResidual()[0]; |
| 341 | |||
| 342 | // lambda for convenient evaluation of the fluxes across scvfs in the neighbors | ||
| 343 | // if the neighbor is a ghost we don't want to add anything to their residual | ||
| 344 | // so we return 0 and omit computing the flux | ||
| 345 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1131388 times.
|
195791726 | const auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf) |
| 346 | { | ||
| 347 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 156344436 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3991788 times.
|
163472524 | if (neighbor.partitionType() == Dune::GhostEntity) |
| 348 | ✗ | return LocalResidualValues(0.0); | |
| 349 | else | ||
| 350 | 169984356 | return this->evalFluxResidual(neighbor, scvf); | |
| 351 | }; | ||
| 352 | |||
| 353 | // get the elements in which we need to evaluate the fluxes | ||
| 354 | // and calculate these in the undeflected state | ||
| 355 | 28110307 | unsigned int j = 1; | |
| 356 |
2/2✓ Branch 0 taken 75347642 times.
✓ Branch 1 taken 15206622 times.
|
170014751 | for (const auto& dataJ : connectivityMap[globalI]) |
| 357 | { | ||
| 358 |
1/2✓ Branch 1 taken 3153050 times.
✗ Branch 2 not taken.
|
148183280 | neighborElements[j-1] = gridGeometry.element(dataJ.globalJ); |
| 359 |
5/5✓ Branch 1 taken 75347642 times.
✓ Branch 2 taken 68391114 times.
✓ Branch 3 taken 13632 times.
✓ Branch 4 taken 13632 times.
✓ Branch 0 taken 10670852 times.
|
291264800 | for (const auto scvfIdx : dataJ.scvfsJ) |
| 360 |
2/3✓ Branch 1 taken 13632 times.
✓ Branch 2 taken 6867374 times.
✗ Branch 3 not taken.
|
163005470 | origResiduals[j] += evalNeighborFlux(neighborElements[j-1], fvGeometry.scvf(scvfIdx)); |
| 361 | |||
| 362 | 141904444 | ++j; | |
| 363 | } | ||
| 364 | |||
| 365 | // reference to the element's scv (needed later) and corresponding vol vars | ||
| 366 |
1/2✓ Branch 0 taken 1549635 times.
✗ Branch 1 not taken.
|
28110307 | const auto& scv = fvGeometry.scv(globalI); |
| 367 | 28110307 | auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv); | |
| 368 | |||
| 369 | // save a copy of the original privars and vol vars in order | ||
| 370 | // to restore the original solution after deflection | ||
| 371 | 28110307 | const auto& curSol = this->curSol()[domainI]; | |
| 372 | 28110307 | const auto origPriVars = curSol[globalI]; | |
| 373 | 28110307 | const auto origVolVars = curVolVars; | |
| 374 | |||
| 375 | // element solution container to be deflected | ||
| 376 | 28110307 | auto elemSol = elementSolution<FVElementGeometry>(origPriVars); | |
| 377 | |||
| 378 | // derivatives in the neighbors with respect to the current elements | ||
| 379 | // in index 0 we save the derivative of the element residual with respect to it's own dofs | ||
| 380 | 28110307 | Residuals partialDerivs(numNeighbors + 1); | |
| 381 | |||
| 382 |
2/2✓ Branch 0 taken 16760666 times.
✓ Branch 1 taken 15206622 times.
|
58589126 | for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) |
| 383 | { | ||
| 384 | 30478819 | partialDerivs = 0.0; | |
| 385 | |||
| 386 | 47263787 | auto evalResiduals = [&](Scalar priVar) | |
| 387 | { | ||
| 388 | 16784968 | Residuals partialDerivsTmp(numNeighbors + 1); | |
| 389 | 16784968 | partialDerivsTmp = 0.0; | |
| 390 | // update the volume variables and the flux var cache | ||
| 391 | 16784968 | elemSol[0][pvIdx] = priVar; | |
| 392 | 16784968 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx); | |
| 393 | 16784968 | curVolVars.update(elemSol, this->problem(), element, scv); | |
| 394 | 16398568 | elemFluxVarsCache.update(element, fvGeometry, curElemVolVars); | |
| 395 | if (enableGridFluxVarsCache) | ||
| 396 | 14373813 | gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars); | |
| 397 | |||
| 398 | // calculate the residual with the deflected primary variables | ||
| 399 | 16784968 | partialDerivsTmp[0] = this->evalLocalResidual()[0]; | |
| 400 | |||
| 401 | // calculate the fluxes in the neighbors with the deflected primary variables | ||
| 402 |
6/6✓ Branch 0 taken 16772134 times.
✓ Branch 1 taken 5071665 times.
✓ Branch 2 taken 65842968 times.
✓ Branch 3 taken 11707522 times.
✓ Branch 4 taken 20892 times.
✓ Branch 5 taken 5781 times.
|
99420962 | for (std::size_t k = 0; k < connectivityMap[globalI].size(); ++k) |
| 403 |
6/6✓ Branch 0 taken 16772134 times.
✓ Branch 1 taken 16772134 times.
✓ Branch 2 taken 74115732 times.
✓ Branch 3 taken 65842968 times.
✓ Branch 4 taken 20892 times.
✓ Branch 5 taken 20892 times.
|
173544752 | for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ) |
| 404 | 110641592 | partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k], fvGeometry.scvf(scvfIdx)); | |
| 405 | |||
| 406 | 16784968 | return partialDerivsTmp; | |
| 407 | }; | ||
| 408 | |||
| 409 | // derive the residuals numerically | ||
| 410 |
5/6✓ Branch 0 taken 264 times.
✓ Branch 1 taken 16760402 times.
✓ Branch 3 taken 231 times.
✓ Branch 4 taken 33 times.
✓ Branch 6 taken 231 times.
✗ Branch 7 not taken.
|
30478819 | static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod"); |
| 411 |
5/6✓ Branch 0 taken 278 times.
✓ Branch 1 taken 16760388 times.
✓ Branch 3 taken 231 times.
✓ Branch 4 taken 47 times.
✓ Branch 6 taken 231 times.
✗ Branch 7 not taken.
|
30478819 | static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()}; |
| 412 |
1/2✓ Branch 1 taken 1168401 times.
✗ Branch 2 not taken.
|
30478819 | NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals, |
| 413 | 30478819 | eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod); | |
| 414 | |||
| 415 | // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the | ||
| 416 | // current primary variable and a 0 elsewhere. As we always solve for a delta of the | ||
| 417 | // solution with respect to the initial one, this results in a delta of 0 for ghosts. | ||
| 418 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16760666 times.
|
30478819 | if (this->elementIsGhost()) |
| 419 | { | ||
| 420 | ✗ | partialDerivs[0] = 0.0; | |
| 421 | ✗ | partialDerivs[0][pvIdx] = 1.0; | |
| 422 | } | ||
| 423 | |||
| 424 | // restore the original state of the scv's volume variables | ||
| 425 | 30478819 | curVolVars = origVolVars; | |
| 426 | |||
| 427 | // restore the current element solution | ||
| 428 |
2/2✓ Branch 0 taken 167 times.
✓ Branch 1 taken 231333 times.
|
30478819 | elemSol[0][pvIdx] = origPriVars[pvIdx]; |
| 429 | |||
| 430 | // restore the undeflected state of the coupling context | ||
| 431 |
2/3✓ Branch 0 taken 705218 times.
✓ Branch 1 taken 463360 times.
✗ Branch 2 not taken.
|
30478819 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx); |
| 432 | |||
| 433 | // add the current partial derivatives to the global jacobian matrix | ||
| 434 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 435 | { | ||
| 436 | // check if own element has internal Dirichlet constraint | ||
| 437 |
2/3✓ Branch 0 taken 705218 times.
✓ Branch 1 taken 463360 times.
✗ Branch 2 not taken.
|
1633478 | const auto internalDirichletConstraintsOwnElement = this->problem().hasInternalDirichletConstraint(this->element(), scv); |
| 438 |
2/2✓ Branch 0 taken 1040 times.
✓ Branch 1 taken 1040 times.
|
1634518 | const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv); |
| 439 | |||
| 440 |
3/3✓ Branch 0 taken 2498678 times.
✓ Branch 1 taken 1633478 times.
✓ Branch 2 taken 200 times.
|
4132356 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) |
| 441 | { | ||
| 442 |
2/2✓ Branch 0 taken 6682 times.
✓ Branch 1 taken 1508196 times.
|
2498878 | if (internalDirichletConstraintsOwnElement[eqIdx]) |
| 443 | { | ||
| 444 |
2/2✓ Branch 0 taken 92 times.
✓ Branch 1 taken 6590 times.
|
6682 | origResiduals[0][eqIdx] = origVolVars.priVars()[eqIdx] - dirichletValues[eqIdx]; |
| 445 |
3/4✓ Branch 0 taken 92 times.
✓ Branch 1 taken 6590 times.
✓ Branch 3 taken 912 times.
✗ Branch 4 not taken.
|
6774 | A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0; |
| 446 | } | ||
| 447 | else | ||
| 448 |
1/2✓ Branch 1 taken 6160 times.
✗ Branch 2 not taken.
|
2492196 | A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx]; |
| 449 | } | ||
| 450 | |||
| 451 | // off-diagonal entries | ||
| 452 | 1633478 | j = 1; | |
| 453 |
2/2✓ Branch 0 taken 6387790 times.
✓ Branch 1 taken 1633478 times.
|
8021268 | for (const auto& dataJ : connectivityMap[globalI]) |
| 454 | { | ||
| 455 |
2/3✓ Branch 0 taken 6348878 times.
✓ Branch 1 taken 3112122 times.
✗ Branch 2 not taken.
|
9730430 | const auto& neighborElement = neighborElements[j-1]; |
| 456 |
2/3✓ Branch 0 taken 6348878 times.
✓ Branch 1 taken 3112122 times.
✗ Branch 2 not taken.
|
9730430 | const auto& neighborScv = fvGeometry.scv(dataJ.globalJ); |
| 457 |
2/3✓ Branch 0 taken 6135478 times.
✓ Branch 1 taken 3112122 times.
✗ Branch 2 not taken.
|
9517030 | const auto internalDirichletConstraintsNeighbor = this->problem().hasInternalDirichletConstraint(neighborElement, neighborScv); |
| 458 | |||
| 459 |
3/3✓ Branch 0 taken 9453800 times.
✓ Branch 1 taken 6387790 times.
✓ Branch 2 taken 265670 times.
|
16107260 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) |
| 460 | { | ||
| 461 |
2/2✓ Branch 0 taken 19592 times.
✓ Branch 1 taken 5930438 times.
|
9719470 | if (internalDirichletConstraintsNeighbor[eqIdx]) |
| 462 |
1/2✓ Branch 1 taken 2624 times.
✗ Branch 2 not taken.
|
19592 | A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0; |
| 463 | else | ||
| 464 |
1/2✓ Branch 1 taken 24640 times.
✗ Branch 2 not taken.
|
9699878 | A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx]; |
| 465 | } | ||
| 466 | |||
| 467 | 6387790 | ++j; | |
| 468 | } | ||
| 469 | } | ||
| 470 | else // no internal Dirichlet constraints specified | ||
| 471 | { | ||
| 472 |
2/2✓ Branch 0 taken 18128684 times.
✓ Branch 1 taken 15127188 times.
|
62816770 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) |
| 473 | { | ||
| 474 | // the diagonal entries | ||
| 475 |
1/2✓ Branch 1 taken 2838641 times.
✗ Branch 2 not taken.
|
33971429 | A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx]; |
| 476 | |||
| 477 | // off-diagonal entries | ||
| 478 | 33971429 | j = 1; | |
| 479 |
2/2✓ Branch 0 taken 91014984 times.
✓ Branch 1 taken 18128684 times.
|
207421191 | for (const auto& dataJ : connectivityMap[globalI]) |
| 480 |
1/2✓ Branch 1 taken 16094574 times.
✗ Branch 2 not taken.
|
173449762 | A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx]; |
| 481 | } | ||
| 482 | } | ||
| 483 | } | ||
| 484 | |||
| 485 | // restore original state of the flux vars cache in case of global caching. | ||
| 486 | // This has to be done in order to guarantee that everything is in an undeflected | ||
| 487 | // state before the assembly of another element is called. In the case of local caching | ||
| 488 | // this is obsolete because the elemFluxVarsCache used here goes out of scope after this. | ||
| 489 | // We only have to do this for the last primary variable, for all others the flux var cache | ||
| 490 | // is updated with the correct element volume variables before residual evaluations | ||
| 491 | if (enableGridFluxVarsCache) | ||
| 492 | 25435893 | gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars); | |
| 493 | |||
| 494 | // compute potential additional derivatives causes by the coupling | ||
| 495 | 28110307 | typename ParentType::LocalResidual::ElementResidualVector origElementResidual({origResiduals[0]}); | |
| 496 | 28110307 | this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origElementResidual, A, gridVariables); | |
| 497 | |||
| 498 | // return the original residual | ||
| 499 | 28110307 | return origResiduals[0]; | |
| 500 | 1138722 | } | |
| 501 | |||
| 502 | /*! | ||
| 503 | * \brief Computes the derivatives of the residual of the given element with respect to primary | ||
| 504 | * variables of domainJ and adds them to the global matrix. | ||
| 505 | */ | ||
| 506 | template<std::size_t otherId, class JacobianBlock, class GridVariables> | ||
| 507 | 31069243 | void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A, | |
| 508 | const LocalResidualValues& res, GridVariables& gridVariables) | ||
| 509 | { | ||
| 510 | //////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 511 | // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. // | ||
| 512 | //////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 513 | |||
| 514 | // get some aliases for convenience | ||
| 515 | 31069243 | const auto& element = this->element(); | |
| 516 | 31069243 | const auto& fvGeometry = this->fvGeometry(); | |
| 517 | 31069243 | const auto& gridGeometry = fvGeometry.gridGeometry(); | |
| 518 | 31069243 | auto&& curElemVolVars = this->curElemVolVars(); | |
| 519 | 31069243 | auto&& elemFluxVarsCache = this->elemFluxVarsCache(); | |
| 520 | |||
| 521 | // get stencil information | ||
| 522 | 31069243 | const auto globalI = gridGeometry.elementMapper().index(element); | |
| 523 | 31069243 | const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); | |
| 524 | 31069243 | const auto& curSolJ = this->curSol()[domainJ]; | |
| 525 | |||
| 526 | // make sure the flux variables cache does not contain any artifacts from prior differentiation | ||
| 527 | 4636352 | elemFluxVarsCache.update(element, fvGeometry, curElemVolVars); | |
| 528 | |||
| 529 | // convenience lambda for call to update self | ||
| 530 | 73022636 | auto updateCoupledVariables = [&] () | |
| 531 | { | ||
| 532 | // Update ourself after the context has been modified. Depending on the | ||
| 533 | // type of caching, other objects might have to be updated. All ifs can be optimized away. | ||
| 534 | if (enableGridFluxVarsCache) | ||
| 535 | { | ||
| 536 | if (enableGridVolVarsCache) | ||
| 537 | { | ||
| 538 | 33023075 | this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache()); | |
| 539 | 41217685 | this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache); | |
| 540 | } | ||
| 541 | else | ||
| 542 | { | ||
| 543 | this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache()); | ||
| 544 | this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache); | ||
| 545 | } | ||
| 546 | } | ||
| 547 | else | ||
| 548 | { | ||
| 549 | if (enableGridVolVarsCache) | ||
| 550 | this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache); | ||
| 551 | else | ||
| 552 | 813400 | this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache); | |
| 553 | } | ||
| 554 | }; | ||
| 555 | |||
| 556 |
2/2✓ Branch 0 taken 16154353 times.
✓ Branch 1 taken 16452553 times.
|
55816251 | for (const auto& globalJ : stencil) |
| 557 | { | ||
| 558 | // undeflected privars and privars to be deflected | ||
| 559 |
1/2✓ Branch 1 taken 8069868 times.
✗ Branch 2 not taken.
|
24747008 | const auto origPriVarsJ = curSolJ[globalJ]; |
| 560 | 24747008 | auto priVarsJ = origPriVarsJ; | |
| 561 | |||
| 562 |
1/2✓ Branch 1 taken 8069868 times.
✗ Branch 2 not taken.
|
25324448 | const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0]; |
| 563 | |||
| 564 |
2/2✓ Branch 0 taken 17109177 times.
✓ Branch 1 taken 16154353 times.
|
50737096 | for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx) |
| 565 | { | ||
| 566 | 43196473 | auto evalCouplingResidual = [&](Scalar priVar) | |
| 567 | { | ||
| 568 | // update the volume variables and the flux var cache | ||
| 569 | 16902073 | priVarsJ[pvIdx] = priVar; | |
| 570 | 17206385 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); | |
| 571 | 17123029 | updateCoupledVariables(); | |
| 572 | 17206385 | return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0]; | |
| 573 | }; | ||
| 574 | |||
| 575 | // derive the residuals numerically | ||
| 576 | 25990088 | LocalResidualValues partialDeriv(0.0); | |
| 577 |
2/2✓ Branch 0 taken 299 times.
✓ Branch 1 taken 17108878 times.
|
25990088 | const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); |
| 578 |
5/6✓ Branch 0 taken 299 times.
✓ Branch 1 taken 17108878 times.
✓ Branch 3 taken 285 times.
✓ Branch 4 taken 14 times.
✓ Branch 6 taken 285 times.
✗ Branch 7 not taken.
|
25990088 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 579 |
5/6✓ Branch 0 taken 303 times.
✓ Branch 1 taken 17108874 times.
✓ Branch 3 taken 285 times.
✓ Branch 4 taken 18 times.
✓ Branch 6 taken 285 times.
✗ Branch 7 not taken.
|
25990096 | static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); |
| 580 |
1/2✓ Branch 1 taken 8115676 times.
✗ Branch 2 not taken.
|
25990088 | NumericDifferentiation::partialDerivative(evalCouplingResidual, priVarsJ[pvIdx], partialDeriv, origResidual, |
| 581 | 25990088 | epsCoupl(priVarsJ[pvIdx], pvIdx), numDiffMethod); | |
| 582 | |||
| 583 | // add the current partial derivatives to the global jacobian matrix | ||
| 584 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 585 | { | ||
| 586 | 5339992 | const auto& scv = fvGeometry.scv(globalI); | |
| 587 |
4/4✓ Branch 0 taken 2820504 times.
✓ Branch 1 taken 1585932 times.
✓ Branch 2 taken 4172748 times.
✓ Branch 3 taken 944 times.
|
7830496 | const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); |
| 588 |
2/2✓ Branch 0 taken 4429652 times.
✓ Branch 1 taken 26340 times.
|
5009992 | if (internalDirichletConstraints.none()) |
| 589 | { | ||
| 590 |
2/2✓ Branch 0 taken 6506884 times.
✓ Branch 1 taken 5093652 times.
|
11600536 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) |
| 591 |
1/2✓ Branch 1 taken 6262244 times.
✗ Branch 2 not taken.
|
6506884 | A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx]; |
| 592 | } | ||
| 593 | else | ||
| 594 | { | ||
| 595 |
2/2✓ Branch 0 taken 26708 times.
✓ Branch 1 taken 26340 times.
|
53048 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) |
| 596 | { | ||
| 597 |
2/3✓ Branch 1 taken 1012 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 25696 times.
|
26708 | if (internalDirichletConstraints[eqIdx]) |
| 598 |
1/2✓ Branch 1 taken 22692 times.
✗ Branch 2 not taken.
|
26340 | A[globalI][globalJ][eqIdx][pvIdx] = 0.0; |
| 599 | else | ||
| 600 |
1/2✓ Branch 1 taken 368 times.
✗ Branch 2 not taken.
|
368 | A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx]; |
| 601 | } | ||
| 602 | } | ||
| 603 | } | ||
| 604 | else | ||
| 605 | { | ||
| 606 |
2/3✓ Branch 0 taken 13445819 times.
✓ Branch 1 taken 11989185 times.
✗ Branch 2 not taken.
|
44233460 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) |
| 607 |
1/2✓ Branch 1 taken 4058196 times.
✗ Branch 2 not taken.
|
23363364 | A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx]; |
| 608 | } | ||
| 609 | |||
| 610 | // restore the current element solution | ||
| 611 |
1/2✓ Branch 1 taken 72912 times.
✗ Branch 2 not taken.
|
25990088 | priVarsJ[pvIdx] = origPriVarsJ[pvIdx]; |
| 612 | |||
| 613 | // restore the undeflected state of the coupling context | ||
| 614 |
1/2✓ Branch 1 taken 83356 times.
✗ Branch 2 not taken.
|
25990088 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); |
| 615 | } | ||
| 616 | |||
| 617 | // restore original state of the flux vars cache and/or vol vars in case of global caching. | ||
| 618 | // This has to be done in order to guarantee that the undeflected residual computation done | ||
| 619 | // above is correct when jumping to the next coupled dof of domainJ. | ||
| 620 | 24747008 | updateCoupledVariables(); | |
| 621 | } | ||
| 622 | 31069243 | } | |
| 623 | }; | ||
| 624 | |||
| 625 | /*! | ||
| 626 | * \ingroup Assembly | ||
| 627 | * \ingroup CCDiscretization | ||
| 628 | * \ingroup MultiDomain | ||
| 629 | * \brief Cell-centered scheme multidomain local assembler using numeric differentiation and explicit time discretization | ||
| 630 | */ | ||
| 631 | template<std::size_t id, class TypeTag, class Assembler> | ||
| 632 | 2028000 | class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false> | |
| 633 | : public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, | ||
| 634 | SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, false>, false> | ||
| 635 | { | ||
| 636 | using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>; | ||
| 637 | using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/false>; | ||
| 638 | |||
| 639 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 640 | using LocalResidualValues = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 641 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 642 | |||
| 643 | static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq(); | ||
| 644 | static constexpr auto domainI = Dune::index_constant<id>(); | ||
| 645 | |||
| 646 | public: | ||
| 647 | 2028000 | using ParentType::ParentType; | |
| 648 | |||
| 649 | /*! | ||
| 650 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 651 | * to the global matrix. | ||
| 652 | * | ||
| 653 | * \note In an explicit scheme, only the storage terms need to be differentiated. | ||
| 654 | * Thus, this can be done as in the uncoupled case as the coupling can only | ||
| 655 | * enter sources or fluxes. | ||
| 656 | * | ||
| 657 | * \return The element residual at the current solution. | ||
| 658 | */ | ||
| 659 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 660 | 2028000 | LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | |
| 661 | { | ||
| 662 | 2028000 | if (this->assembler().isStationaryProblem()) | |
| 663 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual"); | |
| 664 | |||
| 665 | // assemble the undeflected residual | ||
| 666 | 2028000 | const auto residual = this->evalLocalResidual()[0]; | |
| 667 | 2028000 | const auto storageResidual = this->evalLocalStorageResidual(); | |
| 668 | |||
| 669 | ////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 670 | // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the // | ||
| 671 | // neighboring elements all derivatives are zero. For the assembled element only the storage // | ||
| 672 | // derivatives are non-zero. // | ||
| 673 | ////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 674 | |||
| 675 | // get some aliases for convenience | ||
| 676 | 2028000 | const auto& element = this->element(); | |
| 677 | 2028000 | const auto& fvGeometry = this->fvGeometry(); | |
| 678 | 2028000 | const auto& gridGeometry = fvGeometry.gridGeometry(); | |
| 679 | 2028000 | auto&& curElemVolVars = this->curElemVolVars(); | |
| 680 | |||
| 681 | // reference to the element's scv (needed later) and corresponding vol vars | ||
| 682 | 2028000 | const auto globalI = gridGeometry.elementMapper().index(element); | |
| 683 | 2028000 | const auto& scv = fvGeometry.scv(globalI); | |
| 684 | 2028000 | auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv); | |
| 685 | |||
| 686 | // save a copy of the original privars and vol vars in order | ||
| 687 | // to restore the original solution after deflection | ||
| 688 | 2028000 | const auto& curSol = this->curSol()[domainI]; | |
| 689 | 2028000 | const auto origPriVars = curSol[globalI]; | |
| 690 | 2028000 | const auto origVolVars = curVolVars; | |
| 691 | |||
| 692 | // element solution container to be deflected | ||
| 693 | 2028000 | auto elemSol = elementSolution(element, curSol, gridGeometry); | |
| 694 | |||
| 695 | // derivatives in the neighbors with respect to the current elements | ||
| 696 | LocalResidualValues partialDeriv; | ||
| 697 | 4056000 | for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) | |
| 698 | { | ||
| 699 | // reset derivatives of element dof with respect to itself | ||
| 700 | 2028000 | partialDeriv = 0.0; | |
| 701 | |||
| 702 | 4056000 | auto evalStorage = [&](Scalar priVar) | |
| 703 | { | ||
| 704 | // update the volume variables and calculate | ||
| 705 | // the residual with the deflected primary variables | ||
| 706 | 2028000 | elemSol[0][pvIdx] = priVar; | |
| 707 | 2028000 | curVolVars.update(elemSol, this->problem(), element, scv); | |
| 708 | 2028000 | return this->evalLocalStorageResidual(); | |
| 709 | }; | ||
| 710 | |||
| 711 | // for non-ghosts compute the derivative numerically | ||
| 712 | 2028000 | if (!this->elementIsGhost()) | |
| 713 | { | ||
| 714 | 2028000 | static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()}; | |
| 715 | 2028000 | static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod"); | |
| 716 | 2028000 | NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, storageResidual, | |
| 717 | 2028000 | eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod); | |
| 718 | } | ||
| 719 | |||
| 720 | // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else | ||
| 721 | // as we always solve for a delta of the solution with respect to the initial solution this | ||
| 722 | // results in a delta of zero for ghosts | ||
| 723 | ✗ | else partialDeriv[pvIdx] = 1.0; | |
| 724 | |||
| 725 | // restore the original state of the scv's volume variables | ||
| 726 | 2028000 | curVolVars = origVolVars; | |
| 727 | |||
| 728 | // restore the current element solution | ||
| 729 | 2028000 | elemSol[0][pvIdx] = origPriVars[pvIdx]; | |
| 730 | |||
| 731 | // add the current partial derivatives to the global jacobian matrix | ||
| 732 | // only diagonal entries for explicit jacobians | ||
| 733 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 734 | { | ||
| 735 | // check if own element has internal Dirichlet constraint | ||
| 736 | const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); | ||
| 737 | const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv); | ||
| 738 | |||
| 739 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) | ||
| 740 | { | ||
| 741 | if (internalDirichletConstraints[eqIdx]) | ||
| 742 | { | ||
| 743 | residual[eqIdx] = origVolVars.priVars()[eqIdx] - dirichletValues[eqIdx]; | ||
| 744 | A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0; | ||
| 745 | } | ||
| 746 | else | ||
| 747 | A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx]; | ||
| 748 | } | ||
| 749 | } | ||
| 750 | else | ||
| 751 | { | ||
| 752 | 4056000 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) | |
| 753 | 2028000 | A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx]; | |
| 754 | } | ||
| 755 | } | ||
| 756 | |||
| 757 | // return the original residual | ||
| 758 | 2028000 | return residual; | |
| 759 | } | ||
| 760 | |||
| 761 | /*! | ||
| 762 | * \brief Computes the derivatives of the residual of the given element with respect to primary | ||
| 763 | * variables of domainJ and adds them to the global matrix. | ||
| 764 | * \note Since the coupling can only enter sources or fluxes and these are evaluated on | ||
| 765 | * the old time level (explicit scheme), the coupling blocks are empty. | ||
| 766 | */ | ||
| 767 | template<std::size_t otherId, class JacobianBlock, class GridVariables> | ||
| 768 | void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A, | ||
| 769 | const LocalResidualValues& res, GridVariables& gridVariables) | ||
| 770 | {} | ||
| 771 | }; | ||
| 772 | |||
| 773 | /*! | ||
| 774 | * \ingroup Assembly | ||
| 775 | * \ingroup CCDiscretization | ||
| 776 | * \ingroup MultiDomain | ||
| 777 | * \brief Cell-centered scheme local assembler using analytic differentiation and implicit time discretization | ||
| 778 | */ | ||
| 779 | template<std::size_t id, class TypeTag, class Assembler> | ||
| 780 | class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true> | ||
| 781 | : public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, | ||
| 782 | SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, true>, true> | ||
| 783 | { | ||
| 784 | using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>; | ||
| 785 | using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/true>; | ||
| 786 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 787 | using LocalResidualValues = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 788 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 789 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 790 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 791 | |||
| 792 | enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() }; | ||
| 793 | enum { dim = GridView::dimension }; | ||
| 794 | |||
| 795 | static constexpr auto domainI = Dune::index_constant<id>(); | ||
| 796 | |||
| 797 | public: | ||
| 798 | using ParentType::ParentType; | ||
| 799 | |||
| 800 | /*! | ||
| 801 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 802 | * to the global matrix. | ||
| 803 | * | ||
| 804 | * \return The element residual at the current solution. | ||
| 805 | */ | ||
| 806 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 807 | LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | ||
| 808 | { | ||
| 809 | // treat ghost separately, we always want zero update for ghosts | ||
| 810 | if (this->elementIsGhost()) | ||
| 811 | { | ||
| 812 | const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(this->element()); | ||
| 813 | for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) | ||
| 814 | A[globalI][globalI][pvIdx][pvIdx] = 1.0; | ||
| 815 | |||
| 816 | // return zero residual | ||
| 817 | return LocalResidualValues(0.0); | ||
| 818 | } | ||
| 819 | |||
| 820 | // get some aliases for convenience | ||
| 821 | const auto& problem = this->problem(); | ||
| 822 | const auto& element = this->element(); | ||
| 823 | const auto& fvGeometry = this->fvGeometry(); | ||
| 824 | const auto& curElemVolVars = this->curElemVolVars(); | ||
| 825 | const auto& elemFluxVarsCache = this->elemFluxVarsCache(); | ||
| 826 | |||
| 827 | // get reference to the element's current vol vars | ||
| 828 | const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(element); | ||
| 829 | const auto& scv = fvGeometry.scv(globalI); | ||
| 830 | const auto& volVars = curElemVolVars[scv]; | ||
| 831 | |||
| 832 | // if the problem is instationary, add derivative of storage term | ||
| 833 | if (!this->assembler().isStationaryProblem()) | ||
| 834 | this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); | ||
| 835 | |||
| 836 | // add source term derivatives | ||
| 837 | this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); | ||
| 838 | |||
| 839 | // add flux derivatives for each scvf | ||
| 840 | for (const auto& scvf : scvfs(fvGeometry)) | ||
| 841 | { | ||
| 842 | // inner faces | ||
| 843 | if (!scvf.boundary()) | ||
| 844 | this->localResidual().addFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); | ||
| 845 | |||
| 846 | // boundary faces | ||
| 847 | else | ||
| 848 | { | ||
| 849 | const auto& bcTypes = problem.boundaryTypes(element, scvf); | ||
| 850 | |||
| 851 | // add Dirichlet boundary flux derivatives | ||
| 852 | if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) | ||
| 853 | this->localResidual().addCCDirichletFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); | ||
| 854 | |||
| 855 | // add Robin ("solution dependent Neumann") boundary flux derivatives | ||
| 856 | else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) | ||
| 857 | this->localResidual().addRobinFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); | ||
| 858 | |||
| 859 | else | ||
| 860 | DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs"); | ||
| 861 | } | ||
| 862 | } | ||
| 863 | |||
| 864 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 865 | { | ||
| 866 | // check if own element has internal Dirichlet constraint | ||
| 867 | const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); | ||
| 868 | const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv); | ||
| 869 | |||
| 870 | auto residual = this->evalLocalResidual()[0]; | ||
| 871 | |||
| 872 | for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) | ||
| 873 | { | ||
| 874 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) | ||
| 875 | { | ||
| 876 | if (internalDirichletConstraints[eqIdx]) | ||
| 877 | { | ||
| 878 | residual[eqIdx] = volVars.priVars()[eqIdx] - dirichletValues[eqIdx]; | ||
| 879 | A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0; | ||
| 880 | |||
| 881 | // inner faces | ||
| 882 | for (const auto& scvf : scvfs(fvGeometry)) | ||
| 883 | if (!scvf.boundary()) | ||
| 884 | A[globalI][fvGeometry.scv(scvf.outsideScvIdx()).dofIndex()][eqIdx][pvIdx] = 0.0; | ||
| 885 | } | ||
| 886 | } | ||
| 887 | } | ||
| 888 | return residual; | ||
| 889 | } | ||
| 890 | else | ||
| 891 | // return element residual | ||
| 892 | return this->evalLocalResidual()[0]; | ||
| 893 | } | ||
| 894 | |||
| 895 | /*! | ||
| 896 | * \brief Computes the derivatives of the residual of the given element with respect to primary | ||
| 897 | * variables of domainJ and adds them to the global matrix. | ||
| 898 | */ | ||
| 899 | template<std::size_t otherId, class JacobianBlock, class GridVariables> | ||
| 900 | void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A, | ||
| 901 | const LocalResidualValues& res, GridVariables& gridVariables) | ||
| 902 | { | ||
| 903 | //////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 904 | // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. // | ||
| 905 | //////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 906 | |||
| 907 | // get some aliases for convenience | ||
| 908 | const auto& element = this->element(); | ||
| 909 | const auto& fvGeometry = this->fvGeometry(); | ||
| 910 | const auto& gridGeometry = fvGeometry.gridGeometry(); | ||
| 911 | auto&& curElemVolVars = this->curElemVolVars(); | ||
| 912 | // auto&& elemFluxVarsCache = this->elemFluxVarsCache(); | ||
| 913 | |||
| 914 | // get stencil information | ||
| 915 | const auto globalI = gridGeometry.elementMapper().index(element); | ||
| 916 | const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); | ||
| 917 | |||
| 918 | for (const auto globalJ : stencil) | ||
| 919 | { | ||
| 920 | const auto& elementJ = this->assembler().gridGeometry(domainJ).element(globalJ); | ||
| 921 | this->couplingManager().addCouplingDerivatives(A[globalI][globalJ], domainI, element, fvGeometry, curElemVolVars, domainJ, elementJ); | ||
| 922 | |||
| 923 | // add the current partial derivatives to the global jacobian matrix | ||
| 924 | if constexpr (Problem::enableInternalDirichletConstraints()) | ||
| 925 | { | ||
| 926 | const auto& scv = fvGeometry.scv(globalI); | ||
| 927 | const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); | ||
| 928 | if (internalDirichletConstraints.any()) | ||
| 929 | { | ||
| 930 | for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) | ||
| 931 | for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) | ||
| 932 | if (internalDirichletConstraints[eqIdx]) | ||
| 933 | A[globalI][globalJ][eqIdx][pvIdx] = 0.0; | ||
| 934 | } | ||
| 935 | } | ||
| 936 | } | ||
| 937 | } | ||
| 938 | }; | ||
| 939 | |||
| 940 | } // end namespace Dumux | ||
| 941 | |||
| 942 | #endif | ||
| 943 |