| 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 Linear | ||
| 10 | * \brief Scalar Helmholtz operator to be used as a solver component | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_LINEAR_HELMHOLTZ_OPERATOR_HH | ||
| 13 | #define DUMUX_LINEAR_HELMHOLTZ_OPERATOR_HH | ||
| 14 | |||
| 15 | #include <dumux/common/math.hh> | ||
| 16 | #include <dumux/common/properties.hh> | ||
| 17 | #include <dumux/common/numeqvector.hh> | ||
| 18 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 19 | #include <dumux/common/boundarytypes.hh> | ||
| 20 | #include <dumux/common/fvproblem.hh> | ||
| 21 | #include <dumux/assembly/fvassembler.hh> | ||
| 22 | |||
| 23 | namespace Dumux::Detail::HelmholtzOperator { | ||
| 24 | |||
| 25 | template <class Traits> | ||
| 26 | 360000 | class HelmholtzModelVolumeVariables | |
| 27 | { | ||
| 28 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 29 | public: | ||
| 30 | using PrimaryVariables = typename Traits::PrimaryVariables; | ||
| 31 | |||
| 32 | template<class ElementSolution, class Problem, class Element, class SubControlVolume> | ||
| 33 | ✗ | void update(const ElementSolution& elemSol, const Problem&, const Element&, const SubControlVolume& scv) | |
| 34 | 360000 | { priVars_ = elemSol[scv.indexInElement()]; } | |
| 35 | |||
| 36 | 5740800 | Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; } | |
| 37 | const PrimaryVariables& priVars() const { return priVars_; } | ||
| 38 | ✗ | Scalar extrusionFactor() const { return 1.0; } | |
| 39 | |||
| 40 | private: | ||
| 41 | PrimaryVariables priVars_; | ||
| 42 | }; | ||
| 43 | |||
| 44 | template<class TypeTag> | ||
| 45 | class HelmholtzModelLocalResidual | ||
| 46 | : public GetPropType<TypeTag, Properties::BaseLocalResidual> | ||
| 47 | { | ||
| 48 | using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>; | ||
| 49 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 50 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 51 | using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 52 | using VolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 53 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 54 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 55 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 56 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 57 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 58 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 59 | using GridView = typename GridGeometry::GridView; | ||
| 60 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 61 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 62 | public: | ||
| 63 |
2/4✓ Branch 1 taken 180000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180000 times.
✗ Branch 5 not taken.
|
360000 | using ParentType::ParentType; |
| 64 | |||
| 65 | ✗ | NumEqVector computeStorage(const Problem&, | |
| 66 | const SubControlVolume&, | ||
| 67 | const VolumeVariables&) const | ||
| 68 | ✗ | { return NumEqVector(0.0); } | |
| 69 | |||
| 70 | 2870400 | NumEqVector computeFlux(const Problem& problem, | |
| 71 | const Element&, const FVElementGeometry& fvGeometry, | ||
| 72 | const ElementVolumeVariables& elemVolVars, | ||
| 73 | const SubControlVolumeFace& scvf, | ||
| 74 | const ElementFluxVariablesCache& elemFluxVarsCache) const | ||
| 75 | { | ||
| 76 | 2870400 | NumEqVector flux(0.0); | |
| 77 | |||
| 78 | // CVFE schemes | ||
| 79 | if constexpr (GridGeometry::discMethod == DiscretizationMethods::box | ||
| 80 | || GridGeometry::discMethod == DiscretizationMethods::fcdiamond | ||
| 81 | || GridGeometry::discMethod == DiscretizationMethods::pq1bubble) | ||
| 82 | { | ||
| 83 | const auto& fluxVarCache = elemFluxVarsCache[scvf]; | ||
| 84 | Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0); | ||
| 85 | for (const auto& scv : scvs(fvGeometry)) | ||
| 86 | { | ||
| 87 | const auto& volVars = elemVolVars[scv]; | ||
| 88 | gradConcentration.axpy(volVars.concentration(), fluxVarCache.gradN(scv.indexInElement())); | ||
| 89 | } | ||
| 90 | |||
| 91 | NumEqVector flux; | ||
| 92 | flux[0] = -1.0*vtmv(scvf.unitOuterNormal(), problem.a(), gradConcentration)*scvf.area(); | ||
| 93 | } | ||
| 94 | |||
| 95 | // CCTpfa | ||
| 96 | else if constexpr (GridGeometry::discMethod == DiscretizationMethods::cctpfa) | ||
| 97 | { | ||
| 98 | 5740800 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 99 | 2870400 | const auto& insideV = elemVolVars[insideScv]; | |
| 100 | 5740800 | const auto& outsideV = elemVolVars[scvf.outsideScvIdx()]; | |
| 101 | 2870400 | const auto ti = computeTpfaTransmissibility( | |
| 102 | fvGeometry, scvf, insideScv, problem.a(), insideV.extrusionFactor() | ||
| 103 | ); | ||
| 104 | |||
| 105 | 2870400 | const auto tij = [&]{ | |
| 106 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2870400 times.
|
2870400 | if (scvf.boundary()) |
| 107 | ✗ | return scvf.area()*ti; | |
| 108 | else | ||
| 109 | { | ||
| 110 | 5740800 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); | |
| 111 | 8611200 | const Scalar tj = -computeTpfaTransmissibility( | |
| 112 | 2870400 | fvGeometry, scvf, outsideScv, problem.a(), outsideV.extrusionFactor() | |
| 113 | ); | ||
| 114 | |||
| 115 | 2870400 | return scvf.area()*(ti * tj)/(ti + tj); | |
| 116 | } | ||
| 117 | 2870400 | }(); | |
| 118 | |||
| 119 | 5740800 | const auto u = insideV.priVar(0); | |
| 120 | 5740800 | const auto v = outsideV.priVar(0); | |
| 121 | 2870400 | flux[0] = tij*(u - v); | |
| 122 | } | ||
| 123 | else | ||
| 124 | DUNE_THROW(Dune::NotImplemented, "Discretization method " << GridGeometry::discMethod); | ||
| 125 | |||
| 126 | 2870400 | return flux; | |
| 127 | } | ||
| 128 | |||
| 129 | ✗ | NumEqVector computeSource(const Problem& problem, | |
| 130 | const Element&, const FVElementGeometry&, | ||
| 131 | const ElementVolumeVariables& elemVolVars, | ||
| 132 | const SubControlVolume& scv) const | ||
| 133 | { | ||
| 134 | ✗ | NumEqVector source(0.0); | |
| 135 | ✗ | source[0] = -problem.b()*elemVolVars[scv].priVar(0); | |
| 136 | ✗ | return source; | |
| 137 | } | ||
| 138 | }; | ||
| 139 | |||
| 140 | template<class TypeTag> | ||
| 141 | 2 | class HelmholtzModelHomogeneousNeumannProblem : public FVProblem<TypeTag> | |
| 142 | { | ||
| 143 | using ParentType = FVProblem<TypeTag>; | ||
| 144 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 145 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 146 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 147 | public: | ||
| 148 | 2 | HelmholtzModelHomogeneousNeumannProblem(std::shared_ptr<const GridGeometry> gridGeometry, Scalar a, Scalar b) | |
| 149 | : ParentType(gridGeometry) | ||
| 150 |
5/14✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
6 | , a_(a), b_(b) |
| 151 | 2 | {} | |
| 152 | |||
| 153 | template<class GlobalPosition> | ||
| 154 | ✗ | BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const | |
| 155 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 4800 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4800 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4800 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4800 times.
|
19200 | { BoundaryTypes values; values.setAllNeumann(); return values; } |
| 156 | |||
| 157 | ✗ | Scalar a() const { return a_; } | |
| 158 | ✗ | Scalar b() const { return b_; } | |
| 159 | private: | ||
| 160 | Scalar a_, b_; | ||
| 161 | }; | ||
| 162 | |||
| 163 | template<class G, class D, class S> | ||
| 164 | struct Policy | ||
| 165 | { | ||
| 166 | using Scalar = S; | ||
| 167 | using Grid = G; | ||
| 168 | using Discretization = D; | ||
| 169 | }; | ||
| 170 | |||
| 171 | template<class P> | ||
| 172 | struct TTag { | ||
| 173 | using InheritsFrom = std::tuple<typename P::Discretization>; | ||
| 174 | }; | ||
| 175 | |||
| 176 | } // end namespace Dumux::Detail::HelmholtzOperator | ||
| 177 | |||
| 178 | namespace Dumux::Properties { | ||
| 179 | |||
| 180 | //! Set the default type of scalar values to double | ||
| 181 | template<class TypeTag, class P> | ||
| 182 | struct Scalar<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 183 | { using type = typename P::Scalar; }; | ||
| 184 | |||
| 185 | //! Set the default primary variable vector to a vector of size of number of equations | ||
| 186 | template<class TypeTag, class P> | ||
| 187 | struct PrimaryVariables<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 188 | { | ||
| 189 | using type = Dune::FieldVector< | ||
| 190 | GetPropType<TypeTag, Properties::Scalar>, | ||
| 191 | GetPropType<TypeTag, Properties::ModelTraits>::numEq() | ||
| 192 | >; | ||
| 193 | }; | ||
| 194 | |||
| 195 | template<class TypeTag, class P> | ||
| 196 | struct ModelTraits<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 197 | { | ||
| 198 | struct Traits { static constexpr int numEq() { return 1; } }; | ||
| 199 | using type = Traits; | ||
| 200 | }; | ||
| 201 | |||
| 202 | template<class TypeTag, class P> | ||
| 203 | struct LocalResidual<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 204 | { using type = Dumux::Detail::HelmholtzOperator::HelmholtzModelLocalResidual<TypeTag>; }; | ||
| 205 | |||
| 206 | template<class TypeTag, class P> | ||
| 207 | struct Problem<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 208 | { using type = Dumux::Detail::HelmholtzOperator::HelmholtzModelHomogeneousNeumannProblem<TypeTag>; }; | ||
| 209 | |||
| 210 | template<class TypeTag, class P> | ||
| 211 | struct VolumeVariables<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 212 | { | ||
| 213 | struct Traits { using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; }; | ||
| 214 | using type = Dumux::Detail::HelmholtzOperator::HelmholtzModelVolumeVariables<Traits>; | ||
| 215 | }; | ||
| 216 | |||
| 217 | template<class TypeTag, class P> | ||
| 218 | struct EnableGridVolumeVariablesCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 219 | { static constexpr bool value = true; }; | ||
| 220 | |||
| 221 | template<class TypeTag, class P> | ||
| 222 | struct EnableGridFluxVariablesCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 223 | { static constexpr bool value = true; }; | ||
| 224 | |||
| 225 | template<class TypeTag, class P> | ||
| 226 | struct EnableGridGeometryCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 227 | { static constexpr bool value = true; }; | ||
| 228 | |||
| 229 | template<class TypeTag, class P> | ||
| 230 | struct Grid<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>> | ||
| 231 | { using type = typename P::Grid; }; | ||
| 232 | |||
| 233 | } // end namespace Dumux::Properties | ||
| 234 | |||
| 235 | namespace Dumux { | ||
| 236 | |||
| 237 | /*! | ||
| 238 | * \brief make a Helmholtz matrix operator (aΔ + bI) | ||
| 239 | * \tparam Discretization the discretization model type tag | ||
| 240 | * \param gridView the gridView on which to assemble the discrete operator | ||
| 241 | * \param a the weight for the Laplacian | ||
| 242 | * \param b the weight for the mass matrix | ||
| 243 | * | ||
| 244 | * For example: a=-1.0 and b=0.0 will yield a negative Laplacian operator | ||
| 245 | * For example: a=0.0 and b=1.0 will yield the mass matrix | ||
| 246 | * For example: a=1.0 and b=1.0 will yield the unweighted Helmholtz operator | ||
| 247 | * | ||
| 248 | * Useful as a preconditioner component or for testing solvers | ||
| 249 | */ | ||
| 250 | template<class Discretization, class GridView, class Scalar> | ||
| 251 | 2 | auto makeHelmholtzMatrix(const GridView& gridView, const Scalar a = 1.0, const Scalar b = 1.0) | |
| 252 | { | ||
| 253 | using TypeTag = Detail::HelmholtzOperator::TTag< | ||
| 254 | Detail::HelmholtzOperator::Policy<typename GridView::Grid, Discretization, Scalar> | ||
| 255 | >; | ||
| 256 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 257 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
2 | auto gridGeometry = std::make_shared<GridGeometry>(gridView); |
| 258 | |||
| 259 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 260 |
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.
|
4 | auto problem = std::make_shared<Problem>(gridGeometry, a, b); |
| 261 | |||
| 262 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 263 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 264 |
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.
|
4 | auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry); |
| 265 |
4/10✓ 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.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
6 | SolutionVector x(gridGeometry->numDofs()); |
| 266 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | gridVariables->init(x); |
| 267 | |||
| 268 | using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>; | ||
| 269 |
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.
|
4 | auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables); |
| 270 | using A = typename Assembler::JacobianMatrix; | ||
| 271 | using V = typename Assembler::ResidualType; | ||
| 272 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | auto jacobian = std::make_shared<A>(); |
| 273 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
|
6 | auto residual = std::make_shared<V>(); |
| 274 |
3/12✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
6 | assembler->setLinearSystem(jacobian, residual); |
| 275 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | assembler->assembleJacobian(x); |
| 276 | 2 | return jacobian; | |
| 277 | }; | ||
| 278 | |||
| 279 | /*! | ||
| 280 | * \brief make a Laplace matrix operator aΔ | ||
| 281 | */ | ||
| 282 | template<class Discretization, class GridView, class Scalar> | ||
| 283 | auto makeLaplaceMatrix(const GridView& gridView, const Scalar a = 1.0) | ||
| 284 | { return createHelmholtzMatrix(gridView, a, 0.0); }; | ||
| 285 | |||
| 286 | template<class LinearOperator, class Discretization, class GridView, class Scalar> | ||
| 287 | 1 | std::shared_ptr<LinearOperator> makeHelmholtzLinearOperator(const GridView& gridView, const Scalar a, const Scalar b) | |
| 288 |
2/6✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | { return std::make_shared<LinearOperator>(makeHelmholtzMatrix<Discretization>(gridView, a, b)); } |
| 289 | |||
| 290 | template<class LinearOperator, class Discretization, class GridView, class Comm, class Scalar> | ||
| 291 | std::shared_ptr<LinearOperator> makeHelmholtzLinearOperator(const GridView& gridView, const Comm& comm, const Scalar a, const Scalar b) | ||
| 292 | { return std::make_shared<LinearOperator>(makeHelmholtzMatrix<Discretization>(gridView, a, b), comm); } | ||
| 293 | |||
| 294 | template<class LinearOperator, class Discretization, class GridView, class Scalar> | ||
| 295 | std::shared_ptr<LinearOperator> makeLaplaceLinearOperator(const GridView& gridView, const Scalar a) | ||
| 296 | { return std::make_shared<LinearOperator>(makeLaplaceMatrix<Discretization>(gridView, a)); } | ||
| 297 | |||
| 298 | template<class LinearOperator, class Discretization, class GridView, class Comm, class Scalar> | ||
| 299 | std::shared_ptr<LinearOperator> makeLaplaceLinearOperator(const GridView& gridView, const Comm& comm, const Scalar a) | ||
| 300 | { return std::make_shared<LinearOperator>(makeLaplaceMatrix<Discretization>(gridView, a), comm); } | ||
| 301 | |||
| 302 | } // end namespace Dumux | ||
| 303 | |||
| 304 | #endif | ||
| 305 |