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