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 "curlcurloperator.hpp"
5 :
6 : #include <set>
7 : #include "fem/bilinearform.hpp"
8 : #include "fem/coefficient.hpp"
9 : #include "fem/integrator.hpp"
10 : #include "fem/mesh.hpp"
11 : #include "fem/multigrid.hpp"
12 : #include "linalg/hypre.hpp"
13 : #include "linalg/rap.hpp"
14 : #include "utils/communication.hpp"
15 : #include "utils/geodata.hpp"
16 : #include "utils/iodata.hpp"
17 : #include "utils/prettyprint.hpp"
18 :
19 : namespace palace
20 : {
21 :
22 3 : CurlCurlOperator::CurlCurlOperator(const IoData &iodata,
23 3 : const std::vector<std::unique_ptr<Mesh>> &mesh)
24 3 : : print_hdr(true), dbc_attr(SetUpBoundaryProperties(iodata, *mesh.back())),
25 3 : nd_fecs(fem::ConstructFECollections<mfem::ND_FECollection>(
26 3 : iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels,
27 3 : iodata.solver.linear.mg_coarsening, false)),
28 3 : h1_fecs(fem::ConstructFECollections<mfem::H1_FECollection>(
29 3 : iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels,
30 3 : iodata.solver.linear.mg_coarsening, false)),
31 6 : rt_fec(std::make_unique<mfem::RT_FECollection>(iodata.solver.order - 1,
32 3 : mesh.back()->Dimension())),
33 3 : nd_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::ND_FECollection>(
34 3 : iodata.solver.linear.mg_max_levels, mesh, nd_fecs, &dbc_attr, &dbc_tdof_lists)),
35 3 : h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
36 3 : iodata.solver.linear.mg_max_levels, mesh, h1_fecs)),
37 3 : rt_fespace(*mesh.back(), rt_fec.get()), mat_op(iodata, *mesh.back()),
38 3 : surf_j_op(iodata, *mesh.back())
39 : {
40 : // Finalize setup.
41 3 : CheckBoundaryProperties();
42 :
43 : // Print essential BC information.
44 3 : if (dbc_attr.Size())
45 : {
46 3 : Mpi::Print("\nConfiguring Dirichlet BC at attributes:\n");
47 6 : utils::PrettyPrint(dbc_attr);
48 : }
49 3 : }
50 :
51 3 : mfem::Array<int> CurlCurlOperator::SetUpBoundaryProperties(const IoData &iodata,
52 : const mfem::ParMesh &mesh)
53 : {
54 : // Check that boundary attributes have been specified correctly.
55 3 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
56 : mfem::Array<int> bdr_attr_marker;
57 3 : if (!iodata.boundaries.pec.empty())
58 : {
59 : bdr_attr_marker.SetSize(bdr_attr_max);
60 : bdr_attr_marker = 0;
61 21 : for (auto attr : mesh.bdr_attributes)
62 : {
63 18 : bdr_attr_marker[attr - 1] = 1;
64 : }
65 : std::set<int> bdr_warn_list;
66 6 : for (auto attr : iodata.boundaries.pec.attributes)
67 : {
68 : // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
69 : // "PEC boundary attribute tags must be non-negative and correspond to "
70 : // "attributes in the mesh!");
71 : // MFEM_VERIFY(bdr_attr_marker[attr - 1],
72 : // "Unknown PEC boundary attribute " << attr << "!");
73 3 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
74 : {
75 : bdr_warn_list.insert(attr);
76 : }
77 : }
78 3 : if (!bdr_warn_list.empty())
79 : {
80 0 : Mpi::Print("\n");
81 0 : Mpi::Warning("Unknown PEC boundary attributes!\nSolver will just ignore them!");
82 0 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
83 0 : Mpi::Print("\n");
84 : }
85 : }
86 :
87 : // Mark selected boundary attributes from the mesh as essential (Dirichlet).
88 : mfem::Array<int> dbc_bcs;
89 3 : dbc_bcs.Reserve(static_cast<int>(iodata.boundaries.pec.attributes.size()));
90 6 : for (auto attr : iodata.boundaries.pec.attributes)
91 : {
92 3 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
93 : {
94 0 : continue; // Can just ignore if wrong
95 : }
96 3 : dbc_bcs.Append(attr);
97 : }
98 3 : return dbc_bcs;
99 : }
100 :
101 3 : void CurlCurlOperator::CheckBoundaryProperties()
102 : {
103 : // A final check that no boundary attribute is assigned multiple boundary conditions.
104 : const mfem::ParMesh &mesh = GetMesh();
105 3 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
106 3 : const auto dbc_marker = mesh::AttrToMarker(bdr_attr_max, dbc_attr);
107 3 : const auto surf_j_marker = mesh::AttrToMarker(bdr_attr_max, surf_j_op.GetAttrList());
108 21 : for (int i = 0; i < dbc_marker.Size(); i++)
109 : {
110 18 : MFEM_VERIFY(dbc_marker[i] + surf_j_marker[i] <= 1,
111 : "Boundary attributes should not be specified with multiple BC!");
112 : }
113 3 : }
114 :
115 : namespace
116 : {
117 :
118 0 : void PrintHeader(const mfem::ParFiniteElementSpace &h1_fespace,
119 : const mfem::ParFiniteElementSpace &nd_fespace,
120 : const mfem::ParFiniteElementSpace &rt_fespace, bool &print_hdr)
121 : {
122 0 : if (print_hdr)
123 : {
124 0 : Mpi::Print("\nAssembling system matrices, number of global unknowns:\n"
125 : " H1 (p = {:d}): {:d}, ND (p = {:d}): {:d}, RT (p = {:d}): {:d}\n Operator "
126 : "assembly level: {}\n",
127 0 : h1_fespace.GetMaxElementOrder(), h1_fespace.GlobalTrueVSize(),
128 0 : nd_fespace.GetMaxElementOrder(), nd_fespace.GlobalTrueVSize(),
129 0 : rt_fespace.GetMaxElementOrder(), rt_fespace.GlobalTrueVSize(),
130 0 : (nd_fespace.GetMaxElementOrder() >= BilinearForm::pa_order_threshold)
131 0 : ? "Partial"
132 : : "Full");
133 :
134 : const auto &mesh = *nd_fespace.GetParMesh();
135 0 : const auto geom_types = mesh::CheckElements(mesh).GetGeomTypes();
136 0 : Mpi::Print(" Mesh geometries:\n");
137 0 : for (auto geom : geom_types)
138 : {
139 0 : const auto *fe = nd_fespace.FEColl()->FiniteElementForGeometry(geom);
140 0 : MFEM_VERIFY(fe, "MFEM does not support ND spaces on geometry = "
141 : << mfem::Geometry::Name[geom] << "!");
142 0 : const int q_order = fem::DefaultIntegrationOrder::Get(mesh, geom);
143 0 : Mpi::Print(" {}: P = {:d}, Q = {:d} (quadrature order = {:d}){}\n",
144 0 : mfem::Geometry::Name[geom], fe->GetDof(),
145 0 : mfem::IntRules.Get(geom, q_order).GetNPoints(), q_order,
146 0 : (geom == geom_types.back()) ? "" : ",");
147 : }
148 :
149 0 : Mpi::Print("\nAssembling multigrid hierarchy:\n");
150 : }
151 0 : }
152 :
153 : } // namespace
154 :
155 0 : std::unique_ptr<Operator> CurlCurlOperator::GetStiffnessMatrix()
156 : {
157 : // When partially assembled, the coarse operators can reuse the fine operator quadrature
158 : // data if the spaces correspond to the same mesh.
159 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
160 :
161 : constexpr bool skip_zeros = false;
162 : MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
163 0 : mat_op.GetInvPermeability());
164 : BilinearForm k(GetNDSpace());
165 0 : k.AddDomainIntegrator<CurlCurlIntegrator>(muinv_func);
166 : // k.AssembleQuadratureData();
167 0 : auto k_vec = k.Assemble(GetNDSpaces(), skip_zeros);
168 0 : auto K = std::make_unique<MultigridOperator>(GetNDSpaces().GetNumLevels());
169 0 : for (std::size_t l = 0; l < GetNDSpaces().GetNumLevels(); l++)
170 : {
171 : const auto &nd_fespace_l = GetNDSpaces().GetFESpaceAtLevel(l);
172 0 : if (print_hdr)
173 : {
174 0 : Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l,
175 0 : nd_fespace_l.GetMaxElementOrder(), nd_fespace_l.GlobalTrueVSize());
176 0 : if (const auto *k_spm = dynamic_cast<const hypre::HypreCSRMatrix *>(k_vec[l].get()))
177 : {
178 0 : HYPRE_BigInt nnz = k_spm->NNZ();
179 : Mpi::GlobalSum(1, &nnz, nd_fespace_l.GetComm());
180 0 : Mpi::Print(", {:d} NNZ\n", nnz);
181 : }
182 : else
183 : {
184 0 : Mpi::Print("\n");
185 : }
186 : }
187 0 : auto K_l = std::make_unique<ParOperator>(std::move(k_vec[l]), nd_fespace_l);
188 0 : K_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE);
189 0 : K->AddOperator(std::move(K_l));
190 : }
191 :
192 0 : print_hdr = false;
193 0 : return K;
194 0 : }
195 :
196 0 : void CurlCurlOperator::GetExcitationVector(int idx, Vector &RHS)
197 : {
198 : // Assemble the surface current excitation +J. The SurfaceCurrentOperator assembles -J
199 : // (meant for time or frequency domain Maxwell discretization, so we multiply by -1 to
200 : // retrieve +J).
201 : SumVectorCoefficient fb(GetMesh().SpaceDimension());
202 0 : surf_j_op.AddExcitationBdrCoefficients(idx, fb);
203 0 : RHS.SetSize(GetNDSpace().GetTrueVSize());
204 0 : RHS.UseDevice(true);
205 0 : RHS = 0.0;
206 0 : int empty = (fb.empty());
207 : Mpi::GlobalMin(1, &empty, GetComm());
208 0 : if (empty)
209 : {
210 : return;
211 : }
212 0 : mfem::LinearForm rhs(&GetNDSpace().Get());
213 0 : rhs.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb));
214 0 : rhs.UseFastAssembly(false);
215 : rhs.UseDevice(false);
216 0 : rhs.Assemble();
217 : rhs.UseDevice(true);
218 0 : GetNDSpace().GetProlongationMatrix()->AddMultTranspose(rhs, RHS, -1.0);
219 0 : linalg::SetSubVector(RHS, dbc_tdof_lists.back(), 0.0);
220 0 : }
221 :
222 : } // namespace palace
|
