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 : 666 : 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 : 666 : DiagMerger = CkReduction::addReducer( mergeDiag );
51 : 666 : }
52 : :
53 : : bool
54 : 24860 : 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 ) const
59 : : // *****************************************************************************
60 : : // Compute diagnostics, e.g., residuals, norms of errors, etc.
61 : : //! \param[in] d Discretization base class to read from
62 : : //! \param[in] nchGhost Number of chare boundary ghost elements
63 : : //! \param[in] geoElem Element geometry
64 : : //! \param[in] ndofel Vector of local number of degrees of freedom
65 : : //! \param[in] u Current solution vector
66 : : //! \return True if diagnostics have been computed
67 : : //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity,
68 : : //! computed in mesh elements, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ],
69 : : //! where the sum is taken over all mesh elements and V_i is the cell volume.
70 : : //! We send multiple sets of quantities to the host for aggregation across
71 : : //! the whole mesh. The final aggregated solution will end up in
72 : : //! Transporter::diagnostics(). Aggregation of the partially computed
73 : : //! diagnostics is done via potentially different policies for each field.
74 : : //! \see inciter::mergeDiag(), src/Inciter/Diagnostics.h
75 : : // *****************************************************************************
76 : : {
77 : : // Optionally collect diagnostics and send for aggregation across all workers
78 : :
79 : : // Query after how many time steps user wants to dump diagnostics
80 : 24860 : auto diagfreq = g_inputdeck.get< tag::output, tag::iter, tag::diag >();
81 : :
82 [ + + ]: 24860 : if ( !((d.It()+1) % diagfreq) ) { // if remainder, don't compute diagnostics
83 : :
84 : : // Query number of degrees of freedom from user's setting
85 : 13108 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
86 : :
87 : : // Diagnostics vector (of vectors) during aggregation. See
88 : : // Inciter/Diagnostics.h.
89 : : std::vector< std::vector< tk::real > >
90 [ + - ][ + - ]: 39324 : diag( NUMDIAG, std::vector< tk::real >( u.nprop()/rdof, 0.0 ) );
91 : :
92 : : // Compute diagnostics for DG
93 [ + - ]: 13108 : compute_diag(d, rdof, nchGhost, geoElem, ndofel, u, diag);
94 : :
95 : : // Append diagnostics vector with metadata on the current time step
96 : : // ITER: Current iteration count (only the first entry is used)
97 : : // TIME: Current physical time (only the first entry is used)
98 : : // DT: Current physical time step size (only the first entry is used)
99 : 13108 : diag[ITER][0] = static_cast< tk::real >( d.It()+1 );
100 : 13108 : diag[TIME][0] = d.T() + d.Dt();
101 : 13108 : diag[DT][0] = d.Dt();
102 : :
103 : : // Contribute to diagnostics
104 [ + - ]: 13108 : auto stream = serialize( d.MeshId(), diag );
105 [ + - ]: 13108 : d.contribute( stream.first, stream.second.get(), DiagMerger,
106 [ + - ][ + - ]: 26216 : CkCallback(CkIndex_Transporter::diagnostics(nullptr), d.Tr()) );
107 : :
108 : 13108 : return true; // diagnostics have been computed
109 : :
110 : : }
111 : :
112 : 11752 : return false; // diagnostics have not been computed
113 : : }
114 : :
115 : : void
116 : 13108 : ElemDiagnostics::compute_diag( const Discretization& d,
117 : : const std::size_t rdof,
118 : : const std::size_t nchGhost,
119 : : const tk::Fields& geoElem,
120 : : const std::vector< std::size_t >& ndofel,
121 : : const tk::Fields& u,
122 : : std::vector< std::vector< tk::real > >& diag )
123 : : const
124 : : // *****************************************************************************
125 : : // Compute diagnostics, e.g., residuals, norms of errors, etc. for DG
126 : : //! \param[in] d Discretization base class to read from
127 : : //! \param[in] rdof Number of reconstructed degrees of freedom
128 : : //! \param[in] nchGhost Number of chare boundary ghost elements
129 : : //! \param[in] geoElem Element geometry
130 : : //! \param[in] ndofel Vector of local number of degrees of freedom
131 : : //! \param[in] u Current solution vector
132 : : //! \param[in,out] diag Diagnostics vector
133 : : // *****************************************************************************
134 : : {
135 : 13108 : const auto& inpoel = d.Inpoel();
136 : 13108 : const auto& coord = d.Coord();
137 : :
138 : 13108 : const auto& cx = coord[0];
139 : 13108 : const auto& cy = coord[1];
140 : 13108 : const auto& cz = coord[2];
141 : :
142 [ + + ]: 8930738 : for (std::size_t e=0; e<u.nunk()-nchGhost; ++e)
143 : : {
144 : : // Number of quadrature points for volume integration
145 [ + - ]: 8917630 : auto ng = tk::NGdiag(ndofel[e]);
146 : :
147 : : // arrays for quadrature points
148 : 17835260 : std::array< std::vector< tk::real >, 3 > coordgp;
149 : 17835260 : std::vector< tk::real > wgp;
150 : :
151 [ + - ]: 8917630 : coordgp[0].resize( ng );
152 [ + - ]: 8917630 : coordgp[1].resize( ng );
153 [ + - ]: 8917630 : coordgp[2].resize( ng );
154 [ + - ]: 8917630 : wgp.resize( ng );
155 : :
156 [ + - ]: 8917630 : tk::GaussQuadratureTet( ng, coordgp, wgp );
157 : :
158 : : // Extract the element coordinates
159 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
160 : 26752890 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
161 : 26752890 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
162 : 26752890 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
163 : 80258670 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
164 : :
165 [ + + ]: 20647532 : for (std::size_t igp=0; igp<ng; ++igp)
166 : : {
167 : : // Compute the coordinates of quadrature point at physical domain
168 [ + - ]: 11729902 : auto gp = tk::eval_gp( igp, coordel, coordgp );
169 : :
170 : : // Compute the basis function
171 : 23459804 : auto B = tk::eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp],
172 [ + - ]: 46919608 : coordgp[2][igp]);
173 : :
174 [ + - ]: 11729902 : auto wt = wgp[igp] * geoElem(e, 0, 0);
175 : :
176 : 23459804 : std::vector< tk::real > s;
177 : :
178 [ + + ]: 23459804 : for (const auto& eq : g_dgpde)
179 : : // cppcheck-suppress useStlAlgorithm
180 [ + - ]: 11729902 : s = eq.solution( gp[0], gp[1], gp[2], d.T()+d.Dt() );
181 : :
182 [ + + ]: 45176020 : for (std::size_t c=0; c<u.nprop()/rdof; ++c)
183 : : {
184 : 33446118 : auto mark = c*rdof;
185 [ + - ]: 33446118 : auto ugp = u(e, mark, 0);
186 : :
187 [ + + ]: 33446118 : if(ndofel[e] > 1)
188 : : {
189 [ + - ]: 15829400 : ugp += u(e, mark+1, 0) * B[1]
190 [ + - ]: 15829400 : + u(e, mark+2, 0) * B[2]
191 [ + - ]: 15829400 : + u(e, mark+3, 0) * B[3];
192 : :
193 [ + + ]: 15829400 : if(ndofel[e] > 4)
194 : : {
195 [ + - ]: 5719000 : ugp += u(e, mark+4, 0) * B[4]
196 [ + - ]: 5719000 : + u(e, mark+5, 0) * B[5]
197 [ + - ]: 5719000 : + u(e, mark+6, 0) * B[6]
198 [ + - ]: 5719000 : + u(e, mark+7, 0) * B[7]
199 [ + - ]: 5719000 : + u(e, mark+8, 0) * B[8]
200 [ + - ]: 5719000 : + u(e, mark+9, 0) * B[9];
201 : : }
202 : : }
203 : :
204 : : // Compute sum for L2 norm of the numerical solution
205 : 33446118 : diag[L2SOL][c] += wt * ugp * ugp;
206 : :
207 : : // Compute sum for L2 norm of the numerical-analytic solution
208 : 33446118 : diag[L2ERR][c] += wt * (ugp-s[c]) * (ugp-s[c]);
209 : :
210 : : // Compute max for Linf norm of the numerical-analytic solution
211 : 33446118 : auto err = std::abs( ugp - s[c] );
212 [ + + ]: 33446118 : if (err > diag[LINFERR][c]) diag[LINFERR][c] = err;
213 : :
214 : : // Compute sum of the total energy over the entire domain (only the
215 : : // first entry is used)
216 [ + + ]: 33446118 : if (c == u.nprop()/rdof-1) diag[TOTALSOL][0] += wt * ugp;
217 : : }
218 : : }
219 : : }
220 : 13108 : }
|
