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