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