Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "domainpostoperator.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "fem/bilinearform.hpp"
8 : #include "fem/fespace.hpp"
9 : #include "fem/gridfunction.hpp"
10 : #include "fem/integrator.hpp"
11 : #include "models/materialoperator.hpp"
12 : #include "utils/communication.hpp"
13 : #include "utils/iodata.hpp"
14 :
15 : namespace palace
16 : {
17 :
18 9 : DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOperator &mat_op,
19 : const FiniteElementSpace &nd_fespace,
20 9 : const FiniteElementSpace &rt_fespace)
21 : {
22 : // Mass operators are always partially assembled.
23 9 : MFEM_VERIFY(nd_fespace.GetFEColl().GetMapType(nd_fespace.Dimension()) ==
24 : mfem::FiniteElement::H_CURL &&
25 : rt_fespace.GetFEColl().GetMapType(nd_fespace.Dimension()) ==
26 : mfem::FiniteElement::H_DIV,
27 : "Unexpected finite element space types for domain energy postprocessing!");
28 : {
29 : // Construct ND mass matrix to compute the electric field energy integral as:
30 : // E_elec = 1/2 Re{∫_Ω Dᴴ E dV} as (M_eps * e)ᴴ e.
31 : // Only the real part of the permeability contributes to the energy (imaginary part
32 : // cancels out in the inner product due to symmetry).
33 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
34 9 : mat_op.GetPermittivityReal());
35 : BilinearForm m(nd_fespace);
36 9 : m.AddDomainIntegrator<VectorFEMassIntegrator>(epsilon_func);
37 0 : M_elec = m.PartialAssemble();
38 9 : D.SetSize(M_elec->Height());
39 : D.UseDevice(true);
40 9 : }
41 : {
42 : // Construct RT mass matrix to compute the magnetic field energy integral as:
43 : // E_mag = 1/2 Re{∫_Ω Hᴴ B dV} as (M_muinv * b)ᴴ b.
44 : MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
45 9 : mat_op.GetInvPermeability());
46 : BilinearForm m(rt_fespace);
47 9 : m.AddDomainIntegrator<VectorFEMassIntegrator>(muinv_func);
48 0 : M_mag = m.PartialAssemble();
49 9 : H.SetSize(M_mag->Height());
50 : H.UseDevice(true);
51 9 : }
52 :
53 : // Use the provided domain postprocessing indices for postprocessing the electric and
54 : // magnetic field energy in specific regions of the domain.
55 9 : for (const auto &[idx, data] : iodata.domains.postpro.energy)
56 : {
57 : std::unique_ptr<Operator> M_elec_i, M_mag_i;
58 : {
59 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
60 0 : mat_op.GetPermittivityReal());
61 0 : epsilon_func.RestrictCoefficient(mat_op.GetCeedAttributes(data.attributes));
62 : BilinearForm m(nd_fespace);
63 0 : m.AddDomainIntegrator<VectorFEMassIntegrator>(epsilon_func);
64 0 : M_elec_i = m.PartialAssemble();
65 0 : }
66 : {
67 : MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
68 0 : mat_op.GetInvPermeability());
69 0 : muinv_func.RestrictCoefficient(mat_op.GetCeedAttributes(data.attributes));
70 : BilinearForm m(rt_fespace);
71 0 : m.AddDomainIntegrator<VectorFEMassIntegrator>(muinv_func);
72 0 : M_mag_i = m.PartialAssemble();
73 0 : }
74 0 : M_i.emplace(idx, std::make_pair(std::move(M_elec_i), std::move(M_mag_i)));
75 : }
76 9 : }
77 :
78 6 : DomainPostOperator::DomainPostOperator(const IoData &iodata, const MaterialOperator &mat_op,
79 6 : const FiniteElementSpace &fespace)
80 : {
81 6 : const auto map_type = fespace.GetFEColl().GetMapType(fespace.Dimension());
82 6 : if (map_type == mfem::FiniteElement::VALUE)
83 : {
84 : // H1 space for voltage and electric field energy.
85 : {
86 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
87 3 : mat_op.GetPermittivityReal());
88 : BilinearForm m(fespace);
89 3 : m.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
90 0 : M_elec = m.PartialAssemble();
91 3 : D.SetSize(M_elec->Height());
92 : D.UseDevice(true);
93 3 : }
94 :
95 3 : for (const auto &[idx, data] : iodata.domains.postpro.energy)
96 : {
97 : std::unique_ptr<Operator> M_elec_i;
98 : {
99 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
100 0 : mat_op.GetPermittivityReal());
101 0 : epsilon_func.RestrictCoefficient(mat_op.GetCeedAttributes(data.attributes));
102 : BilinearForm m(fespace);
103 0 : m.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
104 0 : M_elec_i = m.PartialAssemble();
105 0 : }
106 0 : M_i.emplace(idx, std::make_pair(std::move(M_elec_i), nullptr));
107 : }
108 : }
109 3 : else if (map_type == mfem::FiniteElement::H_CURL)
110 : {
111 : // H(curl) space for magnetic vector potential and magnetic field energy.
112 : {
113 : MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
114 3 : mat_op.GetInvPermeability());
115 : BilinearForm m(fespace);
116 3 : m.AddDomainIntegrator<CurlCurlIntegrator>(muinv_func);
117 0 : M_mag = m.PartialAssemble();
118 3 : H.SetSize(M_mag->Height());
119 : H.UseDevice(true);
120 3 : }
121 :
122 3 : for (const auto &[idx, data] : iodata.domains.postpro.energy)
123 : {
124 : std::unique_ptr<Operator> M_mag_i;
125 : {
126 : MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
127 0 : mat_op.GetInvPermeability());
128 0 : muinv_func.RestrictCoefficient(mat_op.GetCeedAttributes(data.attributes));
129 : BilinearForm m(fespace);
130 0 : m.AddDomainIntegrator<CurlCurlIntegrator>(muinv_func);
131 0 : M_mag_i = m.PartialAssemble();
132 0 : }
133 0 : M_i.emplace(idx, std::make_pair(nullptr, std::move(M_mag_i)));
134 : }
135 : }
136 : else
137 : {
138 0 : MFEM_ABORT("Unexpected finite element space type for domain energy postprocessing!");
139 : }
140 6 : }
141 :
142 12 : double DomainPostOperator::GetElectricFieldEnergy(const GridFunction &E) const
143 : {
144 12 : if (M_elec)
145 : {
146 12 : M_elec->Mult(E.Real(), D);
147 12 : double dot = linalg::LocalDot(E.Real(), D);
148 12 : if (E.HasImag())
149 : {
150 6 : M_elec->Mult(E.Imag(), D);
151 6 : dot += linalg::LocalDot(E.Imag(), D);
152 : }
153 : Mpi::GlobalSum(1, &dot, E.GetComm());
154 12 : return 0.5 * dot;
155 : }
156 0 : MFEM_ABORT(
157 : "Domain postprocessing is not configured for electric field energy calculation!");
158 : return 0.0;
159 : }
160 :
161 12 : double DomainPostOperator::GetMagneticFieldEnergy(const GridFunction &B) const
162 : {
163 12 : if (M_mag)
164 : {
165 12 : M_mag->Mult(B.Real(), H);
166 12 : double dot = linalg::LocalDot(B.Real(), H);
167 12 : if (B.HasImag())
168 : {
169 6 : M_mag->Mult(B.Imag(), H);
170 6 : dot += linalg::LocalDot(B.Imag(), H);
171 : }
172 : Mpi::GlobalSum(1, &dot, B.GetComm());
173 12 : return 0.5 * dot;
174 : }
175 0 : MFEM_ABORT(
176 : "Domain postprocessing is not configured for magnetic field energy calculation!");
177 : return 0.0;
178 : }
179 :
180 0 : double DomainPostOperator::GetDomainElectricFieldEnergy(int idx,
181 : const GridFunction &E) const
182 : {
183 : // Compute the electric field energy integral for only a portion of the domain.
184 : auto it = M_i.find(idx);
185 0 : MFEM_VERIFY(it != M_i.end(),
186 : "Invalid domain index when postprocessing domain electric field energy!");
187 0 : if (!it->second.first)
188 : {
189 : return 0.0;
190 : }
191 0 : it->second.first->Mult(E.Real(), D);
192 0 : double dot = linalg::LocalDot(E.Real(), D);
193 0 : if (E.HasImag())
194 : {
195 0 : it->second.first->Mult(E.Imag(), D);
196 0 : dot += linalg::LocalDot(E.Imag(), D);
197 : }
198 : Mpi::GlobalSum(1, &dot, E.GetComm());
199 0 : return 0.5 * dot;
200 : }
201 :
202 0 : double DomainPostOperator::GetDomainMagneticFieldEnergy(int idx,
203 : const GridFunction &B) const
204 : {
205 : // Compute the magnetic field energy integral for only a portion of the domain.
206 : auto it = M_i.find(idx);
207 0 : MFEM_VERIFY(it != M_i.end(),
208 : "Invalid domain index when postprocessing domain magnetic field energy!");
209 0 : if (!it->second.second)
210 : {
211 : return 0.0;
212 : }
213 0 : it->second.second->Mult(B.Real(), H);
214 0 : double dot = linalg::LocalDot(B.Real(), H);
215 0 : if (B.HasImag())
216 : {
217 0 : it->second.second->Mult(B.Imag(), H);
218 0 : dot += linalg::LocalDot(B.Imag(), H);
219 : }
220 : Mpi::GlobalSum(1, &dot, B.GetComm());
221 0 : return 0.5 * dot;
222 : }
223 :
224 : } // namespace palace
|
