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