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