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 : 1156 : 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 : 1156 : DiagMerger = CkReduction::addReducer( mergeDiag );
46 : 1156 : }
47 : :
48 : : bool
49 : 16667 : 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,
57 : : const std::unordered_set< std::size_t >& slipwallbcnodes ) const
58 : : // *****************************************************************************
59 : : // Compute diagnostics, e.g., residuals, norms of errors, etc.
60 : : //! \param[in] d Discretization proxy to read from
61 : : //! \param[in] u Current solution vector
62 : : //! \param[in] un Previous solution vector
63 : : //! \param[in] bnorm Face normals in boundary points, key local node id,
64 : : //! first 3 reals of value: unit normal, outer key: side set id
65 : : //! \param[in] symbcnodes Unique set of node ids at which to set symmetry BCs
66 : : //! \param[in] farfieldbcnodes Unique set of node ids at which to set farfield
67 : : //! BCs
68 : : //! \param[in] slipwallbcnodes Unique set of node ids at which to set slip BCs
69 : : //! \return True if diagnostics have been computed
70 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
71 : : //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
72 : : //! where the sum is taken over all mesh nodes and V_i is the nodal volume.
73 : : //! We send multiple sets of quantities to the host for aggregation across
74 : : //! the whole mesh. The final aggregated solution will end up in
75 : : //! Transporter::diagnostics(). Aggregation of the partially computed
76 : : //! diagnostics is done via potentially different policies for each field.
77 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.hpp
78 : : // *****************************************************************************
79 : : {
80 : : // Optionally collect diagnostics and send for aggregation across all workers
81 : :
82 : : // Query after how many time steps user wants to dump diagnostics
83 : 16667 : auto diagfreq = g_inputdeck.get< tag::diagnostics, tag::interval >();
84 : :
85 [ + + ]: 16667 : if ( !((d.It()+1) % diagfreq) ) { // if remainder, don't dump
86 : :
87 : : // Diagnostics vector (of vectors) during aggregation. See
88 : : // Inciter/Diagnostics.h.
89 : : std::vector< std::vector< tk::real > >
90 [ + - ][ + - ]: 48663 : diag( NUMDIAG, std::vector< tk::real >( u.nprop(), 0.0 ) );
[ + - ]
91 : :
92 : : const auto& coord = d.Coord();
93 : : const auto& x = coord[0];
94 : : const auto& y = coord[1];
95 : : const auto& z = coord[2];
96 : : const auto& v = d.V(); // nodal volumes without contributions from others
97 : :
98 : : // Evaluate analytic solution (if exist, if not, IC)
99 : : auto an = u;
100 : : auto mv = d.MeshVel();
101 [ + + ]: 1472431 : for (std::size_t i=0; i<an.nunk(); ++i) {
102 : : // Query analytic solution for all components of all PDEs integrated
103 : : std::vector< tk::real > a;
104 [ + - ]: 1456210 : auto s = g_cgpde[d.MeshId()].solution( x[i], y[i], z[i], d.T()+d.Dt() );
105 : : std::move( begin(s), end(s), std::back_inserter(a) );
106 : : Assert( a.size() == u.nprop(), "Size mismatch" );
107 [ + + ]: 5985644 : for (std::size_t c=0; c<an.nprop(); ++c) an(i,c) = a[c];
108 : : }
109 : : // Apply symmetry BCs on analytic solution (if exist, if not, IC)
110 [ + - ]: 16221 : g_cgpde[d.MeshId()].symbc( an, coord, bnorm, symbcnodes );
111 : : // Apply farfield BCs on analytic solution (if exist, if not, IC)
112 [ + - ]: 16221 : g_cgpde[d.MeshId()].farfieldbc( an, coord, bnorm, farfieldbcnodes );
113 : : // Apply slip wall BCs on analytic solution (if exist, if not, IC)
114 [ + - ]: 16221 : g_cgpde[d.MeshId()].slipwallbc( an, mv, coord, bnorm, slipwallbcnodes );
115 : :
116 : : // Put in norms sweeping our mesh chunk
117 [ + + ]: 1472431 : for (std::size_t i=0; i<u.nunk(); ++i) {
118 : : // Compute sum for L2 norm of the numerical solution
119 [ + + ]: 5985644 : for (std::size_t c=0; c<u.nprop(); ++c)
120 : 4529434 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
121 : : // Compute sum for L2 norm of the numerical-analytic solution
122 [ + + ]: 5985644 : for (std::size_t c=0; c<u.nprop(); ++c)
123 : 4529434 : diag[L2ERR][c] += (u(i,c)-an(i,c)) * (u(i,c)-an(i,c)) * v[i];
124 : : // Compute sum for L2 norm of the residual
125 [ + + ]: 5985644 : for (std::size_t c=0; c<u.nprop(); ++c)
126 : 4529434 : diag[L2RES][c] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
127 : : // Compute max for Linf norm of the numerical-analytic solution
128 [ + + ]: 5985644 : for (std::size_t c=0; c<u.nprop(); ++c) {
129 [ + + ]: 4529434 : auto err = std::abs( u(i,c) - an(i,c) );
130 [ + + ]: 4529434 : if (err > diag[LINFERR][c]) diag[LINFERR][c] = err;
131 : : }
132 : : // Compute sum of the total energy over the entire domain (only the first
133 : : // entry is used)
134 : 1456210 : diag[TOTALSOL][0] += u(i,u.nprop()-1) * v[i];
135 : : }
136 : :
137 : : // Append diagnostics vector with metadata on the current time step
138 : : // ITER:: Current iteration count (only the first entry is used)
139 : : // TIME: Current physical time (only the first entry is used)
140 : : // DT: Current physical time step size (only the first entry is used)
141 : 16221 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
142 : 16221 : diag[TIME][0] = d.T() + d.Dt();
143 : 16221 : diag[DT][0] = d.Dt();
144 : :
145 : : // Contribute to diagnostics
146 [ + - ]: 16221 : auto stream = serialize( d.MeshId(), u.nprop(), diag );
147 [ + - ][ + - ]: 32442 : d.contribute( stream.first, stream.second.get(), DiagMerger,
148 [ + - ][ + - ]: 32442 : CkCallback(CkIndex_Transporter::diagnostics(nullptr), d.Tr()) );
[ - - ]
149 : :
150 : : return true; // diagnostics have been computed
151 : : }
152 : :
153 : : return false; // diagnostics have not been computed
154 : : }
|
