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 : 1144 : 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 : 1144 : DiagMerger = CkReduction::addReducer( mergeDiag );
46 : 1144 : }
47 : :
48 : : bool
49 : 17779 : 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 : 17779 : auto diagfreq = g_inputdeck.get< tag::diagnostics, tag::interval >();
82 : :
83 [ + + ]: 17779 : 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 [ + - ][ + - ]: 50919 : diag( NUMDIAG, std::vector< tk::real >( u.nprop(), 0.0 ) );
89 : :
90 : 16973 : const auto& coord = d.Coord();
91 : 16973 : const auto& x = coord[0];
92 : 16973 : const auto& y = coord[1];
93 : 16973 : const auto& z = coord[2];
94 : 16973 : const auto& v = d.V(); // nodal volumes without contributions from others
95 : :
96 : : // Evaluate analytic solution (if exist, if not, IC)
97 [ + - ]: 33946 : auto an = u;
98 [ + + ]: 1537069 : for (std::size_t i=0; i<an.nunk(); ++i) {
99 : : // Query analytic solution for all components of all PDEs integrated
100 : 3040192 : std::vector< tk::real > a;
101 [ + - ]: 3040192 : auto s = g_cgpde[d.MeshId()].solution( x[i], y[i], z[i], d.T()+d.Dt() );
102 [ + - ][ + - ]: 1520096 : std::move( begin(s), end(s), std::back_inserter(a) );
103 [ - + ][ - - ]: 1520096 : Assert( a.size() == u.nprop(), "Size mismatch" );
[ - - ][ - - ]
104 [ + + ][ + - ]: 6368960 : for (std::size_t c=0; c<an.nprop(); ++c) an(i,c) = a[c];
105 : : }
106 : : // Apply symmetry BCs on analytic solution (if exist, if not, IC)
107 [ + - ]: 16973 : g_cgpde[d.MeshId()].symbc( an, coord, bnorm, symbcnodes );
108 : : // Apply farfield BCs on analytic solution (if exist, if not, IC)
109 [ + - ]: 16973 : g_cgpde[d.MeshId()].farfieldbc( an, coord, bnorm, farfieldbcnodes );
110 : :
111 : : // Put in norms sweeping our mesh chunk
112 [ + + ]: 1537069 : for (std::size_t i=0; i<u.nunk(); ++i) {
113 : : // Compute sum for L2 norm of the numerical solution
114 [ + + ]: 6368960 : for (std::size_t c=0; c<u.nprop(); ++c)
115 [ + - ][ + - ]: 4848864 : diag[L2SOL][c] += u(i,c) * u(i,c) * v[i];
116 : : // Compute sum for L2 norm of the numerical-analytic solution
117 [ + + ]: 6368960 : for (std::size_t c=0; c<u.nprop(); ++c)
118 [ + - ][ + - ]: 4848864 : diag[L2ERR][c] += (u(i,c)-an(i,c)) * (u(i,c)-an(i,c)) * v[i];
[ + - ][ + - ]
119 : : // Compute sum for L2 norm of the residual
120 [ + + ]: 6368960 : for (std::size_t c=0; c<u.nprop(); ++c)
121 [ + - ][ + - ]: 4848864 : diag[L2RES][c] += (u(i,c)-un(i,c)) * (u(i,c)-un(i,c)) * v[i];
[ + - ][ + - ]
122 : : // Compute max for Linf norm of the numerical-analytic solution
123 [ + + ]: 6368960 : for (std::size_t c=0; c<u.nprop(); ++c) {
124 [ + - ][ + - ]: 4848864 : auto err = std::abs( u(i,c) - an(i,c) );
125 [ + + ]: 4848864 : if (err > diag[LINFERR][c]) diag[LINFERR][c] = err;
126 : : }
127 : : // Compute sum of the total energy over the entire domain (only the first
128 : : // entry is used)
129 [ + - ]: 1520096 : diag[TOTALSOL][0] += u(i,u.nprop()-1) * v[i];
130 : : }
131 : :
132 : : // Append diagnostics vector with metadata on the current time step
133 : : // ITER:: Current iteration count (only the first entry is used)
134 : : // TIME: Current physical time (only the first entry is used)
135 : : // DT: Current physical time step size (only the first entry is used)
136 : 16973 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
137 : 16973 : diag[TIME][0] = d.T() + d.Dt();
138 : 16973 : diag[DT][0] = d.Dt();
139 : :
140 : : // Contribute to diagnostics
141 [ + - ]: 16973 : auto stream = serialize( d.MeshId(), u.nprop(), diag );
142 [ + - ]: 16973 : d.contribute( stream.first, stream.second.get(), DiagMerger,
143 [ + - ][ + - ]: 33946 : CkCallback(CkIndex_Transporter::diagnostics(nullptr), d.Tr()) );
144 : :
145 : 16973 : return true; // diagnostics have been computed
146 : : }
147 : :
148 : 806 : return false; // diagnostics have not been computed
149 : : }
|
