TpsaPreconditioner_impl.hpp
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 Copyright 2025 NORCE AS
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20*/
21#ifndef OPM_TPSA_PRECONDITIONER_IMPL_HPP
22#define OPM_TPSA_PRECONDITIONER_IMPL_HPP
23
24#ifndef OPM_TPSA_PRECONDITIONER_HPP
25#include <config.h>
26#include <opm/simulators/linalg/tpsa/TpsaPreconditioner.hpp>
27#endif
28
29namespace Opm
30{
31
32template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
33TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
34TpsaPreconditioner(const Linear::TpsaMatrixView<Scalar>& S, const PropertyTree& prm)
35 requires (!isParallel)
36 : S_(S)
37{
38 initSubSolvers_(prm);
39}
40
41template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
42TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
43TpsaPreconditioner(const Linear::TpsaMatrixView<Scalar>& S,
44 const PropertyTree& prm,
45 const Comm& comm)
46 requires (isParallel)
47 : S_(S)
48 , comm_(&comm)
49{
50 initSubSolvers_(prm);
51}
52
53template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
54void
55TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
56pre(MultiVector&, MultiVector&)
57{
58}
59
60template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
61void
62TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
63post(MultiVector&)
64{
65}
66
67template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
68Dune::SolverCategory::Category
69TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
70category() const
71{
72 if constexpr (isParallel) {
73 return Dune::SolverCategory::overlapping;
74 } else {
75 return Dune::SolverCategory::sequential;
76 }
77}
78
79template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
80void
81TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
82update()
83{
84 dispSolver0_->preconditioner().update();
85 dispSolver1_->preconditioner().update();
86 dispSolver2_->preconditioner().update();
87 rotSolver_->preconditioner().update();
88 sPresSolver_->preconditioner().update();
89}
90
91template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
92bool
93TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
94hasPerfectUpdate() const
95{
96 return true;
97}
98
99template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
100void
101TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
102apply(MultiVector& v, const MultiVector& d)
103{
104 Dune::InverseOperatorResult result;
105
106 // The defects of the coupled fields are updated as the sweep proceeds,
107 // so they need their own copies.
108 auto d0 = d[_0];
109 auto d1 = d[_1];
110 auto d2 = d[_2];
111 auto d3 = d[_3];
112 auto d4 = d[_4];
113
114 // Ensure that v is zero-initialized for apply()
115 v = 0.0;
116
117 dispSolver0_->apply(v[_0], d0, result);
118 dispSolver1_->apply(v[_1], d1, result);
119 dispSolver2_->apply(v[_2], d2, result);
120
121 S_[_3][_0].mmv(v[_0], d3);
122 S_[_3][_1].mmv(v[_1], d3);
123 S_[_3][_2].mmv(v[_2], d3);
124 rotSolver_->apply(v[_3], d3, result);
125
126 S_[_4][_0].mmv(v[_0], d4);
127 S_[_4][_1].mmv(v[_1], d4);
128 S_[_4][_2].mmv(v[_2], d4);
129 sPresSolver_->apply(v[_4], d4, result);
130}
131
132template <class Scalar, class DispOp, class RotOp, class SPresOp, class Comm>
133void
134TpsaPreconditioner<Scalar, DispOp, RotOp, SPresOp, Comm>::
135initSubSolvers_(const PropertyTree& prm)
136{
137 const auto dispPrm = prm.get_child("disp_disp_solver");
138 const auto rotPrm = prm.get_child("rot_rot_solver");
139 const auto sPresPrm = prm.get_child("spres_spres_solver");
140
141 std::function<Linear::DispVector0T<Scalar>()> dispWeightCalc;
142 std::function<Linear::RotVectorT<Scalar>()> rotWeightCalc;
143 std::function<Linear::SPresVectorT<Scalar>()> sPresWeightCalc;
144
145 if constexpr (isParallel) {
146 dispOp0_ = std::make_unique<DispOp>(S_[_0][_0], *comm_);
147 dispSolver0_ = std::make_unique<DispSolver>(*dispOp0_,
148 *comm_,
149 dispPrm,
150 dispWeightCalc,
151 pressureIdx);
152
153 dispOp1_ = std::make_unique<DispOp>(S_[_1][_1], *comm_);
154 dispSolver1_ = std::make_unique<DispSolver>(*dispOp1_,
155 *comm_,
156 dispPrm,
157 dispWeightCalc,
158 pressureIdx);
159
160 dispOp2_ = std::make_unique<DispOp>(S_[_2][_2], *comm_);
161 dispSolver2_ = std::make_unique<DispSolver>(*dispOp2_,
162 *comm_,
163 dispPrm,
164 dispWeightCalc,
165 pressureIdx);
166
167 rotOp_ = std::make_unique<RotOp>(S_[_3][_3], *comm_);
168 rotSolver_ = std::make_unique<RotSolver>(*rotOp_,
169 *comm_,
170 rotPrm,
171 rotWeightCalc,
172 pressureIdx);
173
174 sPresOp_ = std::make_unique<SPresOp>(S_[_4][_4], *comm_);
175 sPresSolver_ = std::make_unique<SPresSolver>(*sPresOp_,
176 *comm_,
177 sPresPrm,
178 sPresWeightCalc,
179 pressureIdx);
180 } else {
181 dispOp0_ = std::make_unique<DispOp>(S_[_0][_0]);
182 dispSolver0_ = std::make_unique<DispSolver>(*dispOp0_,
183 dispPrm,
184 dispWeightCalc,
185 pressureIdx);
186
187 dispOp1_ = std::make_unique<DispOp>(S_[_1][_1]);
188 dispSolver1_ = std::make_unique<DispSolver>(*dispOp1_,
189 dispPrm,
190 dispWeightCalc,
191 pressureIdx);
192
193 dispOp2_ = std::make_unique<DispOp>(S_[_2][_2]);
194 dispSolver2_ = std::make_unique<DispSolver>(*dispOp2_,
195 dispPrm,
196 dispWeightCalc,
197 pressureIdx);
198
199 rotOp_ = std::make_unique<RotOp>(S_[_3][_3]);
200 rotSolver_ = std::make_unique<RotSolver>(*rotOp_,
201 rotPrm,
202 rotWeightCalc,
203 pressureIdx);
204
205 sPresOp_ = std::make_unique<SPresOp>(S_[_4][_4]);
206 sPresSolver_ = std::make_unique<SPresSolver>(*sPresOp_,
207 sPresPrm,
208 sPresWeightCalc,
209 pressureIdx);
210 }
211}
212
213} // namespace Opm
214
215#endif // OPM_TPSA_PRECONDITIONER_IMPL_HPP
Comm
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:394
Opm::Linear::TpsaMatrixView
Lightweight, non-owning 5x5 view over the sub-matrices owned by TpsaMatrix. Provides the operator int...
Definition: TpsaTypes.hpp:132
Opm::PropertyTree
Hierarchical collection of key/value pairs.
Definition: PropertyTree.hpp:39
Opm::PropertyTree::get_child
PropertyTree get_child(const std::string &key) const
Opm::TpsaPreconditioner
Block lower-triangular preconditioner for the field-split TPSA system.
Definition: TpsaPreconditioner.hpp:106
Opm::TpsaPreconditioner::pre
void pre(MultiVector &, MultiVector &) override
Nothing to prepare, the sub-solvers set themselves up.
Definition: TpsaPreconditioner_impl.hpp:56
Opm::TpsaPreconditioner::apply
void apply(MultiVector &v, const MultiVector &d) override
Apply the preconditioner, i.e. approximately solve S v = d.
Definition: TpsaPreconditioner_impl.hpp:102
Opm::TpsaPreconditioner::category
Dune::SolverCategory::Category category() const override
Solver category, overlapping in parallel and sequential otherwise.
Definition: TpsaPreconditioner_impl.hpp:70
Opm::TpsaPreconditioner::hasPerfectUpdate
bool hasPerfectUpdate() const override
The block solvers are rebuilt from the current matrix on update().
Definition: TpsaPreconditioner_impl.hpp:94
Opm::TpsaPreconditioner::post
void post(MultiVector &) override
Nothing to clean up after the last apply().
Definition: TpsaPreconditioner_impl.hpp:63
Opm::TpsaPreconditioner::TpsaPreconditioner
TpsaPreconditioner(const Linear::TpsaMatrixView< Scalar > &S, const PropertyTree &prm)
Sequential constructor.
Definition: TpsaPreconditioner_impl.hpp:34
Opm::TpsaPreconditioner::update
void update() override
Recompute the preconditioners of all five diagonal block solvers.
Definition: TpsaPreconditioner_impl.hpp:82
Opm
Definition: blackoilbioeffectsmodules.hh:45
Opm::InverseOperatorResult
Dune::InverseOperatorResult InverseOperatorResult
Definition: GpuBridge.hpp:32