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 "laplaceoperator.hpp"
5 :
6 : #include <set>
7 : #include "fem/bilinearform.hpp"
8 : #include "fem/integrator.hpp"
9 : #include "fem/mesh.hpp"
10 : #include "fem/multigrid.hpp"
11 : #include "linalg/hypre.hpp"
12 : #include "linalg/rap.hpp"
13 : #include "utils/communication.hpp"
14 : #include "utils/geodata.hpp"
15 : #include "utils/iodata.hpp"
16 : #include "utils/prettyprint.hpp"
17 :
18 : namespace palace
19 : {
20 :
21 3 : LaplaceOperator::LaplaceOperator(const IoData &iodata,
22 3 : const std::vector<std::unique_ptr<Mesh>> &mesh)
23 3 : : print_hdr(true), dbc_attr(SetUpBoundaryProperties(iodata, *mesh.back())),
24 3 : h1_fecs(fem::ConstructFECollections<mfem::H1_FECollection>(
25 3 : iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels,
26 3 : iodata.solver.linear.mg_coarsening, false)),
27 3 : nd_fec(std::make_unique<mfem::ND_FECollection>(iodata.solver.order,
28 3 : mesh.back()->Dimension())),
29 0 : rt_fecs(fem::ConstructFECollections<mfem::RT_FECollection>(
30 3 : iodata.solver.order - 1, mesh.back()->Dimension(),
31 3 : iodata.solver.linear.estimator_mg ? iodata.solver.linear.mg_max_levels : 1,
32 3 : iodata.solver.linear.mg_coarsening, false)),
33 3 : h1_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::H1_FECollection>(
34 3 : iodata.solver.linear.mg_max_levels, mesh, h1_fecs, &dbc_attr, &dbc_tdof_lists)),
35 3 : nd_fespace(*mesh.back(), nd_fec.get()),
36 3 : rt_fespaces(fem::ConstructFiniteElementSpaceHierarchy<mfem::RT_FECollection>(
37 3 : iodata.solver.linear.estimator_mg ? iodata.solver.linear.mg_max_levels : 1, mesh,
38 3 : rt_fecs)),
39 3 : mat_op(iodata, *mesh.back()), source_attr_lists(ConstructSources(iodata))
40 : {
41 : // Print essential BC information.
42 3 : if (dbc_attr.Size())
43 : {
44 3 : Mpi::Print("\nConfiguring Dirichlet BC at attributes:\n");
45 6 : utils::PrettyPrint(dbc_attr);
46 : }
47 3 : }
48 :
49 3 : mfem::Array<int> LaplaceOperator::SetUpBoundaryProperties(const IoData &iodata,
50 : const mfem::ParMesh &mesh)
51 : {
52 : // Check that boundary attributes have been specified correctly.
53 3 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
54 : mfem::Array<int> bdr_attr_marker;
55 3 : if (!iodata.boundaries.pec.empty() || !iodata.boundaries.lumpedport.empty())
56 : {
57 : bdr_attr_marker.SetSize(bdr_attr_max);
58 : bdr_attr_marker = 0;
59 21 : for (auto attr : mesh.bdr_attributes)
60 : {
61 18 : bdr_attr_marker[attr - 1] = 1;
62 : }
63 : std::set<int> bdr_warn_list;
64 6 : for (auto attr : iodata.boundaries.pec.attributes)
65 : {
66 : // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
67 : // "Ground boundary attribute tags must be non-negative and correspond to
68 : // " attributes in the mesh!");
69 : // MFEM_VERIFY(bdr_attr_marker[attr - 1],
70 : // "Unknown ground boundary attribute " << attr << "!");
71 3 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
72 : {
73 : bdr_warn_list.insert(attr);
74 : }
75 : }
76 3 : if (!bdr_warn_list.empty())
77 : {
78 0 : Mpi::Print("\n");
79 0 : Mpi::Warning("Unknown ground boundary attributes!\nSolver will just ignore them!");
80 0 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
81 0 : Mpi::Print("\n");
82 : }
83 3 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
84 : {
85 0 : for (const auto &elem : data.elements)
86 : {
87 0 : for (auto attr : elem.attributes)
88 : {
89 0 : MFEM_VERIFY(
90 : attr > 0 && attr <= bdr_attr_max,
91 : "Terminal boundary attribute tags must be non-negative and correspond to "
92 : "attributes in the mesh!");
93 0 : MFEM_VERIFY(bdr_attr_marker[attr - 1] > 0,
94 : "Unknown terminal boundary attribute " << attr << "!");
95 : }
96 : }
97 : }
98 : }
99 :
100 : // Mark selected boundary attributes from the mesh as essential (Dirichlet).
101 : mfem::Array<int> dbc_bcs;
102 3 : dbc_bcs.Reserve(static_cast<int>(iodata.boundaries.pec.attributes.size()) +
103 : static_cast<int>(iodata.boundaries.lumpedport.size()));
104 6 : for (auto attr : iodata.boundaries.pec.attributes)
105 : {
106 3 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
107 : {
108 0 : continue; // Can just ignore if wrong
109 : }
110 3 : dbc_bcs.Append(attr);
111 : }
112 3 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
113 : {
114 0 : for (const auto &elem : data.elements)
115 : {
116 0 : for (auto attr : elem.attributes)
117 : {
118 0 : dbc_bcs.Append(attr);
119 : }
120 : }
121 : }
122 3 : MFEM_VERIFY(dbc_bcs.Size() > 0,
123 : "Electrostatic problem is ill-posed without any Dirichlet boundaries!");
124 3 : return dbc_bcs;
125 : }
126 :
127 3 : std::map<int, mfem::Array<int>> LaplaceOperator::ConstructSources(const IoData &iodata)
128 : {
129 : // Construct mapping from terminal index to list of associated attributes.
130 : std::map<int, mfem::Array<int>> attr_lists;
131 3 : for (const auto &[idx, data] : iodata.boundaries.lumpedport)
132 : {
133 0 : mfem::Array<int> &attr_list = attr_lists[idx];
134 0 : attr_list.Reserve(
135 : static_cast<int>(data.elements.size())); // Average one attribute per element
136 0 : for (const auto &elem : data.elements)
137 : {
138 0 : for (auto attr : elem.attributes)
139 : {
140 0 : attr_list.Append(attr);
141 : }
142 : }
143 : }
144 3 : return attr_lists;
145 : }
146 :
147 : namespace
148 : {
149 :
150 0 : void PrintHeader(const mfem::ParFiniteElementSpace &h1_fespace,
151 : const mfem::ParFiniteElementSpace &nd_fespace,
152 : const mfem::ParFiniteElementSpace &rt_fespace, bool &print_hdr)
153 : {
154 0 : if (print_hdr)
155 : {
156 0 : Mpi::Print("\nAssembling system matrices, number of global unknowns:\n"
157 : " H1 (p = {:d}): {:d}, ND (p = {:d}): {:d}, RT (p = {:d}): {:d}\n Operator "
158 : "assembly level: {}\n",
159 0 : h1_fespace.GetMaxElementOrder(), h1_fespace.GlobalTrueVSize(),
160 0 : nd_fespace.GetMaxElementOrder(), nd_fespace.GlobalTrueVSize(),
161 0 : rt_fespace.GetMaxElementOrder(), rt_fespace.GlobalTrueVSize(),
162 0 : (h1_fespace.GetMaxElementOrder() >= BilinearForm::pa_order_threshold)
163 0 : ? "Partial"
164 : : "Full");
165 :
166 : const auto &mesh = *h1_fespace.GetParMesh();
167 0 : const auto geom_types = mesh::CheckElements(mesh).GetGeomTypes();
168 0 : Mpi::Print(" Mesh geometries:\n");
169 0 : for (auto geom : geom_types)
170 : {
171 0 : const auto *fe = nd_fespace.FEColl()->FiniteElementForGeometry(geom);
172 0 : MFEM_VERIFY(fe, "MFEM does not support ND spaces on geometry = "
173 : << mfem::Geometry::Name[geom] << "!");
174 0 : const int q_order = fem::DefaultIntegrationOrder::Get(mesh, geom);
175 0 : Mpi::Print(" {}: P = {:d}, Q = {:d} (quadrature order = {:d}){}\n",
176 0 : mfem::Geometry::Name[geom], fe->GetDof(),
177 0 : mfem::IntRules.Get(geom, q_order).GetNPoints(), q_order,
178 0 : (geom == geom_types.back()) ? "" : ",");
179 : }
180 :
181 0 : Mpi::Print("\nAssembling multigrid hierarchy:\n");
182 : }
183 0 : }
184 :
185 : } // namespace
186 :
187 0 : std::unique_ptr<Operator> LaplaceOperator::GetStiffnessMatrix()
188 : {
189 : // When partially assembled, the coarse operators can reuse the fine operator quadrature
190 : // data if the spaces correspond to the same mesh.
191 0 : PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), print_hdr);
192 :
193 : constexpr bool skip_zeros = false;
194 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
195 0 : mat_op.GetPermittivityReal());
196 : BilinearForm k(GetH1Space());
197 0 : k.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
198 : // k.AssembleQuadratureData();
199 0 : auto k_vec = k.Assemble(GetH1Spaces(), skip_zeros);
200 0 : auto K = std::make_unique<MultigridOperator>(GetH1Spaces().GetNumLevels());
201 0 : for (std::size_t l = 0; l < GetH1Spaces().GetNumLevels(); l++)
202 : {
203 : const auto &h1_fespace_l = GetH1Spaces().GetFESpaceAtLevel(l);
204 0 : if (print_hdr)
205 : {
206 0 : Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l,
207 0 : h1_fespace_l.GetMaxElementOrder(), h1_fespace_l.GlobalTrueVSize());
208 0 : if (const auto *k_spm = dynamic_cast<const hypre::HypreCSRMatrix *>(k_vec[l].get()))
209 : {
210 0 : HYPRE_BigInt nnz = k_spm->NNZ();
211 : Mpi::GlobalSum(1, &nnz, h1_fespace_l.GetComm());
212 0 : Mpi::Print(", {:d} NNZ\n", nnz);
213 : }
214 : else
215 : {
216 0 : Mpi::Print("\n");
217 : }
218 : }
219 0 : auto K_l = std::make_unique<ParOperator>(std::move(k_vec[l]), h1_fespace_l);
220 0 : K_l->SetEssentialTrueDofs(dbc_tdof_lists[l], Operator::DiagonalPolicy::DIAG_ONE);
221 0 : K->AddOperator(std::move(K_l));
222 : }
223 :
224 0 : print_hdr = false;
225 0 : return K;
226 0 : }
227 :
228 0 : void LaplaceOperator::GetExcitationVector(int idx, const Operator &K, Vector &X,
229 : Vector &RHS)
230 : {
231 : // Apply the Dirichlet BCs to the solution vector: V = 1 on terminal boundaries with the
232 : // given index, V = 0 on all ground and other terminal boundaries.
233 0 : mfem::ParGridFunction x(&GetH1Space().Get());
234 : x = 0.0;
235 :
236 : // Get a marker of all boundary attributes with the given source surface index.
237 : const mfem::ParMesh &mesh = GetMesh();
238 0 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
239 0 : mfem::Array<int> source_marker = mesh::AttrToMarker(bdr_attr_max, source_attr_lists[idx]);
240 : mfem::ConstantCoefficient one(1.0);
241 : x.ProjectBdrCoefficient(one, source_marker); // Values are only correct on master
242 :
243 : // Eliminate the essential BC to get the RHS vector.
244 0 : X.SetSize(GetH1Space().GetTrueVSize());
245 0 : RHS.SetSize(GetH1Space().GetTrueVSize());
246 0 : X.UseDevice(true);
247 0 : RHS.UseDevice(true);
248 0 : X = 0.0;
249 0 : RHS = 0.0;
250 0 : x.ParallelProject(X); // Restrict to the true dofs
251 0 : const auto *mg_K = dynamic_cast<const MultigridOperator *>(&K);
252 0 : const auto *PtAP_K = mg_K ? dynamic_cast<const ParOperator *>(&mg_K->GetFinestOperator())
253 0 : : dynamic_cast<const ParOperator *>(&K);
254 0 : MFEM_VERIFY(PtAP_K, "LaplaceOperator requires ParOperator for RHS elimination!");
255 0 : PtAP_K->EliminateRHS(X, RHS);
256 0 : }
257 :
258 : } // namespace palace
|
