Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/Problem/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 Field outputs for multi-material equation solver
9 : : \details This file defines functions for field quantites to be output to
10 : : files for compressible multi-material equations.
11 : : */
12 : : // *****************************************************************************
13 : : #include "FieldOutput.hpp"
14 : : #include "MultiMat/MultiMatIndexing.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Inciter/InputDeck/InputDeck.hpp"
17 : : #include "ConfigureMultiMat.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : : extern ctr::InputDeck g_inputdeck;
22 : :
23 [ + - ]: 394 : std::map< std::string, tk::GetVarFn > MultiMatOutVarFn()
24 : : // *****************************************************************************
25 : : // Return a map that associates user-specified strings to functions
26 : : //! \return Map that associates user-specified strings to functions that compute
27 : : //! relevant quantities to be output to file
28 : : // *****************************************************************************
29 : : {
30 : : std::map< std::string, tk::GetVarFn > OutFnMap;
31 : :
32 : : // Allowed strings for user-def field output vars
33 [ + - ][ + - ]: 394 : OutFnMap["density"] = multimat::bulkDensityOutVar;
34 [ + - ][ + - ]: 394 : OutFnMap["pressure"] = multimat::bulkPressureOutVar;
35 [ + - ][ + - ]: 394 : OutFnMap["specific_total_energy"] = multimat::bulkSpecificTotalEnergyOutVar;
36 [ + - ][ + - ]: 394 : OutFnMap["x-velocity"] = multimat::velocityOutVar<0>;
37 [ + - ][ + - ]: 394 : OutFnMap["y-velocity"] = multimat::velocityOutVar<1>;
38 [ + - ][ + - ]: 394 : OutFnMap["z-velocity"] = multimat::velocityOutVar<2>;
39 [ + - ][ + - ]: 394 : OutFnMap["material_indicator"] = multimat::matIndicatorOutVar;
40 : :
41 : 394 : return OutFnMap;
42 : : }
43 : :
44 : : std::vector< std::string >
45 : 0 : MultiMatFieldNames( std::size_t nmat )
46 : : // *****************************************************************************
47 : : // Return multi-material field names to be output to file
48 : : //! \param[in] nmat Number of materials in system
49 : : //! \return Vector of strings labelling fields output in file
50 : : // *****************************************************************************
51 : : {
52 : : std::vector< std::string > n;
53 : :
54 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
55 [ - - ][ - - ]: 0 : n.push_back( "volfrac"+std::to_string(k+1)+"_numerical" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
56 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
57 [ - - ][ - - ]: 0 : n.push_back( "density"+std::to_string(k+1)+"_numerical" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
58 [ - - ]: 0 : n.push_back( "density_numerical" );
59 [ - - ]: 0 : n.push_back( "x-velocity_numerical" );
60 [ - - ]: 0 : n.push_back( "y-velocity_numerical" );
61 [ - - ]: 0 : n.push_back( "z-velocity_numerical" );
62 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
63 [ - - ][ - - ]: 0 : n.push_back( "pressure"+std::to_string(k+1)+"_numerical" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
64 [ - - ]: 0 : n.push_back( "pressure_numerical" );
65 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
66 [ - - ][ - - ]: 0 : n.push_back( "soundspeed"+std::to_string(k+1) );
[ - - ][ - - ]
67 [ - - ]: 0 : n.push_back( "total_energy_density_numerical" );
68 [ - - ]: 0 : n.push_back( "material_indicator" );
69 [ - - ]: 0 : n.push_back( "timestep" );
70 : :
71 : 0 : return n;
72 : : }
73 : :
74 : : std::vector< std::vector< tk::real > >
75 : 0 : MultiMatFieldOutput(
76 : : ncomp_t,
77 : : std::size_t nmat,
78 : : const std::vector< EOS >& mat_blk,
79 : : std::size_t nunk,
80 : : std::size_t rdof,
81 : : const std::vector< tk::real >&,
82 : : const std::array< std::vector< tk::real >, 3 >&,
83 : : const tk::Fields& U,
84 : : const tk::Fields& P )
85 : : // *****************************************************************************
86 : : // Return field output going to file
87 : : //! \param[in] nmat Number of materials in systen
88 : : //! \param[in] nunk Number of unknowns to extract
89 : : //! \param[in] rdof Number of reconstructed degrees of freedom
90 : : //! \param[in] U Solution vector at recent time step
91 : : //! \param[in] P Vector of primitive quantities at recent time step
92 : : //! \return Vector of vectors to be output to file
93 : : // *****************************************************************************
94 : : {
95 : : // field-output vector with:
96 : : // - nmat volume fractions
97 : : // - nmat material densities
98 : : // - 1 bulk density
99 : : // - 3 velocity components
100 : : // - nmat material pressures
101 : : // - 1 bulk pressure
102 : : // - nmat material sound-speeds
103 : : // - 1 bulk total energy density
104 : : // - 1 material indicator
105 : : // - 1 time-step
106 : : // leading to a size of 4*nmat+8
107 : : std::vector< std::vector< tk::real > >
108 [ - - ][ - - ]: 0 : out( 4*nmat+8, std::vector< tk::real >( nunk ) );
109 : :
110 : : //// mesh node coordinates
111 : : //const auto& x = coord[0];
112 : : //const auto& y = coord[1];
113 : : //const auto& z = coord[2];
114 : :
115 : : // material volume-fractions
116 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
117 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i)
118 : 0 : out[k][i] = U(i, volfracDofIdx(nmat, k, rdof, 0));
119 : : }
120 : :
121 : : // material densities
122 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
123 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
124 [ - - ]: 0 : out[nmat+k][i] = U(i, densityDofIdx(nmat, k, rdof, 0)) /
125 [ - - ]: 0 : std::max(1e-16, U(i, volfracDofIdx(nmat, k, rdof, 0)));
126 : : }
127 : : }
128 : :
129 : : // bulk density
130 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
131 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
132 : 0 : out[2*nmat][i] += U(i, densityDofIdx(nmat, k, rdof, 0));
133 : : }
134 : :
135 : : // velocity components
136 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
137 : 0 : out[2*nmat+1][i] = P(i, velocityDofIdx(nmat, 0, rdof, 0));
138 : 0 : out[2*nmat+2][i] = P(i, velocityDofIdx(nmat, 1, rdof, 0));
139 : 0 : out[2*nmat+3][i] = P(i, velocityDofIdx(nmat, 2, rdof, 0));
140 : : }
141 : :
142 : : // material pressures
143 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
144 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
145 : 0 : out[2*nmat+4+k][i] =
146 : 0 : P(i, pressureDofIdx(nmat, k, rdof, 0)) /
147 [ - - ]: 0 : std::max(1e-16, U(i, volfracDofIdx(nmat, k, rdof, 0)));
148 : : }
149 : : }
150 : :
151 : : // bulk pressure
152 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
153 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
154 : 0 : out[3*nmat+4][i] += P(i, pressureDofIdx(nmat, k, rdof, 0));
155 : : }
156 : :
157 : : // material sound speeds
158 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
159 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
160 : 0 : out[3*nmat+5+k][i] =
161 [ - - ]: 0 : mat_blk[k].compute< EOS::soundspeed >(
162 [ - - ]: 0 : std::max(1e-16, U(i, densityDofIdx(nmat,k,rdof,0))),
163 : : P(i, pressureDofIdx(nmat,k,rdof,0)),
164 : : U(i, volfracDofIdx(nmat,k,rdof,0)), k );
165 : : }
166 : : }
167 : :
168 : : // bulk total energy density
169 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
170 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
171 : 0 : out[4*nmat+5][i] += U(i, energyDofIdx(nmat, k, rdof, 0));
172 : : }
173 : :
174 : : // material indicator
175 [ - - ]: 0 : for (std::size_t i=0; i<nunk; ++i) {
176 : 0 : out[4*nmat+6][i] = 0.0;
177 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
178 : 0 : out[4*nmat+6][i] += U(i, volfracDofIdx(nmat,k,rdof,0))
179 : 0 : * static_cast< tk::real >(k+1);
180 : : }
181 : : }
182 : :
183 : : // time-step
184 : : //for (std::size_t i=0; i<nunk; ++i) {
185 : : // // advection velocity
186 : : // auto u = U(i, velocityDofIdx(nmat,0,rdof,0));
187 : : // auto v = U(i, velocityDofIdx(nmat,1,rdof,0));
188 : : // auto w = U(i, velocityDofIdx(nmat,2,rdof,0));
189 : :
190 : : // auto vn = std::sqrt(tk::dot({{u, v, w}}, {{u, v, w}}));
191 : :
192 : : // // acoustic speed
193 : : // auto a = 0.0;
194 : : // for (std::size_t k=0; k<nmat; ++k)
195 : : // {
196 : : // if (U(i, volfracDofIdx(nmat,k,rdof,0)) > 1.0e-04) {
197 : : // a = std::max( a, eos_soundspeed< tag::multimat >( 0,
198 : : // U(i, densityDofIdx(nmat,k,rdof,0)),
199 : : // P(i, pressureDofIdx(nmat,k,rdof,0)),
200 : : // U(i, volfracDofIdx(nmat,k,rdof,0)), k ) );
201 : : // }
202 : : // }
203 : :
204 : : // out[4*nmat+7][i] = geoElem(i,4) / (std::fabs(vn) + a);
205 : : //}
206 : :
207 : 0 : return out;
208 : : }
209 : :
210 [ + - ]: 26 : std::vector< std::string > MultiMatSurfNames()
211 : : // *****************************************************************************
212 : : // Return surface field names to be output to file
213 : : //! \note Every surface will output these fields.
214 : : //! \return Vector of strings labelling surface fields output in file
215 : : // *****************************************************************************
216 : : {
217 : : std::vector< std::string > n;
218 : :
219 [ + - ]: 26 : n.push_back( "density" );
220 [ + - ]: 26 : n.push_back( "x-velocity" );
221 [ + - ]: 26 : n.push_back( "y-velocity" );
222 [ + - ]: 26 : n.push_back( "z-velocity" );
223 [ + - ]: 26 : n.push_back( "specific_total_energy" );
224 [ + - ]: 26 : n.push_back( "pressure" );
225 : :
226 : 26 : return n;
227 : : }
228 : :
229 : : std::vector< std::vector< tk::real > >
230 : 26 : MultiMatSurfOutput(
231 : : const std::size_t nmat,
232 : : const std::size_t rdof,
233 : : const FaceData& fd,
234 : : const tk::Fields& U,
235 : : const tk::Fields& P )
236 : : // *****************************************************************************
237 : : // Return element surface field output (on triangle faces) going to file
238 : : //! \param[in] nmat Number of materials in this PDE system
239 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
240 : : //! \param[in] fd Face connectivity and boundary conditions object
241 : : //! \param[in] U Solution vector at recent time step
242 : : //! \param[in] P Vector of primitives at recent time step
243 : : //! \return Vector of vectors of solution on side set faces to be output to file
244 : : // *****************************************************************************
245 : : {
246 : : std::vector< std::vector< tk::real > > out;
247 : :
248 : : const auto& bface = fd.Bface();
249 : : const auto& esuf = fd.Esuf();
250 : :
251 : : // extract field output along side sets requested
252 [ + + ]: 50 : for (auto s : g_inputdeck.get< tag::field_output, tag::sideset >()) {
253 : : // get face list for side set requested
254 : 24 : auto b = bface.find(static_cast<int>(s));
255 [ - + ]: 24 : if (b == end(bface)) continue;
256 : : const auto& faces = b->second;
257 [ + - ]: 24 : std::vector< tk::real > surfaceSol( faces.size() );
258 [ + - ]: 24 : auto i = out.size();
259 : : out.insert( end(out), 6, surfaceSol );
260 : : std::size_t j = 0;
261 [ + + ]: 3372 : for (auto f : faces) {
262 : : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
263 : 3348 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
264 : :
265 : : // access solutions at boundary element
266 : : tk::real rhob(0.0), rhoE(0.0), pb(0.0);
267 [ + + ]: 10044 : for (std::size_t k=0; k<nmat; ++k) {
268 : 6696 : rhob += U(el, densityDofIdx(nmat,k,rdof,0));
269 : 6696 : rhoE += U(el, energyDofIdx(nmat,k,rdof,0));
270 : 6696 : pb += P(el, pressureDofIdx(nmat,k,rdof,0));
271 : : }
272 : :
273 : 3348 : out[i+0][j] = rhob;
274 : 3348 : out[i+1][j] = P(el, velocityDofIdx(nmat,0,rdof,0));
275 : 3348 : out[i+2][j] = P(el, velocityDofIdx(nmat,1,rdof,0));
276 : 3348 : out[i+3][j] = P(el, velocityDofIdx(nmat,2,rdof,0));
277 : 3348 : out[i+4][j] = rhoE;
278 : 3348 : out[i+5][j] = pb;
279 : 3348 : ++j;
280 : : }
281 : : }
282 : :
283 : 26 : return out;
284 : : }
285 : :
286 [ - - ]: 0 : std::vector< std::string > MultiMatHistNames()
287 : : // *****************************************************************************
288 : : // Return time history field names to be output to file
289 : : //! \note Every time history point will output these fields.
290 : : //! \return Vector of strings labelling time history fields output in file
291 : : // *****************************************************************************
292 : : {
293 : : std::vector< std::string > n;
294 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
295 : :
296 [ - - ]: 0 : n.push_back( "density" );
297 [ - - ]: 0 : n.push_back( "x-velocity" );
298 [ - - ]: 0 : n.push_back( "y-velocity" );
299 [ - - ]: 0 : n.push_back( "z-velocity" );
300 [ - - ]: 0 : n.push_back( "energy" );
301 [ - - ]: 0 : n.push_back( "pressure" );
302 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
303 [ - - ][ - - ]: 0 : n.push_back( "volfrac"+std::to_string(k+1) );
[ - - ][ - - ]
304 : :
305 : 0 : return n;
306 : : }
307 : :
308 : 35 : std::vector< std::string > MultiMatDiagNames(std::size_t nmat)
309 : : // *****************************************************************************
310 : : // Return diagnostic var names to be output to file
311 : : //! \param[in] nmat Number of materials in systen
312 : : //! \return Vector of strings labelling diagnostic fields output in file
313 : : // *****************************************************************************
314 : : {
315 : : std::vector< std::string > n;
316 : :
317 [ + + ]: 125 : for (std::size_t k=0; k<nmat; ++k)
318 [ + - ][ + - ]: 180 : n.push_back( "f"+std::to_string(k+1) );
[ - + ][ - - ]
319 [ + + ]: 125 : for (std::size_t k=0; k<nmat; ++k)
320 [ + - ][ + - ]: 180 : n.push_back( "fr"+std::to_string(k+1) );
[ - + ][ - - ]
321 [ + - ]: 35 : n.push_back( "fru" );
322 [ + - ]: 35 : n.push_back( "frv" );
323 [ + - ]: 70 : n.push_back( "frw" );
324 [ + + ]: 125 : for (std::size_t k=0; k<nmat; ++k)
325 [ + - ][ + - ]: 180 : n.push_back( "fre"+std::to_string(k+1) );
[ - + ][ - - ]
326 : :
327 : 35 : return n;
328 : : }
329 : :
330 : : } //inciter::
|
