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_OPERATOR
13 #define MFEM_OPERATOR
14
15 #include "vector.hpp"
16
17 namespace mfem
18 {
19
20 class ConstrainedOperator;
21 class RectangularConstrainedOperator;
22
23 /// Abstract operator
24 class Operator
25 {
26 protected:
27 int height; ///< Dimension of the output / number of rows in the matrix.
28 int width; ///< Dimension of the input / number of columns in the matrix.
29
30 /// see FormSystemOperator()
31 void FormConstrainedSystemOperator(
32 const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout);
33
34 /// see FormRectangularSystemOperator()
35 void FormRectangularConstrainedSystemOperator(
36 const Array<int> &trial_tdof_list,
37 const Array<int> &test_tdof_list,
38 RectangularConstrainedOperator* &Aout);
39
40 /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
41 @a Pi corresponds to "P", @a Po corresponds to "Rt" */
42 Operator *SetupRAP(const Operator *Pi, const Operator *Po);
43
44 public:
45 /// Defines operator diagonal policy upon elimination of rows and/or columns.
46 enum DiagonalPolicy
47 {
48 DIAG_ZERO, ///< Set the diagonal value to zero
49 DIAG_ONE, ///< Set the diagonal value to one
50 DIAG_KEEP ///< Keep the diagonal value
51 };
52
53 /// Initializes memory for true vectors of linear system
54 void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi,
55 Vector &x, Vector &b,
56 Vector &X, Vector &B) const;
57
58 /// Construct a square Operator with given size s (default 0).
Operator(int s=0)59 explicit Operator(int s = 0) { height = width = s; }
60
61 /** @brief Construct an Operator with the given height (output size) and
62 width (input size). */
Operator(int h,int w)63 Operator(int h, int w) { height = h; width = w; }
64
65 /// Get the height (size of output) of the Operator. Synonym with NumRows().
Height() const66 inline int Height() const { return height; }
67 /** @brief Get the number of rows (size of output) of the Operator. Synonym
68 with Height(). */
NumRows() const69 inline int NumRows() const { return height; }
70
71 /// Get the width (size of input) of the Operator. Synonym with NumCols().
Width() const72 inline int Width() const { return width; }
73 /** @brief Get the number of columns (size of input) of the Operator. Synonym
74 with Width(). */
NumCols() const75 inline int NumCols() const { return width; }
76
77 /// Return the MemoryClass preferred by the Operator.
78 /** This is the MemoryClass that will be used to access the input and output
79 vectors in the Mult() and MultTranspose() methods.
80
81 For example, classes using the MFEM_FORALL macro for implementation can
82 return the value returned by Device::GetMemoryClass().
83
84 The default implementation of this method in class Operator returns
85 MemoryClass::HOST. */
GetMemoryClass() const86 virtual MemoryClass GetMemoryClass() const { return MemoryClass::HOST; }
87
88 /// Operator application: `y=A(x)`.
89 virtual void Mult(const Vector &x, Vector &y) const = 0;
90
91 /** @brief Action of the transpose operator: `y=A^t(x)`. The default behavior
92 in class Operator is to generate an error. */
MultTranspose(const Vector & x,Vector & y) const93 virtual void MultTranspose(const Vector &x, Vector &y) const
94 { mfem_error("Operator::MultTranspose() is not overloaded!"); }
95
96 /** @brief Evaluate the gradient operator at the point @a x. The default
97 behavior in class Operator is to generate an error. */
GetGradient(const Vector & x) const98 virtual Operator &GetGradient(const Vector &x) const
99 {
100 mfem_error("Operator::GetGradient() is not overloaded!");
101 return const_cast<Operator &>(*this);
102 }
103
104 /** @brief Computes the diagonal entries into @a diag. Typically, this
105 operation only makes sense for linear Operator%s. In some cases, only an
106 approximation of the diagonal is computed. */
AssembleDiagonal(Vector & diag) const107 virtual void AssembleDiagonal(Vector &diag) const
108 {
109 MFEM_CONTRACT_VAR(diag);
110 MFEM_ABORT("Not relevant or not implemented for this Operator.");
111 }
112
113 /** @brief Prolongation operator from linear algebra (linear system) vectors,
114 to input vectors for the operator. `NULL` means identity. */
GetProlongation() const115 virtual const Operator *GetProlongation() const { return NULL; }
116 /** @brief Restriction operator from input vectors for the operator to linear
117 algebra (linear system) vectors. `NULL` means identity. */
GetRestriction() const118 virtual const Operator *GetRestriction() const { return NULL; }
119 /** @brief Prolongation operator from linear algebra (linear system) vectors,
120 to output vectors for the operator. `NULL` means identity. */
GetOutputProlongation() const121 virtual const Operator *GetOutputProlongation() const
122 {
123 return GetProlongation(); // Assume square unless specialized
124 }
125 /** @brief Transpose of GetOutputRestriction, directly available in this
126 form to facilitate matrix-free RAP-type operators.
127
128 `NULL` means identity. */
GetOutputRestrictionTranspose() const129 virtual const Operator *GetOutputRestrictionTranspose() const { return NULL; }
130 /** @brief Restriction operator from output vectors for the operator to linear
131 algebra (linear system) vectors. `NULL` means identity. */
GetOutputRestriction() const132 virtual const Operator *GetOutputRestriction() const
133 {
134 return GetRestriction(); // Assume square unless specialized
135 }
136
137 /** @brief Form a constrained linear system using a matrix-free approach.
138
139 Assuming square operator, form the operator linear system `A(X)=B`,
140 corresponding to it and the right-hand side @a b, by applying any
141 necessary transformations such as: parallel assembly, conforming
142 constraints for non-conforming AMR and eliminating boundary conditions.
143 @note Static condensation and hybridization are not supported for general
144 operators (cf. the analogous methods BilinearForm::FormLinearSystem() and
145 ParBilinearForm::FormLinearSystem()).
146
147 The constraints are specified through the prolongation P from
148 GetProlongation(), and restriction R from GetRestriction() methods, which
149 are e.g. available through the (parallel) finite element space of any
150 (parallel) bilinear form operator. We assume that the operator is square,
151 using the same input and output space, so we have: `A(X)=[P^t (*this)
152 P](X)`, `B=P^t(b)`, and `X=R(x)`.
153
154 The vector @a x must contain the essential boundary condition values.
155 These are eliminated through the ConstrainedOperator class and the vector
156 @a X is initialized by setting its essential entries to the boundary
157 conditions and all other entries to zero (@a copy_interior == 0) or
158 copied from @a x (@a copy_interior != 0).
159
160 After solving the system `A(X)=B`, the (finite element) solution @a x can
161 be recovered by calling Operator::RecoverFEMSolution() with the same
162 vectors @a X, @a b, and @a x.
163
164 @note The caller is responsible for destroying the output operator @a A!
165 @note If there are no transformations, @a X simply reuses the data of @a
166 x. */
167 void FormLinearSystem(const Array<int> &ess_tdof_list,
168 Vector &x, Vector &b,
169 Operator* &A, Vector &X, Vector &B,
170 int copy_interior = 0);
171
172 /** @brief Form a column-constrained linear system using a matrix-free approach.
173
174 Form the operator linear system `A(X)=B` corresponding to the operator
175 and the right-hand side @a b, by applying any necessary transformations
176 such as: parallel assembly, conforming constraints for non-conforming AMR
177 and eliminating boundary conditions. @note Static condensation and
178 hybridization are not supported for general operators (cf. the method
179 MixedBilinearForm::FormRectangularLinearSystem())
180
181 The constraints are specified through the input prolongation Pi from
182 GetProlongation(), and output restriction Ro from GetOutputRestriction()
183 methods, which are e.g. available through the (parallel) finite element
184 spaces of any (parallel) mixed bilinear form operator. So we have:
185 `A(X)=[Ro (*this) Pi](X)`, `B=Ro(b)`, and `X=Pi^T(x)`.
186
187 The vector @a x must contain the essential boundary condition values.
188 The "columns" in this operator corresponding to these values are
189 eliminated through the RectangularConstrainedOperator class.
190
191 After solving the system `A(X)=B`, the (finite element) solution @a x can
192 be recovered by calling Operator::RecoverFEMSolution() with the same
193 vectors @a X, @a b, and @a x.
194
195 @note The caller is responsible for destroying the output operator @a A!
196 @note If there are no transformations, @a X simply reuses the data of @a
197 x. */
198 void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
199 const Array<int> &test_tdof_list,
200 Vector &x, Vector &b,
201 Operator* &A, Vector &X, Vector &B);
202
203 /** @brief Reconstruct a solution vector @a x (e.g. a GridFunction) from the
204 solution @a X of a constrained linear system obtained from
205 Operator::FormLinearSystem() or Operator::FormRectangularLinearSystem().
206
207 Call this method after solving a linear system constructed using
208 Operator::FormLinearSystem() to recover the solution as an input vector,
209 @a x, for this Operator (presumably a finite element grid function). This
210 method has identical signature to the analogous method for bilinear
211 forms, though currently @a b is not used in the implementation. */
212 virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
213
214 /** @brief Return in @a A a parallel (on truedofs) version of this square
215 operator.
216
217 This returns the same operator as FormLinearSystem(), but does without
218 the transformations of the right-hand side and initial guess. */
219 void FormSystemOperator(const Array<int> &ess_tdof_list,
220 Operator* &A);
221
222 /** @brief Return in @a A a parallel (on truedofs) version of this
223 rectangular operator (including constraints).
224
225 This returns the same operator as FormRectangularLinearSystem(), but does
226 without the transformations of the right-hand side. */
227 void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
228 const Array<int> &test_tdof_list,
229 Operator* &A);
230
231 /** @brief Return in @a A a parallel (on truedofs) version of this
232 rectangular operator.
233
234 This is similar to FormSystemOperator(), but for dof-to-dof mappings
235 (discrete linear operators), which can also correspond to rectangular
236 matrices. The user should provide specializations of GetProlongation()
237 for the input dofs and GetOutputRestriction() for the output dofs in
238 their Operator implementation that are appropriate for the two spaces the
239 Operator maps between. These are e.g. available through the (parallel)
240 finite element space of any (parallel) bilinear form operator. We have:
241 `A(X)=[Rout (*this) Pin](X)`. */
242 void FormDiscreteOperator(Operator* &A);
243
244 /// Prints operator with input size n and output size m in Matlab format.
245 void PrintMatlab(std::ostream & out, int n = 0, int m = 0) const;
246
247 /// Virtual destructor.
~Operator()248 virtual ~Operator() { }
249
250 /// Enumeration defining IDs for some classes derived from Operator.
251 /** This enumeration is primarily used with class OperatorHandle. */
252 enum Type
253 {
254 ANY_TYPE, ///< ID for the base class Operator, i.e. any type.
255 MFEM_SPARSEMAT, ///< ID for class SparseMatrix.
256 Hypre_ParCSR, ///< ID for class HypreParMatrix.
257 PETSC_MATAIJ, ///< ID for class PetscParMatrix, MATAIJ format.
258 PETSC_MATIS, ///< ID for class PetscParMatrix, MATIS format.
259 PETSC_MATSHELL, ///< ID for class PetscParMatrix, MATSHELL format.
260 PETSC_MATNEST, ///< ID for class PetscParMatrix, MATNEST format.
261 PETSC_MATHYPRE, ///< ID for class PetscParMatrix, MATHYPRE format.
262 PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
263 Complex_Operator, ///< ID for class ComplexOperator.
264 MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
265 Complex_Hypre_ParCSR ///< ID for class ComplexHypreParMatrix.
266 };
267
268 /// Return the type ID of the Operator class.
269 /** This method is intentionally non-virtual, so that it returns the ID of
270 the specific pointer or reference type used when calling this method. If
271 not overridden by derived classes, they will automatically use the type ID
272 of the base Operator class, ANY_TYPE. */
GetType() const273 Type GetType() const { return ANY_TYPE; }
274 };
275
276
277 /// Base abstract class for first order time dependent operators.
278 /** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
279 algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
280 _implicit_ and _explicit_ parts of the operator, respectively. For explicit
281 operators, F(x,k,t) = k, so f(x,t) = G(x,t). */
282 class TimeDependentOperator : public Operator
283 {
284 public:
285 enum Type
286 {
287 EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
288 IMPLICIT, ///< This is the most general type, no assumptions on F and G.
289 HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
290 };
291
292 /// Evaluation mode. See SetEvalMode() for details.
293 enum EvalMode
294 {
295 /** Normal evaluation. */
296 NORMAL,
297 /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
298 first term, f1. */
299 ADDITIVE_TERM_1,
300 /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
301 second term, f2. */
302 ADDITIVE_TERM_2
303 };
304
305 protected:
306 double t; ///< Current time.
307 Type type; ///< Describes the form of the TimeDependentOperator.
308 EvalMode eval_mode; ///< Current evaluation mode.
309
310 public:
311 /** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
312 y have the same dimension @a n. */
TimeDependentOperator(int n=0,double t_=0.0,Type type_=EXPLICIT)313 explicit TimeDependentOperator(int n = 0, double t_ = 0.0,
314 Type type_ = EXPLICIT)
315 : Operator(n) { t = t_; type = type_; eval_mode = NORMAL; }
316
317 /** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
318 dimensions @a w and @a h, respectively. */
TimeDependentOperator(int h,int w,double t_=0.0,Type type_=EXPLICIT)319 TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
320 : Operator(h, w) { t = t_; type = type_; eval_mode = NORMAL; }
321
322 /// Read the currently set time.
GetTime() const323 virtual double GetTime() const { return t; }
324
325 /// Set the current time.
SetTime(const double t_)326 virtual void SetTime(const double t_) { t = t_; }
327
328 /// True if #type is #EXPLICIT.
isExplicit() const329 bool isExplicit() const { return (type == EXPLICIT); }
330 /// True if #type is #IMPLICIT or #HOMOGENEOUS.
isImplicit() const331 bool isImplicit() const { return !isExplicit(); }
332 /// True if #type is #HOMOGENEOUS.
isHomogeneous() const333 bool isHomogeneous() const { return (type == HOMOGENEOUS); }
334
335 /// Return the current evaluation mode. See SetEvalMode() for details.
GetEvalMode() const336 EvalMode GetEvalMode() const { return eval_mode; }
337
338 /// Set the evaluation mode of the time-dependent operator.
339 /** The evaluation mode is a switch that allows time-stepping methods to
340 request evaluation of separate components/terms of the time-dependent
341 operator. For example, IMEX methods typically assume additive split of
342 the operator: f(x,t) = f1(x,t) + f2(x,t) and they rely on the ability to
343 evaluate the two terms separately.
344
345 Generally, setting the evaluation mode should affect the behavior of all
346 evaluation-related methods in the class, such as Mult(), ImplicitSolve(),
347 etc. However, the exact list of methods that need to support a specific
348 mode will depend on the used time-stepping method. */
SetEvalMode(const EvalMode new_eval_mode)349 virtual void SetEvalMode(const EvalMode new_eval_mode)
350 { eval_mode = new_eval_mode; }
351
352 /** @brief Perform the action of the explicit part of the operator, G:
353 @a y = G(@a x, t) where t is the current time.
354
355 Presently, this method is used by some PETSc ODE solvers, for more
356 details, see the PETSc Manual. */
357 virtual void ExplicitMult(const Vector &x, Vector &y) const;
358
359 /** @brief Perform the action of the implicit part of the operator, F:
360 @a y = F(@a x, @a k, t) where t is the current time.
361
362 Presently, this method is used by some PETSc ODE solvers, for more
363 details, see the PETSc Manual.*/
364 virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const;
365
366 /** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
367 k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
368 current time. */
369 virtual void Mult(const Vector &x, Vector &y) const;
370
371 /** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
372 unknown @a k at the current time t.
373
374 For general F and G, the equation for @a k becomes:
375 F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
376
377 The input vector @a x corresponds to time index (or cycle) n, while the
378 currently set time, #t, and the result vector @a k correspond to time
379 index n+1. The time step @a dt corresponds to the time interval between
380 cycles n and n+1.
381
382 This method allows for the abstract implementation of some time
383 integration methods, including diagonal implicit Runge-Kutta (DIRK)
384 methods and the backward Euler method in particular.
385
386 If not re-implemented, this method simply generates an error. */
387 virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
388
389 /** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
390 given @a x, @a k, and the currently set time.
391
392 Presently, this method is used by some PETSc ODE solvers, for more
393 details, see the PETSc Manual. */
394 virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
395 double shift) const;
396
397 /** @brief Return an Operator representing dG/dx at the given point @a x and
398 the currently set time.
399
400 Presently, this method is used by some PETSc ODE solvers, for more
401 details, see the PETSc Manual. */
402 virtual Operator& GetExplicitGradient(const Vector &x) const;
403
404 /** @brief Setup the ODE linear system \f$ A(x,t) = (I - gamma J) \f$ or
405 \f$ A = (M - gamma J) \f$, where \f$ J(x,t) = \frac{df}{dt(x,t)} \f$.
406
407 @param[in] x The state at which \f$A(x,t)\f$ should be evaluated.
408 @param[in] fx The current value of the ODE rhs function, \f$f(x,t)\f$.
409 @param[in] jok Flag indicating if the Jacobian should be updated.
410 @param[out] jcur Flag to signal if the Jacobian was updated.
411 @param[in] gamma The scaled time step value.
412
413 If not re-implemented, this method simply generates an error.
414
415 Presently, this method is used by SUNDIALS ODE solvers, for more
416 details, see the SUNDIALS User Guides. */
417 virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
418 int jok, int *jcur, double gamma);
419
420 /** @brief Solve the ODE linear system \f$ A x = b \f$ as setup by
421 the method SUNImplicitSetup().
422
423 @param[in] b The linear system right-hand side.
424 @param[in,out] x On input, the initial guess. On output, the solution.
425 @param[in] tol Linear solve tolerance.
426
427 If not re-implemented, this method simply generates an error.
428
429 Presently, this method is used by SUNDIALS ODE solvers, for more
430 details, see the SUNDIALS User Guides. */
431 virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
432
433 /** @brief Setup the mass matrix in the ODE system \f$ M y' = f(y,t) \f$ .
434
435 If not re-implemented, this method simply generates an error.
436
437 Presently, this method is used by SUNDIALS ARKStep integrator, for more
438 details, see the ARKode User Guide. */
439 virtual int SUNMassSetup();
440
441 /** @brief Solve the mass matrix linear system \f$ M x = b \f$
442 as setup by the method SUNMassSetup().
443
444 @param[in] b The linear system right-hand side.
445 @param[in,out] x On input, the initial guess. On output, the solution.
446 @param[in] tol Linear solve tolerance.
447
448 If not re-implemented, this method simply generates an error.
449
450 Presently, this method is used by SUNDIALS ARKStep integrator, for more
451 details, see the ARKode User Guide. */
452 virtual int SUNMassSolve(const Vector &b, Vector &x, double tol);
453
454 /** @brief Compute the mass matrix-vector product \f$ v = M x \f$ .
455
456 @param[in] x The vector to multiply.
457 @param[out] v The result of the matrix-vector product.
458
459 If not re-implemented, this method simply generates an error.
460
461 Presently, this method is used by SUNDIALS ARKStep integrator, for more
462 details, see the ARKode User Guide. */
463 virtual int SUNMassMult(const Vector &x, Vector &v);
464
~TimeDependentOperator()465 virtual ~TimeDependentOperator() { }
466 };
467
468
469 /** TimeDependentAdjointOperator is a TimeDependentOperator with Adjoint rate
470 equations to be used with CVODESSolver. */
471 class TimeDependentAdjointOperator : public TimeDependentOperator
472 {
473 public:
474
475 /**
476 \brief The TimedependentAdjointOperator extends the TimeDependentOperator
477 class to use features in SUNDIALS CVODESSolver for computing quadratures
478 and solving adjoint problems.
479
480 To solve adjoint problems one needs to implement the AdjointRateMult
481 method to tell CVODES what the adjoint rate equation is.
482
483 QuadratureIntegration (optional) can be used to compute values over the
484 forward problem
485
486 QuadratureSensitivityMult (optional) can be used to find the sensitivity
487 of the quadrature using the adjoint solution in part.
488
489 SUNImplicitSetupB (optional) can be used to setup custom solvers for the
490 newton solve for the adjoint problem.
491
492 SUNImplicitSolveB (optional) actually uses the solvers from
493 SUNImplicitSetupB to solve the adjoint problem.
494
495 See SUNDIALS user manuals for specifics.
496
497 \param[in] dim Dimension of the forward operator
498 \param[in] adjdim Dimension of the adjoint operator. Typically it is the
499 same size as dim. However, SUNDIALS allows users to specify the size if
500 one wants to perform custom operations.
501 \param[in] t Starting time to set
502 \param[in] type The TimeDependentOperator type
503 */
TimeDependentAdjointOperator(int dim,int adjdim,double t=0.,Type type=EXPLICIT)504 TimeDependentAdjointOperator(int dim, int adjdim, double t = 0.,
505 Type type = EXPLICIT) :
506 TimeDependentOperator(dim, t, type),
507 adjoint_height(adjdim)
508 {}
509
510 /// Destructor
~TimeDependentAdjointOperator()511 virtual ~TimeDependentAdjointOperator() {};
512
513 /**
514 \brief Provide the operator integration of a quadrature equation
515
516 \param[in] y The current value at time t
517 \param[out] qdot The current quadrature rate value at t
518 */
QuadratureIntegration(const Vector & y,Vector & qdot) const519 virtual void QuadratureIntegration(const Vector &y, Vector &qdot) const {};
520
521 /** @brief Perform the action of the operator:
522 @a yBdot = k = f(@a y,@2 yB, t), where
523
524 @param[in] y The primal solution at time t
525 @param[in] yB The adjoint solution at time t
526 @param[out] yBdot the rate at time t
527 */
528 virtual void AdjointRateMult(const Vector &y, Vector & yB,
529 Vector &yBdot) const = 0;
530
531 /**
532 \brief Provides the sensitivity of the quadrature w.r.t to primal and
533 adjoint solutions
534
535 \param[in] y the value of the primal solution at time t
536 \param[in] yB the value of the adjoint solution at time t
537 \param[out] qBdot the value of the sensitivity of the quadrature rate at
538 time t
539 */
QuadratureSensitivityMult(const Vector & y,const Vector & yB,Vector & qBdot) const540 virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
541 Vector &qBdot) const {}
542
543 /** @brief Setup the ODE linear system \f$ A(x,t) = (I - gamma J) \f$ or
544 \f$ A = (M - gamma J) \f$, where \f$ J(x,t) = \frac{df}{dt(x,t)} \f$.
545
546 @param[in] t The current time
547 @param[in] x The state at which \f$A(x,xB,t)\f$ should be evaluated.
548 @param[in] xB The state at which \f$A(x,xB,t)\f$ should be evaluated.
549 @param[in] fxB The current value of the ODE rhs function, \f$f(x,t)\f$.
550 @param[in] jokB Flag indicating if the Jacobian should be updated.
551 @param[out] jcurB Flag to signal if the Jacobian was updated.
552 @param[in] gammaB The scaled time step value.
553
554 If not re-implemented, this method simply generates an error.
555
556 Presently, this method is used by SUNDIALS ODE solvers, for more details,
557 see the SUNDIALS User Guides.
558 */
SUNImplicitSetupB(const double t,const Vector & x,const Vector & xB,const Vector & fxB,int jokB,int * jcurB,double gammaB)559 virtual int SUNImplicitSetupB(const double t, const Vector &x,
560 const Vector &xB, const Vector &fxB,
561 int jokB, int *jcurB, double gammaB)
562 {
563 mfem_error("TimeDependentAdjointOperator::SUNImplicitSetupB() is not "
564 "overridden!");
565 return (-1);
566 }
567
568 /** @brief Solve the ODE linear system \f$ A(x,xB,t) xB = b \f$ as setup by
569 the method SUNImplicitSetup().
570
571 @param[in] b The linear system right-hand side.
572 @param[in,out] x On input, the initial guess. On output, the solution.
573 @param[in] tol Linear solve tolerance.
574
575 If not re-implemented, this method simply generates an error.
576
577 Presently, this method is used by SUNDIALS ODE solvers, for more details,
578 see the SUNDIALS User Guides. */
SUNImplicitSolveB(Vector & x,const Vector & b,double tol)579 virtual int SUNImplicitSolveB(Vector &x, const Vector &b, double tol)
580 {
581 mfem_error("TimeDependentAdjointOperator::SUNImplicitSolveB() is not "
582 "overridden!");
583 return (-1);
584 }
585
586 /// Returns the size of the adjoint problem state space
GetAdjointHeight()587 int GetAdjointHeight() {return adjoint_height;}
588
589 protected:
590 int adjoint_height; /// Size of the adjoint problem
591 };
592
593
594 /// Base abstract class for second order time dependent operators.
595 /** Operator of the form: (x,dxdt,t) -> f(x,dxdt,t), where k = f(x,dxdt,t)
596 generally solves the algebraic equation F(x,dxdt,k,t) = G(x,dxdt,t).
597 The functions F and G represent the_implicit_ and _explicit_ parts of
598 the operator, respectively. For explicit operators,
599 F(x,dxdt,k,t) = k, so f(x,dxdt,t) = G(x,dxdt,t). */
600 class SecondOrderTimeDependentOperator : public TimeDependentOperator
601 {
602 public:
603 /** @brief Construct a "square" SecondOrderTimeDependentOperator
604 y = f(x,dxdt,t), where x, dxdt and y have the same dimension @a n. */
SecondOrderTimeDependentOperator(int n=0,double t_=0.0,Type type_=EXPLICIT)605 explicit SecondOrderTimeDependentOperator(int n = 0, double t_ = 0.0,
606 Type type_ = EXPLICIT)
607 : TimeDependentOperator(n, t_,type_) { }
608
609 /** @brief Construct a SecondOrderTimeDependentOperator y = f(x,dxdt,t),
610 where x, dxdt and y have the same dimension @a n. */
SecondOrderTimeDependentOperator(int h,int w,double t_=0.0,Type type_=EXPLICIT)611 SecondOrderTimeDependentOperator(int h, int w, double t_ = 0.0,
612 Type type_ = EXPLICIT)
613 : TimeDependentOperator(h, w, t_,type_) { }
614
615 using TimeDependentOperator::Mult;
616
617 /** @brief Perform the action of the operator: @a y = k = f(@a x,@ dxdt, t),
618 where k solves the algebraic equation
619 F(@a x,@ dxdt, k, t) = G(@a x,@ dxdt, t) and t is the current time. */
620 virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const;
621
622 using TimeDependentOperator::ImplicitSolve;
623 /** @brief Solve the equation:
624 @a k = f(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t), for the
625 unknown @a k at the current time t.
626
627 For general F and G, the equation for @a k becomes:
628 F(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t)
629 = G(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t).
630
631 The input vectors @a x and @a dxdt corresponds to time index (or cycle) n, while the
632 currently set time, #t, and the result vector @a k correspond to time
633 index n+1.
634
635 This method allows for the abstract implementation of some time
636 integration methods.
637
638 If not re-implemented, this method simply generates an error. */
639 virtual void ImplicitSolve(const double fac0, const double fac1,
640 const Vector &x, const Vector &dxdt, Vector &k);
641
642
~SecondOrderTimeDependentOperator()643 virtual ~SecondOrderTimeDependentOperator() { }
644 };
645
646
647 /// Base class for solvers
648 class Solver : public Operator
649 {
650 public:
651 /// If true, use the second argument of Mult() as an initial guess.
652 bool iterative_mode;
653
654 /** @brief Initialize a square Solver with size @a s.
655
656 @warning Use a Boolean expression for the second parameter (not an int)
657 to distinguish this call from the general rectangular constructor. */
Solver(int s=0,bool iter_mode=false)658 explicit Solver(int s = 0, bool iter_mode = false)
659 : Operator(s) { iterative_mode = iter_mode; }
660
661 /// Initialize a Solver with height @a h and width @a w.
Solver(int h,int w,bool iter_mode=false)662 Solver(int h, int w, bool iter_mode = false)
663 : Operator(h, w) { iterative_mode = iter_mode; }
664
665 /// Set/update the solver for the given operator.
666 virtual void SetOperator(const Operator &op) = 0;
667 };
668
669
670 /// Identity Operator I: x -> x.
671 class IdentityOperator : public Operator
672 {
673 public:
674 /// Create an identity operator of size @a n.
IdentityOperator(int n)675 explicit IdentityOperator(int n) : Operator(n) { }
676
677 /// Operator application
Mult(const Vector & x,Vector & y) const678 virtual void Mult(const Vector &x, Vector &y) const { y = x; }
679
680 /// Application of the transpose
MultTranspose(const Vector & x,Vector & y) const681 virtual void MultTranspose(const Vector &x, Vector &y) const { y = x; }
682 };
683
684 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or
685 /// an IdentityOperator.
IsIdentityProlongation(const Operator * P)686 inline bool IsIdentityProlongation(const Operator *P)
687 {
688 return !P || dynamic_cast<const IdentityOperator*>(P);
689 }
690
691 /// Scaled Operator B: x -> a A(x).
692 class ScaledOperator : public Operator
693 {
694 private:
695 const Operator &A_;
696 double a_;
697
698 public:
699 /// Create an operator which is a scalar multiple of A.
ScaledOperator(const Operator * A,double a)700 explicit ScaledOperator(const Operator *A, double a)
701 : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
702
703 /// Operator application
Mult(const Vector & x,Vector & y) const704 virtual void Mult(const Vector &x, Vector &y) const
705 { A_.Mult(x, y); y *= a_; }
706 };
707
708
709 /** @brief The transpose of a given operator. Switches the roles of the methods
710 Mult() and MultTranspose(). */
711 class TransposeOperator : public Operator
712 {
713 private:
714 const Operator &A;
715
716 public:
717 /// Construct the transpose of a given operator @a *a.
TransposeOperator(const Operator * a)718 TransposeOperator(const Operator *a)
719 : Operator(a->Width(), a->Height()), A(*a) { }
720
721 /// Construct the transpose of a given operator @a a.
TransposeOperator(const Operator & a)722 TransposeOperator(const Operator &a)
723 : Operator(a.Width(), a.Height()), A(a) { }
724
725 /// Operator application. Apply the transpose of the original Operator.
Mult(const Vector & x,Vector & y) const726 virtual void Mult(const Vector &x, Vector &y) const
727 { A.MultTranspose(x, y); }
728
729 /// Application of the transpose. Apply the original Operator.
MultTranspose(const Vector & x,Vector & y) const730 virtual void MultTranspose(const Vector &x, Vector &y) const
731 { A.Mult(x, y); }
732 };
733
734
735 /// General product operator: x -> (A*B)(x) = A(B(x)).
736 class ProductOperator : public Operator
737 {
738 const Operator *A, *B;
739 bool ownA, ownB;
740 mutable Vector z;
741
742 public:
743 ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
744
Mult(const Vector & x,Vector & y) const745 virtual void Mult(const Vector &x, Vector &y) const
746 { B->Mult(x, z); A->Mult(z, y); }
747
MultTranspose(const Vector & x,Vector & y) const748 virtual void MultTranspose(const Vector &x, Vector &y) const
749 { A->MultTranspose(x, z); B->MultTranspose(z, y); }
750
751 virtual ~ProductOperator();
752 };
753
754
755 /// The operator x -> R*A*P*x constructed through the actions of R^T, A and P
756 class RAPOperator : public Operator
757 {
758 private:
759 const Operator & Rt;
760 const Operator & A;
761 const Operator & P;
762 mutable Vector Px;
763 mutable Vector APx;
764 MemoryClass mem_class;
765
766 public:
767 /// Construct the RAP operator given R^T, A and P.
768 RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_);
769
GetMemoryClass() const770 virtual MemoryClass GetMemoryClass() const { return mem_class; }
771
772 /// Operator application.
Mult(const Vector & x,Vector & y) const773 virtual void Mult(const Vector & x, Vector & y) const
774 { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
775
776 /// Approximate diagonal of the RAP Operator.
777 /** Returns the diagonal of A, as returned by its AssembleDiagonal method,
778 multiplied be P^T.
779
780 When P is the FE space prolongation operator on a mesh without hanging
781 nodes and Rt = P, the returned diagonal is exact, as long as the diagonal
782 of A is also exact. */
AssembleDiagonal(Vector & diag) const783 virtual void AssembleDiagonal(Vector &diag) const
784 {
785 A.AssembleDiagonal(APx);
786 P.MultTranspose(APx, diag);
787
788 // TODO: For an AMR mesh, a convergent diagonal can be assembled with
789 // |P^T| APx, where |P^T| has entry-wise absolute values of the conforming
790 // prolongation transpose operator. See BilinearForm::AssembleDiagonal.
791 }
792
793 /// Application of the transpose.
MultTranspose(const Vector & x,Vector & y) const794 virtual void MultTranspose(const Vector & x, Vector & y) const
795 { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
796 };
797
798
799 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
800 class TripleProductOperator : public Operator
801 {
802 const Operator *A;
803 const Operator *B;
804 const Operator *C;
805 bool ownA, ownB, ownC;
806 mutable Vector t1, t2;
807 MemoryClass mem_class;
808
809 public:
810 TripleProductOperator(const Operator *A, const Operator *B,
811 const Operator *C, bool ownA, bool ownB, bool ownC);
812
GetMemoryClass() const813 virtual MemoryClass GetMemoryClass() const { return mem_class; }
814
Mult(const Vector & x,Vector & y) const815 virtual void Mult(const Vector &x, Vector &y) const
816 { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
817
MultTranspose(const Vector & x,Vector & y) const818 virtual void MultTranspose(const Vector &x, Vector &y) const
819 { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
820
821 virtual ~TripleProductOperator();
822 };
823
824
825 /** @brief Square Operator for imposing essential boundary conditions using only
826 the action, Mult(), of a given unconstrained Operator.
827
828 Square operator constrained by fixing certain entries in the solution to
829 given "essential boundary condition" values. This class is used by the
830 general, matrix-free system formulation of Operator::FormLinearSystem.
831
832 Do not confuse with ConstrainedSolver, which despite the name has very
833 different functionality. */
834 class ConstrainedOperator : public Operator
835 {
836 protected:
837 Array<int> constraint_list; ///< List of constrained indices/dofs.
838 Operator *A; ///< The unconstrained Operator.
839 bool own_A; ///< Ownership flag for A.
840 mutable Vector z, w; ///< Auxiliary vectors.
841 MemoryClass mem_class;
842 DiagonalPolicy diag_policy; ///< Diagonal policy for constrained dofs
843
844 public:
845 /** @brief Constructor from a general Operator and a list of essential
846 indices/dofs.
847
848 Specify the unconstrained operator @a *A and a @a list of indices to
849 constrain, i.e. each entry @a list[i] represents an essential dof. If the
850 ownership flag @a own_A is true, the operator @a *A will be destroyed
851 when this object is destroyed. The @a diag_policy determines how the
852 operator sets entries corresponding to essential dofs. */
853 ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false,
854 DiagonalPolicy diag_policy = DIAG_ONE);
855
856 /// Returns the type of memory in which the solution and temporaries are stored.
GetMemoryClass() const857 virtual MemoryClass GetMemoryClass() const { return mem_class; }
858
859 /// Set the diagonal policy for the constrained operator.
SetDiagonalPolicy(const DiagonalPolicy diag_policy_)860 void SetDiagonalPolicy(const DiagonalPolicy diag_policy_)
861 { diag_policy = diag_policy_; }
862
863 /// Diagonal of A, modified according to the used DiagonalPolicy.
864 virtual void AssembleDiagonal(Vector &diag) const;
865
866 /** @brief Eliminate "essential boundary condition" values specified in @a x
867 from the given right-hand side @a b.
868
869 Performs the following steps:
870
871 z = A((0,x_b)); b_i -= z_i; b_b = x_b;
872
873 where the "_b" subscripts denote the essential (boundary) indices/dofs of
874 the vectors, and "_i" -- the rest of the entries. */
875 void EliminateRHS(const Vector &x, Vector &b) const;
876
877 /** @brief Constrained operator action.
878
879 Performs the following steps:
880
881 z = A((x_i,0)); y_i = z_i; y_b = x_b;
882
883 where the "_b" subscripts denote the essential (boundary) indices/dofs of
884 the vectors, and "_i" -- the rest of the entries. */
885 virtual void Mult(const Vector &x, Vector &y) const;
886
887 /// Destructor: destroys the unconstrained Operator, if owned.
~ConstrainedOperator()888 virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
889 };
890
891 /** @brief Rectangular Operator for imposing essential boundary conditions on
892 the input space using only the action, Mult(), of a given unconstrained
893 Operator.
894
895 Rectangular operator constrained by fixing certain entries in the solution
896 to given "essential boundary condition" values. This class is used by the
897 general matrix-free formulation of Operator::FormRectangularLinearSystem. */
898 class RectangularConstrainedOperator : public Operator
899 {
900 protected:
901 Array<int> trial_constraints, test_constraints;
902 Operator *A;
903 bool own_A;
904 mutable Vector z, w;
905 MemoryClass mem_class;
906
907 public:
908 /** @brief Constructor from a general Operator and a list of essential
909 indices/dofs.
910
911 Specify the unconstrained operator @a *A and two lists of indices to
912 constrain, i.e. each entry @a trial_list[i] represents an essential trial
913 dof. If the ownership flag @a own_A is true, the operator @a *A will be
914 destroyed when this object is destroyed. */
915 RectangularConstrainedOperator(Operator *A, const Array<int> &trial_list,
916 const Array<int> &test_list, bool own_A = false);
917 /// Returns the type of memory in which the solution and temporaries are stored.
GetMemoryClass() const918 virtual MemoryClass GetMemoryClass() const { return mem_class; }
919 /** @brief Eliminate columns corresponding to "essential boundary condition"
920 values specified in @a x from the given right-hand side @a b.
921
922 Performs the following steps:
923
924 b -= A((0,x_b));
925 b_j = 0
926
927 where the "_b" subscripts denote the essential (boundary) indices and the
928 "_j" subscript denotes the essential test indices */
929 void EliminateRHS(const Vector &x, Vector &b) const;
930 /** @brief Rectangular-constrained operator action.
931
932 Performs the following steps:
933
934 y = A((x_i,0));
935 y_j = 0
936
937 where the "_i" subscripts denote all the nonessential (boundary) trial
938 indices and the "_j" subscript denotes the essential test indices */
939 virtual void Mult(const Vector &x, Vector &y) const;
940 virtual void MultTranspose(const Vector &x, Vector &y) const;
~RectangularConstrainedOperator()941 virtual ~RectangularConstrainedOperator() { if (own_A) { delete A; } }
942 };
943
944 /** @brief PowerMethod helper class to estimate the largest eigenvalue of an
945 operator using the iterative power method. */
946 class PowerMethod
947 {
948 Vector v1;
949 #ifdef MFEM_USE_MPI
950 MPI_Comm comm;
951 #endif
952
953 public:
954
955 #ifdef MFEM_USE_MPI
PowerMethod()956 PowerMethod() : comm(MPI_COMM_NULL) {}
957 #else
958 PowerMethod() {}
959 #endif
960
961 #ifdef MFEM_USE_MPI
PowerMethod(MPI_Comm comm_)962 PowerMethod(MPI_Comm comm_) : comm(comm_) {}
963 #endif
964
965 /// @brief Returns an estimate of the largest eigenvalue of the operator \p opr
966 /// using the iterative power method.
967 /** \p v0 is being used as the vector for the iterative process and will contain
968 the eigenvector corresponding to the largest eigenvalue after convergence.
969 The maximum number of iterations may set with \p numSteps, the relative
970 tolerance with \p tolerance and the seed of the random initialization of
971 \p v0 with \p seed. If \p seed is 0 \p v0 will not be random-initialized. */
972 double EstimateLargestEigenvalue(Operator& opr, Vector& v0,
973 int numSteps = 10, double tolerance = 1e-8,
974 int seed = 12345);
975 };
976
977 }
978
979 #endif
980