Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_LINALG_ERROR_ESTIMATOR_HPP
5 : #define PALACE_LINALG_ERROR_ESTIMATOR_HPP
6 :
7 : #include <memory>
8 : #include <mfem.hpp>
9 : #include "fem/errorindicator.hpp"
10 : #include "fem/fespace.hpp"
11 : #include "fem/libceed/operator.hpp"
12 : #include "linalg/ksp.hpp"
13 : #include "linalg/operator.hpp"
14 : #include "linalg/vector.hpp"
15 :
16 : namespace palace
17 : {
18 :
19 : class MaterialPropertyCoefficient;
20 : class MaterialOperator;
21 :
22 : //
23 : // Classes used in the estimation of element-wise solution errors via a global L2 projection
24 : // of a discontinuous flux onto a smooth space (flux recovery).
25 : //
26 :
27 : template <typename VecType>
28 : class TimeDependentFluxErrorEstimator;
29 :
30 : // This solver computes a smooth recovery of a discontinuous flux. The difference between
31 : // this resulting smooth flux and the original non-smooth flux provides a localizable error
32 : // estimate.
33 : template <typename VecType>
34 : class FluxProjector
35 : {
36 : using OperType = typename std::conditional<std::is_same<VecType, ComplexVector>::value,
37 : ComplexOperator, Operator>::type;
38 :
39 : private:
40 : // Operator for the mass matrix inversion.
41 : std::unique_ptr<OperType> Flux, M;
42 :
43 : // Linear solver and preconditioner for the projected linear system.
44 : std::unique_ptr<BaseKspSolver<OperType>> ksp;
45 :
46 : // Workspace object for solver application.
47 : mutable VecType rhs;
48 :
49 : public:
50 : FluxProjector(const MaterialPropertyCoefficient &coeff,
51 : const FiniteElementSpaceHierarchy &smooth_fespaces,
52 : const FiniteElementSpace &rhs_fespace, double tol, int max_it, int print,
53 : bool use_mg);
54 :
55 : void Mult(const VecType &x, VecType &y) const;
56 : };
57 :
58 : // Class used for computing gradient flux error estimate, η_K = || ε Eₕ - D ||_K, where D
59 : // denotes a smooth reconstruction of ε Eₕ = ε ∇Vₕ with continuous normal component.
60 : template <typename VecType = Vector>
61 : class GradFluxErrorEstimator
62 : {
63 : friend class TimeDependentFluxErrorEstimator<VecType>;
64 :
65 : private:
66 : // Finite element spaces used to represent E and the recovered D.
67 : const FiniteElementSpace &nd_fespace, &rt_fespace;
68 :
69 : // Global L2 projection solver.
70 : FluxProjector<VecType> projector;
71 :
72 : // Operator which performs the integration of the flux error on each element.
73 : ceed::Operator integ_op;
74 :
75 : // Temporary vectors for error estimation.
76 : mutable VecType E_gf, D, D_gf;
77 :
78 : public:
79 : GradFluxErrorEstimator(const MaterialOperator &mat_op, FiniteElementSpace &nd_fespace,
80 : FiniteElementSpaceHierarchy &rt_fespaces, double tol, int max_it,
81 : int print, bool use_mg);
82 :
83 : // Compute elemental error indicators given the electric field as a vector of true dofs,
84 : // and fold into an existing indicator. The indicators are nondimensionalized using the
85 : // total field energy.
86 : void AddErrorIndicator(const VecType &E, double Et, ErrorIndicator &indicator) const;
87 : };
88 :
89 : // Class used for computing curl flux error estimate, η_K = || μ⁻¹ Bₕ - H ||_K where H
90 : // denotes a smooth reconstruction of μ⁻¹ Bₕ = μ⁻¹ ∇ × Eₕ with continuous tangential
91 : // component.
92 : template <typename VecType = Vector>
93 : class CurlFluxErrorEstimator
94 : {
95 : friend class TimeDependentFluxErrorEstimator<VecType>;
96 :
97 : private:
98 : // Finite element space used to represent B and the recovered H.
99 : const FiniteElementSpace &rt_fespace, &nd_fespace;
100 :
101 : // Global L2 projection solver.
102 : FluxProjector<VecType> projector;
103 :
104 : // Operator which performs the integration of the flux error on each element.
105 : ceed::Operator integ_op;
106 :
107 : // Temporary vectors for error estimation.
108 : mutable VecType B_gf, H, H_gf;
109 :
110 : public:
111 : CurlFluxErrorEstimator(const MaterialOperator &mat_op, FiniteElementSpace &rt_fespace,
112 : FiniteElementSpaceHierarchy &nd_fespaces, double tol, int max_it,
113 : int print, bool use_mg);
114 :
115 : // Compute elemental error indicators given the magnetic flux density as a vector of true
116 : // dofs, and fold into an existing indicator. The indicators are nondimensionalized using
117 : // the total field energy.
118 : void AddErrorIndicator(const VecType &B, double Et, ErrorIndicator &indicator) const;
119 : };
120 :
121 : // Class used for computing sum of the gradient flux and curl flux error estimates,
122 : // η²_K = || ε Eₕ - D ||²_K + || μ⁻¹ Bₕ - H ||²_K, where D and H denote a smooth
123 : // reconstructions of ε Eₕ = ε ∇Vₕ with continuous normal component and μ⁻¹ Bₕ = μ⁻¹ ∇ × Eₕ
124 : // with continuous tangential component.
125 : template <typename VecType>
126 0 : class TimeDependentFluxErrorEstimator
127 : {
128 : private:
129 : GradFluxErrorEstimator<VecType> grad_estimator;
130 : CurlFluxErrorEstimator<VecType> curl_estimator;
131 :
132 : public:
133 : TimeDependentFluxErrorEstimator(const MaterialOperator &mat_op,
134 : FiniteElementSpaceHierarchy &nd_fespaces,
135 : FiniteElementSpaceHierarchy &rt_fespaces, double tol,
136 : int max_it, int print, bool use_mg);
137 :
138 : // Compute elemental error indicators given the electric field and magnetic flux density
139 : // as a vectors of true dofs, and fold into an existing indicator. The indicators are
140 : // nondimensionalized using the total field energy.
141 : void AddErrorIndicator(const VecType &E, const VecType &B, double Et,
142 : ErrorIndicator &indicator) const;
143 : };
144 :
145 : } // namespace palace
146 :
147 : #endif // PALACE_LINALG_ERROR_ESTIMATOR_HPP
|
