Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ElemDiagnostics.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 ElemDiagnostics class for collecting element diagnostics
9 : : \details ElemDiagnostics class for collecting element diagnostics, e.g.,
10 : : residuals, and various norms of errors while solving partial differential
11 : : equations.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <array>
16 : : #include <vector>
17 : : #include <cmath>
18 : :
19 : : #include "DGPDE.hpp"
20 : : #include "ElemDiagnostics.hpp"
21 : : #include "DiagReducer.hpp"
22 : : #include "Discretization.hpp"
23 : : #include "Integrate/Basis.hpp"
24 : : #include "Integrate/Quadrature.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : :
27 : : namespace inciter {
28 : :
29 : : extern ctr::InputDeck g_inputdeck;
30 : : extern std::vector< DGPDE > g_dgpde;
31 : :
32 : : static CkReduction::reducerType DiagMerger;
33 : :
34 : : } // inciter::
35 : :
36 : : using inciter::ElemDiagnostics;
37 : :
38 : : void
39 : 1174 : ElemDiagnostics::registerReducers()
40 : : // *****************************************************************************
41 : : // Configure Charm++ reduction types
42 : : //! \details This routine is supposed to be called from a Charm++ initnode
43 : : //! routine. Since the runtime system executes initnode routines exactly once
44 : : //! on every logical node early on in the Charm++ init sequence, they must be
45 : : //! static as they are called without an object. See also: Section
46 : : //! "Initializations at Program Startup" at in the Charm++ manual
47 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
48 : : // *****************************************************************************
49 : : {
50 : 1174 : DiagMerger = CkReduction::addReducer( mergeDiag );
51 : 1174 : }
52 : :
53 : : bool
54 : 24969 : ElemDiagnostics::compute( Discretization& d,
55 : : const std::size_t nchGhost,
56 : : const tk::Fields& geoElem,
57 : : const std::vector< std::size_t >& ndofel,
58 : : const tk::Fields& u,
59 : : const tk::Fields& un ) const
60 : : // *****************************************************************************
61 : : // Compute diagnostics, e.g., residuals, norms of errors, etc.
62 : : //! \param[in] d Discretization base class to read from
63 : : //! \param[in] nchGhost Number of chare boundary ghost elements
64 : : //! \param[in] geoElem Element geometry
65 : : //! \param[in] ndofel Vector of local number of degrees of freedom
66 : : //! \param[in] u Current solution vector
67 : : //! \param[in] un Previous time-step solution vector
68 : : //! \return True if diagnostics have been computed
69 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
70 : : //! computed in mesh elements, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
71 : : //! where the sum is taken over all mesh elements and V_i is the cell volume.
72 : : //! We send multiple sets of quantities to the host for aggregation across
73 : : //! the whole mesh. The final aggregated solution will end up in
74 : : //! Transporter::diagnostics(). Aggregation of the partially computed
75 : : //! diagnostics is done via potentially different policies for each field.
76 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.h
77 : : // *****************************************************************************
78 : : {
79 : : // Optionally collect diagnostics and send for aggregation across all workers
80 : :
81 : : // Query after how many time steps user wants to dump diagnostics
82 : 24969 : auto diagfreq = g_inputdeck.get< tag::diagnostics, tag::interval >();
83 : :
84 [ + + ][ + + ]: 24969 : if ( !((d.It()+1) % diagfreq) || d.finished() ) { // if remainder, don't compute diagnostics
85 : :
86 : : // Query number of degrees of freedom from user's setting
87 : 13773 : const auto rdof = g_inputdeck.get< tag::rdof >();
88 : :
89 : : // Diagnostics vector (of vectors) during aggregation. See
90 : : // Inciter/Diagnostics.h.
91 : : std::vector< std::vector< tk::real > >
92 [ + - ][ + - ]: 27546 : diag( NUMDIAG, std::vector< tk::real >( u.nprop()/rdof, 0.0 ) );
93 : :
94 : : // Compute diagnostics for DG
95 [ + - ]: 13773 : compute_diag(d, rdof, nchGhost, geoElem, ndofel, u, un, diag);
96 : :
97 : : // Append diagnostics vector with metadata on the current time step
98 : : // ITER: Current iteration count (only the first entry is used)
99 : : // TIME: Current physical time (only the first entry is used)
100 : : // DT: Current physical time step size (only the first entry is used)
101 : 13773 : diag[ITER][0] = static_cast< tk::real >( d.It() );
102 : 13773 : diag[TIME][0] = d.T();
103 : 13773 : diag[DT][0] = d.Dt();
104 : :
105 : : // Contribute to diagnostics
106 [ + - ]: 13773 : auto stream = serialize( d.MeshId(), u.nprop()/rdof, diag );
107 [ + - ][ + - ]: 27546 : d.contribute( stream.first, stream.second.get(), DiagMerger,
108 [ + - ][ + - ]: 27546 : CkCallback(CkIndex_Transporter::diagnostics(nullptr), d.Tr()) );
[ - - ]
109 : :
110 : : return true; // diagnostics have been computed
111 : :
112 : : }
113 : :
114 : : return false; // diagnostics have not been computed
115 : : }
116 : :
117 : : void
118 : 13773 : ElemDiagnostics::compute_diag( const Discretization& d,
119 : : const std::size_t rdof,
120 : : const std::size_t nchGhost,
121 : : const tk::Fields& geoElem,
122 : : const std::vector< std::size_t >& ndofel,
123 : : const tk::Fields& u,
124 : : const tk::Fields& un,
125 : : std::vector< std::vector< tk::real > >& diag )
126 : : const
127 : : // *****************************************************************************
128 : : // Compute diagnostics, e.g., residuals, norms of errors, etc. for DG
129 : : //! \param[in] d Discretization base class to read from
130 : : //! \param[in] rdof Number of reconstructed degrees of freedom
131 : : //! \param[in] nchGhost Number of chare boundary ghost elements
132 : : //! \param[in] geoElem Element geometry
133 : : //! \param[in] ndofel Vector of local number of degrees of freedom
134 : : //! \param[in] u Current solution vector
135 : : //! \param[in] un Previous time-step solution vector
136 : : //! \param[in,out] diag Diagnostics vector
137 : : // *****************************************************************************
138 : : {
139 : 13773 : const auto& inpoel = d.Inpoel();
140 : : const auto& coord = d.Coord();
141 : :
142 : : const auto& cx = coord[0];
143 : : const auto& cy = coord[1];
144 : : const auto& cz = coord[2];
145 : :
146 [ + + ]: 3622784 : for (std::size_t e=0; e<u.nunk()-nchGhost; ++e)
147 : : {
148 : : std::size_t dofe(1);
149 [ + + ]: 3609011 : if (!ndofel.empty()) {
150 : 3534667 : dofe = ndofel[e];
151 : : }
152 : : // Number of quadrature points for volume integration
153 : 3609011 : auto ng = tk::NGdiag(dofe);
154 : :
155 : : // arrays for quadrature points
156 : : std::array< std::vector< tk::real >, 3 > coordgp;
157 : : std::vector< tk::real > wgp;
158 : :
159 [ + - ]: 3609011 : coordgp[0].resize( ng );
160 [ + - ]: 3609011 : coordgp[1].resize( ng );
161 [ + - ]: 3609011 : coordgp[2].resize( ng );
162 [ + - ]: 3609011 : wgp.resize( ng );
163 : :
164 [ + - ]: 3609011 : tk::GaussQuadratureTet( ng, coordgp, wgp );
165 : :
166 : : // Extract the element coordinates
167 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
168 : 3609011 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
169 : 3609011 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
170 : 3609011 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
171 : 3609011 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
172 : :
173 [ + + ]: 10069376 : for (std::size_t igp=0; igp<ng; ++igp)
174 : : {
175 : : // Compute the coordinates of quadrature point at physical domain
176 [ + - ]: 6460365 : auto gp = tk::eval_gp( igp, coordel, coordgp );
177 : :
178 : : // Compute the basis function
179 [ + - ]: 6460365 : auto B = tk::eval_basis( dofe, coordgp[0][igp], coordgp[1][igp],
180 [ + - ]: 6460365 : coordgp[2][igp]);
181 : :
182 [ + - ]: 6460365 : auto wt = wgp[igp] * geoElem(e, 0);
183 : :
184 : : std::vector< tk::real > s;
185 : :
186 : : // cppcheck-suppress useStlAlgorithm
187 [ + - ]: 6460365 : s = g_dgpde[d.MeshId()].solution( gp[0], gp[1], gp[2], d.T() );
188 : :
189 [ + + ]: 35460626 : for (std::size_t c=0; c<u.nprop()/rdof; ++c)
190 : : {
191 [ + + ]: 29000261 : auto mark = c*rdof;
192 : 29000261 : auto ugp = u(e, mark);
193 : :
194 [ + + ]: 29000261 : if(dofe > 1)
195 : : {
196 [ + + ]: 15750434 : ugp += u(e, mark+1) * B[1]
197 : 15750434 : + u(e, mark+2) * B[2]
198 : 15750434 : + u(e, mark+3) * B[3];
199 : :
200 [ + + ]: 15750434 : if(dofe > 4)
201 : : {
202 : 4817974 : ugp += u(e, mark+4) * B[4]
203 : 4817974 : + u(e, mark+5) * B[5]
204 : 4817974 : + u(e, mark+6) * B[6]
205 : 4817974 : + u(e, mark+7) * B[7]
206 : 4817974 : + u(e, mark+8) * B[8]
207 : 4817974 : + u(e, mark+9) * B[9];
208 : : }
209 : : }
210 : :
211 : : // Compute sum for L2 norm of the numerical solution
212 [ + + ]: 29000261 : diag[L2SOL][c] += wt * ugp * ugp;
213 : :
214 : : // Compute sum for L2 norm of the numerical-analytic solution
215 : 29000261 : diag[L2ERR][c] += wt * (ugp-s[c]) * (ugp-s[c]);
216 : :
217 : : // Compute max for Linf norm of the numerical-analytic solution
218 [ + + ]: 29000261 : auto err = std::abs( ugp - s[c] );
219 [ + + ]: 29000261 : if (err > diag[LINFERR][c]) diag[LINFERR][c] = err;
220 : : }
221 : : }
222 : : tk::real sp_te(0.0);
223 : : // cppcheck-suppress useStlAlgorithm
224 [ + - ]: 3609011 : sp_te += g_dgpde[d.MeshId()].sp_totalenergy(e, u);
225 : :
226 : : // Compute sum of the total energy over the entire domain (only the
227 : : // first entry is used)
228 : 3609011 : diag[TOTALSOL][0] += geoElem(e,0) * sp_te;
229 : :
230 : : // Compute sum for L2 norm of the cell-avg residual
231 [ + + ]: 19936094 : for (std::size_t c=0; c<u.nprop()/rdof; ++c)
232 : 16327083 : diag[L2RES][c] += geoElem(e, 0) *
233 : 16327083 : ( (u(e,c*rdof)-un(e,c*rdof)) * (u(e,c*rdof)-un(e,c*rdof)) );
234 : : }
235 : 13773 : }
|
