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-2024 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 : 4608 : 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 : 4608 : DiagMerger = CkReduction::addReducer( diagnostics::mergeDiag );
43 : 4608 : }
44 : :
45 : : bool
46 : 31126 : 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 [ + + ]: 31126 : if ( (d.It()+1) % diag_iter ) return false;
71 : :
72 : 30062 : 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 [ + - ][ + - ]: 30062 : diag( NUMDIAG, std::vector< tk::real >( ncomp, 0.0 ) );
78 : :
79 : 30062 : const auto& v = d.V(); // nodal volumes without contributions from others
80 : :
81 : : // query function to evaluate analytic solution (if defined)
82 [ + - ]: 30062 : auto sol = problems::SOL();
83 : :
84 : : // Evaluate analytic solution (if defined)
85 [ + - ]: 30062 : auto an = u;
86 [ + + ]: 30062 : 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 [ + + ]: 4579430 : for (std::size_t i=0; i<u.nunk(); ++i) {
103 : : // Compute sum for L2 norm of the numerical solution
104 [ + + ]: 27706248 : for (std::size_t c=0; c<ncomp; ++c)
105 [ + - ][ + - ]: 23156880 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
106 : : // Compute sum for L2 norm of the residual
107 [ + + ]: 27706248 : for (std::size_t c=0; c<ncomp; ++c)
108 [ + - ][ + - ]: 23156880 : 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 [ + - ]: 4549368 : diag[TOTALEN][0] += u(i,4) * v[i];
111 : : // Compute sum for L2 norm of the numerical-analytic solution
112 [ + + ]: 4549368 : 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 : 30062 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
136 : 30062 : diag[TIME][0] = d.T() + d.Dt();
137 : 30062 : diag[DT][0] = d.Dt();
138 : :
139 : : // Contribute to diagnostics
140 [ + - ]: 30062 : auto stream = serialize( d.MeshId(), ncomp, diag );
141 [ + - ]: 30062 : d.contribute( stream.first, stream.second.get(), DiagMerger,
142 [ + - ][ + - ]: 60124 : CkCallback(CkIndex_Transporter::rhodiagnostics(nullptr), d.Tr()) );
143 : :
144 : 30062 : return true; // diagnostics have been computed
145 : 30062 : }
146 : :
147 : : bool
148 : 3662 : 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 [ - + ]: 3662 : if ( (d.It()+1) % diag_iter ) return false;
177 : :
178 [ - + ][ - - ]: 3662 : Assert( p.size() == u.nunk(), "Size mismatch" );
[ - - ][ - - ]
179 [ - + ][ - - ]: 3662 : Assert( p.size() == dp.size(), "Size mismatch" );
[ - - ][ - - ]
180 [ - + ][ - - ]: 3662 : Assert( u.nunk() == un.nunk(), "Size mismatch" );
[ - - ][ - - ]
181 : :
182 : 3662 : auto ncomp = u.nprop();
183 : :
184 : 3662 : const auto& v = d.V(); // nodal volumes without contributions from others
185 : :
186 : : // query function to evaluate analytic solution (if defined)
187 [ + - ]: 3662 : auto pressure_sol = problems::PRESSURE_SOL();
188 : :
189 : : // Evaluate analytic solution (if defined)
190 [ + - ]: 3662 : auto an = p;
191 [ + + ]: 3662 : 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 : an[i] = pressure_sol( x[i], y[i], z[i] );
199 : : }
200 : : }
201 : :
202 : : // Diagnostics vector (of vectors) during aggregation. See
203 : : // Inciter/Diagnostics.h.
204 : : std::vector< std::vector< tk::real > >
205 [ + - ][ + - ]: 3662 : diag( NUMDIAG, std::vector< tk::real >( ncomp+1, 0.0 ) );
206 : :
207 : : // Put in norms sweeping our mesh chunk
208 [ + + ]: 361004 : for (std::size_t i=0; i<u.nunk(); ++i) {
209 : : // Compute sum for L2 norm of the numerical solution
210 : 357342 : diag[L2SOL][0] += p[i] * p[i] * v[i];
211 [ + + ]: 1422942 : for (std::size_t c=0; c<ncomp; ++c)
212 [ + - ][ + - ]: 1065600 : diag[L2SOL][c+1] += u(i,c) * u(i,c) * v[i];
213 : : // Compute sum for L2 norm of the residual
214 : 357342 : diag[L2RES][0] += dp[i] * dp[i] * v[i];
215 [ + + ]: 1422942 : for (std::size_t c=0; c<ncomp; ++c)
216 [ + - ][ + - ]: 1065600 : diag[L2RES][c+1] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
217 : : // Compute sum for the total energy over the entire domain
218 : 357342 : diag[TOTALEN][0] += 0.0 * v[i];
219 : : // Compute sum for L2 norm of the numerical-analytic solution
220 [ + + ]: 357342 : if (pressure_sol) {
221 : 2142 : auto pd = p[i] - an[i];
222 : 2142 : diag[L2ERR][0] += pd * pd * v[i];
223 : 2142 : diag[L1ERR][0] += std::abs( pd ) * v[i];
224 : : }
225 : : }
226 : :
227 : : // Append diagnostics vector with metadata on the current time step
228 : : // ITER:: Current iteration count (only the first entry is used)
229 : : // TIME: Current physical time (only the first entry is used)
230 : : // DT: Current physical time step size (only the first entry is used)
231 : 3662 : diag[ITER][0] = static_cast< tk::real >( d.It() );
232 : 3662 : diag[TIME][0] = d.T();
233 : 3662 : diag[DT][0] = d.Dt();
234 : :
235 : : // Contribute to diagnostics
236 [ + - ]: 3662 : auto stream = serialize( d.MeshId(), ncomp+1, diag );
237 [ + - ]: 3662 : d.contribute( stream.first, stream.second.get(), DiagMerger,
238 [ + - ][ + - ]: 7324 : CkCallback(CkIndex_Transporter::prediagnostics(nullptr), d.Tr()) );
239 : :
240 : 3662 : return true; // diagnostics have been computed
241 : 3662 : }
242 : :
243 : : bool
244 : 3530 : NodeDiagnostics::accompute( Discretization& d,
245 : : const tk::Fields& u,
246 : : const tk::Fields& un,
247 : : uint64_t diag_iter ) const
248 : : // *****************************************************************************
249 : : // Compute diagnostics for artificial compressibility solvers
250 : : //! \param[in] d Discretization proxy to read from
251 : : //! \param[in] u Current solution vector
252 : : //! \param[in] un Previous solution vector
253 : : //! \param[in] diag_iter Diagnostics output frequency
254 : : //! \return True if diagnostics have been computed
255 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
256 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
257 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
258 : : //! We send multiple sets of quantities to the host for aggregation across
259 : : //! the whole mesh. The final aggregated solution will end up in
260 : : //! Transporter::diagnostics(). Aggregation of the partially computed
261 : : //! diagnostics is done via potentially different policies for each field.
262 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
263 : : // *****************************************************************************
264 : : {
265 : : using namespace diagnostics;
266 : :
267 : : // Only compute diagnostics if needed in this time step
268 [ - + ]: 3530 : if ( (d.It()+1) % diag_iter ) return false;
269 : :
270 : 3530 : auto ncomp = u.nprop();
271 : :
272 : : // Diagnostics vector (of vectors) during aggregation. See
273 : : // Inciter/Diagnostics.h.
274 : : std::vector< std::vector< tk::real > >
275 [ + - ][ + - ]: 3530 : diag( NUMDIAG, std::vector< tk::real >( ncomp, 0.0 ) );
276 : :
277 : 3530 : const auto& v = d.V(); // nodal volumes without contributions from others
278 : :
279 : : // Put in norms sweeping our mesh chunk
280 [ + + ]: 246420 : for (std::size_t i=0; i<u.nunk(); ++i) {
281 : : // Compute sum for L2 norm of the numerical solution
282 [ + + ]: 1214450 : for (std::size_t c=0; c<ncomp; ++c)
283 [ + - ][ + - ]: 971560 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
284 : : // Compute sum for L2 norm of the residual
285 [ + + ]: 1214450 : for (std::size_t c=0; c<ncomp; ++c)
286 [ + - ][ + - ]: 971560 : diag[L2RES][c] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
287 : : // Compute sum for the total energy over the entire domain
288 : 242890 : diag[TOTALEN][0] += 0.0 * v[i];
289 : : }
290 : :
291 : : // Append diagnostics vector with metadata on the current time step
292 : : // ITER:: Current iteration count (only the first entry is used)
293 : : // TIME: Current physical time (only the first entry is used)
294 : : // DT: Current physical time step size (only the first entry is used)
295 : 3530 : diag[ITER][0] = static_cast< tk::real >( d.It() );
296 : 3530 : diag[TIME][0] = d.T();
297 : 3530 : diag[DT][0] = d.Dt();
298 : :
299 : : // Contribute to diagnostics
300 [ + - ]: 3530 : auto stream = serialize( d.MeshId(), ncomp, diag );
301 [ + - ]: 3530 : d.contribute( stream.first, stream.second.get(), DiagMerger,
302 [ + - ][ + - ]: 7060 : CkCallback(CkIndex_Transporter::acdiagnostics(nullptr), d.Tr()) );
303 : :
304 : 3530 : return true; // diagnostics have been computed
305 : 3530 : }
|
