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 : : const auto& matprop = g_inputdeck.get< tag::material >();
36 : : 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 : 237 : std::vector< tk::real > rho0mat( nmat, 0.0 );
45 : : const auto& ic = g_inputdeck.get< tag::ic >();
46 : : const auto& icbox = ic.get< tag::box >();
47 : : 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 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
128 : 8073 : std::size_t ncomp = U.nprop()/rdof;
129 : : 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 [ + - ]: 2762480 : std::vector< tk::real > B(rdof, 0.0);
150 : 2762480 : B[0] = 1.0;
151 [ + - ][ + - ]: 5524960 : 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 : : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
160 : 2762480 : U(e, energyDofIdx(nmat, kmax, rdof, 0)), almax, gmax );
161 : :
162 : : 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 : : bool ctm_element(false);
185 [ + + ]: 6670420 : auto alk = U(e, volfracDofIdx(nmat, k, rdof, 0));
186 : 6670420 : auto pk = P(e, pressureDofIdx(nmat, k, rdof, 0)) / alk;
187 : :
188 : : // Determine if element-e needs trace-material cleanup. This decision is
189 : : // based on specific scenarios.
190 : : // WARNING: Changing this decision-making logic can adversely affect
191 : : // code stability.
192 : : if (
193 : : // 1. if volume fraction is lesser than threshold
194 [ + + ][ + - ]: 6670420 : (alk < volfracPRelaxLim()) ||
195 : : // 2. if current and majority material is fluid, AND the material
196 : : // (effective) pressure is negative
197 [ + - ][ + - ]: 5629092 : (solidx[k] == 0 && solidx[kmax] == 0 && matExists(alk) &&
[ + - ][ - + ]
198 [ + - ]: 5629092 : pk < mat_blk[k].compute< EOS::min_eff_pressure >(1e-12,
199 : : U(e, densityDofIdx(nmat, k, rdof, 0)), alk))
200 : : )
201 : : ctm_element = true;
202 : :
203 : : if (ctm_element) {
204 : 3855874 : tk::real prelax(0.0);
205 : 3855874 : std::array< std::array< tk::real, 3 >, 3 > gmat {{}};
206 [ - + ]: 3855874 : if (solidx[k] > 0) {
207 : : // for solids, reset deformation gradient and stress
208 [ - - ]: 0 : resetSolidTensors(nmat, k, e, U, P);
209 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
210 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
211 : 0 : gmat[i][j] = U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0));
212 : : }
213 : : // determine target relaxation pressure
214 [ + - ][ + - ]: 3855874 : prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
215 : : U(e, densityDofIdx(nmat, k, rdof, 0)), alk);
216 : 3855874 : prelax = std::max(prelax, p_target);
217 : :
218 : : // energy change
219 [ + - ]: 3855874 : auto arhomat = U(e, densityDofIdx(nmat, k, rdof, 0));
220 : 3855874 : auto arhoEmat = mat_blk[k].compute< EOS::totalenergy >(arhomat, u, v, w,
221 [ + - ][ - - ]: 3855874 : alk*prelax, alk, gmat);
222 : :
223 : : // total energy flux into majority material
224 : 3855874 : d_arE += (U(e, energyDofIdx(nmat, k, rdof, 0))
225 : 3855874 : - arhoEmat);
226 : :
227 : : // update state of trace material
228 : 3855874 : U(e, volfracDofIdx(nmat, k, rdof, 0)) = alk;
229 : 3855874 : U(e, energyDofIdx(nmat, k, rdof, 0)) = arhoEmat;
230 : 3855874 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk*prelax;
231 : :
232 [ + + ]: 7255330 : for (std::size_t i=1; i<rdof; ++i) {
233 : 3399456 : U(e, energyDofIdx(nmat, k, rdof, i)) = 0.0;
234 : : }
235 : : }
236 : : }
237 : :
238 [ + - ]: 2762480 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)) += d_al;
239 : :
240 : : // 2. Flux energy change into majority material
241 : 2762480 : U(e, energyDofIdx(nmat, kmax, rdof, 0)) += d_arE;
242 : 2762480 : P(e, pressureDofIdx(nmat, kmax, rdof, 0)) =
243 [ + - ]: 2762480 : mat_blk[kmax].compute< EOS::pressure >(
244 : : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
245 : : U(e, energyDofIdx(nmat, kmax, rdof, 0)),
246 : : U(e, volfracDofIdx(nmat, kmax, rdof, 0)), kmax, gmax );
247 : :
248 : : // 3. enforce unit sum of volume fractions
249 : : auto alsum = 0.0;
250 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
251 : 6670420 : alsum += U(e, volfracDofIdx(nmat, k, rdof, 0));
252 : :
253 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k) {
254 [ - + ]: 6670420 : U(e, volfracDofIdx(nmat, k, rdof, 0)) /= alsum;
255 [ - + ]: 6670420 : U(e, densityDofIdx(nmat, k, rdof, 0)) /= alsum;
256 : 6670420 : U(e, energyDofIdx(nmat, k, rdof, 0)) /= alsum;
257 : 6670420 : P(e, pressureDofIdx(nmat, k, rdof, 0)) /= alsum;
258 [ - + ]: 6670420 : if (solidx[k] > 0) {
259 [ - - ]: 0 : for (std::size_t i=0; i<6; ++i)
260 : 0 : P(e, stressDofIdx(nmat, solidx[k], i, rdof, 0)) /= alsum;
261 : : }
262 : : }
263 : :
264 : 2762480 : pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0)) /
265 : 2762480 : U(e, volfracDofIdx(nmat, kmax, rdof, 0));
266 : :
267 : : // check for unphysical state
268 [ + + ]: 9432900 : for (std::size_t k=0; k<nmat; ++k)
269 : : {
270 [ - + ]: 6670420 : auto alpha = U(e, volfracDofIdx(nmat, k, rdof, 0));
271 : 6670420 : auto arho = U(e, densityDofIdx(nmat, k, rdof, 0));
272 : 6670420 : auto apr = P(e, pressureDofIdx(nmat, k, rdof, 0));
273 : :
274 : : // lambda for screen outputs
275 : 0 : auto screenout = [&]()
276 : : {
277 : 0 : std::cout << "Physical time: " << t << std::endl;
278 : 0 : std::cout << "Element centroid: " << geoElem(e,1) << ", "
279 : 0 : << geoElem(e,2) << ", " << geoElem(e,3) << std::endl;
280 : 0 : std::cout << "Material-id: " << k << std::endl;
281 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
282 : 0 : std::cout << "Partial density: " << arho << std::endl;
283 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
284 : 0 : std::cout << "Major pressure: " << pmax << " (mat " << kmax << ")"
285 : : << std::endl;
286 : 0 : std::cout << "Major temperature: " << tmax << " (mat " << kmax << ")"
287 : : << std::endl;
288 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
289 : : << std::endl;
290 : 6670420 : };
291 : :
292 [ - + ]: 6670420 : if (arho < 0.0)
293 : : {
294 : : neg_density = true;
295 [ - - ]: 0 : screenout();
296 : : }
297 : : }
298 : : }
299 : 8073 : return neg_density;
300 : : }
301 : :
302 : : tk::real
303 : 760 : timeStepSizeMultiMat(
304 : : const std::vector< EOS >& mat_blk,
305 : : const std::vector< int >& esuf,
306 : : const tk::Fields& geoFace,
307 : : const tk::Fields& geoElem,
308 : : const std::size_t nelem,
309 : : std::size_t nmat,
310 : : const tk::Fields& U,
311 : : const tk::Fields& P )
312 : : // *****************************************************************************
313 : : // Time step restriction for multi material cell-centered schemes
314 : : //! \param[in] mat_blk EOS material block
315 : : //! \param[in] esuf Elements surrounding elements array
316 : : //! \param[in] geoFace Face geometry array
317 : : //! \param[in] geoElem Element geometry array
318 : : //! \param[in] nelem Number of elements
319 : : //! \param[in] nmat Number of materials in this PDE system
320 : : //! \param[in] U High-order solution vector
321 : : //! \param[in] P High-order vector of primitives
322 : : //! \return Maximum allowable time step based on cfl criterion
323 : : // *****************************************************************************
324 : : {
325 : 760 : const auto ndof = g_inputdeck.get< tag::ndof >();
326 : 760 : const auto rdof = g_inputdeck.get< tag::rdof >();
327 : 760 : std::size_t ncomp = U.nprop()/rdof;
328 : 760 : std::size_t nprim = P.nprop()/rdof;
329 : :
330 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
331 : 760 : std::vector< tk::real > delt(U.nunk(), 0.0);
332 [ + - ][ + - ]: 760 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
[ - - ][ - - ]
333 : :
334 : : // compute maximum characteristic speed at all internal element faces
335 [ + + ]: 674860 : for (std::size_t f=0; f<esuf.size()/2; ++f)
336 : : {
337 [ + - ]: 674100 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
338 [ + - ]: 674100 : auto er = esuf[2*f+1];
339 : :
340 : : std::array< tk::real, 3 > fn;
341 : 674100 : fn[0] = geoFace(f,1);
342 : 674100 : fn[1] = geoFace(f,2);
343 : 674100 : fn[2] = geoFace(f,3);
344 : :
345 : : // left element
346 : :
347 : : // Compute the basis function for the left element
348 [ + - ][ - - ]: 674100 : std::vector< tk::real > B_l(rdof, 0.0);
349 : 674100 : B_l[0] = 1.0;
350 : :
351 : : // get conserved quantities
352 [ + - ]: 674100 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
353 : : // get primitive quantities
354 [ + - ]: 674100 : pgp = eval_state(nprim, rdof, ndof, el, P, B_l);
355 : :
356 : : // advection velocity
357 : 674100 : u = pgp[velocityIdx(nmat, 0)];
358 : 674100 : v = pgp[velocityIdx(nmat, 1)];
359 : 674100 : w = pgp[velocityIdx(nmat, 2)];
360 : :
361 : 674100 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
362 : :
363 : : // acoustic speed
364 : 674100 : a = 0.0;
365 [ + + ]: 2022300 : for (std::size_t k=0; k<nmat; ++k)
366 : : {
367 [ + + ]: 1348200 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
368 [ + - ]: 708796 : auto gk = getDeformGrad(nmat, k, ugp);
369 [ + - ]: 708796 : gk = tk::rotateTensor(gk, fn);
370 [ + - ][ + + ]: 1386765 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
371 : : ugp[densityIdx(nmat, k)],
372 [ + - ]: 708796 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
373 : : }
374 : : }
375 : :
376 : 674100 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
377 : :
378 : : // right element
379 [ + + ]: 674100 : if (er > -1) {
380 : 521540 : std::size_t eR = static_cast< std::size_t >( er );
381 : :
382 : : // Compute the basis function for the right element
383 [ + - ][ - - ]: 521540 : std::vector< tk::real > B_r(rdof, 0.0);
384 : 521540 : B_r[0] = 1.0;
385 : :
386 : : // get conserved quantities
387 [ + - ]: 521540 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
388 : : // get primitive quantities
389 [ + - ]: 521540 : pgp = eval_state( nprim, rdof, ndof, eR, P, B_r);
390 : :
391 : : // advection velocity
392 : 521540 : u = pgp[velocityIdx(nmat, 0)];
393 : 521540 : v = pgp[velocityIdx(nmat, 1)];
394 : 521540 : w = pgp[velocityIdx(nmat, 2)];
395 : :
396 : 521540 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
397 : :
398 : : // acoustic speed
399 : 521540 : a = 0.0;
400 [ + + ]: 1564620 : for (std::size_t k=0; k<nmat; ++k)
401 : : {
402 [ + + ]: 1043080 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
403 [ + - ]: 548883 : auto gk = getDeformGrad(nmat, k, ugp);
404 [ + - ]: 548883 : gk = tk::rotateTensor(gk, fn);
405 [ + - ][ + + ]: 1073618 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
406 : : ugp[densityIdx(nmat, k)],
407 [ + - ]: 548883 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
408 : : }
409 : : }
410 : :
411 [ + + ]: 521540 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
412 : :
413 [ + - ]: 521540 : delt[eR] += std::max( dSV_l, dSV_r );
414 : : } else {
415 : 152560 : dSV_r = dSV_l;
416 : : }
417 : :
418 [ + - ]: 674100 : delt[el] += std::max( dSV_l, dSV_r );
419 : : }
420 : :
421 : 760 : tk::real mindt = std::numeric_limits< tk::real >::max();
422 : :
423 : : // compute allowable dt
424 [ + + ]: 279520 : for (std::size_t e=0; e<nelem; ++e)
425 : : {
426 [ + + ]: 284379 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
427 : : }
428 : :
429 [ + - ]: 1520 : return mindt;
430 : : }
431 : :
432 : : tk::real
433 : 100 : timeStepSizeMultiMatFV(
434 : : const std::vector< EOS >& mat_blk,
435 : : const tk::Fields& geoElem,
436 : : std::size_t nelem,
437 : : std::size_t nmat,
438 : : const tk::Fields& U,
439 : : const tk::Fields& P,
440 : : std::vector< tk::real >& local_dte )
441 : : // *****************************************************************************
442 : : // Time step restriction for multi material cell-centered FV scheme
443 : : //! \param[in] mat_blk Material EOS block
444 : : //! \param[in] geoElem Element geometry array
445 : : //! \param[in] nelem Number of elements
446 : : //! \param[in] nmat Number of materials in this PDE system
447 : : //! \param[in] U High-order solution vector
448 : : //! \param[in] P High-order vector of primitives
449 : : //! \param[in,out] local_dte Time step size for each element (for local
450 : : //! time stepping)
451 : : //! \return Maximum allowable time step based on cfl criterion
452 : : // *****************************************************************************
453 : : {
454 : 100 : const auto ndof = g_inputdeck.get< tag::ndof >();
455 : 100 : const auto rdof = g_inputdeck.get< tag::rdof >();
456 : : const auto use_mass_avg =
457 : 100 : g_inputdeck.get< tag::multimat, tag::dt_sos_massavg >();
458 : 100 : std::size_t ncomp = U.nprop()/rdof;
459 : 100 : std::size_t nprim = P.nprop()/rdof;
460 : :
461 [ + - ][ - - ]: 100 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
462 : 100 : tk::real mindt = std::numeric_limits< tk::real >::max();
463 : :
464 : : // loop over all elements
465 [ + + ]: 106220 : for (std::size_t e=0; e<nelem; ++e)
466 : : {
467 : : // basis function at centroid
468 [ + - ][ - - ]: 106120 : std::vector< tk::real > B(rdof, 0.0);
469 : 106120 : B[0] = 1.0;
470 : :
471 : : // get conserved quantities
472 [ + - ]: 106120 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
473 : : // get primitive quantities
474 [ + - ]: 106120 : pgp = eval_state(nprim, rdof, ndof, e, P, B);
475 : :
476 : : // magnitude of advection velocity
477 : 106120 : auto u = pgp[velocityIdx(nmat, 0)];
478 : 106120 : auto v = pgp[velocityIdx(nmat, 1)];
479 : 106120 : auto w = pgp[velocityIdx(nmat, 2)];
480 : 106120 : auto vmag = std::sqrt(tk::dot({u, v, w}, {u, v, w}));
481 : :
482 : : // acoustic speed
483 : 106120 : tk::real a = 0.0;
484 : : tk::real mixtureDensity = 0.0;
485 [ + + ]: 318360 : for (std::size_t k=0; k<nmat; ++k)
486 : : {
487 [ - + ]: 212240 : if (use_mass_avg > 0)
488 : : {
489 : : // mass averaging SoS
490 [ - - ]: 0 : a += ugp[densityIdx(nmat,k)]*mat_blk[k].compute< EOS::soundspeed >(
491 [ - - ]: 0 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
492 : : ugp[volfracIdx(nmat, k)], k );
493 : :
494 : 0 : mixtureDensity += ugp[densityIdx(nmat,k)];
495 : : }
496 : : else
497 : : {
498 [ + + ]: 212240 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04)
499 : : {
500 [ + - ][ + + ]: 215579 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
501 [ + - ]: 108055 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
502 : : ugp[volfracIdx(nmat, k)], k ) );
503 : : }
504 : : }
505 : : }
506 [ - + ]: 106120 : if (use_mass_avg > 0) a /= mixtureDensity;
507 : :
508 : : // characteristic wave speed
509 [ - + ]: 106120 : auto v_char = vmag + a;
510 : :
511 : : // characteristic length (radius of insphere)
512 [ - + ]: 106120 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
513 : 106120 : /std::sqrt(24.0);
514 : :
515 : : // element dt
516 [ + + ]: 106120 : local_dte[e] = dx/v_char;
517 [ + - ]: 106120 : mindt = std::min(mindt, local_dte[e]);
518 : : }
519 : :
520 [ + - ]: 200 : return mindt;
521 : : }
522 : :
523 : : tk::real
524 : 0 : timeStepSizeViscousFV(
525 : : const tk::Fields& geoElem,
526 : : std::size_t nelem,
527 : : std::size_t nmat,
528 : : const tk::Fields& U )
529 : : // *****************************************************************************
530 : : // Compute the time step size restriction based on this physics
531 : : //! \param[in] geoElem Element geometry array
532 : : //! \param[in] nelem Number of elements
533 : : //! \param[in] nmat Number of materials
534 : : //! \param[in] U Solution vector
535 : : //! \return Maximum allowable time step based on viscosity
536 : : // *****************************************************************************
537 : : {
538 : : const auto& ndof = g_inputdeck.get< tag::ndof >();
539 : : const auto& rdof = g_inputdeck.get< tag::rdof >();
540 : 0 : std::size_t ncomp = U.nprop()/rdof;
541 : :
542 : 0 : auto mindt = std::numeric_limits< tk::real >::max();
543 : :
544 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
545 : : {
546 : : // get conserved quantities at centroid
547 [ - - ]: 0 : std::vector< tk::real > B(rdof, 0.0);
548 : 0 : B[0] = 1.0;
549 [ - - ]: 0 : auto ugp = eval_state(ncomp, rdof, ndof, e, U, B);
550 : :
551 : : // Kinematic viscosity
552 : 0 : tk::real nu(0.0);
553 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
554 : : {
555 [ - - ]: 0 : auto alk = ugp[volfracIdx(nmat, k)];
556 [ - - ]: 0 : auto mu = getmatprop< tag::mu >(k);
557 [ - - ][ - - ]: 0 : if (alk > 1.0e-04) nu = std::max(nu, alk*mu/ugp[densityIdx(nmat,k)]);
558 : : }
559 : :
560 : : // characteristic length (radius of insphere)
561 [ - - ]: 0 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
562 : 0 : /std::sqrt(24.0);
563 : :
564 : : // element dt
565 [ - - ][ - - ]: 0 : mindt = std::min(mindt, dx*dx/nu);
566 : : }
567 : :
568 : 0 : return mindt;
569 : : }
570 : :
571 : : void
572 : 0 : resetSolidTensors(
573 : : std::size_t nmat,
574 : : std::size_t k,
575 : : std::size_t e,
576 : : tk::Fields& U,
577 : : tk::Fields& P )
578 : : // *****************************************************************************
579 : : // Reset the solid tensors
580 : : //! \param[in] nmat Number of materials in this PDE system
581 : : //! \param[in] k Material id whose deformation gradient is required
582 : : //! \param[in] e Id of element whose solution is to be limited
583 : : //! \param[in/out] U High-order solution vector which gets modified
584 : : //! \param[in/out] P High-order vector of primitives which gets modified
585 : : // *****************************************************************************
586 : : {
587 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
588 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
589 : :
590 [ - - ]: 0 : if (solidx[k] == 0) return; // fluid: nothing to reset here
591 : :
592 : : // Load g and symmetrize once to reduce roundoff
593 : 0 : std::array<std::array<tk::real,3>,3> g{};
594 [ - - ]: 0 : for (size_t i=0;i<3;++i)
595 [ - - ]: 0 : for (size_t j=0;j<3;++j)
596 : 0 : g[i][j] = U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0));
597 : : // symmetric part
598 [ - - ]: 0 : for (size_t i=0;i<3;++i)
599 [ - - ]: 0 : for (size_t j=i+1;j<3;++j) {
600 : 0 : tk::real sij = 0.5*(g[i][j] + g[j][i]);
601 : 0 : g[i][j] = g[j][i] = sij;
602 : : }
603 : :
604 : : // Robust determinant and spherical replacement
605 : : // Clamp to avoid NaNs; eps should be small but > 0
606 [ - - ]: 0 : auto detg = std::max(1e-18, tk::determinant(g));
607 : 0 : const tk::real new_gii = std::pow(detg, 1.0/3.0); // = J^{-1/3} > 0
608 : :
609 : : // Set g and zero elastic (deviatoric) Cauchy stress DOFs ONLY (not pressure)
610 [ - - ]: 0 : for (size_t i=0;i<3;++i)
611 [ - - ]: 0 : for (size_t j=0;j<3;++j) {
612 [ - - ]: 0 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = (i==j) ? new_gii : 0.0;
613 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, 0)) = 0.0;
614 : : // Clear higher DOFs for g and elastic stress
615 [ - - ]: 0 : for (size_t l=1; l<rdof; ++l) {
616 : 0 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, l)) = 0.0;
617 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, l)) = 0.0;
618 : : }
619 : : }
620 : : }
621 : :
622 : : std::array< std::array< tk::real, 3 >, 3 >
623 : 19617777 : getDeformGrad(
624 : : std::size_t nmat,
625 : : std::size_t k,
626 : : const std::vector< tk::real >& state )
627 : : // *****************************************************************************
628 : : // Get the inverse deformation gradient tensor for a material at given location
629 : : //! \param[in] nmat Number of materials in this PDE system
630 : : //! \param[in] k Material id whose deformation gradient is required
631 : : //! \param[in] state State vector at location
632 : : //! \return Inverse deformation gradient tensor (g_k) of material k
633 : : // *****************************************************************************
634 : : {
635 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
636 : : std::array< std::array< tk::real, 3 >, 3 > gk;
637 : :
638 [ - + ]: 19617777 : if (solidx[k] > 0) {
639 : : // deformation gradient for solids
640 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
641 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
642 : 0 : gk[i][j] = state[deformIdx(nmat,solidx[k],i,j)];
643 : : }
644 : : }
645 : : else {
646 : : // empty vector for fluids
647 : 19617777 : gk = {{}};
648 : : }
649 : :
650 : 19617777 : return gk;
651 : : }
652 : :
653 : : std::array< std::array< tk::real, 3 >, 3 >
654 : 4160000 : getCauchyStress(
655 : : std::size_t nmat,
656 : : std::size_t k,
657 : : std::size_t ncomp,
658 : : const std::vector< tk::real >& state )
659 : : // *****************************************************************************
660 : : // Get the elastic Cauchy stress tensor for a material at given location
661 : : //! \param[in] nmat Number of materials in this PDE system
662 : : //! \param[in] k Material id whose deformation gradient is required
663 : : //! \param[in] ncomp Number of components in the PDE system
664 : : //! \param[in] state State vector at location
665 : : //! \return Elastic Cauchy stress tensor (alpha * \sigma_ij) of material k
666 : : // *****************************************************************************
667 : : {
668 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
669 : : std::array< std::array< tk::real, 3 >, 3 >
670 : 4160000 : asigk{{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
671 : :
672 : : // elastic Cauchy stress for solids
673 [ - + ]: 4160000 : if (solidx[k] > 0) {
674 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
675 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
676 : 0 : asigk[i][j] = state[ncomp +
677 : 0 : stressIdx(nmat, solidx[k], stressCmp[i][j])];
678 : : }
679 : : }
680 : :
681 : 4160000 : return asigk;
682 : : }
683 : :
684 : : bool
685 : 7380705 : haveSolid(
686 : : std::size_t nmat,
687 : : const std::vector< std::size_t >& solidx )
688 : : // *****************************************************************************
689 : : // Check whether we have solid materials in our problem
690 : : //! \param[in] nmat Number of materials in this PDE system
691 : : //! \param[in] solidx Material index indicator
692 : : //! \return true if we have at least one solid, false otherwise.
693 : : // *****************************************************************************
694 : : {
695 : : bool haveSolid = false;
696 [ + + ]: 22142865 : for (std::size_t k=0; k<nmat; ++k)
697 [ - + ]: 14762160 : if (solidx[k] > 0) haveSolid = true;
698 : :
699 : 7380705 : return haveSolid;
700 : : }
701 : :
702 : 11366424 : std::size_t numSolids(
703 : : std::size_t nmat,
704 : : const std::vector< std::size_t >& solidx )
705 : : // *****************************************************************************
706 : : // Count total number of solid materials in the problem
707 : : //! \param[in] nmat Number of materials in this PDE system
708 : : //! \param[in] solidx Material index indicator
709 : : //! \return Total number of solid materials in the problem
710 : : // *****************************************************************************
711 : : {
712 : : // count number of solid materials
713 : : std::size_t nsld(0);
714 [ + + ]: 37333302 : for (std::size_t k=0; k<nmat; ++k)
715 [ - + ]: 25966878 : if (solidx[k] > 0) ++nsld;
716 : :
717 : 11366424 : return nsld;
718 : : }
719 : :
720 : : } //inciter::
|
