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 "basesolver.hpp"
5 :
6 : #include <array>
7 : #include <complex>
8 : #include <numeric>
9 : #include <mfem.hpp>
10 : #include <nlohmann/json.hpp>
11 : #include "drivers/transientsolver.hpp"
12 : #include "fem/errorindicator.hpp"
13 : #include "fem/fespace.hpp"
14 : #include "fem/mesh.hpp"
15 : #include "linalg/ksp.hpp"
16 : #include "models/domainpostoperator.hpp"
17 : #include "models/portexcitations.hpp"
18 : #include "models/postoperator.hpp"
19 : #include "models/surfacepostoperator.hpp"
20 : #include "utils/communication.hpp"
21 : #include "utils/dorfler.hpp"
22 : #include "utils/filesystem.hpp"
23 : #include "utils/geodata.hpp"
24 : #include "utils/iodata.hpp"
25 : #include "utils/timer.hpp"
26 :
27 : namespace palace
28 : {
29 :
30 : using json = nlohmann::json;
31 :
32 : namespace
33 : {
34 :
35 0 : void SaveIteration(MPI_Comm comm, const fs::path &output_dir, int step, int width)
36 : {
37 0 : BlockTimer bt(Timer::IO);
38 : Mpi::Barrier(comm); // Wait for all processes to write postprocessing files
39 0 : if (Mpi::Root(comm))
40 : {
41 : // Create a subfolder for the results of this adaptation.
42 0 : auto step_output = output_dir / fmt::format("iteration{:0{}d}", step, width);
43 : if (!fs::exists(step_output))
44 : {
45 0 : fs::create_directories(step_output);
46 : }
47 : constexpr auto options =
48 : fs::copy_options::recursive | fs::copy_options::overwrite_existing;
49 0 : for (const auto &f : fs::directory_iterator(output_dir))
50 : {
51 0 : if (f.path().filename().string().rfind("iteration") == 0)
52 : {
53 0 : continue;
54 : }
55 0 : fs::copy(f, step_output / f.path().filename(), options);
56 : }
57 0 : }
58 : Mpi::Barrier(comm);
59 0 : }
60 :
61 0 : json LoadMetadata(const fs::path &post_dir)
62 : {
63 0 : std::string path = post_dir / "palace.json";
64 0 : std::ifstream fi(path);
65 0 : if (!fi.is_open())
66 : {
67 0 : MFEM_ABORT("Unable to open metadata file \"" << path << "\"!");
68 : }
69 0 : return json::parse(fi);
70 0 : }
71 :
72 0 : void WriteMetadata(const fs::path &post_dir, const json &meta)
73 : {
74 0 : std::string path = post_dir / "palace.json";
75 0 : std::ofstream fo(path);
76 0 : if (!fo.is_open())
77 : {
78 0 : MFEM_ABORT("Unable to open metadata file \"" << path << "\"!");
79 : }
80 0 : fo << meta.dump(2) << '\n';
81 0 : }
82 :
83 : // Returns an array of indices corresponding to marked elements.
84 0 : mfem::Array<int> MarkedElements(const Vector &e, double threshold)
85 : {
86 : mfem::Array<int> ind;
87 : ind.Reserve(e.Size());
88 0 : for (int i = 0; i < e.Size(); i++)
89 : {
90 0 : if (e[i] >= threshold)
91 : {
92 0 : ind.Append(i);
93 : }
94 : }
95 0 : return ind;
96 : }
97 :
98 : } // namespace
99 :
100 0 : BaseSolver::BaseSolver(const IoData &iodata, bool root, int size, int num_thread,
101 0 : const char *git_tag)
102 0 : : iodata(iodata), post_dir(iodata.problem.output), root(root)
103 : {
104 : // Initialize simulation metadata for this simulation.
105 0 : if (root)
106 : {
107 : json meta;
108 0 : if (git_tag)
109 : {
110 0 : meta["GitTag"] = std::string(git_tag);
111 : }
112 0 : if (size > 0)
113 : {
114 0 : meta["Problem"]["MPISize"] = size;
115 : }
116 0 : if (num_thread > 0)
117 : {
118 0 : meta["Problem"]["OpenMPThreads"] = num_thread;
119 : }
120 0 : WriteMetadata(post_dir, meta);
121 : }
122 0 : }
123 :
124 0 : void BaseSolver::SolveEstimateMarkRefine(std::vector<std::unique_ptr<Mesh>> &mesh) const
125 : {
126 0 : const auto &refinement = iodata.model.refinement;
127 0 : const bool use_amr = [&]()
128 : {
129 0 : if (refinement.max_it > 0 && dynamic_cast<const TransientSolver *>(this) != nullptr)
130 : {
131 0 : Mpi::Warning("AMR is not currently supported for transient simulations!\n");
132 0 : return false;
133 : }
134 0 : return (refinement.max_it > 0);
135 0 : }();
136 0 : if (use_amr && mesh.size() > 1)
137 : {
138 0 : Mpi::Print("\nFlattening mesh sequence:\n AMR will start from the final mesh in "
139 : "the sequence of a priori refinements\n");
140 : mesh.erase(mesh.begin(), mesh.end() - 1);
141 : }
142 0 : MPI_Comm comm = mesh.back()->GetComm();
143 :
144 : // Perform initial solve and estimation.
145 0 : auto [indicators, ntdof] = Solve(mesh);
146 0 : double err = indicators.Norml2(comm);
147 :
148 : // Collection of all tests that might exhaust resources.
149 : auto ExhaustedResources = [&refinement](auto it, auto ntdof)
150 : {
151 : bool ret = false;
152 : // Run out of iterations.
153 0 : ret |= (it >= refinement.max_it);
154 : // Run out of DOFs if a limit was set.
155 0 : ret |= (refinement.max_size > 0 && ntdof > refinement.max_size);
156 : return ret;
157 : };
158 :
159 : // Main AMR loop.
160 0 : int it = 0;
161 0 : while (!ExhaustedResources(it, ntdof) && err >= refinement.tol)
162 : {
163 : // Print timing summary.
164 0 : Mpi::Print("\nCumulative timing statistics:\n");
165 0 : BlockTimer::Print(comm);
166 0 : SaveMetadata(BlockTimer::GlobalTimer());
167 :
168 0 : BlockTimer bt(Timer::ADAPTATION);
169 0 : Mpi::Print("\nAdaptive mesh refinement (AMR) iteration {:d}:\n"
170 : " Indicator norm = {:.3e}, global unknowns = {:d}\n"
171 : " Max. iterations = {:d}, tol. = {:.3e}{}\n",
172 0 : ++it, err, ntdof, refinement.max_it, refinement.tol,
173 0 : (refinement.max_size > 0
174 0 : ? ", max. size = " + std::to_string(refinement.max_size)
175 0 : : ""));
176 :
177 : // Optionally save off the previous solution.
178 0 : if (refinement.save_adapt_iterations)
179 : {
180 0 : SaveIteration(comm, post_dir, it,
181 0 : 1 + static_cast<int>(std::log10(refinement.max_it)));
182 : }
183 :
184 : // Mark.
185 0 : const auto marked_elements = [&comm, &refinement](const auto &indicators)
186 : {
187 0 : const auto [threshold, marked_error] = utils::ComputeDorflerThreshold(
188 0 : comm, indicators.Local(), refinement.update_fraction);
189 0 : const auto marked_elements = MarkedElements(indicators.Local(), threshold);
190 : const auto [glob_marked_elements, glob_elements] =
191 0 : linalg::GlobalSize2(comm, marked_elements, indicators.Local());
192 0 : Mpi::Print(
193 : " Marked {:d}/{:d} elements for refinement ({:.2f}% of the error, θ = {:.2f})\n",
194 0 : glob_marked_elements, glob_elements, 100 * marked_error,
195 0 : refinement.update_fraction);
196 0 : return marked_elements;
197 0 : }(indicators);
198 :
199 : // Refine.
200 : {
201 : mfem::ParMesh &fine_mesh = *mesh.back();
202 0 : const auto initial_elem_count = fine_mesh.GetGlobalNE();
203 0 : fine_mesh.GeneralRefinement(marked_elements, -1, refinement.max_nc_levels);
204 0 : const auto final_elem_count = fine_mesh.GetGlobalNE();
205 0 : Mpi::Print(" {} mesh refinement added {:d} elements (initial = {:d}, final = {:d})\n",
206 0 : fine_mesh.Nonconforming() ? "Nonconforming" : "Conforming",
207 0 : final_elem_count - initial_elem_count, initial_elem_count,
208 : final_elem_count);
209 : }
210 :
211 : // Optionally rebalance and write the adapted mesh to file.
212 : {
213 0 : const auto ratio_pre = mesh::RebalanceMesh(iodata, *mesh.back());
214 0 : if (ratio_pre > refinement.maximum_imbalance)
215 : {
216 : int min_elem, max_elem;
217 0 : min_elem = max_elem = mesh.back()->GetNE();
218 0 : Mpi::GlobalMin(1, &min_elem, comm);
219 0 : Mpi::GlobalMax(1, &max_elem, comm);
220 0 : const auto ratio_post = double(max_elem) / min_elem;
221 0 : Mpi::Print(" Rebalanced mesh: Ratio {:.3f} exceeded max. allowed value {:.3f} "
222 : "(new ratio = {:.3f})\n",
223 0 : ratio_pre, refinement.maximum_imbalance, ratio_post);
224 : }
225 0 : mesh.back()->Update();
226 : }
227 :
228 : // Solve + estimate.
229 0 : Mpi::Print("\nProceeding with solve/estimate iteration {}...\n", it + 1);
230 0 : std::tie(indicators, ntdof) = Solve(mesh);
231 0 : err = indicators.Norml2(comm);
232 0 : }
233 0 : Mpi::Print("\nCompleted {:d} iteration{} of adaptive mesh refinement (AMR):\n"
234 : " Indicator norm = {:.3e}, global unknowns = {:d}\n"
235 : " Max. iterations = {:d}, tol. = {:.3e}{}\n",
236 0 : it, (it == 1 ? "" : "s"), err, ntdof, refinement.max_it, refinement.tol,
237 : (refinement.max_size > 0
238 0 : ? ", max. size = " + std::to_string(refinement.max_size)
239 0 : : ""));
240 0 : }
241 :
242 0 : void BaseSolver::SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const
243 : {
244 : const auto &fespace = fespaces.GetFinestFESpace();
245 0 : HYPRE_BigInt ne = fespace.GetParMesh().GetNE();
246 : Mpi::GlobalSum(1, &ne, fespace.GetComm());
247 0 : std::vector<HYPRE_BigInt> ndofs(fespaces.GetNumLevels());
248 0 : for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++)
249 : {
250 0 : ndofs[l] = fespaces.GetFESpaceAtLevel(l).GlobalTrueVSize();
251 : }
252 0 : if (root)
253 : {
254 0 : json meta = LoadMetadata(post_dir);
255 0 : meta["Problem"]["MeshElements"] = ne;
256 0 : meta["Problem"]["DegreesOfFreedom"] = ndofs.back();
257 0 : meta["Problem"]["MultigridDegreesOfFreedom"] = ndofs;
258 0 : WriteMetadata(post_dir, meta);
259 : }
260 0 : }
261 :
262 : template <typename SolverType>
263 0 : void BaseSolver::SaveMetadata(const SolverType &ksp) const
264 : {
265 0 : if (root)
266 : {
267 0 : json meta = LoadMetadata(post_dir);
268 0 : meta["LinearSolver"]["TotalSolves"] = ksp.NumTotalMult();
269 0 : meta["LinearSolver"]["TotalIts"] = ksp.NumTotalMultIterations();
270 0 : WriteMetadata(post_dir, meta);
271 : }
272 0 : }
273 :
274 0 : void BaseSolver::SaveMetadata(const Timer &timer) const
275 : {
276 0 : if (root)
277 : {
278 0 : json meta = LoadMetadata(post_dir);
279 0 : for (int i = Timer::INIT; i < Timer::NUM_TIMINGS; i++)
280 : {
281 0 : auto key = Timer::descriptions[i];
282 0 : key.erase(std::remove_if(key.begin(), key.end(), isspace), key.end());
283 0 : meta["ElapsedTime"]["Durations"][key] = timer.Data((Timer::Index)i);
284 0 : meta["ElapsedTime"]["Counts"][key] = timer.Counts((Timer::Index)i);
285 : }
286 0 : WriteMetadata(post_dir, meta);
287 : }
288 0 : }
289 :
290 0 : void BaseSolver::SaveMetadata(const PortExcitations &excitation_helper) const
291 : {
292 0 : if (root)
293 : {
294 0 : nlohmann::json meta = LoadMetadata(post_dir);
295 0 : meta["Excitations"] = excitation_helper;
296 0 : WriteMetadata(post_dir, meta);
297 : }
298 0 : }
299 :
300 : template void BaseSolver::SaveMetadata<KspSolver>(const KspSolver &) const;
301 : template void BaseSolver::SaveMetadata<ComplexKspSolver>(const ComplexKspSolver &) const;
302 :
303 : } // namespace palace
|
