| 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-FileCopyrightInfo: 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 |
2/14✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
|
20 | struct MomentumCouplingContext |
| 90 | { | ||
| 91 | FVElementGeometry<freeFlowMassIndex> fvGeometry; | ||
| 92 | ElementVolumeVariables<freeFlowMassIndex> curElemVolVars; | ||
| 93 | ElementVolumeVariables<freeFlowMassIndex> prevElemVolVars; | ||
| 94 | std::size_t eIdx; | ||
| 95 | }; | ||
| 96 | |||
| 97 |
1/12✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
|
4 | struct MassAndEnergyCouplingContext |
| 98 | { | ||
| 99 | 9 | MassAndEnergyCouplingContext(FVElementGeometry<freeFlowMomentumIndex>&& f, const std::size_t i) | |
| 100 | 9 | : fvGeometry(std::move(f)) | |
| 101 |
1/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
16 | , 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 | 9 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | |
| 122 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 123 | GridVariablesTuple&& gridVariables, | ||
| 124 | const SolutionVector& curSol) | ||
| 125 | { | ||
| 126 | 18 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | |
| 127 | 9 | gridVariables_ = gridVariables; | |
| 128 | 9 | this->updateSolution(curSol); | |
| 129 | |||
| 130 | 9 | computeCouplingStencils_(); | |
| 131 | 9 | } | |
| 132 | |||
| 133 | //! use as regular coupling manager in a transient setting | ||
| 134 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | ||
| 135 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 136 | GridVariablesTuple&& gridVariables, | ||
| 137 | const SolutionVector& curSol, | ||
| 138 | const SolutionVector& prevSol) | ||
| 139 | { | ||
| 140 | init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol); | ||
| 141 | prevSol_ = &prevSol; | ||
| 142 | isTransient_ = true; | ||
| 143 | } | ||
| 144 | |||
| 145 | //! use as binary coupling manager in multi model context | ||
| 146 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | ||
| 147 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 148 | GridVariablesTuple&& gridVariables, | ||
| 149 | typename ParentType::SolutionVectorStorage& curSol) | ||
| 150 | { | ||
| 151 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | ||
| 152 | gridVariables_ = gridVariables; | ||
| 153 | this->attachSolution(curSol); | ||
| 154 | |||
| 155 | computeCouplingStencils_(); | ||
| 156 | } | ||
| 157 | |||
| 158 | // \} | ||
| 159 | |||
| 160 | /*! | ||
| 161 | * \name member functions concerning the coupling stencils | ||
| 162 | */ | ||
| 163 | // \{ | ||
| 164 | |||
| 165 | /*! | ||
| 166 | * \brief Returns the pressure at a given sub control volume face | ||
| 167 | */ | ||
| 168 | 36226506 | Scalar pressure(const Element<freeFlowMomentumIndex>& element, | |
| 169 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 170 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 171 | const bool considerPreviousTimeStep = false) const | ||
| 172 | { | ||
| 173 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 36226506 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
36226506 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 174 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 36226506 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 36226506 times.
|
36226506 | const auto& gg = this->problem(freeFlowMassIndex).gridGeometry(); |
| 175 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 36226506 times.
|
36226506 | const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex] |
| 176 | : this->curSol(freeFlowMassIndex); | ||
| 177 | 36226506 | const auto elemSol = elementSolution(element, sol, gg); | |
| 178 |
1/2✓ Branch 3 taken 1545600 times.
✗ Branch 4 not taken.
|
113519918 | return evalSolution(element, element.geometry(), gg, elemSol, scvf.ipGlobal())[pressureIdx]; |
| 179 | } | ||
| 180 | |||
| 181 | /*! | ||
| 182 | * \brief Returns the pressure at a given sub control volume | ||
| 183 | */ | ||
| 184 | Scalar pressure(const Element<freeFlowMomentumIndex>& element, | ||
| 185 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 186 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 187 | const bool considerPreviousTimeStep = false) const | ||
| 188 | { | ||
| 189 | assert(!(considerPreviousTimeStep && !this->isTransient_)); | ||
| 190 | const auto& gg = this->problem(freeFlowMassIndex).gridGeometry(); | ||
| 191 | const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex] | ||
| 192 | : this->curSol(freeFlowMassIndex); | ||
| 193 | const auto elemSol = elementSolution(element, sol, gg); | ||
| 194 | return evalSolution(element, element.geometry(), gg, elemSol, scv.dofPosition())[pressureIdx]; | ||
| 195 | } | ||
| 196 | |||
| 197 | /*! | ||
| 198 | * \brief Returns the density at a given sub control volume face. | ||
| 199 | */ | ||
| 200 | ✗ | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 201 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 202 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 203 | const bool considerPreviousTimeStep = false) const | ||
| 204 | { | ||
| 205 | ✗ | assert(!(considerPreviousTimeStep && !this->isTransient_)); | |
| 206 | ✗ | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 207 | |||
| 208 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 209 | { | ||
| 210 | const auto eIdx = fvGeometry.elementIndex(); | ||
| 211 | const auto& scv = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | ||
| 212 | |||
| 213 | const auto& volVars = considerPreviousTimeStep ? | ||
| 214 | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | ||
| 215 | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | ||
| 216 | |||
| 217 | return volVars.density(); | ||
| 218 | } | ||
| 219 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 220 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 221 | { | ||
| 222 | // TODO: cache the shape values when Box method is used | ||
| 223 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 224 | ✗ | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 225 | ✗ | std::vector<ShapeValue> shapeValues; | |
| 226 | ✗ | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); | |
| 227 | ✗ | localBasis.evaluateFunction(ipLocal, shapeValues); | |
| 228 | |||
| 229 | ✗ | Scalar rho = 0.0; | |
| 230 | ✗ | for (const auto& scv : scvs(this->momentumCouplingContext_()[0].fvGeometry)) | |
| 231 | { | ||
| 232 | ✗ | const auto& volVars = considerPreviousTimeStep ? | |
| 233 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 234 | ✗ | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 235 | ✗ | rho += volVars.density()*shapeValues[scv.indexInElement()][0]; | |
| 236 | } | ||
| 237 | |||
| 238 | ✗ | return rho; | |
| 239 | } | ||
| 240 | else | ||
| 241 | DUNE_THROW(Dune::NotImplemented, | ||
| 242 | "Density interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 243 | ); | ||
| 244 | } | ||
| 245 | |||
| 246 | /*! | ||
| 247 | * \brief Returns the density at a given sub control volume. | ||
| 248 | */ | ||
| 249 | 1113600 | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 250 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 251 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 252 | const bool considerPreviousTimeStep = false) const | ||
| 253 | { | ||
| 254 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1113600 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1113600 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 255 | 1113600 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, scv.elementIndex()); | |
| 256 | |||
| 257 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 258 | { | ||
| 259 | const auto eIdx = scv.elementIndex(); | ||
| 260 | const auto& scvI = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | ||
| 261 | |||
| 262 | const auto& volVars = considerPreviousTimeStep ? | ||
| 263 | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | ||
| 264 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | ||
| 265 | |||
| 266 | return volVars.density(); | ||
| 267 | } | ||
| 268 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 269 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 270 | { | ||
| 271 | // TODO: cache the shape values when Box method is used | ||
| 272 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 273 | 1113600 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 274 |
1/4✓ Branch 1 taken 1113600 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1113600 | std::vector<ShapeValue> shapeValues; |
| 275 |
1/2✓ Branch 1 taken 1113600 times.
✗ Branch 2 not taken.
|
1113600 | const auto ipLocal = element.geometry().local(scv.dofPosition()); |
| 276 |
1/2✓ Branch 1 taken 1113600 times.
✗ Branch 2 not taken.
|
1113600 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 277 | |||
| 278 | 1113600 | Scalar rho = 0.0; | |
| 279 |
4/4✓ Branch 1 taken 3340800 times.
✓ Branch 2 taken 1113600 times.
✓ Branch 3 taken 3340800 times.
✓ Branch 4 taken 1113600 times.
|
7795200 | for (const auto& scvI : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 280 | { | ||
| 281 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3340800 times.
|
3340800 | const auto& volVars = considerPreviousTimeStep ? |
| 282 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | |
| 283 | 3340800 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | |
| 284 | 10022400 | rho += volVars.density()*shapeValues[scvI.indexInElement()][0]; | |
| 285 | } | ||
| 286 |
1/2✓ Branch 0 taken 1113600 times.
✗ Branch 1 not taken.
|
2227200 | return rho; |
| 287 | } | ||
| 288 | else | ||
| 289 | DUNE_THROW(Dune::NotImplemented, | ||
| 290 | "Density interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 291 | ); | ||
| 292 | } | ||
| 293 | |||
| 294 | /*! | ||
| 295 | * \brief Returns the effective viscosity at a given sub control volume face. | ||
| 296 | */ | ||
| 297 | 1550280 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 298 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 299 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 300 | const bool considerPreviousTimeStep = false) const | ||
| 301 | { | ||
| 302 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1550280 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1550280 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 303 | 2850960 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 304 | |||
| 305 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 306 | { | ||
| 307 | const auto eIdx = fvGeometry.elementIndex(); | ||
| 308 | const auto& scv = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | ||
| 309 | const auto& volVars = considerPreviousTimeStep ? | ||
| 310 | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | ||
| 311 | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | ||
| 312 | return volVars.viscosity(); | ||
| 313 | } | ||
| 314 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 315 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 316 | { | ||
| 317 | // TODO: cache the shape values when Box method is used | ||
| 318 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 319 | 1550280 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 320 |
1/4✓ Branch 1 taken 1550280 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1550280 | std::vector<ShapeValue> shapeValues; |
| 321 |
1/2✓ Branch 1 taken 1550280 times.
✗ Branch 2 not taken.
|
1550280 | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); |
| 322 |
1/2✓ Branch 1 taken 1550280 times.
✗ Branch 2 not taken.
|
1550280 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 323 | |||
| 324 | 1550280 | Scalar mu = 0.0; | |
| 325 |
4/4✓ Branch 1 taken 4650840 times.
✓ Branch 2 taken 1550280 times.
✓ Branch 3 taken 4650840 times.
✓ Branch 4 taken 1550280 times.
|
10851960 | for (const auto& scv : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 326 | { | ||
| 327 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4650840 times.
|
4650840 | const auto& volVars = considerPreviousTimeStep ? |
| 328 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scv] | |
| 329 | 4650840 | : this->momentumCouplingContext_()[0].curElemVolVars[scv]; | |
| 330 | 13952520 | mu += volVars.viscosity()*shapeValues[scv.indexInElement()][0]; | |
| 331 | } | ||
| 332 | |||
| 333 |
1/2✓ Branch 0 taken 1550280 times.
✗ Branch 1 not taken.
|
3100560 | return mu; |
| 334 | } | ||
| 335 | else | ||
| 336 | DUNE_THROW(Dune::NotImplemented, | ||
| 337 | "Viscosity interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 338 | ); | ||
| 339 | } | ||
| 340 | |||
| 341 | /*! | ||
| 342 | * \brief Returns the effective viscosity at a given sub control volume. | ||
| 343 | */ | ||
| 344 | 244500 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 345 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 346 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 347 | const bool considerPreviousTimeStep = false) const | ||
| 348 | { | ||
| 349 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
244500 | assert(!(considerPreviousTimeStep && !this->isTransient_)); |
| 350 | 489000 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 351 | |||
| 352 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 353 | { | ||
| 354 | const auto eIdx = fvGeometry.elementIndex(); | ||
| 355 | const auto& scvI = this->momentumCouplingContext_()[0].fvGeometry.scv(eIdx); | ||
| 356 | const auto& volVars = considerPreviousTimeStep ? | ||
| 357 | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | ||
| 358 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | ||
| 359 | return volVars.viscosity(); | ||
| 360 | } | ||
| 361 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 362 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 363 | { | ||
| 364 | // TODO: cache the shape values when Box method is used | ||
| 365 | using ShapeValue = typename Dune::FieldVector<Scalar, 1>; | ||
| 366 | 244500 | const auto& localBasis = this->momentumCouplingContext_()[0].fvGeometry.feLocalBasis(); | |
| 367 |
1/4✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
244500 | std::vector<ShapeValue> shapeValues; |
| 368 |
1/2✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
|
244500 | const auto ipLocal = element.geometry().local(scv.dofPosition()); |
| 369 |
1/2✓ Branch 1 taken 244500 times.
✗ Branch 2 not taken.
|
244500 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 370 | |||
| 371 | 244500 | Scalar mu = 0.0; | |
| 372 |
4/4✓ Branch 1 taken 733500 times.
✓ Branch 2 taken 244500 times.
✓ Branch 3 taken 733500 times.
✓ Branch 4 taken 244500 times.
|
1711500 | for (const auto& scvI : scvs(this->momentumCouplingContext_()[0].fvGeometry)) |
| 373 | { | ||
| 374 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 733500 times.
|
733500 | const auto& volVars = considerPreviousTimeStep ? |
| 375 | ✗ | this->momentumCouplingContext_()[0].prevElemVolVars[scvI] | |
| 376 | 733500 | : this->momentumCouplingContext_()[0].curElemVolVars[scvI]; | |
| 377 | 2200500 | mu += volVars.viscosity()*shapeValues[scvI.indexInElement()][0]; | |
| 378 | } | ||
| 379 | |||
| 380 |
1/2✓ Branch 0 taken 244500 times.
✗ Branch 1 not taken.
|
489000 | return mu; |
| 381 | } | ||
| 382 | else | ||
| 383 | DUNE_THROW(Dune::NotImplemented, | ||
| 384 | "Viscosity interpolation for discretization scheme " << MassDiscretizationMethod{} | ||
| 385 | ); | ||
| 386 | } | ||
| 387 | |||
| 388 | /*! | ||
| 389 | * \brief Returns the velocity at a given sub control volume face. | ||
| 390 | */ | ||
| 391 | 23581296 | VelocityVector faceVelocity(const Element<freeFlowMassIndex>& element, | |
| 392 | const SubControlVolumeFace<freeFlowMassIndex>& scvf) const | ||
| 393 | { | ||
| 394 | // TODO: optimize this function for tpfa where the scvf ip coincides with the dof location | ||
| 395 | |||
| 396 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23581296 times.
|
23581296 | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(element); |
| 397 | 23581296 | bindCouplingContext_(Dune::index_constant<freeFlowMassIndex>(), element, eIdx); | |
| 398 | |||
| 399 | 23581296 | const auto& fvGeometry = this->massAndEnergyCouplingContext_()[0].fvGeometry; | |
| 400 | 23581296 | const auto& localBasis = fvGeometry.feLocalBasis(); | |
| 401 | |||
| 402 |
2/4✓ Branch 1 taken 23581296 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1991108 times.
✗ Branch 4 not taken.
|
25572404 | std::vector<ShapeValue> shapeValues; |
| 403 |
3/8✓ Branch 1 taken 23581296 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22615536 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22615536 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
68812368 | const auto ipLocal = element.geometry().local(scvf.ipGlobal()); |
| 404 |
1/2✓ Branch 1 taken 23581296 times.
✗ Branch 2 not taken.
|
23581296 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 405 | |||
| 406 | // interpolate velocity at scvf | ||
| 407 | 23581296 | VelocityVector velocity(0.0); | |
| 408 |
4/4✓ Branch 0 taken 97699920 times.
✓ Branch 1 taken 23581296 times.
✓ Branch 2 taken 97699920 times.
✓ Branch 3 taken 23581296 times.
|
242562432 | for (const auto& scv : scvs(fvGeometry)) |
| 409 | 586199520 | velocity.axpy(shapeValues[scv.localDofIndex()][0], this->curSol(freeFlowMomentumIndex)[scv.dofIndex()]); | |
| 410 | |||
| 411 |
1/2✓ Branch 0 taken 21590188 times.
✗ Branch 1 not taken.
|
45171484 | return velocity; |
| 412 | } | ||
| 413 | |||
| 414 | /*! | ||
| 415 | * \brief Returns the velocity at the element center. | ||
| 416 | */ | ||
| 417 | 214494 | VelocityVector elementVelocity(const FVElementGeometry<freeFlowMassIndex>& fvGeometry) const | |
| 418 | { | ||
| 419 | 428988 | bindCouplingContext_(Dune::index_constant<freeFlowMassIndex>(), fvGeometry.element()); | |
| 420 | |||
| 421 | 214494 | const auto& momentumFvGeometry = this->massAndEnergyCouplingContext_()[0].fvGeometry; | |
| 422 | 214494 | const auto& localBasis = momentumFvGeometry.feLocalBasis(); | |
| 423 | |||
| 424 | // interpolate velocity at scvf | ||
| 425 | 214494 | VelocityVector velocity(0.0); | |
| 426 |
2/4✓ Branch 1 taken 214494 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 29798 times.
✗ Branch 4 not taken.
|
244292 | std::vector<ShapeValue> shapeValues; |
| 427 |
4/10✓ Branch 1 taken 214494 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 214494 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 214494 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 214494 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
428988 | localBasis.evaluateFunction(referenceElement(fvGeometry.element()).position(0,0), shapeValues); |
| 428 | |||
| 429 |
4/4✓ Branch 0 taken 882108 times.
✓ Branch 1 taken 214494 times.
✓ Branch 2 taken 882108 times.
✓ Branch 3 taken 214494 times.
|
2193204 | for (const auto& scv : scvs(momentumFvGeometry)) |
| 430 | 5292648 | velocity.axpy(shapeValues[scv.localDofIndex()][0], this->curSol(freeFlowMomentumIndex)[scv.dofIndex()]); | |
| 431 | |||
| 432 |
1/2✓ Branch 0 taken 184696 times.
✗ Branch 1 not taken.
|
399190 | return velocity; |
| 433 | } | ||
| 434 | |||
| 435 | /*! | ||
| 436 | * \brief The coupling stencil of domain I, i.e. which domain J DOFs | ||
| 437 | * the given domain I element's residual depends on. | ||
| 438 | */ | ||
| 439 | template<std::size_t j> | ||
| 440 | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | ||
| 441 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 442 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 443 | Dune::index_constant<j> domainJ) const | ||
| 444 | { return emptyStencil_; } | ||
| 445 | |||
| 446 | /*! | ||
| 447 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 448 | * that couple with / influence the element residual of the given element of domain i | ||
| 449 | * | ||
| 450 | * \param domainI the domain index of domain i | ||
| 451 | * \param elementI the coupled element of domain à | ||
| 452 | * \param domainJ the domain index of domain j | ||
| 453 | * | ||
| 454 | * \note The element residual definition depends on the discretization scheme of domain i | ||
| 455 | * box: a container of the residuals of all sub control volumes | ||
| 456 | * cc : the residual of the (sub) control volume | ||
| 457 | * fem: the residual of the element | ||
| 458 | * \note This function has to be implemented by all coupling managers for all combinations of i and j | ||
| 459 | */ | ||
| 460 | ✗ | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 461 | const Element<freeFlowMassIndex>& elementI, | ||
| 462 | Dune::index_constant<freeFlowMomentumIndex> domainJ) const | ||
| 463 | { | ||
| 464 | ✗ | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI); | |
| 465 | ✗ | return massAndEnergyToMomentumStencils_[eIdx]; | |
| 466 | } | ||
| 467 | |||
| 468 | /*! | ||
| 469 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 470 | * that couple with / influence the element residual of the given element of domain i | ||
| 471 | * | ||
| 472 | * \param domainI the domain index of domain i | ||
| 473 | * \param elementI the coupled element of domain à | ||
| 474 | * \param domainJ the domain index of domain j | ||
| 475 | */ | ||
| 476 | ✗ | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 477 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 478 | Dune::index_constant<freeFlowMassIndex> domainJ) const | ||
| 479 | { | ||
| 480 | ✗ | const auto eIdx = this->problem(freeFlowMomentumIndex).gridGeometry().elementMapper().index(elementI); | |
| 481 | ✗ | return momentumToMassAndEnergyStencils_[eIdx]; | |
| 482 | } | ||
| 483 | |||
| 484 | // \} | ||
| 485 | |||
| 486 | /*! | ||
| 487 | * \name member functions concerning variable caching for element residual evaluations | ||
| 488 | */ | ||
| 489 | // \{ | ||
| 490 | |||
| 491 | /*! | ||
| 492 | * \ingroup MultiDomain | ||
| 493 | * \brief updates all data and variables that are necessary to evaluate the residual of the element of domain i | ||
| 494 | * this is called whenever one of the primary variables that the element residual depends on changes in domain j | ||
| 495 | * | ||
| 496 | * \param domainI the domain index of domain i | ||
| 497 | * \param localAssemblerI the local assembler assembling the element residual of an element of domain i | ||
| 498 | * \param domainJ the domain index of domain j | ||
| 499 | * \param dofIdxGlobalJ the index of the degree of freedom of domain j whose solution changed | ||
| 500 | * \param priVarsJ the new solution at the degree of freedom of domain j with index dofIdxGlobalJ | ||
| 501 | * \param pvIdxJ the index of the primary variable of domain j which has been updated | ||
| 502 | * | ||
| 503 | * \note this concerns all data that is used in the evaluation of the element residual and depends on | ||
| 504 | * the primary variables at the degree of freedom location with index dofIdxGlobalJ | ||
| 505 | * \note the element whose residual is to be evaluated can be retrieved from the local assembler | ||
| 506 | * as localAssemblerI.element() | ||
| 507 | * \note per default, we update the solution vector, if the element residual of domain i depends on more than | ||
| 508 | * the primary variables of domain j update the other dependent data here by overloading this function | ||
| 509 | */ | ||
| 510 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 511 | 4005276 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 512 | const LocalAssemblerI& localAssemblerI, | ||
| 513 | Dune::index_constant<j> domainJ, | ||
| 514 | std::size_t dofIdxGlobalJ, | ||
| 515 | const PrimaryVariables<j>& priVarsJ, | ||
| 516 | int pvIdxJ) | ||
| 517 | { | ||
| 518 | 78778346 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; | |
| 519 | |||
| 520 | if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa) | ||
| 521 | { | ||
| 522 | if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex) | ||
| 523 | { | ||
| 524 | 620908 | bindCouplingContext_(domainI, localAssemblerI.element()); | |
| 525 | |||
| 526 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 310454 times.
|
620908 | const auto& problem = this->problem(domainJ); |
| 527 | 1241816 | const auto& deflectedElement = problem.gridGeometry().element(dofIdxGlobalJ); | |
| 528 | 1862724 | const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), problem.gridGeometry()); | |
| 529 | 620908 | const auto& fvGeometry = momentumCouplingContext_()[0].fvGeometry; | |
| 530 | 620908 | const auto& scv = fvGeometry.scv(dofIdxGlobalJ); | |
| 531 | |||
| 532 | if constexpr (ElementVolumeVariables<freeFlowMassIndex>::GridVolumeVariables::cachingEnabled) | ||
| 533 | 620908 | gridVars_(freeFlowMassIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), problem, deflectedElement, scv); | |
| 534 | else | ||
| 535 | momentumCouplingContext_()[0].curElemVolVars[scv].update(std::move(elemSol), problem, deflectedElement, scv); | ||
| 536 | } | ||
| 537 | } | ||
| 538 | else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box | ||
| 539 | || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 540 | { | ||
| 541 | if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex) | ||
| 542 | { | ||
| 543 | 3384368 | bindCouplingContext_(domainI, localAssemblerI.element()); | |
| 544 | |||
| 545 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1692184 times.
|
3384368 | const auto& problem = this->problem(domainJ); |
| 546 | 7018336 | const auto& deflectedElement = problem.gridGeometry().element(this->momentumCouplingContext_()[0].eIdx); | |
| 547 | 10153104 | const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), problem.gridGeometry()); | |
| 548 | 3384368 | const auto& fvGeometry = this->momentumCouplingContext_()[0].fvGeometry; | |
| 549 | |||
| 550 |
4/4✓ Branch 0 taken 5775496 times.
✓ Branch 1 taken 1692184 times.
✓ Branch 2 taken 5775496 times.
✓ Branch 3 taken 1692184 times.
|
18319728 | for (const auto& scv : scvs(fvGeometry)) |
| 551 | { | ||
| 552 |
2/2✓ Branch 0 taken 1692184 times.
✓ Branch 1 taken 4083312 times.
|
11550992 | if(scv.dofIndex() == dofIdxGlobalJ) |
| 553 | { | ||
| 554 | if constexpr (ElementVolumeVariables<freeFlowMassIndex>::GridVolumeVariables::cachingEnabled) | ||
| 555 |
4/8✓ Branch 1 taken 124800 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 124800 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 124800 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 124800 times.
✗ Branch 11 not taken.
|
3384368 | this->gridVars_(freeFlowMassIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), problem, deflectedElement, scv); |
| 556 | else | ||
| 557 | this->momentumCouplingContext_()[0].curElemVolVars[scv].update(std::move(elemSol), problem, deflectedElement, scv); | ||
| 558 | } | ||
| 559 | } | ||
| 560 | } | ||
| 561 | } | ||
| 562 | else | ||
| 563 | DUNE_THROW(Dune::NotImplemented, | ||
| 564 | "Context update for discretization scheme " << MassDiscretizationMethod{} | ||
| 565 | ); | ||
| 566 | 4005276 | } | |
| 567 | |||
| 568 | // \} | ||
| 569 | |||
| 570 | /*! | ||
| 571 | * \brief Compute colors for multithreaded assembly | ||
| 572 | */ | ||
| 573 | void computeColorsForAssembly() | ||
| 574 | { | ||
| 575 | if constexpr (MomentumDiscretizationMethod{} == DiscretizationMethods::fcdiamond) | ||
| 576 | { | ||
| 577 | // use coloring of the mass discretization for both domains | ||
| 578 | // the diamond coloring is a subset (minimum amount of colors) of cctpfa/box coloring | ||
| 579 | elementSets_ = computeColoring(this->problem(freeFlowMassIndex).gridGeometry()).sets; | ||
| 580 | } | ||
| 581 | else | ||
| 582 | { | ||
| 583 | // use coloring of the momentum discretization for both domains | ||
| 584 | elementSets_ = computeColoring(this->problem(freeFlowMomentumIndex).gridGeometry()).sets; | ||
| 585 | } | ||
| 586 | } | ||
| 587 | |||
| 588 | /*! | ||
| 589 | * \brief Execute assembly kernel in parallel | ||
| 590 | * | ||
| 591 | * \param domainId the domain index of domain i | ||
| 592 | * \param assembleElement kernel function to execute for one element | ||
| 593 | */ | ||
| 594 | template<std::size_t i, class AssembleElementFunc> | ||
| 595 | void assembleMultithreaded(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const | ||
| 596 | { | ||
| 597 | if (elementSets_.empty()) | ||
| 598 | DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!"); | ||
| 599 | |||
| 600 | // make this element loop run in parallel | ||
| 601 | // for this we have to color the elements so that we don't get | ||
| 602 | // race conditions when writing into the global matrix | ||
| 603 | // each color can be assembled using multiple threads | ||
| 604 | const auto& grid = this->problem(freeFlowMassIndex).gridGeometry().gridView().grid(); | ||
| 605 | for (const auto& elements : elementSets_) | ||
| 606 | { | ||
| 607 | Dumux::parallelFor(elements.size(), [&](const std::size_t eIdx) | ||
| 608 | { | ||
| 609 | const auto element = grid.entity(elements[eIdx]); | ||
| 610 | assembleElement(element); | ||
| 611 | }); | ||
| 612 | } | ||
| 613 | } | ||
| 614 | |||
| 615 | private: | ||
| 616 | 2002638 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 617 | const Element<freeFlowMomentumIndex>& elementI) const | ||
| 618 | { | ||
| 619 | // The call to this->problem() is expensive because of std::weak_ptr (see base class). Here we try to avoid it if possible. | ||
| 620 |
4/4✓ Branch 1 taken 7 times.
✓ Branch 2 taken 2002631 times.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 2002631 times.
|
2002638 | if (momentumCouplingContext_().empty()) |
| 621 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | bindCouplingContext_(domainI, elementI, this->problem(freeFlowMomentumIndex).gridGeometry().elementMapper().index(elementI)); |
| 622 | else | ||
| 623 | 2002631 | bindCouplingContext_(domainI, elementI, momentumCouplingContext_()[0].fvGeometry.gridGeometry().elementMapper().index(elementI)); | |
| 624 | 2002638 | } | |
| 625 | |||
| 626 | 4911018 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 627 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 628 | const std::size_t eIdx) const | ||
| 629 | { | ||
| 630 |
4/4✓ Branch 1 taken 9 times.
✓ Branch 2 taken 4911009 times.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 4911009 times.
|
4911018 | if (momentumCouplingContext_().empty()) |
| 631 | { | ||
| 632 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
11 | auto fvGeometry = localView(this->problem(freeFlowMassIndex).gridGeometry()); |
| 633 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
9 | fvGeometry.bind(elementI); |
| 634 | |||
| 635 |
5/8✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
12 | auto curElemVolVars = localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 636 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
|
18 | curElemVolVars.bind(elementI, fvGeometry, this->curSol(freeFlowMassIndex)); |
| 637 | |||
| 638 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
9 | auto prevElemVolVars = isTransient_ ? localView(gridVars_(freeFlowMassIndex).prevGridVolVars()) |
| 639 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
9 | : localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 640 | |||
| 641 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (isTransient_) |
| 642 | ✗ | prevElemVolVars.bindElement(elementI, fvGeometry, (*prevSol_)[freeFlowMassIndex]); | |
| 643 | |||
| 644 |
6/12✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
27 | momentumCouplingContext_().emplace_back(MomentumCouplingContext{std::move(fvGeometry), std::move(curElemVolVars), std::move(prevElemVolVars), eIdx}); |
| 645 | } | ||
| 646 |
2/2✓ Branch 1 taken 408126 times.
✓ Branch 2 taken 4502883 times.
|
4911009 | else if (eIdx != momentumCouplingContext_()[0].eIdx) |
| 647 | { | ||
| 648 | 408126 | momentumCouplingContext_()[0].eIdx = eIdx; | |
| 649 | 408126 | momentumCouplingContext_()[0].fvGeometry.bind(elementI); | |
| 650 | 408126 | momentumCouplingContext_()[0].curElemVolVars.bind(elementI, momentumCouplingContext_()[0].fvGeometry, this->curSol(freeFlowMassIndex)); | |
| 651 | |||
| 652 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 408126 times.
|
408126 | if (isTransient_) |
| 653 | ✗ | momentumCouplingContext_()[0].prevElemVolVars.bindElement(elementI, momentumCouplingContext_()[0].fvGeometry, (*prevSol_)[freeFlowMassIndex]); | |
| 654 | } | ||
| 655 | 4911018 | } | |
| 656 | |||
| 657 | 214494 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 658 | const Element<freeFlowMassIndex>& elementI) const | ||
| 659 | { | ||
| 660 | // The call to this->problem() is expensive because of std::weak_ptr (see base class). Here we try to avoid it if possible. | ||
| 661 |
4/4✓ Branch 1 taken 8 times.
✓ Branch 2 taken 214486 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 214486 times.
|
214494 | if (massAndEnergyCouplingContext_().empty()) |
| 662 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | bindCouplingContext_(domainI, elementI, this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI)); |
| 663 | else | ||
| 664 | 214486 | bindCouplingContext_(domainI, elementI, massAndEnergyCouplingContext_()[0].fvGeometry.gridGeometry().elementMapper().index(elementI)); | |
| 665 | 214494 | } | |
| 666 | |||
| 667 | 23795790 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 668 | const Element<freeFlowMassIndex>& elementI, | ||
| 669 | const std::size_t eIdx) const | ||
| 670 | { | ||
| 671 |
4/4✓ Branch 1 taken 9 times.
✓ Branch 2 taken 23795781 times.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 23795781 times.
|
23795790 | if (massAndEnergyCouplingContext_().empty()) |
| 672 | { | ||
| 673 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
|
9 | const auto& gridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 674 |
2/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
13 | auto fvGeometry = localView(gridGeometry); |
| 675 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
9 | fvGeometry.bindElement(elementI); |
| 676 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
9 | massAndEnergyCouplingContext_().emplace_back(std::move(fvGeometry), eIdx); |
| 677 | } | ||
| 678 |
2/2✓ Branch 1 taken 2279473 times.
✓ Branch 2 taken 21516308 times.
|
23795781 | else if (eIdx != massAndEnergyCouplingContext_()[0].eIdx) |
| 679 | { | ||
| 680 | 2279473 | massAndEnergyCouplingContext_()[0].eIdx = eIdx; | |
| 681 | 2279473 | massAndEnergyCouplingContext_()[0].fvGeometry.bindElement(elementI); | |
| 682 | } | ||
| 683 | 23795790 | } | |
| 684 | |||
| 685 | /*! | ||
| 686 | * \brief Return a reference to the grid variables of a sub problem | ||
| 687 | * \param domainIdx The domain index | ||
| 688 | */ | ||
| 689 | template<std::size_t i> | ||
| 690 | 18 | const GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) const | |
| 691 | { | ||
| 692 |
2/4✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
|
36 | if (std::get<i>(gridVariables_)) |
| 693 | 54 | return *std::get<i>(gridVariables_); | |
| 694 | else | ||
| 695 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 696 | } | ||
| 697 | |||
| 698 | /*! | ||
| 699 | * \brief Return a reference to the grid variables of a sub problem | ||
| 700 | * \param domainIdx The domain index | ||
| 701 | */ | ||
| 702 | template<std::size_t i> | ||
| 703 | 2002638 | GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) | |
| 704 | { | ||
| 705 |
2/4✓ Branch 0 taken 2002638 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2002638 times.
✗ Branch 3 not taken.
|
4005276 | if (std::get<i>(gridVariables_)) |
| 706 | 6007914 | return *std::get<i>(gridVariables_); | |
| 707 | else | ||
| 708 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 709 | } | ||
| 710 | |||
| 711 | |||
| 712 | 9 | void computeCouplingStencils_() | |
| 713 | { | ||
| 714 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
|
9 | const auto& momentumGridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 715 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | const auto& massGridGeometry = this->problem(freeFlowMassIndex).gridGeometry(); |
| 716 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
9 | auto momentumFvGeometry = localView(momentumGridGeometry); |
| 717 |
1/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
11 | auto massFvGeometry = localView(massGridGeometry); |
| 718 | |||
| 719 | 9 | massAndEnergyToMomentumStencils_.clear(); | |
| 720 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
18 | massAndEnergyToMomentumStencils_.resize(massGridGeometry.gridView().size(0)); |
| 721 | |||
| 722 | 9 | momentumToMassAndEnergyStencils_.clear(); | |
| 723 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
18 | momentumToMassAndEnergyStencils_.resize(momentumGridGeometry.gridView().size(0)); |
| 724 | |||
| 725 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 9 times.
|
27 | assert(massAndEnergyToMomentumStencils_.size() == momentumToMassAndEnergyStencils_.size()); |
| 726 | |||
| 727 |
14/19✓ Branch 1 taken 2 times.
✓ Branch 2 taken 100447 times.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 100449 times.
✓ Branch 5 taken 7 times.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 10400 times.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 10400 times.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 10400 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 10400 times.
✗ Branch 21 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
|
131671 | for (const auto& element : elements(momentumGridGeometry.gridView())) |
| 728 | { | ||
| 729 |
1/2✓ Branch 1 taken 10400 times.
✗ Branch 2 not taken.
|
110847 | momentumFvGeometry.bindElement(element); |
| 730 |
1/2✓ Branch 1 taken 10400 times.
✗ Branch 2 not taken.
|
110847 | massFvGeometry.bindElement(element); |
| 731 | 110847 | const auto eIdx = momentumFvGeometry.elementIndex(); | |
| 732 | |||
| 733 |
4/4✓ Branch 0 taken 455454 times.
✓ Branch 1 taken 110847 times.
✓ Branch 2 taken 455454 times.
✓ Branch 3 taken 110847 times.
|
1132602 | for (const auto& scv : scvs(momentumFvGeometry)) |
| 734 |
2/4✓ Branch 1 taken 38400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 38400 times.
✗ Branch 5 not taken.
|
910908 | massAndEnergyToMomentumStencils_[eIdx].push_back(scv.dofIndex()); |
| 735 | |||
| 736 |
4/4✓ Branch 0 taken 285735 times.
✓ Branch 1 taken 110847 times.
✓ Branch 2 taken 285735 times.
✓ Branch 3 taken 110847 times.
|
793164 | for (const auto& scv : scvs(massFvGeometry)) |
| 737 |
2/4✓ Branch 1 taken 31200 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31200 times.
✗ Branch 5 not taken.
|
610905 | momentumToMassAndEnergyStencils_[eIdx].push_back(scv.dofIndex()); |
| 738 | } | ||
| 739 | 9 | } | |
| 740 | |||
| 741 | CouplingStencilType emptyStencil_; | ||
| 742 | std::vector<CouplingStencilType> momentumToMassAndEnergyStencils_; | ||
| 743 | std::vector<CouplingStencilType> massAndEnergyToMomentumStencils_; | ||
| 744 | |||
| 745 | // the coupling context exists for each thread | ||
| 746 | // TODO this is a bad pattern, just like mutable caches | ||
| 747 | // we should really construct and pass the context and not store it globally | ||
| 748 | ✗ | std::vector<MomentumCouplingContext>& momentumCouplingContext_() const | |
| 749 | { | ||
| 750 | ✗ | thread_local static std::vector<MomentumCouplingContext> c; | |
| 751 | ✗ | return c; | |
| 752 | } | ||
| 753 | |||
| 754 | // the coupling context exists for each thread | ||
| 755 | ✗ | std::vector<MassAndEnergyCouplingContext>& massAndEnergyCouplingContext_() const | |
| 756 | { | ||
| 757 | ✗ | thread_local static std::vector<MassAndEnergyCouplingContext> c; | |
| 758 | ✗ | return c; | |
| 759 | } | ||
| 760 | |||
| 761 | //! A tuple of std::shared_ptrs to the grid variables of the sub problems | ||
| 762 | GridVariablesTuple gridVariables_; | ||
| 763 | |||
| 764 | const SolutionVector* prevSol_; | ||
| 765 | bool isTransient_; | ||
| 766 | |||
| 767 | std::deque<std::vector<ElementSeed<freeFlowMomentumIndex>>> elementSets_; | ||
| 768 | }; | ||
| 769 | |||
| 770 | namespace Detail { | ||
| 771 | |||
| 772 | // declaration (specialize for different discretization types) | ||
| 773 | template<class Traits, class DiscretizationMethod = typename Detail::MomentumDiscretizationMethod<Traits>::type> | ||
| 774 | struct CouplingManagerSupportsMultithreadedAssemblySelector; | ||
| 775 | |||
| 776 | // disabled for now | ||
| 777 | // the infrastructure for multithreaded assembly is implemented (see code in the class above) | ||
| 778 | // but the current implementation seems to have a bug and may cause race conditions. | ||
| 779 | // The result is different when running in parallel. After this has been fixed activate multithreaded assembly | ||
| 780 | // by removing this specialization | ||
| 781 | template<class Traits, class D> | ||
| 782 | struct CouplingManagerSupportsMultithreadedAssemblySelector<Traits, DiscretizationMethods::CVFE<D>> | ||
| 783 | { using type = std::false_type; }; | ||
| 784 | |||
| 785 | } // end namespace Detail | ||
| 786 | |||
| 787 | //! whether we support multithreaded assembly | ||
| 788 | template<class T> | ||
| 789 | struct CouplingManagerSupportsMultithreadedAssembly<CVFEFreeFlowCouplingManager<T>> | ||
| 790 | : public Detail::CouplingManagerSupportsMultithreadedAssemblySelector<T>::type | ||
| 791 | {}; | ||
| 792 | |||
| 793 | } // end namespace Dumux | ||
| 794 | |||
| 795 | #endif | ||
| 796 |