version 3.10-dev
helmholtzoperator.hh
Go to the documentation of this file.
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//
12#ifndef DUMUX_LINEAR_HELMHOLTZ_OPERATOR_HH
13#define DUMUX_LINEAR_HELMHOLTZ_OPERATOR_HH
14
15#include <dumux/common/math.hh>
22
24
25template <class Traits>
27{
28 using Scalar = typename Traits::PrimaryVariables::value_type;
29public:
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 { priVars_ = elemSol[scv.indexInElement()]; }
35
36 Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; }
37 const PrimaryVariables& priVars() const { return priVars_; }
38 Scalar extrusionFactor() const { return 1.0; }
39
40private:
41 PrimaryVariables priVars_;
42};
43
44template<class TypeTag>
46: public GetPropType<TypeTag, Properties::BaseLocalResidual>
47{
53 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
54 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
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;
62public:
63 using ParentType::ParentType;
64
65 NumEqVector computeStorage(const Problem&,
66 const SubControlVolume&,
67 const VolumeVariables&) const
68 { return NumEqVector(0.0); }
69
70 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 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 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
99 const auto& insideV = elemVolVars[insideScv];
100 const auto& outsideV = elemVolVars[scvf.outsideScvIdx()];
101 const auto ti = computeTpfaTransmissibility(
102 fvGeometry, scvf, insideScv, problem.a(), insideV.extrusionFactor()
103 );
104
105 const auto tij = [&]{
106 if (scvf.boundary())
107 return scvf.area()*ti;
108 else
109 {
110 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
111 const Scalar tj = -computeTpfaTransmissibility(
112 fvGeometry, scvf, outsideScv, problem.a(), outsideV.extrusionFactor()
113 );
114
115 return scvf.area()*(ti * tj)/(ti + tj);
116 }
117 }();
118
119 const auto u = insideV.priVar(0);
120 const auto v = outsideV.priVar(0);
121 flux[0] = tij*(u - v);
122 }
123 else
124 DUNE_THROW(Dune::NotImplemented, "Discretization method " << GridGeometry::discMethod);
125
126 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
140template<class TypeTag>
142{
147public:
148 HelmholtzModelHomogeneousNeumannProblem(std::shared_ptr<const GridGeometry> gridGeometry, Scalar a, Scalar b)
150 , a_(a), b_(b)
151 {}
152
153 template<class GlobalPosition>
154 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
155 { BoundaryTypes values; values.setAllNeumann(); return values; }
156
157 Scalar a() const { return a_; }
158 Scalar b() const { return b_; }
159private:
160 Scalar a_, b_;
161};
162
163template<class G, class D, class S>
164struct Policy
165{
166 using Scalar = S;
167 using Grid = G;
168 using Discretization = D;
169};
170
171template<class P>
172struct TTag {
173 using InheritsFrom = std::tuple<typename P::Discretization>;
174};
175
176} // end namespace Dumux::Detail::HelmholtzOperator
177
178namespace Dumux::Properties {
179
181template<class TypeTag, class P>
182struct Scalar<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
183{ using type = typename P::Scalar; };
184
186template<class TypeTag, class P>
187struct PrimaryVariables<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
188{
189 using type = Dune::FieldVector<
192 >;
193};
194
195template<class TypeTag, class P>
196struct ModelTraits<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
197{
198 struct Traits { static constexpr int numEq() { return 1; } };
199 using type = Traits;
200};
201
202template<class TypeTag, class P>
203struct LocalResidual<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
205
206template<class TypeTag, class P>
207struct Problem<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
209
210template<class TypeTag, class P>
211struct VolumeVariables<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
212{
215};
216
217template<class TypeTag, class P>
218struct EnableGridVolumeVariablesCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
219{ static constexpr bool value = true; };
220
221template<class TypeTag, class P>
222struct EnableGridFluxVariablesCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
223{ static constexpr bool value = true; };
224
225template<class TypeTag, class P>
226struct EnableGridGeometryCache<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
227{ static constexpr bool value = true; };
228
229template<class TypeTag, class P>
230struct Grid<TypeTag, Dumux::Detail::HelmholtzOperator::TTag<P>>
231{ using type = typename P::Grid; };
232
233} // end namespace Dumux::Properties
234
235namespace Dumux {
236
250template<class Discretization, class GridView, class Scalar>
251auto makeHelmholtzMatrix(const GridView& gridView, const Scalar a = 1.0, const Scalar b = 1.0)
252{
253 using TypeTag = Detail::HelmholtzOperator::TTag<
255 >;
257 auto gridGeometry = std::make_shared<GridGeometry>(gridView);
258
260 auto problem = std::make_shared<Problem>(gridGeometry, a, b);
261
264 auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
265 SolutionVector x(gridGeometry->numDofs());
266 gridVariables->init(x);
267
269 auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
270 using A = typename Assembler::JacobianMatrix;
271 using V = typename Assembler::ResidualType;
272 auto jacobian = std::make_shared<A>();
273 auto residual = std::make_shared<V>();
274 assembler->setLinearSystem(jacobian, residual);
275 assembler->assembleJacobian(x);
276 return jacobian;
277};
278
282template<class Discretization, class GridView, class Scalar>
283auto makeLaplaceMatrix(const GridView& gridView, const Scalar a = 1.0)
284{ return makeHelmholtzMatrix<Discretization>(gridView, a, 0.0); };
285
286template<class LinearOperator, class Discretization, class GridView, class Scalar>
287std::shared_ptr<LinearOperator> makeHelmholtzLinearOperator(const GridView& gridView, const Scalar a, const Scalar b)
288{ return std::make_shared<LinearOperator>(makeHelmholtzMatrix<Discretization>(gridView, a, b)); }
289
290template<class LinearOperator, class Discretization, class GridView, class Comm, class Scalar>
291std::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
294template<class LinearOperator, class Discretization, class GridView, class Scalar>
295std::shared_ptr<LinearOperator> makeLaplaceLinearOperator(const GridView& gridView, const Scalar a)
296{ return std::make_shared<LinearOperator>(makeLaplaceMatrix<Discretization>(gridView, a)); }
297
298template<class LinearOperator, class Discretization, class GridView, class Comm, class Scalar>
299std::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
A linear system assembler (residual and Jacobian) for finite volume schemes.
Class to specify the type of a boundary.
Definition: common/boundarytypes.hh:26
void setAllNeumann()
Set all boundary conditions to Neumann.
Definition: common/boundarytypes.hh:90
Scalar b() const
Definition: helmholtzoperator.hh:158
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
Definition: helmholtzoperator.hh:154
Scalar a() const
Definition: helmholtzoperator.hh:157
HelmholtzModelHomogeneousNeumannProblem(std::shared_ptr< const GridGeometry > gridGeometry, Scalar a, Scalar b)
Definition: helmholtzoperator.hh:148
NumEqVector computeStorage(const Problem &, const SubControlVolume &, const VolumeVariables &) const
Definition: helmholtzoperator.hh:65
NumEqVector computeSource(const Problem &problem, const Element &, const FVElementGeometry &, const ElementVolumeVariables &elemVolVars, const SubControlVolume &scv) const
Definition: helmholtzoperator.hh:129
NumEqVector computeFlux(const Problem &problem, const Element &, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const ElementFluxVariablesCache &elemFluxVarsCache) const
Definition: helmholtzoperator.hh:70
typename Traits::PrimaryVariables PrimaryVariables
Definition: helmholtzoperator.hh:30
void update(const ElementSolution &elemSol, const Problem &, const Element &, const SubControlVolume &scv)
Definition: helmholtzoperator.hh:33
const PrimaryVariables & priVars() const
Definition: helmholtzoperator.hh:37
Scalar priVar(const int pvIdx) const
Definition: helmholtzoperator.hh:36
Scalar extrusionFactor() const
Definition: helmholtzoperator.hh:38
A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa,...
Definition: assembly/fvassembler.hh:94
Base class for all finite-volume problems.
Definition: common/fvproblem.hh:43
const GridGeometry & gridGeometry() const
The finite volume grid geometry.
Definition: common/fvproblem.hh:520
Class to specify the type of a boundary.
Base class for all finite volume problems.
Defines all properties used in Dumux.
Tensor::field_type computeTpfaTransmissibility(const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const typename FVElementGeometry::SubControlVolume &scv, const Tensor &T, typename FVElementGeometry::SubControlVolume::Traits::Scalar extrusionFactor)
Free function to evaluate the Tpfa transmissibility associated with the flux (in the form of flux = T...
Definition: tpfa/computetransmissibility.hh:36
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:880
typename NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:34
typename GetProp< TypeTag, Property >::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:296
Define some often used mathematical functions.
Definition: helmholtzoperator.hh:23
constexpr FCDiamond fcdiamond
Definition: method.hh:152
constexpr CCTpfa cctpfa
Definition: method.hh:145
constexpr Box box
Definition: method.hh:147
constexpr PQ1Bubble pq1bubble
Definition: method.hh:148
The energy balance equation for a porous solid.
Definition: common/properties.hh:26
Definition: adapt.hh:17
std::shared_ptr< LinearOperator > makeLaplaceLinearOperator(const GridView &gridView, const Scalar a)
Definition: helmholtzoperator.hh:295
auto makeHelmholtzMatrix(const GridView &gridView, const Scalar a=1.0, const Scalar b=1.0)
make a Helmholtz matrix operator (aΔ + bI)
Definition: helmholtzoperator.hh:251
auto makeLaplaceMatrix(const GridView &gridView, const Scalar a=1.0)
make a Laplace matrix operator aΔ
Definition: helmholtzoperator.hh:283
std::shared_ptr< LinearOperator > makeHelmholtzLinearOperator(const GridView &gridView, const Scalar a, const Scalar b)
Definition: helmholtzoperator.hh:287
A helper to deduce a vector with the same size as numbers of equations.
Definition: helmholtzoperator.hh:165
G Grid
Definition: helmholtzoperator.hh:167
D Discretization
Definition: helmholtzoperator.hh:168
S Scalar
Definition: helmholtzoperator.hh:166
Definition: helmholtzoperator.hh:172
std::tuple< typename P::Discretization > InheritsFrom
Definition: helmholtzoperator.hh:173
typename P::Grid type
Definition: helmholtzoperator.hh:231
Dune::FieldVector< GetPropType< TypeTag, Properties::Scalar >, GetPropType< TypeTag, Properties::ModelTraits >::numEq() > type
Definition: helmholtzoperator.hh:192
typename P::Scalar type
Definition: helmholtzoperator.hh:183
GetPropType< TypeTag, Properties::PrimaryVariables > PrimaryVariables
Definition: helmholtzoperator.hh:213
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...