| 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 MultiDomain | ||
| 10 | * \brief Freeflow coupling managers (Navier-Stokes mass-momentum coupling) | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MULTIDOMAIN_FREEFLOW_COUPLING_MANAGER_CVFE_HH | ||
| 13 | #define DUMUX_MULTIDOMAIN_FREEFLOW_COUPLING_MANAGER_CVFE_HH | ||
| 14 | |||
| 15 | #include <memory> | ||
| 16 | #include <tuple> | ||
| 17 | #include <vector> | ||
| 18 | #include <deque> | ||
| 19 | |||
| 20 | #include <dune/common/exceptions.hh> | ||
| 21 | #include <dune/common/indices.hh> | ||
| 22 | #include <dune/common/float_cmp.hh> | ||
| 23 | #include <dune/geometry/referenceelements.hh> | ||
| 24 | |||
| 25 | #include <dumux/common/properties.hh> | ||
| 26 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 27 | |||
| 28 | #include <dumux/discretization/method.hh> | ||
| 29 | #include <dumux/discretization/evalsolution.hh> | ||
| 30 | #include <dumux/discretization/elementsolution.hh> | ||
| 31 | |||
| 32 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 33 | #include <dumux/multidomain/fvassembler.hh> | ||
| 34 | |||
| 35 | #include <dumux/parallel/parallel_for.hh> | ||
| 36 | #include <dumux/assembly/coloring.hh> | ||
| 37 | |||
| 38 | #include "typetraits.hh" | ||
| 39 | |||
| 40 | namespace Dumux { | ||
| 41 | |||
| 42 | /*! | ||
| 43 | * \ingroup MultiDomain | ||
| 44 | * \brief The interface of the coupling manager for free flow systems | ||
| 45 | * \note coupling manager for control volume finite element schemes | ||
| 46 | */ | ||
| 47 | template<class Traits> | ||
| 48 | class CVFEFreeFlowCouplingManager | ||
| 49 | : public CouplingManager<Traits> | ||
| 50 | { | ||
| 51 | using ParentType = CouplingManager<Traits>; | ||
| 52 | public: | ||
| 53 | static constexpr auto freeFlowMomentumIndex = typename Traits::template SubDomain<0>::Index(); | ||
| 54 | static constexpr auto freeFlowMassIndex = typename Traits::template SubDomain<1>::Index(); | ||
| 55 | |||
| 56 | // this can be used if the coupling manager is used inside a meta-coupling manager (e.g. multi-binary) | ||
| 57 | // to manager the solution vector storage outside this class | ||
| 58 | using SolutionVectorStorage = typename ParentType::SolutionVectorStorage; | ||
| 59 | private: | ||
| 60 | template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomain<id>::TypeTag; | ||
| 61 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 62 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 63 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 64 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 65 | template<std::size_t id> using ElementSeed = typename GridView<id>::Grid::template Codim<0>::EntitySeed; | ||
| 66 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 67 | template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume; | ||
| 68 | template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace; | ||
| 69 | template<std::size_t id> using GridVariables = typename Traits::template SubDomain<id>::GridVariables; | ||
| 70 | template<std::size_t id> using ElementVolumeVariables = typename GridVariables<id>::GridVolumeVariables::LocalView; | ||
| 71 | template<std::size_t id> using GridFluxVariablesCache = typename GridVariables<id>::GridFluxVariablesCache; | ||
| 72 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 73 | template<std::size_t id> using VolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::VolumeVariables>; | ||
| 74 | |||
| 75 | using Scalar = typename Traits::Scalar; | ||
| 76 | using SolutionVector = typename Traits::SolutionVector; | ||
| 77 | |||
| 78 | using CouplingStencilType = std::vector<std::size_t>; | ||
| 79 | |||
| 80 | using GridVariablesTuple = typename Traits::template TupleOfSharedPtr<GridVariables>; | ||
| 81 | |||
| 82 | using FluidSystem = typename VolumeVariables<freeFlowMassIndex>::FluidSystem; | ||
| 83 | |||
| 84 | using VelocityVector = typename SubControlVolumeFace<freeFlowMassIndex>::GlobalPosition; | ||
| 85 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 86 | |||
| 87 | static_assert(std::is_same_v<VelocityVector, typename SubControlVolumeFace<freeFlowMomentumIndex>::GlobalPosition>); | ||
| 88 | |||
| 89 |
3/14✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✗ Branch 5 not taken.
|
18 | struct MomentumCouplingContext |
| 90 | { | ||
| 91 | FVElementGeometry<freeFlowMassIndex> fvGeometry; | ||
| 92 | ElementVolumeVariables<freeFlowMassIndex> curElemVolVars; | ||
| 93 | ElementVolumeVariables<freeFlowMassIndex> prevElemVolVars; | ||
| 94 | std::size_t eIdx; | ||
| 95 | }; | ||
| 96 | |||
| 97 |
1/8✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
4 | struct MassAndEnergyCouplingContext |
| 98 | { | ||
| 99 | 12 | MassAndEnergyCouplingContext(FVElementGeometry<freeFlowMomentumIndex>&& f, const std::size_t i) | |
| 100 |
1/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12 | : fvGeometry(std::move(f)) |
| 101 | 12 | , eIdx(i) | |
| 102 | {} | ||
| 103 | |||
| 104 | FVElementGeometry<freeFlowMomentumIndex> fvGeometry; | ||
| 105 | std::size_t eIdx; | ||
| 106 | }; | ||
| 107 | |||
| 108 | using MomentumDiscretizationMethod = typename GridGeometry<freeFlowMomentumIndex>::DiscretizationMethod; | ||
| 109 | using MassDiscretizationMethod = typename GridGeometry<freeFlowMassIndex>::DiscretizationMethod; | ||
| 110 | |||
| 111 | public: | ||
| 112 | |||
| 113 | static constexpr auto pressureIdx = VolumeVariables<freeFlowMassIndex>::Indices::pressureIdx; | ||
| 114 | |||
| 115 | /*! | ||
| 116 | * \brief Methods to be accessed by main | ||
| 117 | */ | ||
| 118 | // \{ | ||
| 119 | |||
| 120 | //! use as regular coupling manager | ||
| 121 | 12 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | |
| 122 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 123 | GridVariablesTuple&& gridVariables, | ||
| 124 | const SolutionVector& curSol) | ||
| 125 | { | ||
| 126 | 24 | this->momentumCouplingContext_().clear(); | |
| 127 | 24 | this->massAndEnergyCouplingContext_().clear(); | |
| 128 | |||
| 129 | 24 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | |
| 130 | 12 | gridVariables_ = gridVariables; | |
| 131 | 12 | this->updateSolution(curSol); | |
| 132 | |||
| 133 | 12 | computeCouplingStencils_(); | |
| 134 | 12 | } | |
| 135 | |||
| 136 | //! use as regular coupling manager in a transient setting | ||
| 137 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | ||
| 138 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 139 | GridVariablesTuple&& gridVariables, | ||
| 140 | const SolutionVector& curSol, | ||
| 141 | const SolutionVector& prevSol) | ||
| 142 | { | ||
| 143 | init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol); | ||
| 144 | prevSol_ = &prevSol; | ||
| 145 | isTransient_ = true; | ||
| 146 | } | ||
| 147 | |||
| 148 | //! use as binary coupling manager in multi model context | ||
| 149 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | ||
| 150 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 151 | GridVariablesTuple&& gridVariables, | ||
| 152 | const typename ParentType::SolutionVectorStorage& curSol) | ||
| 153 | { | ||
| 154 | this->momentumCouplingContext_().clear(); | ||
| 155 | this->massAndEnergyCouplingContext_().clear(); | ||
| 156 | |||
| 157 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | ||
| 158 | gridVariables_ = gridVariables; | ||
| 159 | this->attachSolution(curSol); | ||
| 160 | |||
| 161 | computeCouplingStencils_(); | ||
| 162 | } | ||
| 163 | |||
| 164 | // \} | ||
| 165 | |||
| 166 | /*! | ||
| 167 | * \name member functions concerning the coupling stencils | ||
| 168 | */ | ||
| 169 | // \{ | ||
| 170 | |||
| 171 | /*! | ||
| 172 | * \brief Returns the pressure at a given sub control volume face | ||
| 173 | */ | ||
| 174 | 40586862 | Scalar pressure(const Element<freeFlowMomentumIndex>& element, | |
| 175 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 176 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 177 | const bool considerPreviousTimeStep = false) const | ||
| 178 | { | ||
| 179 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 40586862 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
40586862 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 180 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 40586862 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40586862 times.
|
40586862 | const auto& gg = this->problem(freeFlowMassIndex).gridGeometry(); |
| 181 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40586862 times.
|
81173724 | const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex] |
| 182 | 40586862 | : this->curSol(freeFlowMassIndex); | |
| 183 | 40586862 | const auto elemSol = elementSolution(element, sol, gg); | |
| 184 |
1/2✓ Branch 2 taken 1545600 times.
✗ Branch 3 not taken.
|
40586862 | return evalSolution(element, element.geometry(), gg, elemSol, scvf.ipGlobal())[pressureIdx]; |
| 185 | } | ||
| 186 | |||
| 187 | /*! | ||
| 188 | * \brief Returns the pressure at a given sub control volume | ||
| 189 | */ | ||
| 190 | Scalar pressure(const Element<freeFlowMomentumIndex>& element, | ||
| 191 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 192 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 193 | const bool considerPreviousTimeStep = false) const | ||
| 194 | { | ||
| 195 | assert(!(considerPreviousTimeStep && !this->isTransient_)); | ||
| 196 | const auto& gg = this->problem(freeFlowMassIndex).gridGeometry(); | ||
| 197 | const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex] | ||
| 198 | : this->curSol(freeFlowMassIndex); | ||
| 199 | const auto elemSol = elementSolution(element, sol, gg); | ||
| 200 | return evalSolution(element, element.geometry(), gg, elemSol, scv.dofPosition())[pressureIdx]; | ||
| 201 | } | ||
| 202 | |||
| 203 | /*! | ||
| 204 | * \brief Returns the density at a given sub control volume face. | ||
| 205 | */ | ||
| 206 | ✗ | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 207 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 208 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 209 | const bool considerPreviousTimeStep = false) const | ||
| 210 | { | ||
| 211 | ✗ | assert(!(considerPreviousTimeStep && !this->isTransient_)); | |
| 212 | ✗ | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 213 | |||
| 214 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 215 | { | ||
| 216 | ✗ | const auto eIdx = fvGeometry.elementIndex(); | |
| 217 | ✗ | const auto& scv = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | |
| 218 | |||
| 219 | ✗ | const auto& volVars = considerPreviousTimeStep ? | |
| 220 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 221 | ✗ | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 222 | |||
| 223 | ✗ | return volVars.density(); | |
| 224 | } | ||
| 225 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 226 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 227 | { | ||
| 228 | // TODO: cache the shape values when Box method is used | ||
| 229 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 230 | ✗ | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 231 | ✗ | std::vector<ShapeValue> shapeValues; | |
| 232 | ✗ | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); | |
| 233 | ✗ | localBasis.evaluateFunction(ipLocal, shapeValues); | |
| 234 | |||
| 235 | ✗ | Scalar rho = 0.0; | |
| 236 | ✗ | for (const auto& scv : scvs(this->momentumCouplingContext_()[0].fvGeometry)) | |
| 237 | { | ||
| 238 | ✗ | const auto& volVars = considerPreviousTimeStep ? | |
| 239 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 240 | ✗ | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 241 | ✗ | rho += volVars.density()*shapeValues[scv.indexInElement()][0]; | |
| 242 | } | ||
| 243 | |||
| 244 | ✗ | return rho; | |
| 245 | ✗ | } | |
| 246 | else | ||
| 247 | DUNE_THROW(Dune::NotImplemented, | ||
| 248 | "Density interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 249 | ); | ||
| 250 | } | ||
| 251 | |||
| 252 | /*! | ||
| 253 | * \brief Returns the density at a given sub control volume. | ||
| 254 | */ | ||
| 255 | 1470054 | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 256 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 257 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 258 | const bool considerPreviousTimeStep = false) const | ||
| 259 | { | ||
| 260 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1470054 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1470054 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 261 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 356454 times.
|
1470054 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, scv.elementIndex()); |
| 262 | |||
| 263 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 264 | { | ||
| 265 | 356454 | const auto eIdx = scv.elementIndex(); | |
| 266 | 356454 | const auto& scvI = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | |
| 267 | |||
| 268 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 356454 times.
|
356454 | const auto& volVars = considerPreviousTimeStep ? |
| 269 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | |
| 270 | 356454 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | |
| 271 | |||
| 272 | 356454 | return volVars.density(); | |
| 273 | } | ||
| 274 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 275 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 276 | { | ||
| 277 | // TODO: cache the shape values when Box method is used | ||
| 278 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 279 | 2227200 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 280 |
1/2✓ Branch 1 taken 1113600 times.
✗ Branch 2 not taken.
|
1113600 | std::vector<ShapeValue> shapeValues; |
| 281 |
1/4✓ Branch 3 taken 1113600 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
1113600 | const auto ipLocal = element.geometry().local(scv.dofPosition()); |
| 282 | 1113600 | localBasis.evaluateFunction(ipLocal, shapeValues); | |
| 283 | |||
| 284 | 1113600 | Scalar rho = 0.0; | |
| 285 |
2/2✓ Branch 0 taken 3340800 times.
✓ Branch 1 taken 1113600 times.
|
4454400 | for (const auto& scvI : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 286 | { | ||
| 287 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3340800 times.
|
6681600 | const auto& volVars = considerPreviousTimeStep ? |
| 288 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | |
| 289 | 3340800 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | |
| 290 | 3340800 | rho += volVars.density()*shapeValues[scvI.indexInElement()][0]; | |
| 291 | } | ||
| 292 |
1/2✓ Branch 0 taken 1113600 times.
✗ Branch 1 not taken.
|
1113600 | return rho; |
| 293 | 1113600 | } | |
| 294 | else | ||
| 295 | DUNE_THROW(Dune::NotImplemented, | ||
| 296 | "Density interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 297 | ); | ||
| 298 | } | ||
| 299 | |||
| 300 | /*! | ||
| 301 | * \brief Returns the effective viscosity at a given sub control volume face. | ||
| 302 | */ | ||
| 303 | 2263188 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 304 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 305 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 306 | const bool considerPreviousTimeStep = false) const | ||
| 307 | { | ||
| 308 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2263188 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2263188 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 309 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 712908 times.
|
2263188 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); |
| 310 | |||
| 311 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 312 | { | ||
| 313 | 712908 | const auto eIdx = fvGeometry.elementIndex(); | |
| 314 | 712908 | const auto& scv = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | |
| 315 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 712908 times.
|
712908 | const auto& volVars = considerPreviousTimeStep ? |
| 316 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 317 | 712908 | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 318 | 712908 | return volVars.viscosity(); | |
| 319 | } | ||
| 320 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 321 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 322 | { | ||
| 323 | // TODO: cache the shape values when Box method is used | ||
| 324 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 325 | 3100560 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 326 |
1/2✓ Branch 1 taken 1550280 times.
✗ Branch 2 not taken.
|
1550280 | std::vector<ShapeValue> shapeValues; |
| 327 |
1/4✓ Branch 3 taken 1550280 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
1550280 | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); |
| 328 | 1550280 | localBasis.evaluateFunction(ipLocal, shapeValues); | |
| 329 | |||
| 330 | 1550280 | Scalar mu = 0.0; | |
| 331 |
2/2✓ Branch 0 taken 4650840 times.
✓ Branch 1 taken 1550280 times.
|
6201120 | for (const auto& scv : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 332 | { | ||
| 333 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4650840 times.
|
9301680 | const auto& volVars = considerPreviousTimeStep ? |
| 334 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 335 | 4650840 | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 336 | 4650840 | mu += volVars.viscosity()*shapeValues[scv.indexInElement()][0]; | |
| 337 | } | ||
| 338 | |||
| 339 |
1/2✓ Branch 0 taken 1550280 times.
✗ Branch 1 not taken.
|
1550280 | return mu; |
| 340 | 1550280 | } | |
| 341 | else | ||
| 342 | DUNE_THROW(Dune::NotImplemented, | ||
| 343 | "Viscosity interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 344 | ); | ||
| 345 | } | ||
| 346 | |||
| 347 | /*! | ||
| 348 | * \brief Returns the effective viscosity at a given sub control volume. | ||
| 349 | */ | ||
| 350 | 244500 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 351 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 352 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 353 | const bool considerPreviousTimeStep = false) const | ||
| 354 | { | ||
| 355 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
244500 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 356 | 244500 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 357 | |||
| 358 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 359 | { | ||
| 360 | const auto eIdx = fvGeometry.elementIndex(); | ||
| 361 | const auto& scvI = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | ||
| 362 | const auto& volVars = considerPreviousTimeStep ? | ||
| 363 | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | ||
| 364 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | ||
| 365 | return volVars.viscosity(); | ||
| 366 | } | ||
| 367 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 368 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 369 | { | ||
| 370 | // TODO: cache the shape values when Box method is used | ||
| 371 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 372 | 489000 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 373 |
1/2✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
|
244500 | std::vector<ShapeValue> shapeValues; |
| 374 |
1/4✓ Branch 3 taken 244500 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
244500 | const auto ipLocal = element.geometry().local(scv.dofPosition()); |
| 375 | 244500 | localBasis.evaluateFunction(ipLocal, shapeValues); | |
| 376 | |||
| 377 | 244500 | Scalar mu = 0.0; | |
| 378 |
2/2✓ Branch 0 taken 733500 times.
✓ Branch 1 taken 244500 times.
|
978000 | for (const auto& scvI : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 379 | { | ||
| 380 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 733500 times.
|
1467000 | const auto& volVars = considerPreviousTimeStep ? |
| 381 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | |
| 382 | 733500 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | |
| 383 | 733500 | mu += volVars.viscosity()*shapeValues[scvI.indexInElement()][0]; | |
| 384 | } | ||
| 385 | |||
| 386 |
1/2✓ Branch 0 taken 244500 times.
✗ Branch 1 not taken.
|
244500 | return mu; |
| 387 | 244500 | } | |
| 388 | else | ||
| 389 | DUNE_THROW(Dune::NotImplemented, | ||
| 390 | "Viscosity interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 391 | ); | ||
| 392 | } | ||
| 393 | |||
| 394 | /*! | ||
| 395 | * \brief Returns the velocity at a given sub control volume face. | ||
| 396 | */ | ||
| 397 | 26761544 | VelocityVector faceVelocity(const Element<freeFlowMassIndex>& element, | |
| 398 | const SubControlVolumeFace<freeFlowMassIndex>& scvf) const | ||
| 399 | { | ||
| 400 | // TODO: optimize this function for tpfa where the scvf ip coincides with the dof location | ||
| 401 | |||
| 402 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 26761544 times.
|
26761544 | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(element); |
| 403 | 26761544 | bindCouplingContext_(Dune::index_constant<freeFlowMassIndex>(), element, eIdx); | |
| 404 | |||
| 405 | 26761544 | const auto& fvGeometry = this->massAndEnergyCouplingContext_()[0].fvGeometry; | |
| 406 | 26761544 | const auto& localBasis = fvGeometry.feLocalBasis(); | |
| 407 | |||
| 408 |
1/2✓ Branch 1 taken 26295284 times.
✗ Branch 2 not taken.
|
26761544 | std::vector<ShapeValue> shapeValues; |
| 409 |
2/6✓ Branch 1 taken 25329524 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 257920 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
26761544 | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); |
| 410 |
1/2✓ Branch 1 taken 26503624 times.
✗ Branch 2 not taken.
|
26761544 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 411 | |||
| 412 | // interpolate velocity at scvf | ||
| 413 | 26761544 | VelocityVector velocity(0.0); | |
| 414 |
2/2✓ Branch 0 taken 110138980 times.
✓ Branch 1 taken 26761544 times.
|
136900524 | for (const auto& scv : scvs(fvGeometry)) |
| 415 | 220277960 | velocity.axpy(shapeValues[scv.localDofIndex()][0], this->curSol(freeFlowMomentumIndex)[scv.dofIndex()]); | |
| 416 | |||
| 417 |
1/2✓ Branch 0 taken 24304176 times.
✗ Branch 1 not taken.
|
26761544 | return velocity; |
| 418 | 26761544 | } | |
| 419 | |||
| 420 | /*! | ||
| 421 | * \brief Returns the velocity at the element center. | ||
| 422 | */ | ||
| 423 | 332172 | VelocityVector elementVelocity(const FVElementGeometry<freeFlowMassIndex>& fvGeometry) const | |
| 424 | { | ||
| 425 | 332172 | bindCouplingContext_(Dune::index_constant<freeFlowMassIndex>(), fvGeometry.element()); | |
| 426 | |||
| 427 | 332172 | const auto& momentumFvGeometry = this->massAndEnergyCouplingContext_()[0].fvGeometry; | |
| 428 | 332172 | const auto& localBasis = momentumFvGeometry.feLocalBasis(); | |
| 429 | |||
| 430 | // interpolate velocity at scvf | ||
| 431 | 332172 | VelocityVector velocity(0.0); | |
| 432 |
1/2✓ Branch 1 taken 6400 times.
✗ Branch 2 not taken.
|
332172 | std::vector<ShapeValue> shapeValues; |
| 433 |
2/6✓ Branch 1 taken 332172 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 332172 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
332172 | localBasis.evaluateFunction(referenceElement(fvGeometry.element()).position(0,0), shapeValues); |
| 434 | |||
| 435 |
2/2✓ Branch 0 taken 1298696 times.
✓ Branch 1 taken 332172 times.
|
1630868 | for (const auto& scv : scvs(momentumFvGeometry)) |
| 436 | 2597392 | velocity.axpy(shapeValues[scv.localDofIndex()][0], this->curSol(freeFlowMomentumIndex)[scv.dofIndex()]); | |
| 437 | |||
| 438 |
1/2✓ Branch 0 taken 300488 times.
✗ Branch 1 not taken.
|
332172 | return velocity; |
| 439 | 332172 | } | |
| 440 | |||
| 441 | /*! | ||
| 442 | * \brief The coupling stencil of domain I, i.e. which domain J DOFs | ||
| 443 | * the given domain I element's residual depends on. | ||
| 444 | */ | ||
| 445 | template<std::size_t j> | ||
| 446 | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | ||
| 447 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 448 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 449 | Dune::index_constant<j> domainJ) const | ||
| 450 | { return emptyStencil_; } | ||
| 451 | |||
| 452 | /*! | ||
| 453 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 454 | * that couple with / influence the element residual of the given element of domain i | ||
| 455 | * | ||
| 456 | * \param domainI the domain index of domain i | ||
| 457 | * \param elementI the coupled element of domain à | ||
| 458 | * \param domainJ the domain index of domain j | ||
| 459 | * | ||
| 460 | * \note The element residual definition depends on the discretization scheme of domain i | ||
| 461 | * box: a container of the residuals of all sub control volumes | ||
| 462 | * cc : the residual of the (sub) control volume | ||
| 463 | * fem: the residual of the element | ||
| 464 | * \note This function has to be implemented by all coupling managers for all combinations of i and j | ||
| 465 | */ | ||
| 466 | 638546 | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 467 | const Element<freeFlowMassIndex>& elementI, | ||
| 468 | Dune::index_constant<freeFlowMomentumIndex> domainJ) const | ||
| 469 | { | ||
| 470 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 638546 times.
|
638546 | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI); |
| 471 | 638546 | return massAndEnergyToMomentumStencils_[eIdx]; | |
| 472 | } | ||
| 473 | |||
| 474 | /*! | ||
| 475 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 476 | * that couple with / influence the element residual of the given element of domain i | ||
| 477 | * | ||
| 478 | * \param domainI the domain index of domain i | ||
| 479 | * \param elementI the coupled element of domain à | ||
| 480 | * \param domainJ the domain index of domain j | ||
| 481 | */ | ||
| 482 | 638546 | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 483 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 484 | Dune::index_constant<freeFlowMassIndex> domainJ) const | ||
| 485 | { | ||
| 486 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 638546 times.
|
638546 | const auto eIdx = this->problem(freeFlowMomentumIndex).gridGeometry().elementMapper().index(elementI); |
| 487 | 638546 | return momentumToMassAndEnergyStencils_[eIdx]; | |
| 488 | } | ||
| 489 | |||
| 490 | // \} | ||
| 491 | |||
| 492 | /*! | ||
| 493 | * \name member functions concerning variable caching for element residual evaluations | ||
| 494 | */ | ||
| 495 | // \{ | ||
| 496 | |||
| 497 | /*! | ||
| 498 | * \ingroup MultiDomain | ||
| 499 | * \brief updates all data and variables that are necessary to evaluate the residual of the element of domain i | ||
| 500 | * this is called whenever one of the primary variables that the element residual depends on changes in domain j | ||
| 501 | * | ||
| 502 | * \param domainI the domain index of domain i | ||
| 503 | * \param localAssemblerI the local assembler assembling the element residual of an element of domain i | ||
| 504 | * \param domainJ the domain index of domain j | ||
| 505 | * \param dofIdxGlobalJ the index of the degree of freedom of domain j whose solution changed | ||
| 506 | * \param priVarsJ the new solution at the degree of freedom of domain j with index dofIdxGlobalJ | ||
| 507 | * \param pvIdxJ the index of the primary variable of domain j which has been updated | ||
| 508 | * | ||
| 509 | * \note this concerns all data that is used in the evaluation of the element residual and depends on | ||
| 510 | * the primary variables at the degree of freedom location with index dofIdxGlobalJ | ||
| 511 | * \note the element whose residual is to be evaluated can be retrieved from the local assembler | ||
| 512 | * as localAssemblerI.element() | ||
| 513 | * \note per default, we update the solution vector, if the element residual of domain i depends on more than | ||
| 514 | * the primary variables of domain j update the other dependent data here by overloading this function | ||
| 515 | */ | ||
| 516 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 517 | 19256512 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 518 | const LocalAssemblerI& localAssemblerI, | ||
| 519 | Dune::index_constant<j> domainJ, | ||
| 520 | std::size_t dofIdxGlobalJ, | ||
| 521 | const PrimaryVariables<j>& priVarsJ, | ||
| 522 | int pvIdxJ) | ||
| 523 | { | ||
| 524 | 19256512 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; | |
| 525 | |||
| 526 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 527 | { | ||
| 528 | if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex) | ||
| 529 | { | ||
| 530 | 374008 | bindCouplingContext_(domainI, localAssemblerI.element()); | |
| 531 | |||
| 532 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 374008 times.
|
374008 | const auto& problem = this->problem(domainJ); |
| 533 | 374008 | const auto& deflectedElement = problem.gridGeometry().element(dofIdxGlobalJ); | |
| 534 | 374008 | const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), problem.gridGeometry()); | |
| 535 | 374008 | const auto& fvGeometry = momentumCouplingContext_()[0].fvGeometry; | |
| 536 | 374008 | const auto& scv = fvGeometry.scv(dofIdxGlobalJ); | |
| 537 | |||
| 538 | if constexpr (ElementVolumeVariables<freeFlowMassIndex>::GridVolumeVariables::cachingEnabled) | ||
| 539 | 374008 | gridVars_(freeFlowMassIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), problem, deflectedElement, scv); | |
| 540 | else | ||
| 541 | momentumCouplingContext_()[0].curElemVolVars[scv].update(std::move(elemSol), problem, deflectedElement, scv); | ||
| 542 | } | ||
| 543 | } | ||
| 544 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 545 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 546 | { | ||
| 547 | if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex) | ||
| 548 | { | ||
| 549 | 1865872 | bindCouplingContext_(domainI, localAssemblerI.element()); | |
| 550 | |||
| 551 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1865872 times.
|
1865872 | const auto& problem = this->problem(domainJ); |
| 552 | 1865872 | const auto& deflectedElement = problem.gridGeometry().element(this->momentumCouplingContext_()[0].eIdx); | |
| 553 | 1865872 | const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), problem.gridGeometry()); | |
| 554 | 1865872 | const auto& fvGeometry = this->momentumCouplingContext_()[0].fvGeometry; | |
| 555 | |||
| 556 |
4/4✓ Branch 0 taken 1865872 times.
✓ Branch 1 taken 4430688 times.
✓ Branch 2 taken 6296560 times.
✓ Branch 3 taken 1865872 times.
|
8162432 | for (const auto& scv : scvs(fvGeometry)) |
| 557 | { | ||
| 558 |
2/2✓ Branch 0 taken 1865872 times.
✓ Branch 1 taken 4430688 times.
|
6296560 | if(scv.dofIndex() == dofIdxGlobalJ) |
| 559 | { | ||
| 560 | if constexpr (ElementVolumeVariables<freeFlowMassIndex>::GridVolumeVariables::cachingEnabled) | ||
| 561 |
2/4✓ Branch 1 taken 124800 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124800 times.
✗ Branch 5 not taken.
|
8162432 | this->gridVars_(freeFlowMassIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), problem, deflectedElement, scv); |
| 562 | else | ||
| 563 | this->momentumCouplingContext_()[0].curElemVolVars[scv].update(std::move(elemSol), problem, deflectedElement, scv); | ||
| 564 | } | ||
| 565 | } | ||
| 566 | 124800 | } | |
| 567 | } | ||
| 568 | else | ||
| 569 | DUNE_THROW(Dune::NotImplemented, | ||
| 570 | "Context update for discretization scheme " << MassDiscretizationMethod{} | ||
| 571 | ); | ||
| 572 | 2239880 | } | |
| 573 | |||
| 574 | // \} | ||
| 575 | |||
| 576 | /*! | ||
| 577 | * \brief Compute colors for multithreaded assembly | ||
| 578 | */ | ||
| 579 | void computeColorsForAssembly() | ||
| 580 | { | ||
| 581 | if constexpr (MomentumDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 582 | { | ||
| 583 | // use coloring of the mass discretization for both domains | ||
| 584 | // the diamond coloring is a subset (minimum amount of colors) of cctpfa/box coloring | ||
| 585 | elementSets_ = computeColoring(this->problem(freeFlowMassIndex).gridGeometry()).sets; | ||
| 586 | } | ||
| 587 | else | ||
| 588 | { | ||
| 589 | // use coloring of the momentum discretization for both domains | ||
| 590 | elementSets_ = computeColoring(this->problem(freeFlowMomentumIndex).gridGeometry()).sets; | ||
| 591 | } | ||
| 592 | } | ||
| 593 | |||
| 594 | /*! | ||
| 595 | * \brief Execute assembly kernel in parallel | ||
| 596 | * | ||
| 597 | * \param domainId the domain index of domain i | ||
| 598 | * \param assembleElement kernel function to execute for one element | ||
| 599 | */ | ||
| 600 | template<std::size_t i, class AssembleElementFunc> | ||
| 601 | void assembleMultithreaded(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const | ||
| 602 | { | ||
| 603 | if (elementSets_.empty()) | ||
| 604 | DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!"); | ||
| 605 | |||
| 606 | // make this element loop run in parallel | ||
| 607 | // for this we have to color the elements so that we don't get | ||
| 608 | // race conditions when writing into the global matrix | ||
| 609 | // each color can be assembled using multiple threads | ||
| 610 | const auto& grid = this->problem(freeFlowMassIndex).gridGeometry().gridView().grid(); | ||
| 611 | for (const auto& elements : elementSets_) | ||
| 612 | { | ||
| 613 | Dumux::parallelFor(elements.size(), [&](const std::size_t eIdx) | ||
| 614 | { | ||
| 615 | const auto element = grid.entity(elements[eIdx]); | ||
| 616 | assembleElement(element); | ||
| 617 | }); | ||
| 618 | } | ||
| 619 | } | ||
| 620 | |||
| 621 | private: | ||
| 622 | 2239880 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 623 | const Element<freeFlowMomentumIndex>& elementI) const | ||
| 624 | { | ||
| 625 | // The call to this->problem() is expensive because of std::weak_ptr (see base class). Here we try to avoid it if possible. | ||
| 626 | 2239880 | if (momentumCouplingContext_().empty()) | |
| 627 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | bindCouplingContext_(domainI, elementI, this->problem(freeFlowMomentumIndex).gridGeometry().elementMapper().index(elementI)); |
| 628 | else | ||
| 629 | 2239871 | bindCouplingContext_(domainI, elementI, momentumCouplingContext_()[0].fvGeometry.gridGeometry().elementMapper().index(elementI)); | |
| 630 | 2239880 | } | |
| 631 | |||
| 632 | 6217622 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 633 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 634 | const std::size_t eIdx) const | ||
| 635 | { | ||
| 636 | 6217622 | if (momentumCouplingContext_().empty()) | |
| 637 | { | ||
| 638 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
|
12 | auto fvGeometry = localView(this->problem(freeFlowMassIndex).gridGeometry()); |
| 639 | 12 | fvGeometry.bind(elementI); | |
| 640 | |||
| 641 |
3/5✓ Branch 1 taken 2 times.
✓ Branch 2 taken 5 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
|
12 | auto curElemVolVars = localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 642 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
12 | curElemVolVars.bind(elementI, fvGeometry, this->curSol(freeFlowMassIndex)); |
| 643 | |||
| 644 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
17 | auto prevElemVolVars = isTransient_ ? localView(gridVars_(freeFlowMassIndex).prevGridVolVars()) |
| 645 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
12 | : localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 646 | |||
| 647 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | if (isTransient_) |
| 648 | ✗ | prevElemVolVars.bindElement(elementI, fvGeometry, (*prevSol_)[freeFlowMassIndex]); | |
| 649 | |||
| 650 |
2/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12 | momentumCouplingContext_().emplace_back(MomentumCouplingContext{std::move(fvGeometry), std::move(curElemVolVars), std::move(prevElemVolVars), eIdx}); |
| 651 | 12 | } | |
| 652 | 6217610 | else if (eIdx != momentumCouplingContext_()[0].eIdx) | |
| 653 | { | ||
| 654 | 468848 | momentumCouplingContext_()[0].eIdx = eIdx; | |
| 655 | 937696 | momentumCouplingContext_()[0].fvGeometry.bind(elementI); | |
| 656 | 750697 | momentumCouplingContext_()[0].curElemVolVars.bind(elementI, momentumCouplingContext_()[0].fvGeometry, this->curSol(freeFlowMassIndex)); | |
| 657 | |||
| 658 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 281849 times.
|
281849 | if (isTransient_) |
| 659 | ✗ | momentumCouplingContext_()[0].prevElemVolVars.bindElement(elementI, momentumCouplingContext_()[0].fvGeometry, (*prevSol_)[freeFlowMassIndex]); | |
| 660 | } | ||
| 661 | 6217629 | } | |
| 662 | |||
| 663 | 332172 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 664 | const Element<freeFlowMassIndex>& elementI) const | ||
| 665 | { | ||
| 666 | // The call to this->problem() is expensive because of std::weak_ptr (see base class). Here we try to avoid it if possible. | ||
| 667 | 332172 | if (massAndEnergyCouplingContext_().empty()) | |
| 668 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | bindCouplingContext_(domainI, elementI, this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI)); |
| 669 | else | ||
| 670 | 332161 | bindCouplingContext_(domainI, elementI, massAndEnergyCouplingContext_()[0].fvGeometry.gridGeometry().elementMapper().index(elementI)); | |
| 671 | 332172 | } | |
| 672 | |||
| 673 | 27093716 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 674 | const Element<freeFlowMassIndex>& elementI, | ||
| 675 | const std::size_t eIdx) const | ||
| 676 | { | ||
| 677 | 27093716 | if (massAndEnergyCouplingContext_().empty()) | |
| 678 | { | ||
| 679 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
|
12 | const auto& gridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 680 | 12 | auto fvGeometry = localView(gridGeometry); | |
| 681 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
12 | fvGeometry.bindElement(elementI); |
| 682 | 12 | massAndEnergyCouplingContext_().emplace_back(std::move(fvGeometry), eIdx); | |
| 683 | 2 | } | |
| 684 | 27093704 | else if (eIdx != massAndEnergyCouplingContext_()[0].eIdx) | |
| 685 | { | ||
| 686 | 2776583 | massAndEnergyCouplingContext_()[0].eIdx = eIdx; | |
| 687 | 2776583 | massAndEnergyCouplingContext_()[0].fvGeometry.bindElement(elementI); | |
| 688 | } | ||
| 689 | 27093716 | } | |
| 690 | |||
| 691 | /*! | ||
| 692 | * \brief Return a reference to the grid variables of a sub problem | ||
| 693 | * \param domainIdx The domain index | ||
| 694 | */ | ||
| 695 | template<std::size_t i> | ||
| 696 | 24 | const GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) const | |
| 697 | { | ||
| 698 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | if (std::get<i>(gridVariables_)) |
| 699 | 24 | return *std::get<i>(gridVariables_); | |
| 700 | else | ||
| 701 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 702 | } | ||
| 703 | |||
| 704 | /*! | ||
| 705 | * \brief Return a reference to the grid variables of a sub problem | ||
| 706 | * \param domainIdx The domain index | ||
| 707 | */ | ||
| 708 | template<std::size_t i> | ||
| 709 | 2239880 | GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) | |
| 710 | { | ||
| 711 |
1/2✓ Branch 0 taken 2239880 times.
✗ Branch 1 not taken.
|
2239880 | if (std::get<i>(gridVariables_)) |
| 712 | 2239880 | return *std::get<i>(gridVariables_); | |
| 713 | else | ||
| 714 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 715 | } | ||
| 716 | |||
| 717 | |||
| 718 | 12 | void computeCouplingStencils_() | |
| 719 | { | ||
| 720 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
|
12 | const auto& momentumGridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 721 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
|
12 | const auto& massGridGeometry = this->problem(freeFlowMassIndex).gridGeometry(); |
| 722 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
12 | auto momentumFvGeometry = localView(momentumGridGeometry); |
| 723 | 12 | auto massFvGeometry = localView(massGridGeometry); | |
| 724 | |||
| 725 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | massAndEnergyToMomentumStencils_.clear(); |
| 726 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
12 | massAndEnergyToMomentumStencils_.resize(massGridGeometry.gridView().size(0)); |
| 727 | |||
| 728 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | momentumToMassAndEnergyStencils_.clear(); |
| 729 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
12 | momentumToMassAndEnergyStencils_.resize(momentumGridGeometry.gridView().size(0)); |
| 730 | |||
| 731 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | assert(massAndEnergyToMomentumStencils_.size() == momentumToMassAndEnergyStencils_.size()); |
| 732 | |||
| 733 |
9/13✓ Branch 1 taken 2 times.
✓ Branch 2 taken 159287 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10400 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 10400 times.
✓ Branch 13 taken 2 times.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 3 taken 9 times.
|
340335 | for (const auto& element : elements(momentumGridGeometry.gridView())) |
| 734 | { | ||
| 735 |
1/2✓ Branch 1 taken 10400 times.
✗ Branch 2 not taken.
|
169686 | momentumFvGeometry.bindElement(element); |
| 736 |
1/2✓ Branch 1 taken 10400 times.
✗ Branch 2 not taken.
|
169686 | massFvGeometry.bindElement(element); |
| 737 | 169686 | const auto eIdx = momentumFvGeometry.elementIndex(); | |
| 738 | |||
| 739 |
2/2✓ Branch 0 taken 663748 times.
✓ Branch 1 taken 169686 times.
|
833434 | for (const auto& scv : scvs(momentumFvGeometry)) |
| 740 |
1/2✓ Branch 1 taken 38400 times.
✗ Branch 2 not taken.
|
663748 | massAndEnergyToMomentumStencils_[eIdx].push_back(scv.dofIndex()); |
| 741 | |||
| 742 |
2/2✓ Branch 0 taken 402470 times.
✓ Branch 1 taken 169686 times.
|
572156 | for (const auto& scv : scvs(massFvGeometry)) |
| 743 |
1/2✓ Branch 1 taken 31200 times.
✗ Branch 2 not taken.
|
402470 | momentumToMassAndEnergyStencils_[eIdx].push_back(scv.dofIndex()); |
| 744 | } | ||
| 745 | 14 | } | |
| 746 | |||
| 747 | CouplingStencilType emptyStencil_; | ||
| 748 | std::vector<CouplingStencilType> momentumToMassAndEnergyStencils_; | ||
| 749 | std::vector<CouplingStencilType> massAndEnergyToMomentumStencils_; | ||
| 750 | |||
| 751 | 12 | std::vector<MomentumCouplingContext>& momentumCouplingContext_() const | |
| 752 |
18/22✓ Branch 0 taken 9 times.
✓ Branch 1 taken 2822321 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 537600 times.
✓ Branch 4 taken 2940801 times.
✓ Branch 5 taken 1670590 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 1443370 times.
✓ Branch 12 taken 41228 times.
✓ Branch 13 taken 1158595 times.
✓ Branch 6 taken 545599 times.
✓ Branch 7 taken 2481180 times.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 356455 times.
✓ Branch 10 taken 1741068 times.
✓ Branch 14 taken 5 times.
✗ Branch 11 not taken.
✗ Branch 15 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 5658 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
|
28174540 | { return momentumCouplingContextImpl_; } |
| 753 | |||
| 754 | 12 | std::vector<MassAndEnergyCouplingContext>& massAndEnergyCouplingContext_() const | |
| 755 |
10/12✓ Branch 0 taken 11 times.
✓ Branch 1 taken 332161 times.
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 27093704 times.
✓ Branch 7 taken 2742187 times.
✓ Branch 8 taken 23372161 times.
✓ Branch 9 taken 34398 times.
✓ Branch 10 taken 944960 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 10 times.
|
54519628 | { return massAndEnergyCouplingContextImpl_; } |
| 756 | |||
| 757 | mutable std::vector<MassAndEnergyCouplingContext> massAndEnergyCouplingContextImpl_; | ||
| 758 | mutable std::vector<MomentumCouplingContext> momentumCouplingContextImpl_; | ||
| 759 | |||
| 760 | //! A tuple of std::shared_ptrs to the grid variables of the sub problems | ||
| 761 | GridVariablesTuple gridVariables_; | ||
| 762 | |||
| 763 | const SolutionVector* prevSol_; | ||
| 764 | bool isTransient_; | ||
| 765 | |||
| 766 | std::deque<std::vector<ElementSeed<freeFlowMomentumIndex>>> elementSets_; | ||
| 767 | }; | ||
| 768 | |||
| 769 | namespace Detail { | ||
| 770 | |||
| 771 | // declaration (specialize for different discretization types) | ||
| 772 | template<class Traits, class DiscretizationMethod = typename Detail::MomentumDiscretizationMethod<Traits>::type> | ||
| 773 | struct CouplingManagerSupportsMultithreadedAssemblySelector; | ||
| 774 | |||
| 775 | // multi-threading is not supported because we have only one coupling context instance and a mutable cache | ||
| 776 | template<class Traits, class D> | ||
| 777 | struct CouplingManagerSupportsMultithreadedAssemblySelector<Traits, DiscretizationMethods::CVFE<D>> | ||
| 778 | { using type = std::false_type; }; | ||
| 779 | |||
| 780 | } // end namespace Detail | ||
| 781 | |||
| 782 | //! whether we support multithreaded assembly | ||
| 783 | template<class T> | ||
| 784 | struct CouplingManagerSupportsMultithreadedAssembly<CVFEFreeFlowCouplingManager<T>> | ||
| 785 | : public Detail::CouplingManagerSupportsMultithreadedAssemblySelector<T>::type | ||
| 786 | {}; | ||
| 787 | |||
| 788 | } // end namespace Dumux | ||
| 789 | |||
| 790 | #endif | ||
| 791 |