| 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_STAGGERED_HH | ||
| 13 | #define DUMUX_MULTIDOMAIN_FREEFLOW_COUPLING_MANAGER_STAGGERED_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 | |||
| 24 | #include <dumux/common/properties.hh> | ||
| 25 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 26 | |||
| 27 | #include <dumux/discretization/method.hh> | ||
| 28 | #include <dumux/discretization/evalsolution.hh> | ||
| 29 | #include <dumux/discretization/elementsolution.hh> | ||
| 30 | |||
| 31 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 32 | #include <dumux/multidomain/fvassembler.hh> | ||
| 33 | #include <dumux/discretization/facecentered/staggered/consistentlyorientedgrid.hh> | ||
| 34 | |||
| 35 | #include <dumux/parallel/parallel_for.hh> | ||
| 36 | #include <dumux/assembly/coloring.hh> | ||
| 37 | |||
| 38 | namespace Dumux { | ||
| 39 | |||
| 40 | /*! | ||
| 41 | * \ingroup MultiDomain | ||
| 42 | * \brief The interface of the coupling manager for free flow systems | ||
| 43 | * \note coupling manager the face-centered staggered discretization scheme | ||
| 44 | */ | ||
| 45 | template<class Traits> | ||
| 46 | class FCStaggeredFreeFlowCouplingManager | ||
| 47 | : public CouplingManager<Traits> | ||
| 48 | { | ||
| 49 | using ParentType = CouplingManager<Traits>; | ||
| 50 | public: | ||
| 51 | static constexpr auto freeFlowMomentumIndex = typename Traits::template SubDomain<0>::Index(); | ||
| 52 | static constexpr auto freeFlowMassIndex = typename Traits::template SubDomain<1>::Index(); | ||
| 53 | |||
| 54 | // this can be used if the coupling manager is used inside a meta-coupling manager (e.g. multi-binary) | ||
| 55 | // to manager the solution vector storage outside this class | ||
| 56 | using SolutionVectorStorage = typename ParentType::SolutionVectorStorage; | ||
| 57 | private: | ||
| 58 | template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomain<id>::TypeTag; | ||
| 59 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 60 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 61 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 62 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 63 | template<std::size_t id> using ElementSeed = typename GridView<id>::Grid::template Codim<0>::EntitySeed; | ||
| 64 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 65 | template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume; | ||
| 66 | template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace; | ||
| 67 | template<std::size_t id> using GridVariables = typename Traits::template SubDomain<id>::GridVariables; | ||
| 68 | template<std::size_t id> using ElementVolumeVariables = typename GridVariables<id>::GridVolumeVariables::LocalView; | ||
| 69 | template<std::size_t id> using GridFluxVariablesCache = typename GridVariables<id>::GridFluxVariablesCache; | ||
| 70 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 71 | template<std::size_t id> using VolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::VolumeVariables>; | ||
| 72 | |||
| 73 | using Scalar = typename Traits::Scalar; | ||
| 74 | using SolutionVector = typename Traits::SolutionVector; | ||
| 75 | |||
| 76 | using CouplingStencilType = std::vector<std::size_t>; | ||
| 77 | |||
| 78 | using GridVariablesTuple = typename Traits::template TupleOfSharedPtr<GridVariables>; | ||
| 79 | |||
| 80 | using FluidSystem = typename VolumeVariables<freeFlowMassIndex>::FluidSystem; | ||
| 81 | |||
| 82 | using VelocityVector = typename SubControlVolumeFace<freeFlowMassIndex>::GlobalPosition; | ||
| 83 | static_assert(std::is_same_v<VelocityVector, typename SubControlVolumeFace<freeFlowMomentumIndex>::GlobalPosition>); | ||
| 84 | |||
| 85 | 96 | struct MomentumCouplingContext | |
| 86 | { | ||
| 87 | FVElementGeometry<freeFlowMassIndex> fvGeometry; | ||
| 88 | ElementVolumeVariables<freeFlowMassIndex> curElemVolVars; | ||
| 89 | ElementVolumeVariables<freeFlowMassIndex> prevElemVolVars; | ||
| 90 | std::size_t eIdx; | ||
| 91 | }; | ||
| 92 | |||
| 93 | ✗ | struct MassAndEnergyCouplingContext | |
| 94 | { | ||
| 95 | 55 | MassAndEnergyCouplingContext(FVElementGeometry<freeFlowMomentumIndex>&& f, const std::size_t i) | |
| 96 | 55 | : fvGeometry(std::move(f)) | |
| 97 | 55 | , eIdx(i) | |
| 98 | {} | ||
| 99 | |||
| 100 | FVElementGeometry<freeFlowMomentumIndex> fvGeometry; | ||
| 101 | std::size_t eIdx; | ||
| 102 | }; | ||
| 103 | |||
| 104 | public: | ||
| 105 | |||
| 106 | static constexpr auto pressureIdx = VolumeVariables<freeFlowMassIndex>::Indices::pressureIdx; | ||
| 107 | |||
| 108 | /*! | ||
| 109 | * \brief Methods to be accessed by main | ||
| 110 | */ | ||
| 111 | // \{ | ||
| 112 | |||
| 113 | //! use as regular coupling manager | ||
| 114 | 39 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | |
| 115 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 116 | GridVariablesTuple&& gridVariables, | ||
| 117 | const SolutionVector& curSol) | ||
| 118 | { | ||
| 119 | 78 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | |
| 120 | 39 | gridVariables_ = gridVariables; | |
| 121 | 39 | this->updateSolution(curSol); | |
| 122 | |||
| 123 | 39 | computeCouplingStencils_(); | |
| 124 | 39 | } | |
| 125 | |||
| 126 | //! use as regular coupling manager in a transient setting | ||
| 127 | 21 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | |
| 128 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 129 | GridVariablesTuple&& gridVariables, | ||
| 130 | const SolutionVector& curSol, | ||
| 131 | const SolutionVector& prevSol) | ||
| 132 | { | ||
| 133 |
3/10✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
42 | init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol); |
| 134 | 21 | prevSol_ = &prevSol; | |
| 135 | 21 | isTransient_ = true; | |
| 136 | 21 | } | |
| 137 | |||
| 138 | //! use as binary coupling manager in multi model context | ||
| 139 | 16 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem, | |
| 140 | std::shared_ptr<Problem<freeFlowMassIndex>> massProblem, | ||
| 141 | GridVariablesTuple&& gridVariables, | ||
| 142 | typename ParentType::SolutionVectorStorage& curSol) | ||
| 143 | { | ||
| 144 | 32 | this->setSubProblems(std::make_tuple(momentumProblem, massProblem)); | |
| 145 | 16 | gridVariables_ = gridVariables; | |
| 146 | 16 | this->attachSolution(curSol); | |
| 147 | |||
| 148 | 16 | computeCouplingStencils_(); | |
| 149 | 16 | } | |
| 150 | |||
| 151 | // \} | ||
| 152 | |||
| 153 | using CouplingManager<Traits>::evalCouplingResidual; | ||
| 154 | |||
| 155 | /*! | ||
| 156 | * \brief evaluates the element residual of a coupled element of domain i which depends on the variables | ||
| 157 | * at the degree of freedom with index dofIdxGlobalJ of domain j | ||
| 158 | * | ||
| 159 | * \param domainI the domain index of domain i | ||
| 160 | * \param localAssemblerI the local assembler assembling the element residual of an element of domain i | ||
| 161 | * \param scvI the sub-control-volume of domain i | ||
| 162 | * \param domainJ the domain index of domain j | ||
| 163 | * \param dofIdxGlobalJ the index of the degree of freedom of domain j which has an influence on the element residual of domain i | ||
| 164 | * | ||
| 165 | * \note the element whose residual is to be evaluated can be retrieved from the local assembler | ||
| 166 | * as localAssemblerI.element() as well as all up-to-date variables and caches. | ||
| 167 | * \note the default implementation evaluates the complete element residual | ||
| 168 | * if only parts (i.e. only certain scvs, or only certain terms of the residual) of the residual are coupled | ||
| 169 | * to dof with index dofIdxGlobalJ the function can be overloaded in the coupling manager | ||
| 170 | * \return the element residual | ||
| 171 | */ | ||
| 172 | template<std::size_t j, class LocalAssemblerI> | ||
| 173 | 23478912 | decltype(auto) evalCouplingResidual(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 174 | const LocalAssemblerI& localAssemblerI, | ||
| 175 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 176 | Dune::index_constant<j> domainJ, | ||
| 177 | std::size_t dofIdxGlobalJ) const | ||
| 178 | { | ||
| 179 | 23478912 | const auto& problem = localAssemblerI.problem(); | |
| 180 | 23478912 | const auto& element = localAssemblerI.element(); | |
| 181 | 23478912 | const auto& fvGeometry = localAssemblerI.fvGeometry(); | |
| 182 | 23478912 | const auto& curElemVolVars = localAssemblerI.curElemVolVars(); | |
| 183 | 23478912 | const auto& prevElemVolVars = localAssemblerI.prevElemVolVars(); | |
| 184 | 46899648 | typename LocalAssemblerI::ElementResidualVector residual(localAssemblerI.element().subEntities(1)); | |
| 185 | 23478912 | const auto& localResidual = localAssemblerI.localResidual(); | |
| 186 | |||
| 187 | 23478912 | localResidual.evalSource(residual, problem, element, fvGeometry, curElemVolVars, scvI); | |
| 188 | |||
| 189 |
4/4✓ Branch 1 taken 71081504 times.
✓ Branch 2 taken 23478912 times.
✓ Branch 3 taken 71081504 times.
✓ Branch 4 taken 23478912 times.
|
94560416 | for (const auto& scvf : scvfs(fvGeometry, scvI)) |
| 190 | 213244512 | localResidual.evalFlux(residual, problem, element, fvGeometry, curElemVolVars, localAssemblerI.elemBcTypes(), localAssemblerI.elemFluxVarsCache(), scvf); | |
| 191 | |||
| 192 |
4/4✓ Branch 0 taken 16207792 times.
✓ Branch 1 taken 7271120 times.
✓ Branch 2 taken 16207792 times.
✓ Branch 3 taken 7271120 times.
|
46957824 | if (!localAssemblerI.assembler().isStationaryProblem()) |
| 193 | { | ||
| 194 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16207792 times.
|
16207792 | assert(isTransient_); |
| 195 | 16207792 | localResidual.evalStorage(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, scvI); | |
| 196 | } | ||
| 197 | |||
| 198 | 23478912 | return residual; | |
| 199 | } | ||
| 200 | |||
| 201 | /*! | ||
| 202 | * \name member functions concerning the coupling stencils | ||
| 203 | */ | ||
| 204 | // \{ | ||
| 205 | |||
| 206 | /*! | ||
| 207 | * \brief Returns the pressure at a given _frontal_ sub control volume face. | ||
| 208 | */ | ||
| 209 | ✗ | Scalar pressure(const Element<freeFlowMomentumIndex>& element, | |
| 210 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 211 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const | ||
| 212 | { | ||
| 213 | ✗ | assert(scvf.isFrontal() && !scvf.isLateral() && !scvf.boundary()); | |
| 214 | ✗ | return this->curSol(freeFlowMassIndex)[fvGeometry.elementIndex()][pressureIdx]; | |
| 215 | } | ||
| 216 | |||
| 217 | /*! | ||
| 218 | * \brief Returns the pressure at the center of a sub control volume corresponding to a given sub control volume face. | ||
| 219 | * This is used for setting a Dirichlet pressure for the mass model when a fixed pressure for the momentum balance is set at another | ||
| 220 | * boundary. Since the the pressure at the given scvf is solution-dependent and thus unknown a priori, we just use the value | ||
| 221 | * of the interior cell here. | ||
| 222 | */ | ||
| 223 | ✗ | Scalar cellPressure(const Element<freeFlowMassIndex>& element, | |
| 224 | const SubControlVolumeFace<freeFlowMassIndex>& scvf) const | ||
| 225 | { | ||
| 226 |
20/20✓ Branch 0 taken 1420 times.
✓ Branch 1 taken 9300 times.
✓ Branch 2 taken 1420 times.
✓ Branch 3 taken 9300 times.
✓ Branch 4 taken 1420 times.
✓ Branch 5 taken 9300 times.
✓ Branch 6 taken 1420 times.
✓ Branch 7 taken 9300 times.
✓ Branch 8 taken 1420 times.
✓ Branch 9 taken 9300 times.
✓ Branch 10 taken 2050 times.
✓ Branch 11 taken 1750 times.
✓ Branch 12 taken 2050 times.
✓ Branch 13 taken 1750 times.
✓ Branch 14 taken 2050 times.
✓ Branch 15 taken 1750 times.
✓ Branch 16 taken 2050 times.
✓ Branch 17 taken 1750 times.
✓ Branch 18 taken 2050 times.
✓ Branch 19 taken 1750 times.
|
83300 | return this->curSol(freeFlowMassIndex)[scvf.insideScvIdx()][pressureIdx]; |
| 227 | } | ||
| 228 | |||
| 229 | /*! | ||
| 230 | * \brief Returns the density at a given sub control volume face. | ||
| 231 | */ | ||
| 232 | 65362606 | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 233 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 234 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 235 | const bool considerPreviousTimeStep = false) const | ||
| 236 | { | ||
| 237 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 65362606 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
65362606 | assert(!(considerPreviousTimeStep && !isTransient_)); |
| 238 | 128813060 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 239 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 63450454 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 63450454 times.
|
130725212 | const auto& insideMomentumScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 240 |
2/2✓ Branch 1 taken 1912152 times.
✓ Branch 2 taken 63450454 times.
|
65362606 | const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(insideMomentumScv.elementIndex()); |
| 241 | |||
| 242 | 196002078 | const auto rho = [&](const auto& elemVolVars) | |
| 243 | { | ||
| 244 |
2/2✓ Branch 0 taken 85740 times.
✓ Branch 1 taken 65276866 times.
|
65362606 | if (scvf.boundary()) |
| 245 | 85740 | return elemVolVars[insideMassScv].density(); | |
| 246 | else | ||
| 247 | { | ||
| 248 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 63373338 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 63373338 times.
|
130553732 | const auto& outsideMomentumScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 249 |
1/2✓ Branch 1 taken 1903528 times.
✗ Branch 2 not taken.
|
65276866 | const auto& outsideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(outsideMomentumScv.elementIndex()); |
| 250 | // TODO distance weighting | ||
| 251 | 65276866 | return 0.5*(elemVolVars[insideMassScv].density() + elemVolVars[outsideMassScv].density()); | |
| 252 | } | ||
| 253 | }; | ||
| 254 | |||
| 255 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 65362606 times.
|
65362606 | return considerPreviousTimeStep ? rho(momentumCouplingContext_()[0].prevElemVolVars) |
| 256 | 65362606 | : rho(momentumCouplingContext_()[0].curElemVolVars); | |
| 257 | } | ||
| 258 | |||
| 259 | 129364612 | auto insideAndOutsideDensity(const Element<freeFlowMomentumIndex>& element, | |
| 260 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 261 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 262 | const bool considerPreviousTimeStep = false) const | ||
| 263 | { | ||
| 264 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 129364612 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
129364612 | assert(!(considerPreviousTimeStep && !isTransient_)); |
| 265 | 255012104 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 266 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 125647492 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 125647492 times.
|
258729224 | const auto& insideMomentumScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 267 |
3/4✓ Branch 1 taken 3717120 times.
✓ Branch 2 taken 125647492 times.
✓ Branch 3 taken 3717120 times.
✗ Branch 4 not taken.
|
129364612 | const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(insideMomentumScv.elementIndex()); |
| 268 | |||
| 269 | 388042948 | const auto result = [&](const auto& elemVolVars) | |
| 270 | { | ||
| 271 |
2/2✓ Branch 0 taken 50888 times.
✓ Branch 1 taken 129313724 times.
|
129364612 | if (scvf.boundary()) |
| 272 | 50888 | return std::make_pair(elemVolVars[insideMassScv].density(), elemVolVars[insideMassScv].density()); | |
| 273 | else | ||
| 274 | { | ||
| 275 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 125597132 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 125597132 times.
|
258627448 | const auto& outsideMomentumScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 276 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3716592 times.
|
129313724 | const auto& outsideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(outsideMomentumScv.elementIndex()); |
| 277 | 129313724 | return std::make_pair(elemVolVars[insideMassScv].density(), elemVolVars[outsideMassScv].density()); | |
| 278 | } | ||
| 279 | }; | ||
| 280 | |||
| 281 | ✗ | return considerPreviousTimeStep ? result(momentumCouplingContext_()[0].prevElemVolVars) | |
| 282 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 129364612 times.
|
129364612 | : result(momentumCouplingContext_()[0].curElemVolVars); |
| 283 | } | ||
| 284 | |||
| 285 | /*! | ||
| 286 | * \brief Returns the density at a given sub control volume. | ||
| 287 | */ | ||
| 288 | 200944950 | Scalar density(const Element<freeFlowMomentumIndex>& element, | |
| 289 | const SubControlVolume<freeFlowMomentumIndex>& scv, | ||
| 290 | const bool considerPreviousTimeStep = false) const | ||
| 291 | { | ||
| 292 |
3/4✓ Branch 0 taken 52411520 times.
✓ Branch 1 taken 148533430 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 52411520 times.
|
200944950 | assert(!(considerPreviousTimeStep && !isTransient_)); |
| 293 | 200944950 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, scv.elementIndex()); | |
| 294 |
4/4✓ Branch 1 taken 52411520 times.
✓ Branch 2 taken 148533430 times.
✓ Branch 3 taken 50464520 times.
✓ Branch 4 taken 144395958 times.
|
200944950 | const auto& massScv = (*scvs(momentumCouplingContext_()[0].fvGeometry).begin()); |
| 295 | |||
| 296 |
2/2✓ Branch 0 taken 52411520 times.
✓ Branch 1 taken 148533430 times.
|
493874338 | return considerPreviousTimeStep ? momentumCouplingContext_()[0].prevElemVolVars[massScv].density() |
| 297 | 148533430 | : momentumCouplingContext_()[0].curElemVolVars[massScv].density(); | |
| 298 | } | ||
| 299 | |||
| 300 | /*! | ||
| 301 | * \brief Returns the pressure at a given sub control volume face. | ||
| 302 | */ | ||
| 303 | 280919120 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 304 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 305 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const | ||
| 306 | { | ||
| 307 | 555680536 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 308 | |||
| 309 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 274761416 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 274761416 times.
|
561838240 | const auto& insideMomentumScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 310 |
2/2✓ Branch 1 taken 8925270 times.
✓ Branch 2 taken 271993850 times.
|
280919120 | const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(insideMomentumScv.elementIndex()); |
| 311 | |||
| 312 |
2/2✓ Branch 0 taken 2867150 times.
✓ Branch 1 taken 278051970 times.
|
280919120 | if (scvf.boundary()) |
| 313 | 2867150 | return momentumCouplingContext_()[0].curElemVolVars[insideMassScv].viscosity(); | |
| 314 | |||
| 315 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 271993850 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 271993850 times.
|
556103940 | const auto& outsideMomentumScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 316 |
2/2✓ Branch 1 taken 2052568 times.
✓ Branch 2 taken 4005552 times.
|
278051970 | const auto& outsideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(outsideMomentumScv.elementIndex()); |
| 317 | |||
| 318 | 834155910 | const auto mu = [&](const auto& elemVolVars) | |
| 319 | { | ||
| 320 | // TODO distance weighting | ||
| 321 | 278051970 | return 0.5*(elemVolVars[insideMassScv].viscosity() + elemVolVars[outsideMassScv].viscosity()); | |
| 322 | }; | ||
| 323 | |||
| 324 | 278051970 | return mu(momentumCouplingContext_()[0].curElemVolVars); | |
| 325 | } | ||
| 326 | |||
| 327 | /*! | ||
| 328 | * \brief Returns the pressure at a given sub control volume. | ||
| 329 | */ | ||
| 330 | 90344 | Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element, | |
| 331 | const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 332 | const SubControlVolume<freeFlowMomentumIndex>& scv) const | ||
| 333 | { | ||
| 334 | 180688 | bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex()); | |
| 335 | 90344 | const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(scv.elementIndex()); | |
| 336 | 90344 | return momentumCouplingContext_()[0].curElemVolVars[insideMassScv].viscosity(); | |
| 337 | } | ||
| 338 | |||
| 339 | /*! | ||
| 340 | * \brief Returns the velocity at a given sub control volume face. | ||
| 341 | */ | ||
| 342 | 127071136 | VelocityVector faceVelocity(const Element<freeFlowMassIndex>& element, | |
| 343 | const SubControlVolumeFace<freeFlowMassIndex>& scvf) const | ||
| 344 | { | ||
| 345 | // TODO: rethink this! Maybe we only need scvJ.dofIndex() | ||
| 346 | 254142272 | bindCouplingContext_(Dune::index_constant<freeFlowMassIndex>(), element, scvf.insideScvIdx()/*eIdx*/); | |
| 347 | |||
| 348 | // the TPFA scvf index corresponds to the staggered scv index (might need mapping) | ||
| 349 | 127071136 | const auto localMomentumScvIdx = massScvfToMomentumScvIdx_(scvf, massAndEnergyCouplingContext_()[0].fvGeometry); | |
| 350 |
2/3✗ Branch 0 not taken.
✓ Branch 1 taken 123609520 times.
✓ Branch 2 taken 192768 times.
|
127071136 | const auto& scvJ = massAndEnergyCouplingContext_()[0].fvGeometry.scv(localMomentumScvIdx); |
| 351 | |||
| 352 | // create a unit normal vector oriented in positive coordinate direction | ||
| 353 | 127071136 | typename SubControlVolumeFace<freeFlowMassIndex>::GlobalPosition velocity; | |
| 354 | 127071136 | velocity[scvJ.dofAxis()] = 1.0; | |
| 355 | |||
| 356 | // create the actual velocity vector | ||
| 357 | 381216016 | velocity *= this->curSol(freeFlowMomentumIndex)[scvJ.dofIndex()]; | |
| 358 | |||
| 359 | 127071136 | return velocity; | |
| 360 | } | ||
| 361 | |||
| 362 | /*! | ||
| 363 | * \brief The coupling stencil of domain I, i.e. which domain J DOFs | ||
| 364 | * the given domain I scv's residual depends on. | ||
| 365 | * | ||
| 366 | * \param domainI the domain index of domain i | ||
| 367 | * \param elementI the coupled element of domain à | ||
| 368 | * \param scvI the sub-control volume of domain i | ||
| 369 | * \param domainJ the domain index of domain j | ||
| 370 | */ | ||
| 371 | template<std::size_t j> | ||
| 372 | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | ||
| 373 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 374 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 375 | Dune::index_constant<j> domainJ) const | ||
| 376 | { return emptyStencil_; } | ||
| 377 | |||
| 378 | /*! | ||
| 379 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 380 | * that couple with / influence the element residual of the given element of domain i | ||
| 381 | * | ||
| 382 | * \param domainI the domain index of domain i | ||
| 383 | * \param elementI the coupled element of domain à | ||
| 384 | * \param domainJ the domain index of domain j | ||
| 385 | * | ||
| 386 | * \note The element residual definition depends on the discretization scheme of domain i | ||
| 387 | * box: a container of the residuals of all sub control volumes | ||
| 388 | * cc : the residual of the (sub) control volume | ||
| 389 | * fem: the residual of the element | ||
| 390 | * \note This function has to be implemented by all coupling managers for all combinations of i and j | ||
| 391 | */ | ||
| 392 | ✗ | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 393 | const Element<freeFlowMassIndex>& elementI, | ||
| 394 | Dune::index_constant<freeFlowMomentumIndex> domainJ) const | ||
| 395 | { | ||
| 396 | ✗ | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI); | |
| 397 | ✗ | return massAndEnergyToMomentumStencils_[eIdx]; | |
| 398 | } | ||
| 399 | |||
| 400 | /*! | ||
| 401 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 402 | * that couple with / influence the residual of the given sub-control volume of domain i | ||
| 403 | * | ||
| 404 | * \param domainI the domain index of domain i | ||
| 405 | * \param elementI the coupled element of domain à | ||
| 406 | * \param scvI the sub-control volume of domain i | ||
| 407 | * \param domainJ the domain index of domain j | ||
| 408 | */ | ||
| 409 | ✗ | const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 410 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 411 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 412 | Dune::index_constant<freeFlowMassIndex> domainJ) const | ||
| 413 | { | ||
| 414 | 17959632 | return momentumToMassAndEnergyStencils_[scvI.index()]; | |
| 415 | } | ||
| 416 | |||
| 417 | // \} | ||
| 418 | |||
| 419 | /*! | ||
| 420 | * \name member functions concerning variable caching for element residual evaluations | ||
| 421 | */ | ||
| 422 | // \{ | ||
| 423 | |||
| 424 | //! \copydoc CouplingManager::updateCouplingContext | ||
| 425 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 426 | 55596704 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 427 | const LocalAssemblerI& localAssemblerI, | ||
| 428 | Dune::index_constant<j> domainJ, | ||
| 429 | std::size_t dofIdxGlobalJ, | ||
| 430 | const PrimaryVariables<j>& priVarsJ, | ||
| 431 | int pvIdxJ) | ||
| 432 | { | ||
| 433 |
14/16✗ Branch 0 not taken.
✓ Branch 1 taken 1199200 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1904300 times.
✓ Branch 4 taken 246504 times.
✓ Branch 5 taken 1904300 times.
✓ Branch 6 taken 246620 times.
✓ Branch 7 taken 897656 times.
✓ Branch 8 taken 246620 times.
✓ Branch 9 taken 185535 times.
✓ Branch 10 taken 53037 times.
✓ Branch 11 taken 185535 times.
✓ Branch 12 taken 45965 times.
✓ Branch 13 taken 192556 times.
✓ Branch 14 taken 116 times.
✓ Branch 15 taken 185484 times.
|
596355424 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; |
| 434 | |||
| 435 | if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex) | ||
| 436 | { | ||
| 437 | 55596704 | bindCouplingContext_(domainI, localAssemblerI.element()); | |
| 438 | |||
| 439 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 27798352 times.
|
55596704 | const auto& problem = this->problem(domainJ); |
| 440 | 111193408 | const auto& deflectedElement = problem.gridGeometry().element(dofIdxGlobalJ); | |
| 441 |
3/6✓ Branch 1 taken 56576 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 56576 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 56576 times.
✗ Branch 8 not taken.
|
166790112 | const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), problem.gridGeometry()); |
| 442 | 55596704 | const auto& fvGeometry = momentumCouplingContext_()[0].fvGeometry; | |
| 443 | 55596704 | const auto scvIdxJ = dofIdxGlobalJ; | |
| 444 |
2/3✓ Branch 0 taken 468288 times.
✓ Branch 1 taken 881888 times.
✗ Branch 2 not taken.
|
55596704 | const auto& scv = fvGeometry.scv(scvIdxJ); |
| 445 | |||
| 446 | if constexpr (ElementVolumeVariables<freeFlowMassIndex>::GridVolumeVariables::cachingEnabled) | ||
| 447 |
3/6✓ Branch 1 taken 28288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28288 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28288 times.
✗ Branch 8 not taken.
|
52952928 | gridVars_(freeFlowMassIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), problem, deflectedElement, scv); |
| 448 | else | ||
| 449 |
1/2✓ Branch 3 taken 28288 times.
✗ Branch 4 not taken.
|
2643776 | momentumCouplingContext_()[0].curElemVolVars[scv].update(std::move(elemSol), problem, deflectedElement, scv); |
| 450 | } | ||
| 451 | 55596704 | } | |
| 452 | |||
| 453 | // \} | ||
| 454 | |||
| 455 | /*! | ||
| 456 | * \brief Compute colors for multithreaded assembly | ||
| 457 | */ | ||
| 458 | void computeColorsForAssembly() | ||
| 459 | { | ||
| 460 | // use coloring of the fc staggered discretization for both domains | ||
| 461 | elementSets_ = computeColoring(this->problem(freeFlowMomentumIndex).gridGeometry()).sets; | ||
| 462 | } | ||
| 463 | |||
| 464 | /*! | ||
| 465 | * \brief Execute assembly kernel in parallel | ||
| 466 | * | ||
| 467 | * \param domainI the domain index of domain i | ||
| 468 | * \param assembleElement kernel function to execute for one element | ||
| 469 | */ | ||
| 470 | template<std::size_t i, class AssembleElementFunc> | ||
| 471 | void assembleMultithreaded(Dune::index_constant<i> domainI, AssembleElementFunc&& assembleElement) const | ||
| 472 | { | ||
| 473 | if (elementSets_.empty()) | ||
| 474 | DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!"); | ||
| 475 | |||
| 476 | // make this element loop run in parallel | ||
| 477 | // for this we have to color the elements so that we don't get | ||
| 478 | // race conditions when writing into the global matrix | ||
| 479 | // each color can be assembled using multiple threads | ||
| 480 | const auto& grid = this->problem(freeFlowMomentumIndex).gridGeometry().gridView().grid(); | ||
| 481 | for (const auto& elements : elementSets_) | ||
| 482 | { | ||
| 483 | Dumux::parallelFor(elements.size(), [&](const std::size_t eIdx) | ||
| 484 | { | ||
| 485 | const auto element = grid.entity(elements[eIdx]); | ||
| 486 | assembleElement(element); | ||
| 487 | }); | ||
| 488 | } | ||
| 489 | } | ||
| 490 | |||
| 491 | private: | ||
| 492 | ✗ | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 493 | const Element<freeFlowMomentumIndex>& elementI) const | ||
| 494 | { | ||
| 495 | ✗ | const auto eIdx = this->problem(freeFlowMomentumIndex).gridGeometry().elementMapper().index(elementI); | |
| 496 | ✗ | bindCouplingContext_(domainI, elementI, eIdx); | |
| 497 | ✗ | } | |
| 498 | |||
| 499 | 704479984 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 500 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 501 | const std::size_t eIdx) const | ||
| 502 | { | ||
| 503 |
4/4✓ Branch 1 taken 55 times.
✓ Branch 2 taken 704479929 times.
✓ Branch 3 taken 55 times.
✓ Branch 4 taken 704479929 times.
|
704479984 | if (momentumCouplingContext_().empty()) |
| 504 | { | ||
| 505 |
3/8✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
62 | auto fvGeometry = localView(this->problem(freeFlowMassIndex).gridGeometry()); |
| 506 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
55 | fvGeometry.bind(elementI); |
| 507 | |||
| 508 |
6/10✓ Branch 1 taken 7 times.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
|
110 | auto curElemVolVars = localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 509 |
2/4✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
|
110 | curElemVolVars.bind(elementI, fvGeometry, this->curSol(freeFlowMassIndex)); |
| 510 | |||
| 511 |
3/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 34 times.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
55 | auto prevElemVolVars = isTransient_ ? localView(gridVars_(freeFlowMassIndex).prevGridVolVars()) |
| 512 |
1/2✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
|
34 | : localView(gridVars_(freeFlowMassIndex).curGridVolVars()); |
| 513 | |||
| 514 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
|
4 | if (isTransient_) |
| 515 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | prevElemVolVars.bindElement(elementI, fvGeometry, (*prevSol_)[freeFlowMassIndex]); |
| 516 | |||
| 517 |
7/11✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
✓ Branch 4 taken 44 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11 times.
✓ Branch 7 taken 44 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 44 times.
✗ Branch 11 not taken.
|
235 | momentumCouplingContext_().emplace_back(MomentumCouplingContext{std::move(fvGeometry), std::move(curElemVolVars), std::move(prevElemVolVars), eIdx}); |
| 518 | } | ||
| 519 |
2/2✓ Branch 1 taken 2599055 times.
✓ Branch 2 taken 701880874 times.
|
704479929 | else if (eIdx != momentumCouplingContext_()[0].eIdx) |
| 520 | { | ||
| 521 | 2599055 | momentumCouplingContext_()[0].eIdx = eIdx; | |
| 522 | 2599055 | momentumCouplingContext_()[0].fvGeometry.bind(elementI); | |
| 523 | 2599055 | momentumCouplingContext_()[0].curElemVolVars.bind(elementI, momentumCouplingContext_()[0].fvGeometry, this->curSol(freeFlowMassIndex)); | |
| 524 | |||
| 525 |
2/2✓ Branch 0 taken 1244783 times.
✓ Branch 1 taken 1354272 times.
|
2599055 | if (isTransient_) |
| 526 | 1244783 | momentumCouplingContext_()[0].prevElemVolVars.bindElement(elementI, momentumCouplingContext_()[0].fvGeometry, (*prevSol_)[freeFlowMassIndex]); | |
| 527 | } | ||
| 528 | 704479984 | } | |
| 529 | |||
| 530 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | ||
| 531 | const Element<freeFlowMassIndex>& elementI) const | ||
| 532 | { | ||
| 533 | const auto eIdx = this->problem(freeFlowMassIndex).gridGeometry().elementMapper().index(elementI); | ||
| 534 | bindCouplingContext_(domainI, elementI, eIdx); | ||
| 535 | } | ||
| 536 | |||
| 537 | 127071136 | void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, | |
| 538 | const Element<freeFlowMassIndex>& elementI, | ||
| 539 | const std::size_t eIdx) const | ||
| 540 | { | ||
| 541 |
4/4✓ Branch 1 taken 55 times.
✓ Branch 2 taken 127071081 times.
✓ Branch 3 taken 55 times.
✓ Branch 4 taken 127071081 times.
|
127071136 | if (massAndEnergyCouplingContext_().empty()) |
| 542 | { | ||
| 543 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
|
55 | const auto& gridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 544 | 55 | auto fvGeometry = localView(gridGeometry); | |
| 545 | 55 | fvGeometry.bindElement(elementI); | |
| 546 | 55 | massAndEnergyCouplingContext_().emplace_back(std::move(fvGeometry), eIdx); | |
| 547 | } | ||
| 548 |
2/2✓ Branch 1 taken 24476096 times.
✓ Branch 2 taken 102594985 times.
|
127071081 | else if (eIdx != massAndEnergyCouplingContext_()[0].eIdx) |
| 549 | { | ||
| 550 | 24476096 | massAndEnergyCouplingContext_()[0].eIdx = eIdx; | |
| 551 | 24476096 | massAndEnergyCouplingContext_()[0].fvGeometry.bindElement(elementI); | |
| 552 | } | ||
| 553 | 127071136 | } | |
| 554 | |||
| 555 | /*! | ||
| 556 | * \brief Return a reference to the grid variables of a sub problem | ||
| 557 | * \param domainIdx The domain index | ||
| 558 | */ | ||
| 559 | template<std::size_t i> | ||
| 560 | 110 | const GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) const | |
| 561 | { | ||
| 562 |
2/4✓ Branch 0 taken 110 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
|
220 | if (std::get<i>(gridVariables_)) |
| 563 | 330 | return *std::get<i>(gridVariables_); | |
| 564 | else | ||
| 565 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 566 | } | ||
| 567 | |||
| 568 | /*! | ||
| 569 | * \brief Return a reference to the grid variables of a sub problem | ||
| 570 | * \param domainIdx The domain index | ||
| 571 | */ | ||
| 572 | template<std::size_t i> | ||
| 573 | 26476464 | GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) | |
| 574 | { | ||
| 575 |
2/4✓ Branch 0 taken 26476464 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 26476464 times.
✗ Branch 3 not taken.
|
52952928 | if (std::get<i>(gridVariables_)) |
| 576 | 79429392 | return *std::get<i>(gridVariables_); | |
| 577 | else | ||
| 578 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 579 | } | ||
| 580 | |||
| 581 | |||
| 582 | 55 | void computeCouplingStencils_() | |
| 583 | { | ||
| 584 | // TODO higher order | ||
| 585 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
|
55 | const auto& momentumGridGeometry = this->problem(freeFlowMomentumIndex).gridGeometry(); |
| 586 | 55 | auto momentumFvGeometry = localView(momentumGridGeometry); | |
| 587 | 55 | massAndEnergyToMomentumStencils_.clear(); | |
| 588 | 111 | massAndEnergyToMomentumStencils_.resize(momentumGridGeometry.gridView().size(0)); | |
| 589 | |||
| 590 | 55 | momentumToMassAndEnergyStencils_.clear(); | |
| 591 | 106 | momentumToMassAndEnergyStencils_.resize(momentumGridGeometry.numScv()); | |
| 592 | |||
| 593 |
11/14✓ Branch 2 taken 327323 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 100 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3536 times.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 3536 times.
✓ Branch 11 taken 6 times.
✓ Branch 13 taken 3536 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3536 times.
✗ Branch 17 not taken.
|
665281 | for (const auto& element : elements(momentumGridGeometry.gridView())) |
| 594 | { | ||
| 595 |
2/4✓ Branch 1 taken 3536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3536 times.
✗ Branch 5 not taken.
|
661620 | const auto eIdx = momentumGridGeometry.elementMapper().index(element); |
| 596 |
1/2✓ Branch 1 taken 3536 times.
✗ Branch 2 not taken.
|
330810 | momentumFvGeometry.bind(element); |
| 597 |
5/5✓ Branch 0 taken 12072 times.
✓ Branch 1 taken 1315986 times.
✓ Branch 2 taken 327792 times.
✓ Branch 3 taken 1312968 times.
✓ Branch 4 taken 327792 times.
|
2946418 | for (const auto& scv : scvs(momentumFvGeometry)) |
| 598 | { | ||
| 599 |
2/4✓ Branch 1 taken 14144 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14144 times.
✗ Branch 5 not taken.
|
2650080 | massAndEnergyToMomentumStencils_[eIdx].push_back(scv.dofIndex()); |
| 600 |
2/4✓ Branch 1 taken 14144 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14144 times.
✗ Branch 5 not taken.
|
2650080 | momentumToMassAndEnergyStencils_[scv.index()].push_back(eIdx); |
| 601 | |||
| 602 | // extend the stencil for fluids with variable viscosity and density, | ||
| 603 | if constexpr (FluidSystem::isCompressible(0/*phaseIdx*/)) | ||
| 604 | // if constexpr (FluidSystem::isCompressible(0/*phaseIdx*/) || !FluidSystem::viscosityIsConstant(0/*phaseIdx*/)) // TODO fix on master | ||
| 605 | { | ||
| 606 |
4/4✓ Branch 1 taken 83170 times.
✓ Branch 2 taken 27400 times.
✓ Branch 3 taken 83170 times.
✓ Branch 4 taken 27400 times.
|
110570 | for (const auto& scvf : scvfs(momentumFvGeometry, scv)) |
| 607 | { | ||
| 608 |
4/4✓ Branch 0 taken 54800 times.
✓ Branch 1 taken 28370 times.
✓ Branch 2 taken 52860 times.
✓ Branch 3 taken 1940 times.
|
83170 | if (scvf.isLateral() && !scvf.boundary()) |
| 609 | { | ||
| 610 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 43160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 43160 times.
|
105720 | const auto& outsideScv = momentumFvGeometry.scv(scvf.outsideScvIdx()); |
| 611 | 105720 | momentumToMassAndEnergyStencils_[scv.index()].push_back(outsideScv.elementIndex()); | |
| 612 | } | ||
| 613 | } | ||
| 614 | } | ||
| 615 | } | ||
| 616 | } | ||
| 617 | 55 | } | |
| 618 | |||
| 619 | 375616 | std::size_t massScvfToMomentumScvIdx_(const SubControlVolumeFace<freeFlowMassIndex>& massScvf, | |
| 620 | [[maybe_unused]] const FVElementGeometry<freeFlowMomentumIndex>& momentumFVGeometry) const | ||
| 621 | { | ||
| 622 | if constexpr (ConsistentlyOrientedGrid<typename GridView<freeFlowMomentumIndex>::Grid>{}) | ||
| 623 | 126695520 | return massScvf.index(); | |
| 624 | else | ||
| 625 | { | ||
| 626 |
4/6✓ Branch 0 taken 7 times.
✓ Branch 1 taken 375609 times.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
|
375616 | static const bool makeConsistentlyOriented = getParam<bool>("Grid.MakeConsistentlyOriented", true); |
| 627 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 375616 times.
|
375616 | if (!makeConsistentlyOriented) |
| 628 | ✗ | return massScvf.index(); | |
| 629 | |||
| 630 |
3/5✓ Branch 0 taken 457120 times.
✓ Branch 1 taken 482000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 482000 times.
✗ Branch 4 not taken.
|
1121968 | for (const auto& momentumScv : scvs(momentumFVGeometry)) |
| 631 | { | ||
| 632 | 939120 | typename SubControlVolumeFace<freeFlowMassIndex>::GlobalPosition momentumUnitOuterNormal(0.0); | |
| 633 | 939120 | momentumUnitOuterNormal[momentumScv.dofAxis()] = momentumScv.directionSign(); | |
| 634 |
4/4✓ Branch 0 taken 375696 times.
✓ Branch 1 taken 563424 times.
✓ Branch 2 taken 375616 times.
✓ Branch 3 taken 563504 times.
|
3193056 | if (Dune::FloatCmp::eq<typename GridView<freeFlowMomentumIndex>::ctype>(massScvf.unitOuterNormal()*momentumUnitOuterNormal, 1.0)) |
| 635 | 375616 | return momentumScv.index(); | |
| 636 | } | ||
| 637 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No Momentum SCV found"); | |
| 638 | } | ||
| 639 | } | ||
| 640 | |||
| 641 | CouplingStencilType emptyStencil_; | ||
| 642 | std::vector<CouplingStencilType> momentumToMassAndEnergyStencils_; | ||
| 643 | std::vector<CouplingStencilType> massAndEnergyToMomentumStencils_; | ||
| 644 | |||
| 645 | // the coupling context exists for each thread | ||
| 646 | // TODO this is a bad pattern, just like mutable caches | ||
| 647 | // we should really construct and pass the context and not store it globally | ||
| 648 | ✗ | std::vector<MomentumCouplingContext>& momentumCouplingContext_() const | |
| 649 | { | ||
| 650 | ✗ | thread_local static std::vector<MomentumCouplingContext> c; | |
| 651 | ✗ | return c; | |
| 652 | } | ||
| 653 | |||
| 654 | // the coupling context exists for each thread | ||
| 655 | ✗ | std::vector<MassAndEnergyCouplingContext>& massAndEnergyCouplingContext_() const | |
| 656 | { | ||
| 657 | ✗ | thread_local static std::vector<MassAndEnergyCouplingContext> c; | |
| 658 | ✗ | return c; | |
| 659 | } | ||
| 660 | |||
| 661 | //! A tuple of std::shared_ptrs to the grid variables of the sub problems | ||
| 662 | GridVariablesTuple gridVariables_; | ||
| 663 | |||
| 664 | const SolutionVector* prevSol_; | ||
| 665 | bool isTransient_; | ||
| 666 | |||
| 667 | std::deque<std::vector<ElementSeed<freeFlowMomentumIndex>>> elementSets_; | ||
| 668 | }; | ||
| 669 | |||
| 670 | //! TODO The infrastructure for multithreaded assembly is implemented (see code in the class above) but the current implementation seems to have a bug and may cause race conditions. The result is different when running in parallel. After this has been fixed activate multithreaded assembly by inheriting from std::true_type here. | ||
| 671 | template<class T> | ||
| 672 | struct CouplingManagerSupportsMultithreadedAssembly<FCStaggeredFreeFlowCouplingManager<T>> | ||
| 673 | : public std::false_type {}; | ||
| 674 | |||
| 675 | } // end namespace Dumux | ||
| 676 | |||
| 677 | #endif | ||
| 678 |