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_MODELS_SPACE_OPERATOR_HPP
5 : #define PALACE_MODELS_SPACE_OPERATOR_HPP
6 :
7 : #include <complex>
8 : #include <memory>
9 : #include <vector>
10 : #include <mfem.hpp>
11 : #include "fem/fespace.hpp"
12 : #include "linalg/operator.hpp"
13 : #include "linalg/vector.hpp"
14 : #include "models/farfieldboundaryoperator.hpp"
15 : #include "models/lumpedportoperator.hpp"
16 : #include "models/materialoperator.hpp"
17 : #include "models/portexcitations.hpp"
18 : #include "models/surfaceconductivityoperator.hpp"
19 : #include "models/surfacecurrentoperator.hpp"
20 : #include "models/surfaceimpedanceoperator.hpp"
21 : #include "models/waveportoperator.hpp"
22 :
23 : namespace palace
24 : {
25 :
26 : class IoData;
27 : class Mesh;
28 :
29 : //
30 : // A class handling spatial discretization of the governing equations.
31 : //
32 : class SpaceOperator
33 : {
34 : private:
35 : const bool pc_mat_real; // Use real-valued matrix for preconditioner
36 : const bool pc_mat_shifted; // Use shifted mass matrix for preconditioner
37 :
38 : // Helper variables for log file printing.
39 : bool print_hdr, print_prec_hdr;
40 :
41 : // Perfect electrical conductor essential boundary condition attributes.
42 : mfem::Array<int> dbc_attr, aux_bdr_attr;
43 : std::vector<mfem::Array<int>> nd_dbc_tdof_lists, h1_dbc_tdof_lists, aux_bdr_tdof_lists;
44 :
45 : // Objects defining the finite element spaces for the electric field (Nedelec) and
46 : // magnetic flux density (Raviart-Thomas) on the given mesh. The H1 spaces are used for
47 : // various purposes throughout the code including postprocessing.
48 : std::vector<std::unique_ptr<mfem::ND_FECollection>> nd_fecs;
49 : std::vector<std::unique_ptr<mfem::H1_FECollection>> h1_fecs;
50 : std::vector<std::unique_ptr<mfem::RT_FECollection>> rt_fecs;
51 : FiniteElementSpaceHierarchy nd_fespaces, h1_fespaces, rt_fespaces;
52 :
53 : // Operator for domain material properties.
54 : MaterialOperator mat_op;
55 :
56 : // Operators for boundary conditions and source excitations.
57 : FarfieldBoundaryOperator farfield_op;
58 : SurfaceConductivityOperator surf_sigma_op;
59 : SurfaceImpedanceOperator surf_z_op;
60 : LumpedPortOperator lumped_port_op;
61 : WavePortOperator wave_port_op;
62 : SurfaceCurrentOperator surf_j_op;
63 :
64 : PortExcitations port_excitation_helper;
65 :
66 : mfem::Array<int> SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh);
67 : void CheckBoundaryProperties();
68 :
69 : // Helper functions for building the bilinear forms corresponding to the discretized
70 : // operators in Maxwell's equations.
71 : void AddStiffnessCoefficients(double coeff, MaterialPropertyCoefficient &df,
72 : MaterialPropertyCoefficient &f);
73 : void AddStiffnessBdrCoefficients(double coeff, MaterialPropertyCoefficient &fb);
74 : void AddDampingCoefficients(double coeff, MaterialPropertyCoefficient &f);
75 : void AddDampingBdrCoefficients(double coeff, MaterialPropertyCoefficient &fb);
76 : void AddRealMassCoefficients(double coeff, MaterialPropertyCoefficient &f);
77 : void AddRealMassBdrCoefficients(double coeff, MaterialPropertyCoefficient &fb);
78 : void AddImagMassCoefficients(double coeff, MaterialPropertyCoefficient &f);
79 : void AddAbsMassCoefficients(double coeff, MaterialPropertyCoefficient &f);
80 : void AddExtraSystemBdrCoefficients(double omega, MaterialPropertyCoefficient &dfbr,
81 : MaterialPropertyCoefficient &dfbi,
82 : MaterialPropertyCoefficient &fbr,
83 : MaterialPropertyCoefficient &fbi);
84 : void AddRealPeriodicCoefficients(double coeff, MaterialPropertyCoefficient &f);
85 : void AddImagPeriodicCoefficients(double coeff, MaterialPropertyCoefficient &f);
86 :
87 : // Helper functions for excitation vector assembly.
88 : bool AddExcitationVector1Internal(int excitation_idx, Vector &RHS);
89 : bool AddExcitationVector2Internal(int excitation_idx, double omega, ComplexVector &RHS);
90 :
91 : // Helper functions to build the preconditioner matrix.
92 : void AssemblePreconditioner(std::complex<double> a0, std::complex<double> a1,
93 : std::complex<double> a2, double a3,
94 : std::vector<std::unique_ptr<Operator>> &br_vec,
95 : std::vector<std::unique_ptr<Operator>> &br_aux_vec,
96 : std::vector<std::unique_ptr<Operator>> &bi_vec,
97 : std::vector<std::unique_ptr<Operator>> &bi_aux_vec);
98 : void AssemblePreconditioner(std::complex<double> a0, std::complex<double> a1,
99 : std::complex<double> a2, double a3,
100 : std::vector<std::unique_ptr<Operator>> &br_vec,
101 : std::vector<std::unique_ptr<Operator>> &br_aux_vec);
102 : void AssemblePreconditioner(double a0, double a1, double a2, double a3,
103 : std::vector<std::unique_ptr<Operator>> &br_vec,
104 : std::vector<std::unique_ptr<Operator>> &br_aux_vec);
105 :
106 : public:
107 : SpaceOperator(const IoData &iodata, const std::vector<std::unique_ptr<Mesh>> &mesh);
108 :
109 : // Return list of all PEC boundary true dofs for all finite element space levels.
110 : const std::vector<mfem::Array<int>> &GetNDDbcTDofLists() const
111 : {
112 : return nd_dbc_tdof_lists;
113 : }
114 : const std::vector<mfem::Array<int>> &GetH1DbcTDofLists() const
115 : {
116 : return h1_dbc_tdof_lists;
117 : }
118 :
119 : // Returns lists of all boundary condition true dofs, PEC included, for the auxiliary
120 : // H1 space hierarchy. These are all boundaries which affect the stiffness and damping
121 : // (K and C) matrices, used for nullspace corrections.
122 : const std::vector<mfem::Array<int>> &GetAuxBdrTDofLists() const
123 : {
124 : return aux_bdr_tdof_lists;
125 : }
126 :
127 : // Return material operator for postprocessing.
128 75 : const MaterialOperator &GetMaterialOp() const { return mat_op; }
129 :
130 : // Access to underlying BC operator objects for postprocessing.
131 11 : auto &GetLumpedPortOp() { return lumped_port_op; }
132 0 : auto &GetWavePortOp() { return wave_port_op; }
133 8 : auto &GetSurfaceCurrentOp() { return surf_j_op; }
134 : const auto &GetLumpedPortOp() const { return lumped_port_op; }
135 : const auto &GetWavePortOp() const { return wave_port_op; }
136 : const auto &GetSurfaceCurrentOp() const { return surf_j_op; }
137 :
138 5 : const auto &GetPortExcitations() const { return port_excitation_helper; }
139 :
140 : // Return the parallel finite element space objects.
141 0 : auto &GetNDSpaces() { return nd_fespaces; }
142 : const auto &GetNDSpaces() const { return nd_fespaces; }
143 : auto &GetNDSpace() { return nd_fespaces.GetFinestFESpace(); }
144 : const auto &GetNDSpace() const { return nd_fespaces.GetFinestFESpace(); }
145 0 : auto &GetH1Spaces() { return h1_fespaces; }
146 : const auto &GetH1Spaces() const { return h1_fespaces; }
147 : auto &GetH1Space() { return h1_fespaces.GetFinestFESpace(); }
148 : const auto &GetH1Space() const { return h1_fespaces.GetFinestFESpace(); }
149 0 : auto &GetRTSpaces() { return rt_fespaces; }
150 : const auto &GetRTSpaces() const { return rt_fespaces; }
151 : auto &GetRTSpace() { return rt_fespaces.GetFinestFESpace(); }
152 : const auto &GetRTSpace() const { return rt_fespaces.GetFinestFESpace(); }
153 :
154 : // Access the underlying mesh object.
155 : const auto &GetMesh() const { return GetNDSpace().GetMesh(); }
156 :
157 : // Return the number of true (conforming) dofs on the finest ND space.
158 : auto GlobalTrueVSize() const { return GetNDSpace().GlobalTrueVSize(); }
159 :
160 : // Construct any part of the frequency-dependent complex linear system matrix:
161 : // A = K + iω C - ω² (Mr + i Mi) + A2(ω).
162 : // For time domain problems, any one of K, C, or M = Mr can be constructed. The argument
163 : // ω is required only for the constructing the "extra" matrix A2(ω).
164 : template <typename OperType>
165 : std::unique_ptr<OperType> GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy);
166 : template <typename OperType>
167 : std::unique_ptr<OperType> GetDampingMatrix(Operator::DiagonalPolicy diag_policy);
168 : template <typename OperType>
169 : std::unique_ptr<OperType> GetMassMatrix(Operator::DiagonalPolicy diag_policy);
170 : template <typename OperType>
171 : std::unique_ptr<OperType> GetExtraSystemMatrix(double omega,
172 : Operator::DiagonalPolicy diag_policy);
173 :
174 : // Construct the complete frequency or time domain system matrix using the provided
175 : // stiffness, damping, mass, and extra matrices:
176 : // A = a0 K + a1 C + a2 (Mr + i Mi) + A2.
177 : // It is assumed that the inputs have been constructed using previous calls to
178 : // GetSystemMatrix() and the returned operator does not inherit ownership of any of them.
179 : template <typename OperType, typename ScalarType>
180 : std::unique_ptr<OperType>
181 : GetSystemMatrix(ScalarType a0, ScalarType a1, ScalarType a2, const OperType *K,
182 : const OperType *C, const OperType *M, const OperType *A2 = nullptr);
183 :
184 : // Construct the real, SPD matrix for weighted L2 or H(curl) inner products:
185 : // B = a0 Kr + a2 Mr .
186 : // It is assumed that the inputs have been constructed using previous calls to
187 : // GetSystemMatrix() and the returned operator does not inherit ownership of any of them.
188 : // If K or M have eliminated boundary conditions, they are not eliminated from the
189 : // returned operator.
190 : std::unique_ptr<Operator> GetInnerProductMatrix(double a0, double a2,
191 : const ComplexOperator *K,
192 : const ComplexOperator *M);
193 :
194 : // Construct the matrix for frequency or time domain linear system preconditioning. If it
195 : // is real-valued (Mr > 0, Mi < 0, |Mr + Mi| is done on the material property coefficient,
196 : // not the matrix entries themselves):
197 : // B = a0 K + a1 C -/+ a2 |Mr + Mi| + A2r(a3) + A2i(a3).
198 : template <typename OperType, typename ScalarType>
199 : std::unique_ptr<OperType> GetPreconditionerMatrix(ScalarType a0, ScalarType a1,
200 : ScalarType a2, double a3);
201 :
202 : // Construct and return the discrete curl or gradient matrices.
203 : const Operator &GetGradMatrix() const
204 : {
205 : return GetNDSpace().GetDiscreteInterpolator(GetH1Space());
206 : }
207 : const Operator &GetCurlMatrix() const
208 : {
209 0 : return GetRTSpace().GetDiscreteInterpolator(GetNDSpace());
210 : }
211 :
212 : // Assemble the right-hand side source term vector for an incident field or current source
213 : // applied on specified excited boundaries. The return value indicates whether or not the
214 : // excitation is nonzero (and thus is true most of the time).
215 : bool GetExcitationVector(int excitation_idx, Vector &RHS);
216 : bool GetExcitationVector(int excitation_idx, double omega, ComplexVector &RHS);
217 :
218 : // Separate out RHS vector as RHS = iω RHS1 + RHS2(ω). The return value indicates whether
219 : // or not the excitation is nonzero (and thus is true most of the time).
220 : bool GetExcitationVector1(int excitation_idx, ComplexVector &RHS1);
221 : bool GetExcitationVector2(int excitation_idx, double omega, ComplexVector &RHS2);
222 :
223 : // Construct a constant or randomly initialized vector which satisfies the PEC essential
224 : // boundary conditions.
225 : void GetRandomInitialVector(ComplexVector &v);
226 : void GetConstantInitialVector(ComplexVector &v);
227 :
228 : // Get the associated MPI communicator.
229 : MPI_Comm GetComm() const { return GetNDSpace().GetComm(); }
230 : };
231 :
232 : } // namespace palace
233 :
234 : #endif // PALACE_MODELS_SPACE_OPERATOR_HPP
|
