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 : : All rights reserved. See the LICENSE file for details.
8 : : \brief NodeDiagnostics class for collecting nodal diagnostics
9 : : \details NodeDiagnostics class for collecting nodal diagnostics, e.g.,
10 : : residuals, and various norms of errors while solving partial differential
11 : : equations.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include "CGPDE.hpp"
16 : : #include "NodeDiagnostics.hpp"
17 : : #include "DiagReducer.hpp"
18 : : #include "Discretization.hpp"
19 : : #include "Inciter/InputDeck/InputDeck.hpp"
20 : : #include "Refiner.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : : extern std::vector< CGPDE > g_cgpde;
26 : :
27 : : static CkReduction::reducerType DiagMerger;
28 : :
29 : : } // inciter::
30 : :
31 : : using inciter::NodeDiagnostics;
32 : :
33 : : void
34 : 1332 : NodeDiagnostics::registerReducers()
35 : : // *****************************************************************************
36 : : // Configure Charm++ reduction types
37 : : //! \details This routine is supposed to be called from a Charm++ initnode
38 : : //! routine. Since the runtime system executes initnode routines exactly once
39 : : //! on every logical node early on in the Charm++ init sequence, they must be
40 : : //! static as they are called without an object. See also: Section
41 : : //! "Initializations at Program Startup" at in the Charm++ manual
42 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
43 : : // *****************************************************************************
44 : : {
45 : 1332 : DiagMerger = CkReduction::addReducer( mergeDiag );
46 : 1332 : }
47 : :
48 : : bool
49 : 33885 : NodeDiagnostics::compute(
50 : : Discretization& d,
51 : : const tk::Fields& u,
52 : : const tk::Fields& un,
53 : : const std::unordered_map< int,
54 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
55 : : const std::unordered_set< std::size_t >& symbcnodes,
56 : : const std::unordered_set< std::size_t >& farfieldbcnodes ) const
57 : : // *****************************************************************************
58 : : // Compute diagnostics, e.g., residuals, norms of errors, etc.
59 : : //! \param[in] d Discretization proxy to read from
60 : : //! \param[in] u Current solution vector
61 : : //! \param[in] un Previous solution vector
62 : : //! \param[in] bnorm Face normals in boundary points, key local node id,
63 : : //! first 3 reals of value: unit normal, outer key: side set id
64 : : //! \param[in] symbcnodes Unique set of node ids at which to set symmetry BCs
65 : : //! \param[in] farfieldbcnodes Unique set of node ids at which to set farfield
66 : : //! BCs
67 : : //! \return True if diagnostics have been computed
68 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
69 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
70 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
71 : : //! We send multiple sets of quantities to the host for aggregation across
72 : : //! the whole mesh. The final aggregated solution will end up in
73 : : //! Transporter::diagnostics(). Aggregation of the partially computed
74 : : //! diagnostics is done via potentially different policies for each field.
75 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
76 : : // *****************************************************************************
77 : : {
78 : : // Optionally collect diagnostics and send for aggregation across all workers
79 : :
80 : : // Query after how many time steps user wants to dump diagnostics
81 : 33885 : auto diagfreq = g_inputdeck.get< tag::output, tag::iter, tag::diag >();
82 : :
83 [ + + ]: 33885 : if ( !((d.It()+1) % diagfreq) ) { // if remainder, don't dump
84 : :
85 : : // Diagnostics vector (of vectors) during aggregation. See
86 : : // Inciter/Diagnostics.h.
87 : : std::vector< std::vector< tk::real > >
88 [ + - ][ + - ]: 98421 : diag( NUMDIAG, std::vector< tk::real >( u.nprop(), 0.0 ) );
89 : :
90 : 32807 : const auto& coord = d.Coord();
91 : 32807 : const auto& x = coord[0];
92 : 32807 : const auto& y = coord[1];
93 : 32807 : const auto& z = coord[2];
94 : 32807 : const auto& v = d.V(); // nodal volumes without contributions from others
95 : :
96 : : // Evaluate analytic solution (if exist, if not, IC)
97 [ + - ]: 65614 : auto an = u;
98 [ + + ]: 5602521 : for (std::size_t i=0; i<an.nunk(); ++i) {
99 : : // Query analytic solution for all components of all PDEs integrated
100 : 11139428 : std::vector< tk::real > a;
101 [ + + ]: 11139428 : for (const auto& eq : g_cgpde) {
102 [ + - ]: 11139428 : auto s = eq.solution( x[i], y[i], z[i], d.T()+d.Dt() );
103 [ + - ][ + - ]: 5569714 : std::move( begin(s), end(s), std::back_inserter(a) );
104 : : }
105 [ - + ][ - - ]: 5569714 : Assert( a.size() == u.nprop(), "Size mismatch" );
[ - - ][ - - ]
106 [ + + ][ + - ]: 16971824 : for (std::size_t c=0; c<an.nprop(); ++c) an(i,c,0) = a[c];
107 : : }
108 : : // Apply symmetry BCs on analytic solution (if exist, if not, IC)
109 [ + + ]: 65614 : for (const auto& eq : g_cgpde)
110 [ + - ]: 32807 : eq.symbc( an, coord, bnorm, symbcnodes );
111 : : // Apply farfield BCs on analytic solution (if exist, if not, IC)
112 [ + + ]: 65614 : for (const auto& eq : g_cgpde)
113 [ + - ]: 32807 : eq.farfieldbc( an, coord, bnorm, farfieldbcnodes );
114 : :
115 : : // Put in norms sweeping our mesh chunk
116 [ + + ]: 5602521 : for (std::size_t i=0; i<u.nunk(); ++i) {
117 : : // Compute sum for L2 norm of the numerical solution
118 [ + + ]: 16971824 : for (std::size_t c=0; c<u.nprop(); ++c)
119 [ + - ][ + - ]: 11402110 : diag[L2SOL][c] += u(i,c,0) * u(i,c,0) * v[i];
120 : : // Compute sum for L2 norm of the numerical-analytic solution
121 [ + + ]: 16971824 : for (std::size_t c=0; c<u.nprop(); ++c)
122 [ + - ][ + - ]: 11402110 : diag[L2ERR][c] += (u(i,c,0)-an(i,c,0)) * (u(i,c,0)-an(i,c,0)) * v[i];
[ + - ][ + - ]
123 : : // Compute sum for L2 norm of the residual
124 [ + + ]: 16971824 : for (std::size_t c=0; c<u.nprop(); ++c)
125 [ + - ][ + - ]: 11402110 : diag[L2RES][c] += (u(i,c,0)-un(i,c,0)) * (u(i,c,0)-un(i,c,0)) * v[i];
[ + - ][ + - ]
126 : : // Compute max for Linf norm of the numerical-analytic solution
127 [ + + ]: 16971824 : for (std::size_t c=0; c<u.nprop(); ++c) {
128 [ + - ][ + - ]: 11402110 : auto err = std::abs( u(i,c,0) - an(i,c,0) );
129 [ + + ]: 11402110 : if (err > diag[LINFERR][c]) diag[LINFERR][c] = err;
130 : : }
131 : : // Compute sum of the total energy over the entire domain (only the first
132 : : // entry is used)
133 [ + - ]: 5569714 : diag[TOTALSOL][0] += u(i,u.nprop()-1,0) * v[i];
134 : : }
135 : :
136 : : // Append diagnostics vector with metadata on the current time step
137 : : // ITER:: Current iteration count (only the first entry is used)
138 : : // TIME: Current physical time (only the first entry is used)
139 : : // DT: Current physical time step size (only the first entry is used)
140 : 32807 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
141 : 32807 : diag[TIME][0] = d.T() + d.Dt();
142 : 32807 : diag[DT][0] = d.Dt();
143 : :
144 : : // Contribute to diagnostics
145 [ + - ]: 32807 : auto stream = serialize( d.MeshId(), diag );
146 [ + - ]: 32807 : d.contribute( stream.first, stream.second.get(), DiagMerger,
147 [ + - ][ + - ]: 65614 : CkCallback(CkIndex_Transporter::diagnostics(nullptr), d.Tr()) );
148 : :
149 : 32807 : return true; // diagnostics have been computed
150 : : }
151 : :
152 : 1078 : return false; // diagnostics have not been computed
153 : : }
|
