MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
constraints.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 
12 #ifndef MFEM_CONSTRAINED
13 #define MFEM_CONSTRAINED
14 
15 #include "solvers.hpp"
16 #include "blockoperator.hpp"
17 #include "sparsemat.hpp"
18 #include "hypre.hpp"
19 
20 namespace mfem
21 {
22 
23 class FiniteElementSpace;
24 #ifdef MFEM_USE_MPI
25 class ParFiniteElementSpace;
26 #endif
27 
28 /** @brief An abstract class to solve the constrained system \f$ Ax = f \f$
29  subject to the constraint \f$ B x = r \f$.
30 
31  Although implementations may not use the below formulation, for
32  understanding some of its methods and notation you can think of
33  it as solving the saddle-point system
34 
35  ( A B^T ) ( x ) ( f )
36  ( B ) ( lambda ) = ( r )
37 
38  Do not confuse with ConstrainedOperator, which handles only simple
39  pointwise constraints and is not a Solver.
40 
41  The height and width of this object as an IterativeSolver are the same as
42  just the unconstrained operator \f$ A \f$, and the Mult() interface just
43  takes \f$ f \f$ as an argument. You can set \f$ r \f$ with
44  SetConstraintRHS() (it defaults to zero) and get the Lagrange multiplier
45  solution with GetMultiplierSolution().
46 
47  Alternatively, you can use LagrangeSystemMult() to solve the block system
48  shown above.
49 
50  This abstract object unifies this interface so that derived classes can
51  solve whatever linear system makes sense and the interface will provide
52  uniform access to solutions, Lagrange multipliers, etc. */
54 {
55 public:
56 #ifdef MFEM_USE_MPI
57  ConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_);
58 #endif
60 
61  virtual ~ConstrainedSolver() { }
62 
63  virtual void SetOperator(const Operator& op) override { }
64 
65  /** @brief Set the right-hand side r for the constraint B x = r
66 
67  (r defaults to zero if you don't call this) */
68  virtual void SetConstraintRHS(const Vector& r);
69 
70  /** @brief Return the Lagrange multiplier solution in lambda
71 
72  Mult() only gives you x, this provides access to lambda
73 
74  Does not make sense unless you've already solved the constrained
75  system with Mult() or LagrangeSystemMult() */
76  void GetMultiplierSolution(Vector& lambda) const { lambda = multiplier_sol; }
77 
78  /** @brief Solve for \f$ x \f$ given \f$ f \f$.
79 
80  If you want to set \f$ r \f$, call SetConstraintRHS() before this.
81 
82  If you want to get \f$ \lambda \f$, call GetMultiplierSolution() after
83  this.
84 
85  The base class implementation calls LagrangeSystemMult(), so derived
86  classes must implement either this or LagrangeSystemMult() */
87  virtual void Mult(const Vector& f, Vector& x) const override;
88 
89  /** @brief Solve for (x, lambda) given (f, r)
90 
91  The base class implementation calls Mult(), so derived classes
92  must implement either this or Mult() */
93  virtual void LagrangeSystemMult(const Vector& f_and_r,
94  Vector& x_and_lambda) const;
95 
96 protected:
99 
102  mutable Vector workb;
103  mutable Vector workx;
104 
105 private:
106  void Initialize();
107 };
108 
109 
110 /** @brief Perform elimination of a single constraint.
111 
112  See EliminationProjection, EliminationCGSolver
113 
114  This keeps track of primary / secondary tdofs and does small dense block
115  solves to eliminate constraints from a global system.
116 
117  \f$ B_s^{-1} \f$ maps the lagrange space into secondary dofs, while
118  \f$ -B_s^{-1} B_p \f$ maps primary dofs to secondary dofs. */
120 {
121 public:
122  Eliminator(const SparseMatrix& B, const Array<int>& lagrange_dofs,
123  const Array<int>& primary_tdofs,
124  const Array<int>& secondary_tdofs);
125 
126  const Array<int>& LagrangeDofs() const { return lagrange_tdofs; }
127  const Array<int>& PrimaryDofs() const { return primary_tdofs; }
128  const Array<int>& SecondaryDofs() const { return secondary_tdofs; }
129 
130  /// Given primary dofs in in, return secondary dofs in out
131  /// This applies \f$ -B_s^{-1} B_p \f$.
132  void Eliminate(const Vector& in, Vector& out) const;
133 
134  /// Transpose of Eliminate(), applies \f$ -B_p^T B_s^{-T} \f$
135  void EliminateTranspose(const Vector& in, Vector& out) const;
136 
137  /// Maps Lagrange multipliers to secondary dofs, applies \f$ B_s^{-1} \f$
138  void LagrangeSecondary(const Vector& in, Vector& out) const;
139 
140  /// Transpose of LagrangeSecondary()
141  void LagrangeSecondaryTranspose(const Vector& in, Vector& out) const;
142 
143  /// Return \f$ -B_s^{-1} B_p \f$ explicitly assembled in mat
144  void ExplicitAssembly(DenseMatrix& mat) const;
145 
146 private:
147  Array<int> lagrange_tdofs;
148  Array<int> primary_tdofs;
149  Array<int> secondary_tdofs;
150 
151  DenseMatrix Bp;
152  DenseMatrix Bs; // gets inverted in place
153  LUFactors Bsinverse;
154  DenseMatrix BsT; // gets inverted in place
155  LUFactors BsTinverse;
156  Array<int> ipiv;
157  Array<int> ipivT;
158 };
159 
160 
161 /** Collects action of several Eliminator objects to perform elimination of
162  constraints.
163 
164  Works in parallel, but each Eliminator must be processor local, and must
165  operate on disjoint degrees of freedom (ie, the primary and secondary dofs
166  for one Eliminator must have no overlap with any dofs from a different
167  Eliminator). */
169 {
170 public:
171  EliminationProjection(const Operator& A, Array<Eliminator*>& eliminators);
172 
173  void Mult(const Vector& x, Vector& y) const override;
174 
175  void MultTranspose(const Vector& x, Vector& y) const override;
176 
177  /** @brief Assemble this projector as a (processor-local) SparseMatrix.
178 
179  Some day we may also want to try approximate variants. */
180  SparseMatrix * AssembleExact() const;
181 
182  /** Given Lagrange multiplier right-hand-side \f$ g \f$, return
183  \f$ \tilde{g} \f$ */
184  void BuildGTilde(const Vector& g, Vector& gtilde) const;
185 
186  /** After a solve, recover the Lagrange multiplier. */
187  void RecoverMultiplier(const Vector& primalrhs,
188  const Vector& primalvars, Vector& lm) const;
189 
190 private:
191  const Operator& Aop;
192  Array<Eliminator*> eliminators;
193 };
194 
195 
196 #ifdef MFEM_USE_MPI
197 /** @brief Solve constrained system by eliminating the constraint; see
198  ConstrainedSolver
199 
200  Solves the system with the operator \f$ P^T A P + Z_P \f$, where P is
201  EliminationProjection and Z_P is the identity on the eliminated dofs. */
203 {
204 public:
205  /** @brief Constructor, with explicit splitting into primary/secondary dofs.
206 
207  This constructor uses a single elimination block (per processor), which
208  provides the most general algorithm but is also not scalable
209 
210  The secondary_dofs are eliminated from the system in this algorithm,
211  as they can be written in terms of the primary_dofs.
212 
213  Both primary_dofs and secondary_dofs are in the local truedof numbering;
214  All elimination has to be done locally on processor, though the global
215  system can be parallel. */
217  Array<int>& primary_dofs,
218  Array<int>& secondary_dofs);
219 
220  /** @brief Constructor, elimination is by blocks.
221 
222  Each block is eliminated independently; if the blocks are reasonably
223  small this can be reasonably efficient.
224 
225  The nonzeros in B are assumed to be in disjoint rows and columns; the
226  rows are identified with the constraint_rowstarts array, the secondary
227  dofs are assumed to be the first nonzeros in the rows. */
229  Array<int>& constraint_rowstarts);
230 
232 
233  void Mult(const Vector& x, Vector& y) const override;
234 
235  void SetOperator(const Operator& op) override
236  { MFEM_ABORT("Operator cannot be reset!"); }
237 
238 protected:
239  /// Internal utility routine; assembles eliminated matrix explicitly
240  void BuildExplicitOperator();
241 
242  /// Build preconditioner for eliminated system
243  virtual Solver* BuildPreconditioner() const = 0;
244  /// Select krylov solver for eliminated system
245  virtual IterativeSolver* BuildKrylov() const = 0;
246 
252  mutable Solver* prec;
253 };
254 
255 
256 /** EliminationSolver using CG and HypreBoomerAMG */
258 {
259 public:
261  Array<int>& constraint_rowstarts,
262  int dimension_=0, bool reorder_=false) :
263  EliminationSolver(A, B, constraint_rowstarts),
264  dimension(dimension_), reorder(reorder_)
265  { }
266 
267 protected:
268  virtual Solver* BuildPreconditioner() const override
269  {
271  h_prec->SetPrintLevel(0);
272  if (dimension > 0) { h_prec->SetSystemsOptions(dimension, reorder); }
273  return h_prec;
274  }
275 
276  virtual IterativeSolver* BuildKrylov() const override
277  { return new CGSolver(GetComm()); }
278 
279 private:
280  int dimension;
281  bool reorder;
282 };
283 
284 
285 /** @brief Solve constrained system with penalty method; see ConstrainedSolver.
286 
287  Only approximates the solution, better approximation with higher penalty,
288  but with higher penalty the preconditioner is less effective. */
290 {
291 public:
293  double penalty_);
294 
296  double penalty_);
297 
299 
300  void Mult(const Vector& x, Vector& y) const override;
301 
302  void SetOperator(const Operator& op) override
303  { MFEM_ABORT("Operator cannot be reset!"); }
304 
305 protected:
306  void Initialize(HypreParMatrix& A, HypreParMatrix& B);
307 
308  /// Build preconditioner for penalized system
309  virtual Solver* BuildPreconditioner() const = 0;
310  /// Select krylov solver for penalized system
311  virtual IterativeSolver* BuildKrylov() const = 0;
312 
313  double penalty;
317  mutable Solver * prec;
318 };
319 
320 
321 /** Uses CG and a HypreBoomerAMG preconditioner for the penalized system. */
323 {
324 public:
326  int dimension=0, bool reorder=false) :
327  PenaltyConstrainedSolver(A, B, penalty_),
328  dimension_(dimension), reorder_(reorder)
329  { }
330 
332  int dimension=0, bool reorder=false) :
333  PenaltyConstrainedSolver(A, B, penalty_),
334  dimension_(dimension), reorder_(reorder)
335  { }
336 
337 protected:
338  virtual Solver* BuildPreconditioner() const override
339  {
341  h_prec->SetPrintLevel(0);
342  if (dimension_ > 0) { h_prec->SetSystemsOptions(dimension_, reorder_); }
343  return h_prec;
344  }
345 
346  virtual IterativeSolver* BuildKrylov() const override
347  { return new CGSolver(GetComm()); }
348 
349 private:
350  int dimension_;
351  bool reorder_;
352 };
353 
354 #endif
355 
356 /** @brief Solve constrained system by solving original mixed system;
357  see ConstrainedSolver.
358 
359  Solves the saddle-point problem with a block-diagonal preconditioner, with
360  user-provided preconditioner in the top-left block and (by default) an
361  identity matrix in the bottom-right.
362 
363  This is the most general ConstrainedSolver, needing only Operator objects
364  to function. But in general it is not very efficient or scalable. */
366 {
367 public:
368  /// Setup constrained system, with primal_pc a user-provided preconditioner
369  /// for the top-left block.
370 #ifdef MFEM_USE_MPI
371  SchurConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_,
372  Solver& primal_pc_);
373 #endif
374  SchurConstrainedSolver(Operator& A_, Operator& B_, Solver& primal_pc_);
375  virtual ~SchurConstrainedSolver();
376 
377  virtual void LagrangeSystemMult(const Vector& x, Vector& y) const override;
378 
379 protected:
380 #ifdef MFEM_USE_MPI
381  SchurConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_);
382 #endif
384 
388  Solver * primal_pc; // NOT owned
390  Solver * dual_pc; // owned
391 
392 private:
393  void Initialize();
394 };
395 
396 
397 #ifdef MFEM_USE_MPI
398 /** @brief Basic saddle-point solver with assembled blocks (ie, the
399  operators are assembled HypreParMatrix objects.)
400 
401  This uses a block-diagonal preconditioner that approximates
402  \f$ [ A^{-1} 0; 0 (B A^{-1} B^T)^{-1} ] \f$.
403 
404  In the top-left block, we approximate \f$ A^{-1} \f$ with HypreBoomerAMG.
405  In the bottom-right, we approximate \f$ A^{-1} \f$ with the inverse of the
406  diagonal of \f$ A \f$, assemble \f$ B diag(A)^{-1} B^T \f$, and use
407  HypreBoomerAMG on that assembled matrix. */
409 {
410 public:
411  SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix& hA_,
412  HypreParMatrix& hB_, int dimension=0,
413  bool reorder=false);
415 
416 private:
417  HypreParMatrix& hA;
418  HypreParMatrix& hB;
419  HypreParMatrix * schur_mat;
420 };
421 #endif
422 
423 /** @brief Build a matrix constraining normal components to zero.
424 
425  Given a vector space fespace, and the array constrained_att that
426  includes the boundary *attributes* that are constrained to have normal
427  component zero, this returns a SparseMatrix representing the
428  constraints that need to be imposed.
429 
430  Each row of the returned matrix corresponds to a node that is
431  constrained. The rows are arranged in (contiguous) blocks corresponding
432  to a physical constraint; in 3D, a one-row constraint means the node
433  is free to move along a plane, a two-row constraint means it is free
434  to move along a line (e.g. the intersection of two normal-constrained
435  planes), and a three-row constraint is fully constrained (equivalent
436  to MFEM's usual essential boundary conditions).
437 
438  The constraint_rowstarts array is filled in to describe the structure of
439  these constraints, so that (block) constraint k is encoded in rows
440  constraint_rowstarts[k] to constraint_rowstarts[k + 1] - 1, inclusive,
441  of the returned matrix.
442 
443  Constraints are imposed on "true" degrees of freedom, which are different
444  in serial and parallel, so we need different numbering systems for the
445  serial and parallel versions of this function.
446 
447  When two attributes intersect, this version will combine constraints,
448  so in 2D the point at the intersection is fully constrained (ie,
449  fixed in both directions). This is the wrong thing to do if the
450  two boundaries are (close to) parallel at that point.
451 
452  @param[in] fespace A vector finite element space
453  @param[in] constrained_att Boundary attributes to constrain
454  @param[out] constraint_rowstarts The rowstarts for separately
455  eliminated constraints, possible
456  input to EliminationCGSolver
457  @param[in] parallel Indicate that fespace is actually a
458  ParFiniteElementSpace and the numbering
459  in the returned matrix should be based
460  on truedofs.
461 
462  @return a constraint matrix
463 */
465  Array<int>& constrained_att,
466  Array<int>& constraint_rowstarts,
467  bool parallel=false);
468 
469 #ifdef MFEM_USE_MPI
470 /// Parallel wrapper for BuildNormalConstraints
472  Array<int>& constrained_att,
473  Array<int>& constraint_rowstarts);
474 #endif
475 
476 }
477 
478 #endif
HypreParMatrix & hA
void Eliminate(const Vector &in, Vector &out) const
Definition: constraints.cpp:51
Conjugate gradient method.
Definition: solvers.hpp:316
void Mult(const Vector &x, Vector &y) const override
Solve for given .
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:120
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: constraints.hpp:63
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
Solve constrained system by solving original mixed system; see ConstrainedSolver. ...
EliminationSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &primary_dofs, Array< int > &secondary_dofs)
Constructor, with explicit splitting into primary/secondary dofs.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
Definition: constraints.cpp:95
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix &hA_, HypreParMatrix &hB_, int dimension=0, bool reorder=false)
EliminationProjection(const Operator &A, Array< Eliminator * > &eliminators)
Definition: constraints.cpp:86
An abstract class to solve the constrained system subject to the constraint .
Definition: constraints.hpp:53
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual Solver * BuildPreconditioner() const override
Build preconditioner for penalized system.
Basic saddle-point solver with assembled blocks (ie, the operators are assembled HypreParMatrix objec...
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for eliminated system.
HypreParMatrix * h_explicit_operator
Eliminator(const SparseMatrix &B, const Array< int > &lagrange_dofs, const Array< int > &primary_tdofs, const Array< int > &secondary_tdofs)
Definition: constraints.cpp:22
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 ...
virtual Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for penalized system.
Array< Eliminator * > eliminators
ConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_)
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Definition: constraints.cpp:66
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
Solve constrained system by eliminating the constraint; see ConstrainedSolver.
void BuildGTilde(const Vector &g, Vector &gtilde) const
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:41
Perform elimination of a single constraint.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for penalized system.
PenaltyPCGSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_, int dimension=0, bool reorder=false)
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1467
const Array< int > & PrimaryDofs() const
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
EliminationCGSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &constraint_rowstarts, int dimension_=0, bool reorder_=false)
virtual void SetConstraintRHS(const Vector &r)
Set the right-hand side r for the constraint B x = r.
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
Definition: constraints.cpp:78
Solve constrained system with penalty method; see ConstrainedSolver.
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
Definition: constraints.cpp:58
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_)
Abstract base class for iterative solver.
Definition: solvers.hpp:66
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
Definition: constraints.hpp:76
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
virtual void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
SchurConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_, Solver &primal_pc_)
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for eliminated system.
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:711
EliminationProjection * projector
void Mult(const Vector &x, Vector &y) const override
Solve for given .
IterativeSolver * krylov
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, double penalty_, int dimension=0, bool reorder=false)
virtual void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
Vector data type.
Definition: vector.hpp:60
const Array< int > & LagrangeDofs() const
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4466
BlockDiagonalPreconditioner * block_pc
Base class for solvers.
Definition: operator.hpp:648
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
A class to handle Block systems in a matrix-free implementation.
TransposeOperator * tr_B
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
Definition: constraints.cpp:72
const Array< int > & SecondaryDofs() const