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