Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/FieldOutput.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 Extract field output for inciter
9 : : \details Extract field output for inciter.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "FieldOutput.hpp"
14 : : #include "ContainerUtil.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Integrate/Basis.hpp"
17 : :
18 : : namespace inciter {
19 : :
20 : : std::vector< std::string >
21 : 5807 : numericFieldNames( tk::Centering c )
22 : : // *****************************************************************************
23 : : // Collect field output names from numerical solution based on user input
24 : : //! \param[in] c Extract variable names only with this centering
25 : : //! \return Output field names requested by user
26 : : // *****************************************************************************
27 : : {
28 : 5807 : std::vector< std::string > f;
29 [ + + ]: 32721 : for (const auto& v : g_inputdeck.get< tag::field_output, tag::outvar >()) {
30 [ + + ][ + + ]: 26914 : if (v.centering == c && !v.analytic()) {
[ + + ]
31 [ + - ]: 14037 : std::stringstream s;
32 [ + + ][ + - ]: 14037 : if (v.alias.empty()) s << v.name; else s << v.alias;
[ + - ]
33 [ + - ][ + - ]: 14037 : f.push_back( s.str() );
34 : : }
35 : : }
36 : :
37 : 5807 : return f;
38 : : }
39 : :
40 : : std::vector< std::vector< tk::real > >
41 : 5807 : numericFieldOutput( const tk::Fields& U,
42 : : tk::Centering c,
43 : : const std::map< std::string, tk::GetVarFn >& outvarfn,
44 : : const tk::Fields& P )
45 : : // *****************************************************************************
46 : : // Collect field output from numerical solution based on user input
47 : : //! \param[in] U Solution data to extract from
48 : : //! \param[in] c Extract variables only with this centering
49 : : //! \param[in] outvarfn Map of outvar functions
50 : : //! \param[in] P Optional primitive variable solution data to extract from
51 : : //! \return Output fields requested by user
52 : : // *****************************************************************************
53 : : {
54 : : // will not use P if empty
55 [ + + ]: 5807 : const auto& p = P.empty() ? U : P;
56 : :
57 : : //auto rdof =
58 : : // c == tk::Centering::NODE ? 1 : g_inputdeck.get< tag::rdof >();
59 : 5807 : std::size_t rdof = 1;
60 : :
61 : 5807 : std::vector< std::vector< tk::real > > f;
62 [ + + ]: 32721 : for (const auto& v : g_inputdeck.get< tag::field_output, tag::outvar >()) {
63 [ + + ][ + + ]: 26914 : if (v.centering == c && !v.analytic()) {
[ + + ]
64 : 14037 : auto iPDE = g_inputdeck.get< tag::pde >();
65 [ + + ]: 14037 : const auto& F = v.primitive(iPDE) ? p : U;
66 [ + + ]: 14037 : if (v.varFnIdx == "null") { // depvar-based direct access
67 [ + - ][ + - ]: 2283 : f.push_back( F.extract_comp( v.field*rdof ) );
68 : : } else { // human-readable non-analytic via custom fn
69 [ + - ][ - + ]: 11754 : Assert(outvarfn.find(v.varFnIdx) != outvarfn.end(),
[ - - ][ - - ]
[ - - ]
70 : : "getvar() not configured for " + v.name );
71 [ + - ][ + - ]: 11754 : f.push_back( outvarfn.at(v.varFnIdx)( F, rdof ) );
[ + - ]
72 : : }
73 : : }
74 : : }
75 : :
76 : 5807 : return f;
77 : : }
78 : :
79 : : void
80 : 2202 : evalSolution(
81 : : const Discretization& D,
82 : : const std::vector< std::size_t >& inpoel,
83 : : const tk::UnsMesh::Coords& coord,
84 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
85 : : const std::vector< std::size_t >& ndofel,
86 : : const tk::Fields& U,
87 : : const tk::Fields& P,
88 : : tk::Fields& uElemfields,
89 : : tk::Fields& pElemfields,
90 : : tk::Fields& uNodefields,
91 : : tk::Fields& pNodefields )
92 : : // *****************************************************************************
93 : : // Evaluate solution on incoming (a potentially refined) mesh
94 : : //! \param[in] D Discretization base class to read from
95 : : //! \param[in] inpoel Incoming (potentially refined field-output) mesh
96 : : //! connectivity
97 : : //! \param[in] coord Incoming (potentially refined Field-output) mesh node
98 : : //! coordinates
99 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
100 : : //! \param[in] ndofel Vector of local number of degrees of freedom
101 : : //! \param[in] U Solution vector
102 : : //! \param[in] P Vector of primitives
103 : : //! \param[in,out] uElemfields Solution elem output fields
104 : : //! \param[in,out] pElemfields Primitive elem output fields
105 : : //! \param[in,out] uNodefields Solution nodal output fields
106 : : //! \param[in,out] pNodefields Primitive nodal output fields
107 : : //! \details This function evaluates the solution on the incoming mesh, and
108 : : //! stores it in uElemfields, pElemfields, uNodefields, and pNodefields
109 : : //! appropriately. The incoming mesh can be refined but can also be just the
110 : : //! mesh the numerical solution is computed on.
111 : : //! \note If the incoming mesh is refined (for field output) compared to the
112 : : //! mesh the numerical solution is computed on, the solution is evaluated in
113 : : //! cells as wells as in nodes. If the solution is not refined, the solution
114 : : //! is evaluated in nodes.
115 : : // *****************************************************************************
116 : : {
117 : : using tk::dot;
118 : : using tk::real;
119 : :
120 : 2202 : const auto nelem = inpoel.size()/4;
121 : 2202 : const auto rdof = g_inputdeck.get< tag::rdof >();
122 : 2202 : const auto uncomp = U.nprop() / rdof;
123 : 2202 : const auto pncomp = P.nprop() / rdof;
124 : :
125 : : // If mesh is not refined for field output, cut off ghosts from element
126 : : // solution. (No need to output ghosts and writer would error.) If mesh is
127 : : // refined for field output, resize element solution fields to refined mesh.
128 : 2202 : uElemfields.resize( nelem );
129 : 2202 : pElemfields.resize( nelem );
130 : :
131 : 2202 : auto npoin = coord[0].size();
132 : 2202 : uNodefields.resize( npoin );
133 : 2202 : pNodefields.resize( npoin );
134 : 2202 : uNodefields.fill(0.0);
135 : 2202 : pNodefields.fill(0.0);
136 : :
137 : 2202 : const auto& x = coord[0];
138 : 2202 : const auto& y = coord[1];
139 : 2202 : const auto& z = coord[2];
140 : :
141 : : // Assign values to element-fields
142 [ + + ]: 1607409 : for (std::size_t e=0; e<U.nunk(); ++e) {
143 [ + + ]: 1605207 : if (e < nelem) {
144 [ + + ]: 4354284 : for (std::size_t i=0; i<uncomp; ++i) {
145 : 3127885 : uElemfields(e,i) = U(e,rdof*i);
146 : : }
147 [ + + ]: 1618299 : for (std::size_t i=0; i<pncomp; ++i) {
148 : 391900 : pElemfields(e,i) = P(e,rdof*i);
149 : : }
150 : : }
151 : : }
152 : :
153 : : // If mesh is not refined for output, evaluate solution in nodes
154 [ + + ]: 2202 : if (addedTets.empty()) {
155 : :
156 [ + + ]: 1204985 : for (std::size_t e=0; e<nelem; ++e) {
157 : 1202816 : std::size_t dofe(1);
158 [ + + ]: 1202816 : if (!ndofel.empty()) {
159 : 1196752 : dofe = ndofel[e];
160 : : }
161 : 1202816 : auto e4 = e*4;
162 : : // Extract element node coordinates
163 : : std::array< std::array< real, 3>, 4 > ce{{
164 : 3608448 : {{ x[inpoel[e4 ]], y[inpoel[e4 ]], z[inpoel[e4 ]] }},
165 : 3608448 : {{ x[inpoel[e4+1]], y[inpoel[e4+1]], z[inpoel[e4+1]] }},
166 : 3608448 : {{ x[inpoel[e4+2]], y[inpoel[e4+2]], z[inpoel[e4+2]] }},
167 : 10825344 : {{ x[inpoel[e4+3]], y[inpoel[e4+3]], z[inpoel[e4+3]] }} }};
168 : : // Compute inverse Jacobian
169 : 1202816 : auto J = tk::inverseJacobian( ce[0], ce[1], ce[2], ce[3] );
170 : : // Evaluate solution in child nodes
171 [ + + ]: 6014080 : for (std::size_t j=0; j<4; ++j) {
172 : : std::array< real, 3 >
173 : 4811264 : h{{ce[j][0]-ce[0][0], ce[j][1]-ce[0][1], ce[j][2]-ce[0][2] }};
174 : : auto Bn = tk::eval_basis( dofe,
175 [ + - ]: 9622528 : dot(J[0],h), dot(J[1],h), dot(J[2],h) );
176 [ + - ]: 9622528 : auto u = eval_state( uncomp, rdof, dofe, e, U, Bn );
177 [ + - ]: 9622528 : auto p = eval_state( pncomp, rdof, dofe, e, P, Bn );
178 : : // Assign child node solution
179 [ + + ][ + - ]: 16851144 : for (std::size_t i=0; i<uncomp; ++i) uNodefields(inpoel[e4+j],i) += u[i];
180 [ + + ][ + - ]: 6378864 : for (std::size_t i=0; i<pncomp; ++i) pNodefields(inpoel[e4+j],i) += p[i];
181 : : }
182 : : }
183 : :
184 : : // If mesh is refed for output, evaluate solution in elements and nodes of
185 : : // refined mesh
186 : : } else {
187 : :
188 : 33 : const auto& pinpoel = D.Inpoel(); // unrefined (parent) mesh
189 : :
190 [ + + ]: 157545 : for ([[maybe_unused]] const auto& [child,parent] : addedTets) {
191 [ - + ][ - - ]: 157512 : Assert( child < nelem, "Indexing out of new solution vector" );
[ - - ][ - - ]
192 [ - + ][ - - ]: 157512 : Assert( parent < pinpoel.size()/4,
[ - - ][ - - ]
193 : : "Indexing out of old solution vector" );
194 : : }
195 : :
196 [ + + ]: 157545 : for (const auto& [child,parent] : addedTets) {
197 : 157512 : std::size_t dofe(1);
198 [ + - ]: 157512 : if (!ndofel.empty()) {
199 : 157512 : dofe = ndofel[parent];
200 : : }
201 : : // Extract parent element's node coordinates
202 : 157512 : auto p4 = 4*parent;
203 : : std::array< std::array< real, 3>, 4 > cp{{
204 : 472536 : {{ x[pinpoel[p4 ]], y[pinpoel[p4 ]], z[pinpoel[p4 ]] }},
205 : 472536 : {{ x[pinpoel[p4+1]], y[pinpoel[p4+1]], z[pinpoel[p4+1]] }},
206 : 472536 : {{ x[pinpoel[p4+2]], y[pinpoel[p4+2]], z[pinpoel[p4+2]] }},
207 : 1417608 : {{ x[pinpoel[p4+3]], y[pinpoel[p4+3]], z[pinpoel[p4+3]] }} }};
208 : : // Evaluate inverse Jacobian of the parent
209 : 157512 : auto Jp = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
210 : : // Compute child cell centroid
211 : 157512 : auto c4 = 4*child;
212 : 157512 : auto cx = (x[inpoel[c4 ]] + x[inpoel[c4+1]] +
213 : 157512 : x[inpoel[c4+2]] + x[inpoel[c4+3]]) / 4.0;
214 : 157512 : auto cy = (y[inpoel[c4 ]] + y[inpoel[c4+1]] +
215 : 157512 : y[inpoel[c4+2]] + y[inpoel[c4+3]]) / 4.0;
216 : 157512 : auto cz = (z[inpoel[c4 ]] + z[inpoel[c4+1]] +
217 : 157512 : z[inpoel[c4+2]] + z[inpoel[c4+3]]) / 4.0;
218 : : // Compute solution in child centroid
219 : 157512 : std::array< real, 3 > h{{cx-cp[0][0], cy-cp[0][1], cz-cp[0][2] }};
220 [ + - ]: 315024 : auto B = tk::eval_basis( dofe, dot(Jp[0],h), dot(Jp[1],h), dot(Jp[2],h) );
221 [ + - ]: 315024 : auto u = eval_state( uncomp, rdof, dofe, parent, U, B );
222 [ + - ]: 315024 : auto p = eval_state( pncomp, rdof, dofe, parent, P, B );
223 : : // Assign cell center solution from parent to child
224 [ + + ][ + - ]: 945072 : for (std::size_t i=0; i<uncomp; ++i) uElemfields(child,i) = u[i];
225 [ - + ][ - - ]: 157512 : for (std::size_t i=0; i<pncomp; ++i) pElemfields(child,i) = p[i];
226 : : // Extract child element's node coordinates
227 : : std::array< std::array< real, 3>, 4 > cc{{
228 : 472536 : {{ x[inpoel[c4 ]], y[inpoel[c4 ]], z[inpoel[c4 ]] }},
229 : 472536 : {{ x[inpoel[c4+1]], y[inpoel[c4+1]], z[inpoel[c4+1]] }},
230 : 472536 : {{ x[inpoel[c4+2]], y[inpoel[c4+2]], z[inpoel[c4+2]] }},
231 : 1417608 : {{ x[inpoel[c4+3]], y[inpoel[c4+3]], z[inpoel[c4+3]] }} }};
232 : : // Evaluate solution in child nodes
233 [ + + ]: 787560 : for (std::size_t j=0; j<4; ++j) {
234 : : std::array< real, 3 >
235 : 630048 : hn{{cc[j][0]-cp[0][0], cc[j][1]-cp[0][1], cc[j][2]-cp[0][2] }};
236 : : auto Bn = tk::eval_basis( dofe,
237 [ + - ]: 1260096 : dot(Jp[0],hn), dot(Jp[1],hn), dot(Jp[2],hn) );
238 [ + - ]: 1260096 : auto cnu = eval_state(uncomp, rdof, dofe, parent, U, Bn);
239 [ + - ]: 1260096 : auto cnp = eval_state(pncomp, rdof, dofe, parent, P, Bn);
240 : : // Assign child node solution
241 [ + + ][ + - ]: 3780288 : for (std::size_t i=0; i<uncomp; ++i) uNodefields(inpoel[c4+j],i) += cnu[i];
242 [ - + ][ - - ]: 630048 : for (std::size_t i=0; i<pncomp; ++i) pNodefields(inpoel[c4+j],i) += cnp[i];
243 : : }
244 : : }
245 : : }
246 : 2202 : }
247 : :
248 : : } // inciter::
|
