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