Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/NodeDiagnostics.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2025 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief NodeDiagnostics class for collecting nodal diagnostics
10 : : \details NodeDiagnostics class for collecting nodal diagnostics, e.g.,
11 : : residuals, and various norms of errors while solving partial differential
12 : : equations.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include "Diagnostics.hpp"
17 : : #include "NodeDiagnostics.hpp"
18 : : #include "DiagReducer.hpp"
19 : : #include "Discretization.hpp"
20 : : #include "Problems.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : static CkReduction::reducerType DiagMerger;
25 : :
26 : : } // inciter::
27 : :
28 : : using inciter::NodeDiagnostics;
29 : :
30 : : void
31 : 4812 : NodeDiagnostics::registerReducers()
32 : : // *****************************************************************************
33 : : // Configure Charm++ reduction types
34 : : //! \details This routine is supposed to be called from a Charm++ initnode
35 : : //! routine. Since the runtime system executes initnode routines exactly once
36 : : //! on every logical node early on in the Charm++ init sequence, they must be
37 : : //! static as they are called without an object. See also: Section
38 : : //! "Initializations at Program Startup" at in the Charm++ manual
39 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
40 : : // *****************************************************************************
41 : : {
42 : 4812 : DiagMerger = CkReduction::addReducer( diagnostics::mergeDiag );
43 : 4812 : }
44 : :
45 : : bool
46 : 31176 : NodeDiagnostics::rhocompute( Discretization& d,
47 : : const tk::Fields& u,
48 : : const tk::Fields& un,
49 : : uint64_t diag_iter ) const
50 : : // *****************************************************************************
51 : : // Compute diagnostics for density-based solvers
52 : : //! \param[in] d Discretization proxy to read from
53 : : //! \param[in] u Current solution vector
54 : : //! \param[in] un Previous solution vector
55 : : //! \param[in] diag_iter Diagnostics output frequency
56 : : //! \return True if diagnostics have been computed
57 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
58 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
59 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
60 : : //! We send multiple sets of quantities to the host for aggregation across
61 : : //! the whole mesh. The final aggregated solution will end up in
62 : : //! Transporter::diagnostics(). Aggregation of the partially computed
63 : : //! diagnostics is done via potentially different policies for each field.
64 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
65 : : // *****************************************************************************
66 : : {
67 : : using namespace diagnostics;
68 : :
69 : : // Only compute diagnostics if needed in this time step
70 [ + + ]: 31176 : if ( (d.It()+1) % diag_iter ) return false;
71 : :
72 : 30112 : auto ncomp = u.nprop();
73 : :
74 : : // Diagnostics vector (of vectors) during aggregation. See
75 : : // Inciter/Diagnostics.h.
76 : : std::vector< std::vector< tk::real > >
77 [ + - ][ + - ]: 30112 : diag( NUMDIAG, std::vector< tk::real >( ncomp, 0.0 ) );
78 : :
79 : 30112 : const auto& v = d.V(); // nodal volumes without contributions from others
80 : :
81 : : // query function to evaluate analytic solution (if defined)
82 [ + - ]: 30112 : auto sol = problems::SOL();
83 : :
84 : : // Evaluate analytic solution (if defined)
85 [ + - ]: 30112 : auto an = u;
86 [ + + ]: 30112 : if (sol) {
87 : 13772 : const auto& coord = d.Coord();
88 : 13772 : const auto& x = coord[0];
89 : 13772 : const auto& y = coord[1];
90 : 13772 : const auto& z = coord[2];
91 [ + + ]: 888450 : for (std::size_t i=0; i<u.nunk(); ++i) {
92 [ + - ]: 874678 : auto s = sol( x[i], y[i], z[i], d.T()+d.Dt() );
93 : 874678 : s[1] /= s[0];
94 : 874678 : s[2] /= s[0];
95 : 874678 : s[3] /= s[0];
96 : 874678 : s[4] = s[4] / s[0] - 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
97 [ + - ][ + + ]: 5648818 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
98 : 874678 : }
99 : : }
100 : :
101 : : // Put in norms sweeping our mesh chunk
102 [ + + ]: 4864400 : for (std::size_t i=0; i<u.nunk(); ++i) {
103 : : // Compute sum for L2 norm of the numerical solution
104 [ + + ]: 29415768 : for (std::size_t c=0; c<ncomp; ++c)
105 [ + - ][ + - ]: 24581480 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
106 : : // Compute sum for L2 norm of the residual
107 [ + + ]: 29415768 : for (std::size_t c=0; c<ncomp; ++c)
108 [ + - ][ + - ]: 24581480 : diag[L2RES][c] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
109 : : // Compute sum for the total energy over the entire domain (first entry)
110 [ + - ]: 4834288 : diag[TOTALEN][0] += u(i,4) * v[i];
111 : : // Compute sum for L2 norm of the numerical-analytic solution
112 [ + + ]: 4834288 : if (sol) {
113 [ + - ]: 874678 : auto nu = u[i];
114 : 874678 : nu[1] /= nu[0];
115 : 874678 : nu[2] /= nu[0];
116 : 874678 : nu[3] /= nu[0];
117 : 874678 : nu[4] = nu[4] / nu[0] - 0.5*(nu[1]*nu[1] + nu[2]*nu[2] + nu[3]*nu[3]);
118 [ + + ]: 5248068 : for (std::size_t c=0; c<5; ++c) {
119 [ + - ]: 4373390 : auto du = nu[c] - an(i,c);
120 : 4373390 : diag[L2ERR][c] += du * du * v[i];
121 : 4373390 : diag[L1ERR][c] += std::abs( du ) * v[i];
122 : : }
123 [ + + ]: 1275428 : for (std::size_t c=5; c<ncomp; ++c) {
124 [ + - ][ + - ]: 400750 : auto du = u(i,c) - an(i,c);
125 : 400750 : diag[L2ERR][c] += du * du * v[i];
126 : 400750 : diag[L1ERR][c] += std::abs( du ) * v[i];
127 : : }
128 : 874678 : }
129 : : }
130 : :
131 : : // Append diagnostics vector with metadata on the current time step
132 : : // ITER:: Current iteration count (only the first entry is used)
133 : : // TIME: Current physical time (only the first entry is used)
134 : : // DT: Current physical time step size (only the first entry is used)
135 : 30112 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
136 : 30112 : diag[TIME][0] = d.T() + d.Dt();
137 : 30112 : diag[DT][0] = d.Dt();
138 : :
139 : : // Contribute to diagnostics
140 [ + - ]: 30112 : auto stream = serialize( d.MeshId(), ncomp, diag );
141 [ + - ]: 30112 : d.contribute( stream.first, stream.second.get(), DiagMerger,
142 [ + - ][ + - ]: 60224 : CkCallback(CkIndex_Transporter::rhodiagnostics(nullptr), d.Tr()) );
143 : :
144 : 30112 : return true; // diagnostics have been computed
145 : 30112 : }
146 : :
147 : : bool
148 : 4662 : NodeDiagnostics::precompute( Discretization& d,
149 : : const tk::Fields& u,
150 : : const tk::Fields& un,
151 : : const std::vector< tk::real >& p,
152 : : const std::vector< tk::real >& dp,
153 : : uint64_t diag_iter ) const
154 : : // *****************************************************************************
155 : : // Compute diagnostics for pressure-based solvers
156 : : //! \param[in] d Discretization proxy to read from
157 : : //! \param[in] u Current solution vector
158 : : //! \param[in] un Previous solution vector
159 : : //! \param[in] p Current pressure solution
160 : : //! \param[in] dp Recent pressure solution increment
161 : : //! \param[in] diag_iter Diagnostics output frequency
162 : : //! \return True if diagnostics have been computed
163 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
164 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
165 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
166 : : //! We send multiple sets of quantities to the host for aggregation across
167 : : //! the whole mesh. The final aggregated solution will end up in
168 : : //! Transporter::diagnostics(). Aggregation of the partially computed
169 : : //! diagnostics is done via potentially different policies for each field.
170 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
171 : : // *****************************************************************************
172 : : {
173 : : using namespace diagnostics;
174 : :
175 : : // Only compute diagnostics if needed in this time step
176 [ - + ]: 4662 : if ( (d.It()+1) % diag_iter ) return false;
177 : :
178 [ - + ][ - - ]: 4662 : Assert( p.size() == u.nunk(), "Size mismatch" );
[ - - ][ - - ]
179 [ - + ][ - - ]: 4662 : Assert( p.size() == dp.size(), "Size mismatch" );
[ - - ][ - - ]
180 [ - + ][ - - ]: 4662 : Assert( u.nunk() == un.nunk(), "Size mismatch" );
[ - - ][ - - ]
181 : :
182 : 4662 : auto ncomp = u.nprop();
183 : :
184 : 4662 : const auto& v = d.V(); // nodal volumes without contributions from others
185 : :
186 : : // query function to evaluate analytic solution (if defined)
187 [ + - ]: 4662 : auto pressure_sol = problems::PRESSURE_SOL();
188 : :
189 : : // Evaluate analytic solution (if defined)
190 [ + - ]: 4662 : auto ap = p;
191 [ + + ]: 4662 : if (pressure_sol) {
192 : 32 : ncomp = 0;
193 : 32 : const auto& coord = d.Coord();
194 : 32 : const auto& x = coord[0];
195 : 32 : const auto& y = coord[1];
196 : 32 : const auto& z = coord[2];
197 [ + + ]: 2174 : for (std::size_t i=0; i<p.size(); ++i) {
198 [ + - ]: 2142 : ap[i] = pressure_sol( x[i], y[i], z[i] );
199 : : }
200 : : }
201 : :
202 : : // query function to evaluate analytic solution (if defined)
203 [ + - ]: 4662 : auto sol = problems::SOL();
204 : :
205 : : // Evaluate analytic solution (if defined)
206 [ + - ]: 4662 : auto an = u;
207 [ + + ]: 4662 : if (sol) {
208 : 1032 : const auto& coord = d.Coord();
209 : 1032 : const auto& x = coord[0];
210 : 1032 : const auto& y = coord[1];
211 : 1032 : const auto& z = coord[2];
212 [ + + ]: 118214 : for (std::size_t i=0; i<u.nunk(); ++i) {
213 [ + - ]: 117182 : auto s = sol( x[i], y[i], z[i], d.T()+d.Dt() );
214 [ + - ][ + + ]: 583768 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
215 : 117182 : }
216 : : }
217 : :
218 : : // Diagnostics vector (of vectors) during aggregation. See
219 : : // Inciter/Diagnostics.h.
220 : : std::vector< std::vector< tk::real > >
221 [ + - ][ + - ]: 4662 : diag( NUMDIAG, std::vector< tk::real >( ncomp+1, 0.0 ) );
222 : :
223 : : // Put in norms sweeping our mesh chunk
224 [ + + ]: 472754 : for (std::size_t i=0; i<u.nunk(); ++i) {
225 : : // Compute sum for L2 norm of the numerical solution
226 : 468092 : diag[L2SOL][0] += p[i] * p[i] * v[i];
227 [ + + ]: 1980982 : for (std::size_t c=0; c<ncomp; ++c) {
228 [ + - ][ + - ]: 1512890 : diag[L2SOL][c+1] += u(i,c) * u(i,c) * v[i];
229 : : }
230 : : // Compute sum for L2 norm of the residual
231 : 468092 : diag[L2RES][0] += dp[i] * dp[i] * v[i];
232 [ + + ]: 1980982 : for (std::size_t c=0; c<ncomp; ++c) {
233 [ + - ][ + - ]: 1512890 : diag[L2RES][c+1] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
234 : : }
235 : : // Compute sum for the total energy over the entire domain
236 : 468092 : diag[TOTALEN][0] += 0.0 * v[i];
237 : : // Compute sum for norms of the numerical-analytic pressure solution
238 [ + + ]: 468092 : if (pressure_sol) {
239 : 2142 : auto pd = p[i] - ap[i];
240 : 2142 : diag[L2ERR][0] += pd * pd * v[i];
241 : 2142 : diag[L1ERR][0] += std::abs( pd ) * v[i];
242 : : }
243 : : // Compute sum for norms of the numerical-analytic adv/diff solution
244 [ + + ]: 468092 : if (sol) {
245 [ + + ]: 577342 : for (std::size_t c=0; c<ncomp; ++c) {
246 [ + - ][ + - ]: 460160 : auto du = u(i,c) - an(i,c);
247 : 460160 : diag[L2ERR][c+1] += du * du * v[i];
248 : 460160 : diag[L1ERR][c+1] += std::abs( du ) * v[i];
249 : : }
250 : : }
251 : : }
252 : :
253 : : // Append diagnostics vector with metadata on the current time step
254 : : // ITER:: Current iteration count (only the first entry is used)
255 : : // TIME: Current physical time (only the first entry is used)
256 : : // DT: Current physical time step size (only the first entry is used)
257 : 4662 : diag[ITER][0] = static_cast< tk::real >( d.It() );
258 : 4662 : diag[TIME][0] = d.T();
259 : 4662 : diag[DT][0] = d.Dt();
260 : :
261 : : // Contribute to diagnostics
262 [ + - ]: 4662 : auto stream = serialize( d.MeshId(), ncomp+1, diag );
263 [ + - ]: 4662 : d.contribute( stream.first, stream.second.get(), DiagMerger,
264 [ + - ][ + - ]: 9324 : CkCallback(CkIndex_Transporter::prediagnostics(nullptr), d.Tr()) );
265 : :
266 : 4662 : return true; // diagnostics have been computed
267 : 4662 : }
268 : :
269 : : bool
270 : 4730 : NodeDiagnostics::accompute( Discretization& d,
271 : : const tk::Fields& u,
272 : : const tk::Fields& un,
273 : : uint64_t diag_iter ) const
274 : : // *****************************************************************************
275 : : // Compute diagnostics for artificial compressibility solvers
276 : : //! \param[in] d Discretization proxy to read from
277 : : //! \param[in] u Current solution vector
278 : : //! \param[in] un Previous solution vector
279 : : //! \param[in] diag_iter Diagnostics output frequency
280 : : //! \return True if diagnostics have been computed
281 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
282 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
283 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
284 : : //! We send multiple sets of quantities to the host for aggregation across
285 : : //! the whole mesh. The final aggregated solution will end up in
286 : : //! Transporter::diagnostics(). Aggregation of the partially computed
287 : : //! diagnostics is done via potentially different policies for each field.
288 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
289 : : // *****************************************************************************
290 : : {
291 : : using namespace diagnostics;
292 : :
293 : : // Only compute diagnostics if needed in this time step
294 [ + + ]: 4730 : if ( (d.It()+1) % diag_iter ) return false;
295 : :
296 : 4532 : auto ncomp = u.nprop();
297 : :
298 : : // Diagnostics vector (of vectors) during aggregation. See
299 : : // Inciter/Diagnostics.h.
300 : : std::vector< std::vector< tk::real > >
301 [ + - ][ + - ]: 4532 : diag( NUMDIAG, std::vector< tk::real >( ncomp, 0.0 ) );
302 : :
303 : 4532 : const auto& v = d.V(); // nodal volumes without contributions from others
304 : :
305 : : // query function to evaluate analytic solution (if defined)
306 [ + - ]: 4532 : auto sol = problems::SOL();
307 : :
308 : : // Evaluate analytic solution (if defined)
309 [ + - ]: 4532 : auto an = u;
310 [ + + ]: 4532 : if (sol) {
311 : 1000 : const auto& coord = d.Coord();
312 : 1000 : const auto& x = coord[0];
313 : 1000 : const auto& y = coord[1];
314 : 1000 : const auto& z = coord[2];
315 [ + + ]: 116040 : for (std::size_t i=0; i<u.nunk(); ++i) {
316 [ + - ]: 115040 : auto s = sol( x[i], y[i], z[i], d.T()+d.Dt() );
317 [ + - ][ + + ]: 690240 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
318 : 115040 : }
319 : : }
320 : :
321 : : // Put in norms sweeping our mesh chunk
322 [ + + ]: 343192 : for (std::size_t i=0; i<u.nunk(); ++i) {
323 : : // Compute sum for L2 norm of the numerical solution
324 [ + + ]: 1808340 : for (std::size_t c=0; c<ncomp; ++c)
325 [ + - ][ + - ]: 1469680 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
326 : : // Compute sum for L2 norm of the residual
327 [ + + ]: 1808340 : for (std::size_t c=0; c<ncomp; ++c)
328 [ + - ][ + - ]: 1469680 : diag[L2RES][c] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
329 : : // Compute sum for the total energy over the entire domain
330 : 338660 : diag[TOTALEN][0] += 0.0 * v[i];
331 : : // Compute sum for norms of the numerical-analytic adv/diff solution
332 [ + + ]: 338660 : if (sol) {
333 [ + + ]: 575200 : for (std::size_t c=1; c<ncomp; ++c) {
334 [ + - ][ + - ]: 460160 : auto du = u(i,c) - an(i,c);
335 : 460160 : diag[L2ERR][c] += du * du * v[i];
336 : 460160 : diag[L1ERR][c] += std::abs( du ) * v[i];
337 : : }
338 : : }
339 : : }
340 : :
341 : : // Append diagnostics vector with metadata on the current time step
342 : : // ITER:: Current iteration count (only the first entry is used)
343 : : // TIME: Current physical time (only the first entry is used)
344 : : // DT: Current physical time step size (only the first entry is used)
345 : 4532 : diag[ITER][0] = static_cast< tk::real >( d.It() );
346 : 4532 : diag[TIME][0] = d.T();
347 : 4532 : diag[DT][0] = d.Dt();
348 : :
349 : : // Contribute to diagnostics
350 [ + - ]: 4532 : auto stream = serialize( d.MeshId(), ncomp, diag );
351 [ + - ]: 4532 : d.contribute( stream.first, stream.second.get(), DiagMerger,
352 [ + - ][ + - ]: 9064 : CkCallback(CkIndex_Transporter::acdiagnostics(nullptr), d.Tr()) );
353 : :
354 : 4532 : return true; // diagnostics have been computed
355 : 4532 : }
|
