MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
hdiv_linear_solver.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
15
16namespace mfem
17{
18
19/// Replace x[i] with 1.0/x[i] for all i.
21{
22 const int n = x.Size();
23 real_t *d_x = x.ReadWrite();
24 MFEM_FORALL(i, n, d_x[i] = 1.0/d_x[i]; );
25}
26
27/// Return a new HypreParMatrix with given diagonal entries
29 const ParFiniteElementSpace &fes)
30{
31 const int n = diag.Size();
32
33 SparseMatrix diag_spmat;
34 diag_spmat.OverrideSize(n, n);
35 diag_spmat.GetMemoryI().New(n+1, Device::GetDeviceMemoryType());
38
39 {
40 int *I = diag_spmat.WriteI();
41 int *J = diag_spmat.WriteJ();
42 real_t *A = diag_spmat.WriteData();
43 const real_t *d_diag = diag.Read();
44 MFEM_FORALL(i, n+1, I[i] = i;);
45 MFEM_FORALL(i, n,
46 {
47 J[i] = i;
48 A[i] = d_diag[i];
49 });
50 }
51
52 HYPRE_BigInt global_size = fes.GlobalTrueVSize();
53 HYPRE_BigInt *row_starts = fes.GetTrueDofOffsets();
54 HypreParMatrix D(fes.GetComm(), global_size, row_starts, &diag_spmat);
55 return new HypreParMatrix(D); // make a deep copy
56}
57
59{
60 Mesh *mesh = fes_l2.GetMesh();
61 const FiniteElement *fe = fes_l2.GetFE(0);
62 return MassIntegrator::GetRule(*fe, *fe, *mesh->GetElementTransformation(0));
63}
64
66 ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
67 Coefficient &L_coeff_, Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_,
68 Mode mode_)
69 : minres(mesh.GetComm()),
70 order(fes_rt_.GetMaxElementOrder()),
71 fec_l2(order - 1, mesh.Dimension(), b2, mt),
72 fes_l2(&mesh, &fec_l2),
73 fec_rt(order - 1, mesh.Dimension(), b1, b2),
74 fes_rt(&mesh, &fec_rt),
75 ess_rt_dofs(ess_rt_dofs_),
76 basis_l2(fes_l2_),
77 basis_rt(fes_rt_),
78 convert_map_type(fes_l2_.GetFE(0)->GetMapType() == FiniteElement::VALUE),
79 mass_l2(&fes_l2),
80 mass_rt(&fes_rt),
81 L_coeff(L_coeff_),
82 R_coeff(R_coeff_),
83 mode(mode_),
84 qs(mesh, GetMassIntRule(fes_l2)),
85 W_coeff_qf(qs),
86 W_mix_coeff_qf(qs),
87 W_coeff(W_coeff_qf),
88 W_mix_coeff(W_mix_coeff_qf)
89{
90 // If the user gives zero L coefficient, switch mode to DARCY_ZERO
91 auto *L_const_coeff = dynamic_cast<ConstantCoefficient*>(&L_coeff);
92 zero_l2_block = (L_const_coeff && L_const_coeff->constant == 0.0);
93
94 if (mode == Mode::GRAD_DIV)
95 {
96 MFEM_VERIFY(!zero_l2_block,
97 "Mode::GRAD_DIV incompatible with zero coefficient.");
98 }
99
100 mass_l2.AddDomainIntegrator(new MassIntegrator(W_coeff));
102
103 mass_rt.AddDomainIntegrator(new VectorFEMassIntegrator(&R_coeff));
105
106 D.reset(FormDiscreteDivergenceMatrix(fes_rt, fes_l2, ess_rt_dofs));
107 Dt.reset(D->Transpose());
108
109 // Versions without BCs needed for elimination
110 D_e.reset(FormDiscreteDivergenceMatrix(fes_rt, fes_l2, empty));
111 mass_rt.FormSystemMatrix(empty, R_e);
112
113 offsets.SetSize(3);
114 offsets[0] = 0;
115 offsets[1] = fes_l2.GetTrueVSize();
116 offsets[2] = offsets[1] + fes_rt.GetTrueVSize();
117
118 minres.SetAbsTol(0.0);
119 minres.SetRelTol(1e-12);
120 minres.SetMaxIter(500);
122 minres.iterative_mode = false;
123
124 R_diag.SetSize(fes_rt.GetTrueVSize());
125 L_diag.SetSize(fes_l2.GetTrueVSize());
126
127 S_inv.SetPrintLevel(0);
128
129 if (mode == Mode::DARCY && !zero_l2_block)
130 {
131 ParBilinearForm mass_l2_unweighted(&fes_l2);
132 QuadratureFunction det_J_qf(qs);
133 QuadratureFunctionCoefficient det_J_coeff(det_J_qf);
134 if (convert_map_type)
135 {
136 const auto flags = GeometricFactors::DETERMINANTS;
137 auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
138 det_J_qf = geom->detJ;
139 mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator(det_J_coeff));
140 }
141 else
142 {
143 mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator);
144 }
145 mass_l2_unweighted.SetAssemblyLevel(AssemblyLevel::PARTIAL);
146 mass_l2_unweighted.Assemble();
147 const int n_l2 = fes_l2.GetTrueVSize();
148 L_diag_unweighted.SetSize(n_l2);
149 mass_l2_unweighted.AssembleDiagonal(L_diag_unweighted);
150 }
151
152 Setup();
153}
154
156 ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
157 Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_)
158 : HdivSaddlePointSolver(mesh, fes_rt_, fes_l2_, zero, R_coeff_,
159 ess_rt_dofs_, Mode::DARCY)
160{ }
161
163{
164 const auto flags = GeometricFactors::DETERMINANTS;
165 auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
166
167 if (!zero_l2_block) { L_coeff.Project(W_coeff_qf); }
168 // In "grad-div mode", the transformation matrix is scaled by the coefficient
169 // of the mass and divergence matrices.
170 // In "Darcy mode", the transformation matrix is unweighted.
171 if (mode == Mode::GRAD_DIV) { W_mix_coeff_qf = W_coeff_qf; }
172 else { W_mix_coeff_qf = 1.0; }
173
174 // The transformation matrix has to be "mixed" value and integral map type,
175 // which means that the coefficient has to be scaled like the Jacobian
176 // determinant.
177 if (convert_map_type)
178 {
179 const int n = W_mix_coeff_qf.Size();
180 const real_t *d_detJ = geom->detJ.Read();
181 real_t *d_w_mix = W_mix_coeff_qf.ReadWrite();
182 real_t *d_w = W_coeff_qf.ReadWrite();
183 const bool zero_l2 = zero_l2_block;
184 MFEM_FORALL(i, n,
185 {
186 const real_t detJ = d_detJ[i];
187 if (!zero_l2) { d_w[i] *= detJ*detJ; }
188 d_w_mix[i] *= detJ;
189 });
190 }
191
192 L_inv.reset(new DGMassInverse(fes_l2, W_mix_coeff));
193
194 if (zero_l2_block)
195 {
196 A_11.reset();
197 }
198 else
199 {
200 mass_l2.Assemble();
201 mass_l2.AssembleDiagonal(L_diag);
202 mass_l2.FormSystemMatrix(empty, L);
203
204 A_11.reset(new RAPOperator(*L_inv, *L, *L_inv));
205
206 if (mode == GRAD_DIV)
207 {
208 L_diag_unweighted.SetSize(L_diag.Size());
209
210 BilinearForm mass_l2_mix(&fes_l2);
211 mass_l2_mix.AddDomainIntegrator(new MassIntegrator(W_mix_coeff));
213 mass_l2_mix.Assemble();
214 mass_l2_mix.AssembleDiagonal(L_diag_unweighted);
215 }
216
217 const real_t *d_L_diag_unweighted = L_diag_unweighted.Read();
218 real_t *d_L_diag = L_diag.ReadWrite();
219 MFEM_FORALL(i, L_diag.Size(),
220 {
221 const real_t d = d_L_diag_unweighted[i];
222 d_L_diag[i] /= d*d;
223 });
224 }
225
226 // Reassemble the RT mass operator with the new coefficient
227 mass_rt.Update();
228 mass_rt.Assemble();
229 mass_rt.FormSystemMatrix(ess_rt_dofs, R);
230
231 // Form the updated approximate Schur complement
232 mass_rt.AssembleDiagonal(R_diag);
233
234 // Update the mass RT diagonal for essential DOFs
235 {
236 const int *d_I = ess_rt_dofs.Read();
237 real_t *d_R_diag = R_diag.ReadWrite();
238 MFEM_FORALL(i, ess_rt_dofs.Size(), d_R_diag[d_I[i]] = 1.0;);
239 }
240
241 // Form the approximate Schur complement
242 {
243 Reciprocal(R_diag);
244 std::unique_ptr<HypreParMatrix> R_diag_inv(MakeDiagonalMatrix(R_diag, fes_rt));
245 if (zero_l2_block)
246 {
247 S.reset(RAP(R_diag_inv.get(), Dt.get()));
248 }
249 else
250 {
251 std::unique_ptr<HypreParMatrix> D_Minv_Dt(RAP(R_diag_inv.get(), Dt.get()));
252 std::unique_ptr<HypreParMatrix> L_diag_inv(MakeDiagonalMatrix(L_diag, fes_l2));
253 S.reset(ParAdd(D_Minv_Dt.get(), L_diag_inv.get()));
254 }
255 }
256
257 // Reassemble the preconditioners
258 R_inv.reset(new OperatorJacobiSmoother(mass_rt, ess_rt_dofs));
259 S_inv.SetOperator(*S);
260
261 // Set up the block operators
262 A_block.reset(new BlockOperator(offsets));
263 // Omit the (1,1)-block when the L coefficient is identically zero.
264 if (A_11) { A_block->SetBlock(0, 0, A_11.get()); }
265 A_block->SetBlock(0, 1, D.get());
266 A_block->SetBlock(1, 0, Dt.get());
267 A_block->SetBlock(1, 1, R.Ptr(), -1.0);
268
269 D_prec.reset(new BlockDiagonalPreconditioner(offsets));
270 D_prec->SetDiagonalBlock(0, &S_inv);
271 D_prec->SetDiagonalBlock(1, R_inv.get());
272
273 minres.SetPreconditioner(*D_prec);
274 minres.SetOperator(*A_block);
275}
276
278{
279 const int n_ess_dofs = ess_rt_dofs.Size();
280 if (fes_l2.GetParMesh()->ReduceInt(n_ess_dofs) == 0) { return; }
281
282 const int n_l2 = offsets[1];
283 const int n_rt = offsets[2]-offsets[1];
284 Vector bE(b, 0, n_l2);
285 Vector bF(b, n_l2, n_rt);
286
287 // SetBC must be called first
288 MFEM_VERIFY(x_bc.Size() == n_rt || n_ess_dofs == 0, "BCs not set");
289
290 // Create a vector z that has the BC values at essential DOFs, zero elsewhere
291 z.SetSize(n_rt);
292 z.UseDevice(true);
293 z = 0.0;
294 const int *d_I = ess_rt_dofs.Read();
295 const real_t *d_x_bc = x_bc.Read();
296 real_t *d_z = z.ReadWrite();
297 MFEM_FORALL(i, n_ess_dofs,
298 {
299 const int j = d_I[i];
300 d_z[j] = d_x_bc[j];
301 });
302
303 // Convert to the IntegratedGLL basis used internally
304 w.SetSize(n_rt);
305 basis_rt.MultInverse(z, w);
306
307 // Eliminate the BCs in the L2 RHS
308 D_e->Mult(-1.0, w, 1.0, bE);
309
310 // Eliminate the BCs in the RT RHS
311 // Flip the sign because the R block appears with multiplier -1
312 z.SetSize(n_rt);
313 R_e->Mult(w, z);
314 bF += z;
315
316 // Insert the RT BCs into the RHS at the essential DOFs.
317 const real_t *d_w = w.Read();
318 real_t *d_bF = bF.ReadWrite(); // Need read-write access to set subvector
319 MFEM_FORALL(i, n_ess_dofs,
320 {
321 const int j = d_I[i];
322 d_bF[j] = -d_w[j];
323 });
324
325 // Make sure the monolithic RHS is updated
326 bE.SyncAliasMemory(b);
327 bF.SyncAliasMemory(b);
328}
329
331{
332 w.SetSize(fes_l2.GetTrueVSize());
333 b_prime.SetSize(b.Size());
334 x_prime.SetSize(x.Size());
335
336 // Transform RHS to the IntegratedGLL basis
337 Vector bE_prime(b_prime, offsets[0], offsets[1]-offsets[0]);
338 Vector bF_prime(b_prime, offsets[1], offsets[2]-offsets[1]);
339
340 const Vector bE(const_cast<Vector&>(b), offsets[0], offsets[1]-offsets[0]);
341 const Vector bF(const_cast<Vector&>(b), offsets[1], offsets[2]-offsets[1]);
342
343 z.SetSize(bE.Size());
344 basis_l2.MultTranspose(bE, z);
345 basis_rt.MultTranspose(bF, bF_prime);
346 // Transform by the inverse of the L2 mass matrix
347 L_inv->Mult(z, bE_prime);
348
349 // Update the monolithic transformed RHS
350 bE_prime.SyncAliasMemory(b_prime);
351 bF_prime.SyncAliasMemory(b_prime);
352
353 // Eliminate the RT essential BCs
354 EliminateBC(b_prime);
355
356 // Solve the transformed system
357 minres.Mult(b_prime, x_prime);
358
359 // Transform the solution back to the user's basis
360 Vector xE_prime(x_prime, offsets[0], offsets[1]-offsets[0]);
361 Vector xF_prime(x_prime, offsets[1], offsets[2]-offsets[1]);
362
363 Vector xE(x, offsets[0], offsets[1]-offsets[0]);
364 Vector xF(x, offsets[1], offsets[2]-offsets[1]);
365
366 z.SetSize(bE.Size()); // Size of z may have changed in EliminateBC
367 L_inv->Mult(xE_prime, z);
368
369 basis_l2.Mult(z, xE);
370 basis_rt.Mult(xF_prime, xF);
371
372 // Update the monolithic solution vector
373 xE.SyncAliasMemory(x);
374 xF.SyncAliasMemory(x);
375}
376
377} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:317
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A class to handle Block diagonal preconditioners in a matrix-free implementation.
A class to handle Block systems in a matrix-free implementation.
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void MultInverse(const Vector &x, Vector &y) const
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
A coefficient that is constant across space and time.
Solver for the discontinuous Galerkin mass matrix.
Definition dgmassinv.hpp:28
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:274
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
Abstract class for all finite elements.
Definition fe_base.hpp:239
Solve the H(div) saddle-point system using MINRES with matrix-free block-diagonal preconditioning.
HdivSaddlePointSolver(ParMesh &mesh_, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_, Coefficient &L_coeff_, Coefficient &R_coeff_, const Array< int > &ess_rt_dofs_, Mode mode_)
Creates a solver for the H(div) saddle-point system.
void Setup()
Build the linear operator and solver. Must be called when the coefficients change.
void EliminateBC(Vector &) const
Eliminates the BCs (called internally, not public interface).
Mode
Which type of saddle-point problem is being solved?
@ DARCY
Darcy/mixed Poisson problem.
void Mult(const Vector &b, Vector &x) const override
Solve the linear system for L2 (scalar) and RT (flux) unknowns.
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition hypre.cpp:5092
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the MINRES method.
Definition solvers.cpp:1616
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:1595
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.hpp:640
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Mesh data type.
Definition mesh.hpp:56
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition mesh.cpp:856
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Class for parallel bilinear form.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a true-dof Vector.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
HYPRE_BigInt * GetTrueDofOffsets() const
Definition pfespace.hpp:282
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
Class for parallel meshes.
Definition pmesh.hpp:34
long long ReduceInt(int value) const override
Utility function: sum integers from all processors (Allreduce).
Definition pmesh.cpp:6307
Quadrature function coefficient which requires that the quadrature rules used for this coefficient be...
Represents values or vectors of values at quadrature points on a mesh.
Definition qfunction.hpp:24
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition qspace.hpp:81
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Definition operator.hpp:817
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Data type sparse matrix.
Definition sparsemat.hpp:51
int * WriteJ(bool on_dev=true)
Memory< int > & GetMemoryI()
int * WriteI(bool on_dev=true)
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Memory< int > & GetMemoryJ()
Memory< real_t > & GetMemoryData()
real_t * WriteData(bool on_dev=true)
Vector data type.
Definition vector.hpp:80
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:490
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition vector.hpp:259
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:136
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
HypreParMatrix * MakeDiagonalMatrix(Vector &diag, const ParFiniteElementSpace &fes)
Return a new HypreParMatrix with given diagonal entries.
HypreParMatrix * FormDiscreteDivergenceMatrix(ParFiniteElementSpace &fes_rt, ParFiniteElementSpace &fes_l2, const Array< int > &ess_dofs)
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A + B.
Definition hypre.cpp:2909
void Reciprocal(Vector &x)
Replace x[i] with 1.0/x[i] for all i.
float real_t
Definition config.hpp:43
const IntegrationRule & GetMassIntRule(FiniteElementSpace &fes_l2)
Settings for the output behavior of the IterativeSolver.
Definition solvers.hpp:79