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. */ 53 class ConstrainedSolver : public IterativeSolver 54 { 55 public: 56 #ifdef MFEM_USE_MPI 57 ConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_); 58 #endif 59 ConstrainedSolver(Operator& A_, Operator& B_); 60 ~ConstrainedSolver()61 virtual ~ConstrainedSolver() { } 62 SetOperator(const Operator & op)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() */ GetMultiplierSolution(Vector & lambda) const76 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: 97 Operator& A; 98 Operator& B; 99 100 mutable Vector constraint_rhs; 101 mutable Vector multiplier_sol; 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. */ 119 class Eliminator 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 LagrangeDofs() const126 const Array<int>& LagrangeDofs() const { return lagrange_tdofs; } PrimaryDofs() const127 const Array<int>& PrimaryDofs() const { return primary_tdofs; } SecondaryDofs() const128 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). */ 168 class EliminationProjection : public Operator 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. */ 202 class EliminationSolver : public ConstrainedSolver 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. */ 216 EliminationSolver(HypreParMatrix& A, SparseMatrix& B, 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. */ 228 EliminationSolver(HypreParMatrix& A, SparseMatrix& B, 229 Array<int>& constraint_rowstarts); 230 231 ~EliminationSolver(); 232 233 void Mult(const Vector& x, Vector& y) const override; 234 SetOperator(const Operator & op)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 247 HypreParMatrix& hA; 248 Array<Eliminator*> eliminators; 249 EliminationProjection * projector; 250 HypreParMatrix * h_explicit_operator; 251 mutable IterativeSolver* krylov; 252 mutable Solver* prec; 253 }; 254 255 256 /** EliminationSolver using CG and HypreBoomerAMG */ 257 class EliminationCGSolver : public EliminationSolver 258 { 259 public: EliminationCGSolver(HypreParMatrix & A,SparseMatrix & B,Array<int> & constraint_rowstarts,int dimension_=0,bool reorder_=false)260 EliminationCGSolver(HypreParMatrix& A, SparseMatrix& B, 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: BuildPreconditioner() const268 virtual Solver* BuildPreconditioner() const override 269 { 270 HypreBoomerAMG * h_prec = new HypreBoomerAMG(*h_explicit_operator); 271 h_prec->SetPrintLevel(0); 272 if (dimension > 0) { h_prec->SetSystemsOptions(dimension, reorder); } 273 return h_prec; 274 } 275 BuildKrylov() const276 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. */ 289 class PenaltyConstrainedSolver : public ConstrainedSolver 290 { 291 public: 292 PenaltyConstrainedSolver(HypreParMatrix& A, SparseMatrix& B, 293 double penalty_); 294 295 PenaltyConstrainedSolver(HypreParMatrix& A, HypreParMatrix& B, 296 double penalty_); 297 298 ~PenaltyConstrainedSolver(); 299 300 void Mult(const Vector& x, Vector& y) const override; 301 SetOperator(const Operator & op)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; 314 Operator& constraintB; 315 HypreParMatrix * penalized_mat; 316 mutable IterativeSolver * krylov; 317 mutable Solver * prec; 318 }; 319 320 321 /** Uses CG and a HypreBoomerAMG preconditioner for the penalized system. */ 322 class PenaltyPCGSolver : public PenaltyConstrainedSolver 323 { 324 public: PenaltyPCGSolver(HypreParMatrix & A,SparseMatrix & B,double penalty_,int dimension=0,bool reorder=false)325 PenaltyPCGSolver(HypreParMatrix& A, SparseMatrix& B, double penalty_, 326 int dimension=0, bool reorder=false) : 327 PenaltyConstrainedSolver(A, B, penalty_), 328 dimension_(dimension), reorder_(reorder) 329 { } 330 PenaltyPCGSolver(HypreParMatrix & A,HypreParMatrix & B,double penalty_,int dimension=0,bool reorder=false)331 PenaltyPCGSolver(HypreParMatrix& A, HypreParMatrix& B, double penalty_, 332 int dimension=0, bool reorder=false) : 333 PenaltyConstrainedSolver(A, B, penalty_), 334 dimension_(dimension), reorder_(reorder) 335 { } 336 337 protected: BuildPreconditioner() const338 virtual Solver* BuildPreconditioner() const override 339 { 340 HypreBoomerAMG* h_prec = new HypreBoomerAMG(*penalized_mat); 341 h_prec->SetPrintLevel(0); 342 if (dimension_ > 0) { h_prec->SetSystemsOptions(dimension_, reorder_); } 343 return h_prec; 344 } 345 BuildKrylov() const346 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. */ 365 class SchurConstrainedSolver : public ConstrainedSolver 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 383 SchurConstrainedSolver(Operator& A_, Operator& B_); 384 385 Array<int> offsets; 386 BlockOperator * block_op; // owned 387 TransposeOperator * tr_B; // owned 388 Solver * primal_pc; // NOT owned 389 BlockDiagonalPreconditioner * block_pc; // 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. */ 408 class SchurConstrainedHypreSolver : public SchurConstrainedSolver 409 { 410 public: 411 SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix& hA_, 412 HypreParMatrix& hB_, int dimension=0, 413 bool reorder=false); 414 virtual ~SchurConstrainedHypreSolver(); 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 */ 464 SparseMatrix * BuildNormalConstraints(FiniteElementSpace& fespace, 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 471 SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace& fespace, 472 Array<int>& constrained_att, 473 Array<int>& constraint_rowstarts); 474 #endif 475 476 } 477 478 #endif 479