Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/MiscMultiMatFns.hpp
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 Misc multi-material system functions
9 : : \details This file defines functions that required for multi-material
10 : : compressible hydrodynamics.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <iostream>
15 : :
16 : : #include "MiscMultiMatFns.hpp"
17 : : #include "Inciter/InputDeck/InputDeck.hpp"
18 : : #include "Integrate/Basis.hpp"
19 : : #include "MultiMat/MultiMatIndexing.hpp"
20 : : #include "EoS/GetMatProp.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : :
26 : 237 : void initializeMaterialEoS( std::vector< EOS >& mat_blk )
27 : : // *****************************************************************************
28 : : // Initialize the material block with configured EOS
29 : : //! \param[in,out] mat_blk Material block that gets initialized
30 : : // *****************************************************************************
31 : : {
32 : : // EoS initialization
33 : : // ---------------------------------------------------------------------------
34 : 237 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
35 : 237 : const auto& matprop = g_inputdeck.get< tag::material >();
36 : 237 : const auto& matidxmap = g_inputdeck.get< tag::matidxmap >();
37 [ + + ]: 862 : for (std::size_t k=0; k<nmat; ++k) {
38 : 625 : auto mateos = matprop[matidxmap.get< tag::eosidx >()[k]].get<tag::eos>();
39 [ + - ]: 625 : mat_blk.emplace_back(mateos, EqType::multimat, k);
40 : : }
41 : :
42 : : // set rho0 for all materials
43 : : // ---------------------------------------------------------------------------
44 [ + - ]: 474 : std::vector< tk::real > rho0mat( nmat, 0.0 );
45 : 237 : const auto& ic = g_inputdeck.get< tag::ic >();
46 : 237 : const auto& icbox = ic.get< tag::box >();
47 : 237 : const auto& icmbk = ic.get< tag::meshblock >();
48 : : // Get background properties
49 : 237 : std::size_t k = ic.get< tag::materialid >() - 1;
50 : 237 : tk::real pr = ic.get< tag::pressure >();
51 : 237 : tk::real tmp = ic.get< tag::temperature >();
52 [ + - ]: 237 : rho0mat[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
53 : :
54 : : // Check inside used defined box
55 [ + + ]: 237 : if (!icbox.empty())
56 : : {
57 [ + + ]: 86 : for (const auto& b : icbox) { // for all boxes
58 : 43 : k = b.get< tag::materialid >() - 1;
59 : 43 : auto boxmas = b.get< tag::mass >();
60 : : // if mass and volume are given, compute density as mass/volume
61 [ - + ]: 43 : if (boxmas > 0.0) {
62 : : std::vector< tk::real > box
63 : 0 : { b.get< tag::xmin >(), b.get< tag::xmax >(),
64 : 0 : b.get< tag::ymin >(), b.get< tag::ymax >(),
65 [ - - ]: 0 : b.get< tag::zmin >(), b.get< tag::zmax >() };
66 : 0 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
67 : 0 : rho0mat[k] = boxmas / V_ex;
68 : : }
69 : : // else compute density from eos
70 : : else {
71 : 43 : pr = b.get< tag::pressure >();
72 : 43 : tmp = b.get< tag::temperature >();
73 [ + - ]: 43 : rho0mat[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
74 : : }
75 : : }
76 : : }
77 : :
78 : : // Check inside user-specified mesh blocks
79 [ - + ]: 237 : if (!icmbk.empty())
80 : : {
81 [ - - ]: 0 : for (const auto& b : icmbk) { // for all blocks
82 : 0 : k = b.get< tag::materialid >() - 1;
83 : 0 : auto boxmas = b.get< tag::mass >();
84 : : // if mass and volume are given, compute density as mass/volume
85 [ - - ]: 0 : if (boxmas > 0.0) {
86 : 0 : auto V_ex = b.get< tag::volume >();
87 : 0 : rho0mat[k] = boxmas / V_ex;
88 : : }
89 : : // else compute density from eos
90 : : else {
91 : 0 : pr = b.get< tag::pressure >();
92 : 0 : tmp = b.get< tag::temperature >();
93 [ - - ]: 0 : rho0mat[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
94 : : }
95 : : }
96 : : }
97 : :
98 : : // Store computed rho0's in the EOS-block
99 [ + + ]: 862 : for (std::size_t i=0; i<nmat; ++i) {
100 [ + - ]: 625 : mat_blk[i].set< EOS::setRho0 >(rho0mat[i]);
101 : : }
102 : 237 : }
103 : :
104 : : bool
105 : 8073 : cleanTraceMultiMat(
106 : : tk::real t,
107 : : std::size_t nelem,
108 : : const std::vector< EOS >& mat_blk,
109 : : const tk::Fields& geoElem,
110 : : std::size_t nmat,
111 : : tk::Fields& U,
112 : : tk::Fields& P )
113 : : // *****************************************************************************
114 : : // Clean up the state of trace materials for multi-material PDE system
115 : : //! \param[in] t Physical time
116 : : //! \param[in] nelem Number of elements
117 : : //! \param[in] mat_blk EOS material block
118 : : //! \param[in] geoElem Element geometry array
119 : : //! \param[in] nmat Number of materials in this PDE system
120 : : //! \param[in/out] U High-order solution vector which gets modified
121 : : //! \param[in/out] P High-order vector of primitives which gets modified
122 : : //! \return Boolean indicating if an unphysical material state was found
123 : : // *****************************************************************************
124 : : {
125 : 8073 : const auto ndof = g_inputdeck.get< tag::ndof >();
126 : 8073 : const auto rdof = g_inputdeck.get< tag::rdof >();
127 : 8073 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
128 : 8073 : std::size_t ncomp = U.nprop()/rdof;
129 : 8073 : auto neg_density = false;
130 : :
131 [ + - ]: 8073 : std::vector< tk::real > ugp(ncomp, 0.0);
132 : :
133 [ + + ]: 2770553 : for (std::size_t e=0; e<nelem; ++e)
134 : : {
135 : : // find material in largest quantity.
136 : 2762480 : auto almax(0.0);
137 : 2762480 : std::size_t kmax = 0;
138 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
139 : : {
140 [ + - ]: 6670420 : auto al = U(e, volfracDofIdx(nmat, k, rdof, 0));
141 [ + + ]: 6670420 : if (al > almax)
142 : : {
143 : 4261454 : almax = al;
144 : 4261454 : kmax = k;
145 : : }
146 : : }
147 : :
148 : : // get conserved quantities
149 [ + - ]: 5524960 : std::vector< tk::real > B(rdof, 0.0);
150 : 2762480 : B[0] = 1.0;
151 [ + - ]: 2762480 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
152 : :
153 [ + - ]: 2762480 : auto u = P(e, velocityDofIdx(nmat, 0, rdof, 0));
154 [ + - ]: 2762480 : auto v = P(e, velocityDofIdx(nmat, 1, rdof, 0));
155 [ + - ]: 2762480 : auto w = P(e, velocityDofIdx(nmat, 2, rdof, 0));
156 [ + - ]: 2762480 : auto pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0))/almax;
157 [ + - ]: 2762480 : auto gmax = getDeformGrad(nmat, kmax, ugp);
158 : 2762480 : auto tmax = mat_blk[kmax].compute< EOS::temperature >(
159 [ + - ][ + - ]: 2762480 : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
160 [ + - ]: 2762480 : U(e, energyDofIdx(nmat, kmax, rdof, 0)), almax, gmax );
161 : :
162 : 2762480 : tk::real p_target(0.0), d_al(0.0), d_arE(0.0);
163 : : //// get equilibrium pressure
164 : : //std::vector< tk::real > kmat(nmat, 0.0);
165 : : //tk::real ratio(0.0);
166 : : //for (std::size_t k=0; k<nmat; ++k)
167 : : //{
168 : : // auto arhok = U(e, densityDofIdx(nmat,k,rdof,0));
169 : : // auto alk = U(e, volfracDofIdx(nmat,k,rdof,0));
170 : : // auto apk = P(e, pressureDofIdx(nmat,k,rdof,0)3);
171 : : // auto ak = mat_blk[k].compute< EOS::soundspeed >(arhok, apk, alk, k);
172 : : // kmat[k] = arhok * ak * ak;
173 : :
174 : : // p_target += alk * apk / kmat[k];
175 : : // ratio += alk * alk / kmat[k];
176 : : //}
177 : : //p_target /= ratio;
178 : : //p_target = std::max(p_target, 1e-14);
179 : 2762480 : p_target = std::max(pmax, 1e-14);
180 : :
181 : : // 1. Correct minority materials and store volume/energy changes
182 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
183 : : {
184 [ + - ]: 6670420 : auto alk = U(e, volfracDofIdx(nmat, k, rdof, 0));
185 [ + - ]: 6670420 : auto pk = P(e, pressureDofIdx(nmat, k, rdof, 0)) / alk;
186 : : // for positive volume fractions
187 [ + - ][ + - ]: 6670420 : if (solidx[k] == 0 && solidx[kmax] == 0 && matExists(alk))
[ + + ][ + + ]
188 : : {
189 : : // check if volume fraction is lesser than threshold (volfracPRelaxLim)
190 : : // and if the material (effective) pressure is negative. If either of
191 : : // these conditions is true, perform pressure relaxation.
192 [ + + ]: 5806671 : if ((alk < volfracPRelaxLim()) ||
193 [ - + ][ + + ]: 8621217 : (pk < mat_blk[k].compute< EOS::min_eff_pressure >(1e-12,
194 [ + - ][ + - ]: 2814546 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk))
195 : : /*&& (std::fabs((pk-pmax)/pmax) > 1e-08)*/)
196 : : {
197 : : // determine target relaxation pressure
198 : 177579 : auto prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
199 [ + - ][ + - ]: 177579 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk);
200 : 177579 : prelax = std::max(prelax, p_target);
201 : :
202 : : // energy change
203 [ + - ]: 177579 : auto arhomat = U(e, densityDofIdx(nmat, k, rdof, 0));
204 [ + - ]: 177579 : auto gmat = getDeformGrad(nmat, k, ugp);
205 : 177579 : auto arhoEmat = mat_blk[k].compute< EOS::totalenergy >(arhomat, u, v, w,
206 [ + - ]: 177579 : alk*prelax, alk, gmat);
207 : :
208 : : // total energy flux into majority material
209 [ + - ]: 177579 : d_arE += (U(e, energyDofIdx(nmat, k, rdof, 0))
210 : 177579 : - arhoEmat);
211 : :
212 : : // update state of trace material
213 [ + - ]: 177579 : U(e, volfracDofIdx(nmat, k, rdof, 0)) = alk;
214 [ + - ]: 177579 : U(e, energyDofIdx(nmat, k, rdof, 0)) = arhoEmat;
215 [ + - ]: 177579 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk*prelax;
216 : : }
217 : : }
218 [ + - ]: 3678295 : else if (!matExists(alk)) { // condition so that else-branch not exec'ed for solids
219 : : // determine target relaxation pressure
220 : 3678295 : auto prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
221 [ + - ][ + - ]: 3678295 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk);
222 : 3678295 : prelax = std::max(prelax, p_target);
223 [ + - ]: 3678295 : auto arhok = U(e, densityDofIdx(nmat, k, rdof, 0));
224 : 3678295 : auto gk = std::array< std::array< tk::real, 3 >, 3 >
225 : : {{ {{1, 0, 0}},
226 : : {{0, 1, 0}},
227 : : {{0, 0, 1}} }};
228 : : // update state of trace material
229 [ + - ]: 3678295 : U(e, energyDofIdx(nmat, k, rdof, 0)) =
230 [ + - ]: 3678295 : mat_blk[k].compute< EOS::totalenergy >( arhok, u, v, w, alk*prelax,
231 : : alk, gk );
232 [ + - ]: 3678295 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk *
233 : : prelax;
234 [ + - ]: 3678295 : resetSolidTensors(nmat, k, e, U, P);
235 [ + + ]: 6976558 : for (std::size_t i=1; i<rdof; ++i) {
236 [ + - ]: 3298263 : U(e, energyDofIdx(nmat, k, rdof, i)) = 0.0;
237 [ + - ]: 3298263 : P(e, pressureDofIdx(nmat, k, rdof, i)) = 0.0;
238 : : }
239 : : }
240 : : }
241 : :
242 [ + - ]: 2762480 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)) += d_al;
243 : :
244 : : // 2. Flux energy change into majority material
245 [ + - ]: 2762480 : U(e, energyDofIdx(nmat, kmax, rdof, 0)) += d_arE;
246 [ + - ]: 2762480 : P(e, pressureDofIdx(nmat, kmax, rdof, 0)) =
247 : 5524960 : mat_blk[kmax].compute< EOS::pressure >(
248 [ + - ][ + - ]: 2762480 : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
249 [ + - ]: 2762480 : U(e, energyDofIdx(nmat, kmax, rdof, 0)),
250 [ + - ]: 2762480 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)), kmax, gmax );
251 : :
252 : : // 3. enforce unit sum of volume fractions
253 : 2762480 : auto alsum = 0.0;
254 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
255 [ + - ]: 6670420 : alsum += U(e, volfracDofIdx(nmat, k, rdof, 0));
256 : :
257 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k) {
258 [ + - ]: 6670420 : U(e, volfracDofIdx(nmat, k, rdof, 0)) /= alsum;
259 [ + - ]: 6670420 : U(e, densityDofIdx(nmat, k, rdof, 0)) /= alsum;
260 [ + - ]: 6670420 : U(e, energyDofIdx(nmat, k, rdof, 0)) /= alsum;
261 [ + - ]: 6670420 : P(e, pressureDofIdx(nmat, k, rdof, 0)) /= alsum;
262 [ - + ]: 6670420 : if (solidx[k] > 0) {
263 [ - - ]: 0 : for (std::size_t i=0; i<6; ++i)
264 [ - - ]: 0 : P(e, stressDofIdx(nmat, solidx[k], i, rdof, 0)) /= alsum;
265 : : }
266 : : }
267 : :
268 [ + - ]: 2762480 : pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0)) /
269 [ + - ]: 2762480 : U(e, volfracDofIdx(nmat, kmax, rdof, 0));
270 : :
271 : : // check for unphysical state
272 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
273 : : {
274 [ + - ]: 6670420 : auto alpha = U(e, volfracDofIdx(nmat, k, rdof, 0));
275 [ + - ]: 6670420 : auto arho = U(e, densityDofIdx(nmat, k, rdof, 0));
276 [ + - ]: 6670420 : auto apr = P(e, pressureDofIdx(nmat, k, rdof, 0));
277 : :
278 : : // lambda for screen outputs
279 : 0 : auto screenout = [&]()
280 : : {
281 : 0 : std::cout << "Physical time: " << t << std::endl;
282 : 0 : std::cout << "Element centroid: " << geoElem(e,1) << ", "
283 : 0 : << geoElem(e,2) << ", " << geoElem(e,3) << std::endl;
284 : 0 : std::cout << "Material-id: " << k << std::endl;
285 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
286 : 0 : std::cout << "Partial density: " << arho << std::endl;
287 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
288 : 0 : std::cout << "Major pressure: " << pmax << " (mat " << kmax << ")"
289 : 0 : << std::endl;
290 : 0 : std::cout << "Major temperature: " << tmax << " (mat " << kmax << ")"
291 : 0 : << std::endl;
292 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
293 : 0 : << std::endl;
294 : 6670420 : };
295 : :
296 [ - + ]: 6670420 : if (arho < 0.0)
297 : : {
298 : 0 : neg_density = true;
299 [ - - ]: 0 : screenout();
300 : : }
301 : : }
302 : : }
303 : 16146 : return neg_density;
304 : : }
305 : :
306 : : tk::real
307 : 760 : timeStepSizeMultiMat(
308 : : const std::vector< EOS >& mat_blk,
309 : : const std::vector< int >& esuf,
310 : : const tk::Fields& geoFace,
311 : : const tk::Fields& geoElem,
312 : : const std::size_t nelem,
313 : : std::size_t nmat,
314 : : const tk::Fields& U,
315 : : const tk::Fields& P )
316 : : // *****************************************************************************
317 : : // Time step restriction for multi material cell-centered schemes
318 : : //! \param[in] mat_blk EOS material block
319 : : //! \param[in] esuf Elements surrounding elements array
320 : : //! \param[in] geoFace Face geometry array
321 : : //! \param[in] geoElem Element geometry array
322 : : //! \param[in] nelem Number of elements
323 : : //! \param[in] nmat Number of materials in this PDE system
324 : : //! \param[in] U High-order solution vector
325 : : //! \param[in] P High-order vector of primitives
326 : : //! \return Maximum allowable time step based on cfl criterion
327 : : // *****************************************************************************
328 : : {
329 : 760 : const auto ndof = g_inputdeck.get< tag::ndof >();
330 : 760 : const auto rdof = g_inputdeck.get< tag::rdof >();
331 : 760 : std::size_t ncomp = U.nprop()/rdof;
332 : 760 : std::size_t nprim = P.nprop()/rdof;
333 : :
334 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
335 [ + - ]: 1520 : std::vector< tk::real > delt(U.nunk(), 0.0);
336 [ + - ][ + - ]: 1520 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
337 : :
338 : : // compute maximum characteristic speed at all internal element faces
339 [ + + ]: 674860 : for (std::size_t f=0; f<esuf.size()/2; ++f)
340 : : {
341 : 674100 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
342 : 674100 : auto er = esuf[2*f+1];
343 : :
344 : : std::array< tk::real, 3 > fn;
345 [ + - ]: 674100 : fn[0] = geoFace(f,1);
346 [ + - ]: 674100 : fn[1] = geoFace(f,2);
347 [ + - ]: 674100 : fn[2] = geoFace(f,3);
348 : :
349 : : // left element
350 : :
351 : : // Compute the basis function for the left element
352 [ + - ]: 674100 : std::vector< tk::real > B_l(rdof, 0.0);
353 : 674100 : B_l[0] = 1.0;
354 : :
355 : : // get conserved quantities
356 [ + - ]: 674100 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
357 : : // get primitive quantities
358 [ + - ]: 674100 : pgp = eval_state(nprim, rdof, ndof, el, P, B_l);
359 : :
360 : : // advection velocity
361 : 674100 : u = pgp[velocityIdx(nmat, 0)];
362 : 674100 : v = pgp[velocityIdx(nmat, 1)];
363 : 674100 : w = pgp[velocityIdx(nmat, 2)];
364 : :
365 [ + - ][ + - ]: 674100 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ + - ]
366 : :
367 : : // acoustic speed
368 : 674100 : a = 0.0;
369 [ + + ]: 2022300 : for (std::size_t k=0; k<nmat; ++k)
370 : : {
371 [ + + ]: 1348200 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
372 [ + - ]: 708787 : auto gk = getDeformGrad(nmat, k, ugp);
373 [ + - ]: 708787 : gk = tk::rotateTensor(gk, fn);
374 [ + - ]: 2126361 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
375 : 708787 : ugp[densityIdx(nmat, k)],
376 : 1417574 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
377 : : }
378 : : }
379 : :
380 [ + - ]: 674100 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
381 : :
382 : : // right element
383 [ + + ]: 674100 : if (er > -1) {
384 : 521540 : std::size_t eR = static_cast< std::size_t >( er );
385 : :
386 : : // Compute the basis function for the right element
387 [ + - ]: 521540 : std::vector< tk::real > B_r(rdof, 0.0);
388 : 521540 : B_r[0] = 1.0;
389 : :
390 : : // get conserved quantities
391 [ + - ]: 521540 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
392 : : // get primitive quantities
393 [ + - ]: 521540 : pgp = eval_state( nprim, rdof, ndof, eR, P, B_r);
394 : :
395 : : // advection velocity
396 : 521540 : u = pgp[velocityIdx(nmat, 0)];
397 : 521540 : v = pgp[velocityIdx(nmat, 1)];
398 : 521540 : w = pgp[velocityIdx(nmat, 2)];
399 : :
400 [ + - ][ + - ]: 521540 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ + - ]
401 : :
402 : : // acoustic speed
403 : 521540 : a = 0.0;
404 [ + + ]: 1564620 : for (std::size_t k=0; k<nmat; ++k)
405 : : {
406 [ + + ]: 1043080 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
407 [ + - ]: 548892 : auto gk = getDeformGrad(nmat, k, ugp);
408 [ + - ]: 548892 : gk = tk::rotateTensor(gk, fn);
409 [ + - ]: 1646676 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
410 : 548892 : ugp[densityIdx(nmat, k)],
411 : 1097784 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
412 : : }
413 : : }
414 : :
415 [ + - ]: 521540 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
416 : :
417 : 521540 : delt[eR] += std::max( dSV_l, dSV_r );
418 : : } else {
419 : 152560 : dSV_r = dSV_l;
420 : : }
421 : :
422 : 674100 : delt[el] += std::max( dSV_l, dSV_r );
423 : : }
424 : :
425 : 760 : tk::real mindt = std::numeric_limits< tk::real >::max();
426 : :
427 : : // compute allowable dt
428 [ + + ]: 279520 : for (std::size_t e=0; e<nelem; ++e)
429 : : {
430 [ + - ]: 278760 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
431 : : }
432 : :
433 : 1520 : return mindt;
434 : : }
435 : :
436 : : tk::real
437 : 100 : timeStepSizeMultiMatFV(
438 : : const std::vector< EOS >& mat_blk,
439 : : const tk::Fields& geoElem,
440 : : std::size_t nelem,
441 : : std::size_t nmat,
442 : : const tk::Fields& U,
443 : : const tk::Fields& P,
444 : : std::vector< tk::real >& local_dte )
445 : : // *****************************************************************************
446 : : // Time step restriction for multi material cell-centered FV scheme
447 : : //! \param[in] mat_blk Material EOS block
448 : : //! \param[in] geoElem Element geometry array
449 : : //! \param[in] nelem Number of elements
450 : : //! \param[in] nmat Number of materials in this PDE system
451 : : //! \param[in] U High-order solution vector
452 : : //! \param[in] P High-order vector of primitives
453 : : //! \param[in,out] local_dte Time step size for each element (for local
454 : : //! time stepping)
455 : : //! \return Maximum allowable time step based on cfl criterion
456 : : // *****************************************************************************
457 : : {
458 : 100 : const auto ndof = g_inputdeck.get< tag::ndof >();
459 : 100 : const auto rdof = g_inputdeck.get< tag::rdof >();
460 : : const auto use_mass_avg =
461 : 100 : g_inputdeck.get< tag::multimat, tag::dt_sos_massavg >();
462 : 100 : std::size_t ncomp = U.nprop()/rdof;
463 : 100 : std::size_t nprim = P.nprop()/rdof;
464 : :
465 [ + - ][ + - ]: 200 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
466 : 100 : tk::real mindt = std::numeric_limits< tk::real >::max();
467 : :
468 : : // loop over all elements
469 [ + + ]: 106220 : for (std::size_t e=0; e<nelem; ++e)
470 : : {
471 : : // basis function at centroid
472 [ + - ]: 106120 : std::vector< tk::real > B(rdof, 0.0);
473 : 106120 : B[0] = 1.0;
474 : :
475 : : // get conserved quantities
476 [ + - ]: 106120 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
477 : : // get primitive quantities
478 [ + - ]: 106120 : pgp = eval_state(nprim, rdof, ndof, e, P, B);
479 : :
480 : : // magnitude of advection velocity
481 : 106120 : auto u = pgp[velocityIdx(nmat, 0)];
482 : 106120 : auto v = pgp[velocityIdx(nmat, 1)];
483 : 106120 : auto w = pgp[velocityIdx(nmat, 2)];
484 : 106120 : auto vmag = std::sqrt(tk::dot({u, v, w}, {u, v, w}));
485 : :
486 : : // acoustic speed
487 : 106120 : tk::real a = 0.0;
488 : 106120 : tk::real mixtureDensity = 0.0;
489 [ + + ]: 318360 : for (std::size_t k=0; k<nmat; ++k)
490 : : {
491 [ - + ]: 212240 : if (use_mass_avg > 0)
492 : : {
493 : : // mass averaging SoS
494 [ - - ]: 0 : a += ugp[densityIdx(nmat,k)]*mat_blk[k].compute< EOS::soundspeed >(
495 : 0 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
496 : 0 : ugp[volfracIdx(nmat, k)], k );
497 : :
498 : 0 : mixtureDensity += ugp[densityIdx(nmat,k)];
499 : : }
500 : : else
501 : : {
502 [ + + ]: 212240 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04)
503 : : {
504 [ + - ]: 324165 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
505 : 108055 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
506 : 216110 : ugp[volfracIdx(nmat, k)], k ) );
507 : : }
508 : : }
509 : : }
510 [ - + ]: 106120 : if (use_mass_avg > 0) a /= mixtureDensity;
511 : :
512 : : // characteristic wave speed
513 : 106120 : auto v_char = vmag + a;
514 : :
515 : : // characteristic length (radius of insphere)
516 [ + - ][ + - ]: 106120 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
517 : 106120 : /std::sqrt(24.0);
518 : :
519 : : // element dt
520 : 106120 : local_dte[e] = dx/v_char;
521 : 106120 : mindt = std::min(mindt, local_dte[e]);
522 : : }
523 : :
524 : 200 : return mindt;
525 : : }
526 : :
527 : : tk::real
528 : 0 : timeStepSizeViscousFV(
529 : : const tk::Fields& geoElem,
530 : : std::size_t nelem,
531 : : std::size_t nmat,
532 : : const tk::Fields& U )
533 : : // *****************************************************************************
534 : : // Compute the time step size restriction based on this physics
535 : : //! \param[in] geoElem Element geometry array
536 : : //! \param[in] nelem Number of elements
537 : : //! \param[in] nmat Number of materials
538 : : //! \param[in] U Solution vector
539 : : //! \return Maximum allowable time step based on viscosity
540 : : // *****************************************************************************
541 : : {
542 : 0 : const auto& ndof = g_inputdeck.get< tag::ndof >();
543 : 0 : const auto& rdof = g_inputdeck.get< tag::rdof >();
544 : 0 : std::size_t ncomp = U.nprop()/rdof;
545 : :
546 : 0 : auto mindt = std::numeric_limits< tk::real >::max();
547 : :
548 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
549 : : {
550 : : // get conserved quantities at centroid
551 [ - - ]: 0 : std::vector< tk::real > B(rdof, 0.0);
552 : 0 : B[0] = 1.0;
553 [ - - ]: 0 : auto ugp = eval_state(ncomp, rdof, ndof, e, U, B);
554 : :
555 : : // Kinematic viscosity
556 : 0 : tk::real nu(0.0);
557 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
558 : : {
559 : 0 : auto alk = ugp[volfracIdx(nmat, k)];
560 [ - - ]: 0 : auto mu = getmatprop< tag::mu >(k);
561 [ - - ]: 0 : if (alk > 1.0e-04) nu = std::max(nu, alk*mu/ugp[densityIdx(nmat,k)]);
562 : : }
563 : :
564 : : // characteristic length (radius of insphere)
565 [ - - ][ - - ]: 0 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
566 : 0 : /std::sqrt(24.0);
567 : :
568 : : // element dt
569 : 0 : mindt = std::min(mindt, dx*dx/nu);
570 : : }
571 : :
572 : 0 : return mindt;
573 : : }
574 : :
575 : : void
576 : 3678295 : resetSolidTensors(
577 : : std::size_t nmat,
578 : : std::size_t k,
579 : : std::size_t e,
580 : : tk::Fields& U,
581 : : tk::Fields& P )
582 : : // *****************************************************************************
583 : : // Reset the solid tensors
584 : : //! \param[in] nmat Number of materials in this PDE system
585 : : //! \param[in] k Material id whose deformation gradient is required
586 : : //! \param[in] e Id of element whose solution is to be limited
587 : : //! \param[in/out] U High-order solution vector which gets modified
588 : : //! \param[in/out] P High-order vector of primitives which gets modified
589 : : // *****************************************************************************
590 : : {
591 : 3678295 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
592 : 3678295 : const auto rdof = g_inputdeck.get< tag::rdof >();
593 : :
594 [ - + ]: 3678295 : if (solidx[k] > 0) {
595 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
596 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
597 : : // deformation gradient reset
598 [ - - ]: 0 : if (i==j) U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = 1.0;
599 : 0 : else U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = 0.0;
600 : :
601 : : // elastic Cauchy-stress reset
602 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, 0)) = 0.0;
603 : :
604 : : // high-order reset
605 [ - - ]: 0 : for (std::size_t l=1; l<rdof; ++l) {
606 : 0 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, l)) = 0.0;
607 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, l)) = 0.0;
608 : : }
609 : : }
610 : : }
611 : : }
612 : 3678295 : }
613 : :
614 : : std::array< std::array< tk::real, 3 >, 3 >
615 : 19795356 : getDeformGrad(
616 : : std::size_t nmat,
617 : : std::size_t k,
618 : : const std::vector< tk::real >& state )
619 : : // *****************************************************************************
620 : : // Get the inverse deformation gradient tensor for a material at given location
621 : : //! \param[in] nmat Number of materials in this PDE system
622 : : //! \param[in] k Material id whose deformation gradient is required
623 : : //! \param[in] state State vector at location
624 : : //! \return Inverse deformation gradient tensor (g_k) of material k
625 : : // *****************************************************************************
626 : : {
627 : 19795356 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
628 : : std::array< std::array< tk::real, 3 >, 3 > gk;
629 : :
630 [ - + ]: 19795356 : if (solidx[k] > 0) {
631 : : // deformation gradient for solids
632 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
633 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
634 : 0 : gk[i][j] = state[deformIdx(nmat,solidx[k],i,j)];
635 : : }
636 : : }
637 : : else {
638 : : // empty vector for fluids
639 : 19795356 : gk = {{}};
640 : : }
641 : :
642 : 19795356 : return gk;
643 : : }
644 : :
645 : : std::array< std::array< tk::real, 3 >, 3 >
646 : 4160000 : getCauchyStress(
647 : : std::size_t nmat,
648 : : std::size_t k,
649 : : std::size_t ncomp,
650 : : const std::vector< tk::real >& state )
651 : : // *****************************************************************************
652 : : // Get the elastic Cauchy stress tensor for a material at given location
653 : : //! \param[in] nmat Number of materials in this PDE system
654 : : //! \param[in] k Material id whose deformation gradient is required
655 : : //! \param[in] ncomp Number of components in the PDE system
656 : : //! \param[in] state State vector at location
657 : : //! \return Elastic Cauchy stress tensor (alpha * \sigma_ij) of material k
658 : : // *****************************************************************************
659 : : {
660 : 4160000 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
661 : : std::array< std::array< tk::real, 3 >, 3 >
662 : 4160000 : asigk{{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
663 : :
664 : : // elastic Cauchy stress for solids
665 [ - + ]: 4160000 : if (solidx[k] > 0) {
666 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
667 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
668 : 0 : asigk[i][j] = state[ncomp +
669 : 0 : stressIdx(nmat, solidx[k], stressCmp[i][j])];
670 : : }
671 : : }
672 : :
673 : 4160000 : return asigk;
674 : : }
675 : :
676 : : bool
677 : 7380705 : haveSolid(
678 : : std::size_t nmat,
679 : : const std::vector< std::size_t >& solidx )
680 : : // *****************************************************************************
681 : : // Check whether we have solid materials in our problem
682 : : //! \param[in] nmat Number of materials in this PDE system
683 : : //! \param[in] solidx Material index indicator
684 : : //! \return true if we have at least one solid, false otherwise.
685 : : // *****************************************************************************
686 : : {
687 : 7380705 : bool haveSolid = false;
688 [ + + ]: 22142865 : for (std::size_t k=0; k<nmat; ++k)
689 [ - + ]: 14762160 : if (solidx[k] > 0) haveSolid = true;
690 : :
691 : 7380705 : return haveSolid;
692 : : }
693 : :
694 : 3079899 : std::size_t numSolids(
695 : : std::size_t nmat,
696 : : const std::vector< std::size_t >& solidx )
697 : : // *****************************************************************************
698 : : // Count total number of solid materials in the problem
699 : : //! \param[in] nmat Number of materials in this PDE system
700 : : //! \param[in] solidx Material index indicator
701 : : //! \return Total number of solid materials in the problem
702 : : // *****************************************************************************
703 : : {
704 : : // count number of solid materials
705 : 3079899 : std::size_t nsld(0);
706 [ + + ]: 9911877 : for (std::size_t k=0; k<nmat; ++k)
707 [ - + ]: 6831978 : if (solidx[k] > 0) ++nsld;
708 : :
709 : 3079899 : return nsld;
710 : : }
711 : :
712 : : } //inciter::
|
