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 : 235 : 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 : 235 : 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 [ + + ]: 856 : for (std::size_t k=0; k<nmat; ++k) {
38 : 621 : auto mateos = matprop[matidxmap.get< tag::eosidx >()[k]].get<tag::eos>();
39 : 621 : mat_blk.emplace_back(mateos, EqType::multimat, k);
40 : : }
41 : :
42 : : // set rho0 for all materials
43 : : // ---------------------------------------------------------------------------
44 : 235 : 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 : 235 : std::size_t k = ic.get< tag::materialid >() - 1;
50 : 235 : tk::real pr = ic.get< tag::pressure >();
51 : 235 : tk::real tmp = ic.get< tag::temperature >();
52 [ + - ][ + + ]: 235 : rho0mat[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
53 : :
54 : : // Check inside used defined box
55 [ + + ]: 235 : if (!icbox.empty())
56 : : {
57 [ + + ]: 168 : for (const auto& b : icbox) { // for all boxes
58 : 84 : k = b.get< tag::materialid >() - 1;
59 : 84 : auto boxmas = b.get< tag::mass >();
60 : : // if mass and volume are given, compute density as mass/volume
61 [ - + ]: 84 : 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 : 84 : pr = b.get< tag::pressure >();
72 : 84 : tmp = b.get< tag::temperature >();
73 [ + - ]: 84 : rho0mat[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
74 : : }
75 : : }
76 : : }
77 : :
78 : : // Check inside user-specified mesh blocks
79 [ - + ]: 235 : 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 [ + + ]: 856 : for (std::size_t i=0; i<nmat; ++i) {
100 [ + - ]: 621 : mat_blk[i].set< EOS::setRho0 >(rho0mat[i]);
101 : : }
102 : 235 : }
103 : :
104 : : bool
105 : 7818 : 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 : 7818 : const auto ndof = g_inputdeck.get< tag::ndof >();
126 : 7818 : const auto rdof = g_inputdeck.get< tag::rdof >();
127 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
128 : 7818 : std::size_t ncomp = U.nprop()/rdof;
129 : : auto neg_density = false;
130 : :
131 : 7818 : std::vector< tk::real > ugp(ncomp, 0.0);
132 : :
133 [ + + ]: 2588378 : for (std::size_t e=0; e<nelem; ++e)
134 : : {
135 : : // find material in largest quantity.
136 : 2580560 : auto almax(0.0);
137 : 2580560 : std::size_t kmax = 0;
138 [ + + ]: 8887140 : for (std::size_t k=0; k<nmat; ++k)
139 : : {
140 : 6306580 : auto al = U(e, volfracDofIdx(nmat, k, rdof, 0));
141 [ + + ]: 6306580 : if (al > almax)
142 : : {
143 : 3989290 : almax = al;
144 : 3989290 : kmax = k;
145 : : }
146 : : }
147 : :
148 : : // get conserved quantities
149 [ + - ]: 2580560 : std::vector< tk::real > B(rdof, 0.0);
150 : 2580560 : B[0] = 1.0;
151 [ + - ][ + - ]: 5161120 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
152 : :
153 [ + - ]: 2580560 : auto u = P(e, velocityDofIdx(nmat, 0, rdof, 0));
154 : 2580560 : auto v = P(e, velocityDofIdx(nmat, 1, rdof, 0));
155 : 2580560 : auto w = P(e, velocityDofIdx(nmat, 2, rdof, 0));
156 [ + - ]: 2580560 : auto pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0))/almax;
157 [ + - ]: 2580560 : auto gmax = getDeformGrad(nmat, kmax, ugp);
158 [ + - ]: 2580560 : auto tmax = mat_blk[kmax].compute< EOS::temperature >(
159 : : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
160 : 2580560 : 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 [ - + ]: 2580560 : p_target = std::max(pmax, 1e-14);
180 : :
181 : : // 1. Correct minority materials and store volume/energy changes
182 [ + + ]: 8887140 : for (std::size_t k=0; k<nmat; ++k)
183 : : {
184 : : bool ctm_element(false);
185 [ + + ]: 6306580 : auto alk = U(e, volfracDofIdx(nmat, k, rdof, 0));
186 : 6306580 : 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 [ + + ][ + - ]: 6306580 : (alk < volfracPRelaxLim()) ||
195 : : // 2. if current and majority material is fluid, AND the material
196 : : // (effective) pressure is negative
197 [ + - ][ + - ]: 5247396 : (solidx[k] == 0 && solidx[kmax] == 0 && matExists(alk) &&
[ + - ][ + + ]
198 [ + - ]: 5247394 : 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 : 3682884 : tk::real prelax(0.0);
205 : 3682884 : std::array< std::array< tk::real, 3 >, 3 > gmat {{}};
206 [ - + ]: 3682884 : 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 : :
214 : : // check for unbounded volume fractions
215 [ + - ][ - + ]: 3682884 : if (alk < 0.0 || !std::isfinite(alk)) {
216 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(p_target,
217 [ - - ]: 0 : std::max(1e-8,tmax));
218 [ - - ]: 0 : if (std::isfinite(alk)) d_al += (alk - 1e-14);
219 : 0 : alk = 1e-14;
220 : 0 : U(e, densityDofIdx(nmat, k, rdof, 0)) = alk * rhok;
221 : : }
222 : :
223 : : // determine target relaxation pressure
224 [ + - ][ + + ]: 3682884 : prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
225 : : U(e, densityDofIdx(nmat, k, rdof, 0)), alk);
226 : 3682884 : prelax = std::max(prelax, p_target);
227 : :
228 : : // energy change
229 [ + - ]: 3682884 : auto arhomat = U(e, densityDofIdx(nmat, k, rdof, 0));
230 : 3682884 : auto arhoEmat = mat_blk[k].compute< EOS::totalenergy >(arhomat, u, v, w,
231 [ + - ][ - - ]: 3682884 : alk*prelax, alk, gmat);
232 : :
233 : : // total energy flux into majority material
234 : 3682884 : d_arE += (U(e, energyDofIdx(nmat, k, rdof, 0))
235 : 3682884 : - arhoEmat);
236 : :
237 : : // update state of trace material
238 : 3682884 : U(e, volfracDofIdx(nmat, k, rdof, 0)) = alk;
239 : 3682884 : U(e, energyDofIdx(nmat, k, rdof, 0)) = arhoEmat;
240 : 3682884 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk*prelax;
241 : :
242 [ + + ]: 6563370 : for (std::size_t i=1; i<rdof; ++i) {
243 : 2880486 : U(e, energyDofIdx(nmat, k, rdof, i)) = 0.0;
244 : : }
245 : : }
246 : : }
247 : :
248 [ + - ]: 2580560 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)) += d_al;
249 : :
250 : : // 2. Flux energy change into majority material
251 : 2580560 : U(e, energyDofIdx(nmat, kmax, rdof, 0)) += d_arE;
252 : 2580560 : P(e, pressureDofIdx(nmat, kmax, rdof, 0)) =
253 [ + - ]: 2580560 : mat_blk[kmax].compute< EOS::pressure >(
254 : : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
255 : : U(e, energyDofIdx(nmat, kmax, rdof, 0)),
256 : : U(e, volfracDofIdx(nmat, kmax, rdof, 0)), kmax, gmax );
257 : :
258 : : // 3. enforce unit sum of volume fractions
259 : : auto alsum = 0.0;
260 [ + + ]: 8887140 : for (std::size_t k=0; k<nmat; ++k)
261 : 6306580 : alsum += U(e, volfracDofIdx(nmat, k, rdof, 0));
262 : :
263 [ + + ]: 8887140 : for (std::size_t k=0; k<nmat; ++k) {
264 [ - + ]: 6306580 : U(e, volfracDofIdx(nmat, k, rdof, 0)) /= alsum;
265 [ - + ]: 6306580 : U(e, densityDofIdx(nmat, k, rdof, 0)) /= alsum;
266 : 6306580 : U(e, energyDofIdx(nmat, k, rdof, 0)) /= alsum;
267 : 6306580 : P(e, pressureDofIdx(nmat, k, rdof, 0)) /= alsum;
268 [ - + ]: 6306580 : if (solidx[k] > 0) {
269 [ - - ]: 0 : for (std::size_t i=0; i<6; ++i)
270 : 0 : P(e, stressDofIdx(nmat, solidx[k], i, rdof, 0)) /= alsum;
271 : : }
272 : : }
273 : :
274 : 2580560 : pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0)) /
275 : 2580560 : U(e, volfracDofIdx(nmat, kmax, rdof, 0));
276 : :
277 : : // check for unphysical state
278 [ + + ]: 8887140 : for (std::size_t k=0; k<nmat; ++k)
279 : : {
280 [ - + ]: 6306580 : auto alpha = U(e, volfracDofIdx(nmat, k, rdof, 0));
281 : 6306580 : auto arho = U(e, densityDofIdx(nmat, k, rdof, 0));
282 : 6306580 : auto apr = P(e, pressureDofIdx(nmat, k, rdof, 0));
283 : :
284 : : // lambda for screen outputs
285 : 0 : auto screenout = [&]()
286 : : {
287 : 0 : std::cout << "Physical time: " << t << std::endl;
288 : 0 : std::cout << "Element centroid: " << geoElem(e,1) << ", "
289 : 0 : << geoElem(e,2) << ", " << geoElem(e,3) << std::endl;
290 : 0 : std::cout << "Material-id: " << k << std::endl;
291 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
292 : 0 : std::cout << "Partial density: " << arho << std::endl;
293 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
294 : 0 : std::cout << "Major pressure: " << pmax << " (mat " << kmax << ")"
295 : : << std::endl;
296 : 0 : std::cout << "Major temperature: " << tmax << " (mat " << kmax << ")"
297 : : << std::endl;
298 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
299 : : << std::endl;
300 : 6306580 : };
301 : :
302 [ - + ]: 6306580 : if (arho < 0.0)
303 : : {
304 : : neg_density = true;
305 [ - - ]: 0 : screenout();
306 : : }
307 : : }
308 : : }
309 : 7818 : return neg_density;
310 : : }
311 : :
312 : : tk::real
313 : 675 : timeStepSizeMultiMat(
314 : : const std::vector< EOS >& mat_blk,
315 : : const std::vector< int >& esuf,
316 : : const tk::Fields& geoFace,
317 : : const tk::Fields& geoElem,
318 : : const std::size_t nelem,
319 : : std::size_t nmat,
320 : : const tk::Fields& U,
321 : : const tk::Fields& P )
322 : : // *****************************************************************************
323 : : // Time step restriction for multi material cell-centered schemes
324 : : //! \param[in] mat_blk EOS material block
325 : : //! \param[in] esuf Elements surrounding elements array
326 : : //! \param[in] geoFace Face geometry array
327 : : //! \param[in] geoElem Element geometry array
328 : : //! \param[in] nelem Number of elements
329 : : //! \param[in] nmat Number of materials in this PDE system
330 : : //! \param[in] U High-order solution vector
331 : : //! \param[in] P High-order vector of primitives
332 : : //! \return Maximum allowable time step based on cfl criterion
333 : : // *****************************************************************************
334 : : {
335 : 675 : const auto ndof = g_inputdeck.get< tag::ndof >();
336 : 675 : const auto rdof = g_inputdeck.get< tag::rdof >();
337 : 675 : std::size_t ncomp = U.nprop()/rdof;
338 : 675 : std::size_t nprim = P.nprop()/rdof;
339 : :
340 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
341 : 675 : std::vector< tk::real > delt(U.nunk(), 0.0);
342 [ + - ][ + - ]: 675 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
[ - - ][ - - ]
343 : :
344 : : // compute maximum characteristic speed at all internal element faces
345 [ + + ]: 535385 : for (std::size_t f=0; f<esuf.size()/2; ++f)
346 : : {
347 [ + - ]: 534710 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
348 [ + - ]: 534710 : auto er = esuf[2*f+1];
349 : :
350 : : std::array< tk::real, 3 > fn;
351 : 534710 : fn[0] = geoFace(f,1);
352 : 534710 : fn[1] = geoFace(f,2);
353 : 534710 : fn[2] = geoFace(f,3);
354 : :
355 : : // left element
356 : :
357 : : // Compute the basis function for the left element
358 [ + - ][ - - ]: 534710 : std::vector< tk::real > B_l(rdof, 0.0);
359 : 534710 : B_l[0] = 1.0;
360 : :
361 : : // get conserved quantities
362 [ + - ]: 534710 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
363 : : // get primitive quantities
364 [ + - ]: 534710 : pgp = eval_state(nprim, rdof, ndof, el, P, B_l);
365 : :
366 : : // advection velocity
367 : 534710 : u = pgp[velocityIdx(nmat, 0)];
368 : 534710 : v = pgp[velocityIdx(nmat, 1)];
369 : 534710 : w = pgp[velocityIdx(nmat, 2)];
370 : :
371 : 534710 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
372 : :
373 : : // acoustic speed
374 : 534710 : a = 0.0;
375 [ + + ]: 1604130 : for (std::size_t k=0; k<nmat; ++k)
376 : : {
377 [ + + ]: 1069420 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
378 [ + - ]: 561186 : auto gk = getDeformGrad(nmat, k, ugp);
379 [ + - ]: 561186 : gk = tk::rotateTensor(gk, fn);
380 [ + - ][ + + ]: 1096504 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
381 : : ugp[densityIdx(nmat, k)],
382 [ + - ]: 561186 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
383 : : }
384 : : }
385 : :
386 : 534710 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
387 : :
388 : : // right element
389 [ + + ]: 534710 : if (er > -1) {
390 : 413990 : std::size_t eR = static_cast< std::size_t >( er );
391 : :
392 : : // Compute the basis function for the right element
393 [ + - ][ - - ]: 413990 : std::vector< tk::real > B_r(rdof, 0.0);
394 : 413990 : B_r[0] = 1.0;
395 : :
396 : : // get conserved quantities
397 [ + - ]: 413990 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
398 : : // get primitive quantities
399 [ + - ]: 413990 : pgp = eval_state( nprim, rdof, ndof, eR, P, B_r);
400 : :
401 : : // advection velocity
402 : 413990 : u = pgp[velocityIdx(nmat, 0)];
403 : 413990 : v = pgp[velocityIdx(nmat, 1)];
404 : 413990 : w = pgp[velocityIdx(nmat, 2)];
405 : :
406 : 413990 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
407 : :
408 : : // acoustic speed
409 : 413990 : a = 0.0;
410 [ + + ]: 1241970 : for (std::size_t k=0; k<nmat; ++k)
411 : : {
412 [ + + ]: 827980 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
413 [ + - ]: 434557 : auto gk = getDeformGrad(nmat, k, ugp);
414 [ + - ]: 434557 : gk = tk::rotateTensor(gk, fn);
415 [ + - ][ + + ]: 849019 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
416 : : ugp[densityIdx(nmat, k)],
417 [ + - ]: 434557 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
418 : : }
419 : : }
420 : :
421 [ + + ]: 413990 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
422 : :
423 [ + - ]: 413990 : delt[eR] += std::max( dSV_l, dSV_r );
424 : : } else {
425 : 120720 : dSV_r = dSV_l;
426 : : }
427 : :
428 [ + - ]: 534710 : delt[el] += std::max( dSV_l, dSV_r );
429 : : }
430 : :
431 : 675 : tk::real mindt = std::numeric_limits< tk::real >::max();
432 : :
433 : : // compute allowable dt
434 [ + + ]: 218795 : for (std::size_t e=0; e<nelem; ++e)
435 : : {
436 [ + + ]: 223222 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
437 : : }
438 : :
439 [ + - ]: 1350 : return mindt;
440 : : }
441 : :
442 : : tk::real
443 : 100 : timeStepSizeMultiMatFV(
444 : : const std::vector< EOS >& mat_blk,
445 : : const tk::Fields& geoElem,
446 : : std::size_t nelem,
447 : : std::size_t nmat,
448 : : const tk::Fields& U,
449 : : const tk::Fields& P,
450 : : std::vector< tk::real >& local_dte )
451 : : // *****************************************************************************
452 : : // Time step restriction for multi material cell-centered FV scheme
453 : : //! \param[in] mat_blk Material EOS block
454 : : //! \param[in] geoElem Element geometry array
455 : : //! \param[in] nelem Number of elements
456 : : //! \param[in] nmat Number of materials in this PDE system
457 : : //! \param[in] U High-order solution vector
458 : : //! \param[in] P High-order vector of primitives
459 : : //! \param[in,out] local_dte Time step size for each element (for local
460 : : //! time stepping)
461 : : //! \return Maximum allowable time step based on cfl criterion
462 : : // *****************************************************************************
463 : : {
464 : 100 : const auto ndof = g_inputdeck.get< tag::ndof >();
465 : 100 : const auto rdof = g_inputdeck.get< tag::rdof >();
466 : : const auto use_mass_avg =
467 : 100 : g_inputdeck.get< tag::multimat, tag::dt_sos_massavg >();
468 : 100 : std::size_t ncomp = U.nprop()/rdof;
469 : 100 : std::size_t nprim = P.nprop()/rdof;
470 : :
471 [ + - ][ - - ]: 100 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
472 : 100 : tk::real mindt = std::numeric_limits< tk::real >::max();
473 : :
474 : : // loop over all elements
475 [ + + ]: 106220 : for (std::size_t e=0; e<nelem; ++e)
476 : : {
477 : : // basis function at centroid
478 [ + - ][ - - ]: 106120 : std::vector< tk::real > B(rdof, 0.0);
479 : 106120 : B[0] = 1.0;
480 : :
481 : : // get conserved quantities
482 [ + - ]: 106120 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
483 : : // get primitive quantities
484 [ + - ]: 106120 : pgp = eval_state(nprim, rdof, ndof, e, P, B);
485 : :
486 : : // magnitude of advection velocity
487 : 106120 : auto u = pgp[velocityIdx(nmat, 0)];
488 : 106120 : auto v = pgp[velocityIdx(nmat, 1)];
489 : 106120 : auto w = pgp[velocityIdx(nmat, 2)];
490 : 106120 : auto vmag = std::sqrt(tk::dot({u, v, w}, {u, v, w}));
491 : :
492 : : // acoustic speed
493 : 106120 : tk::real a = 0.0;
494 : : tk::real mixtureDensity = 0.0;
495 [ + + ]: 318360 : for (std::size_t k=0; k<nmat; ++k)
496 : : {
497 [ - + ]: 212240 : if (use_mass_avg > 0)
498 : : {
499 : : // mass averaging SoS
500 [ - - ]: 0 : a += ugp[densityIdx(nmat,k)]*mat_blk[k].compute< EOS::soundspeed >(
501 [ - - ]: 0 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
502 : : ugp[volfracIdx(nmat, k)], k );
503 : :
504 : 0 : mixtureDensity += ugp[densityIdx(nmat,k)];
505 : : }
506 : : else
507 : : {
508 [ + + ]: 212240 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04)
509 : : {
510 [ + - ][ + + ]: 215579 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
511 [ + - ]: 108055 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
512 : : ugp[volfracIdx(nmat, k)], k ) );
513 : : }
514 : : }
515 : : }
516 [ - + ]: 106120 : if (use_mass_avg > 0) a /= mixtureDensity;
517 : :
518 : : // characteristic wave speed
519 [ - + ]: 106120 : auto v_char = vmag + a;
520 : :
521 : : // characteristic length (radius of insphere)
522 [ - + ]: 106120 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
523 : 106120 : /std::sqrt(24.0);
524 : :
525 : : // element dt
526 [ + + ]: 106120 : local_dte[e] = dx/v_char;
527 [ + - ]: 106120 : mindt = std::min(mindt, local_dte[e]);
528 : : }
529 : :
530 [ + - ]: 200 : return mindt;
531 : : }
532 : :
533 : : tk::real
534 : 0 : timeStepSizeViscousFV(
535 : : const tk::Fields& geoElem,
536 : : std::size_t nelem,
537 : : std::size_t nmat,
538 : : const tk::Fields& U )
539 : : // *****************************************************************************
540 : : // Compute the time step size restriction based on this physics
541 : : //! \param[in] geoElem Element geometry array
542 : : //! \param[in] nelem Number of elements
543 : : //! \param[in] nmat Number of materials
544 : : //! \param[in] U Solution vector
545 : : //! \return Maximum allowable time step based on viscosity
546 : : // *****************************************************************************
547 : : {
548 : : const auto& ndof = g_inputdeck.get< tag::ndof >();
549 : : const auto& rdof = g_inputdeck.get< tag::rdof >();
550 : 0 : std::size_t ncomp = U.nprop()/rdof;
551 : :
552 : 0 : auto mindt = std::numeric_limits< tk::real >::max();
553 : :
554 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
555 : : {
556 : : // get conserved quantities at centroid
557 [ - - ]: 0 : std::vector< tk::real > B(rdof, 0.0);
558 : 0 : B[0] = 1.0;
559 [ - - ]: 0 : auto ugp = eval_state(ncomp, rdof, ndof, e, U, B);
560 : :
561 : : // Kinematic viscosity
562 : 0 : tk::real nu(0.0);
563 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
564 : : {
565 [ - - ]: 0 : auto alk = ugp[volfracIdx(nmat, k)];
566 [ - - ]: 0 : auto mu = getmatprop< tag::mu >(k);
567 [ - - ][ - - ]: 0 : if (alk > 1.0e-04) nu = std::max(nu, alk*mu/ugp[densityIdx(nmat,k)]);
568 : : }
569 : :
570 : : // characteristic length (radius of insphere)
571 [ - - ]: 0 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
572 : 0 : /std::sqrt(24.0);
573 : :
574 : : // element dt
575 [ - - ][ - - ]: 0 : mindt = std::min(mindt, dx*dx/nu);
576 : : }
577 : :
578 : 0 : return mindt;
579 : : }
580 : :
581 : : void
582 : 0 : resetSolidTensors(
583 : : std::size_t nmat,
584 : : std::size_t k,
585 : : std::size_t e,
586 : : tk::Fields& U,
587 : : tk::Fields& P )
588 : : // *****************************************************************************
589 : : // Reset the solid tensors
590 : : //! \param[in] nmat Number of materials in this PDE system
591 : : //! \param[in] k Material id whose deformation gradient is required
592 : : //! \param[in] e Id of element whose solution is to be limited
593 : : //! \param[in/out] U High-order solution vector which gets modified
594 : : //! \param[in/out] P High-order vector of primitives which gets modified
595 : : // *****************************************************************************
596 : : {
597 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
598 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
599 : :
600 [ - - ]: 0 : if (solidx[k] == 0) return; // fluid: nothing to reset here
601 : :
602 : : // Load g and symmetrize once to reduce roundoff
603 : 0 : std::array<std::array<tk::real,3>,3> g{};
604 [ - - ]: 0 : for (size_t i=0;i<3;++i)
605 [ - - ]: 0 : for (size_t j=0;j<3;++j)
606 : 0 : g[i][j] = U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0));
607 : : // symmetric part
608 [ - - ]: 0 : for (size_t i=0;i<3;++i)
609 [ - - ]: 0 : for (size_t j=i+1;j<3;++j) {
610 : 0 : tk::real sij = 0.5*(g[i][j] + g[j][i]);
611 : 0 : g[i][j] = g[j][i] = sij;
612 : : }
613 : :
614 : : // Robust determinant and spherical replacement
615 : : // Clamp to avoid NaNs; eps should be small but > 0
616 [ - - ]: 0 : auto detg = std::max(1e-18, tk::determinant(g));
617 : 0 : const tk::real new_gii = std::pow(detg, 1.0/3.0); // = J^{-1/3} > 0
618 : :
619 : : // Set g and zero elastic (deviatoric) Cauchy stress DOFs ONLY (not pressure)
620 [ - - ]: 0 : for (size_t i=0;i<3;++i)
621 [ - - ]: 0 : for (size_t j=0;j<3;++j) {
622 [ - - ]: 0 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = (i==j) ? new_gii : 0.0;
623 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, 0)) = 0.0;
624 : : // Clear higher DOFs for g and elastic stress
625 [ - - ]: 0 : for (size_t l=1; l<rdof; ++l) {
626 : 0 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, l)) = 0.0;
627 : 0 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, l)) = 0.0;
628 : : }
629 : : }
630 : : }
631 : :
632 : : std::array< std::array< tk::real, 3 >, 3 >
633 : 17794361 : getDeformGrad(
634 : : std::size_t nmat,
635 : : std::size_t k,
636 : : const std::vector< tk::real >& state )
637 : : // *****************************************************************************
638 : : // Get the inverse deformation gradient tensor for a material at given location
639 : : //! \param[in] nmat Number of materials in this PDE system
640 : : //! \param[in] k Material id whose deformation gradient is required
641 : : //! \param[in] state State vector at location
642 : : //! \return Inverse deformation gradient tensor (g_k) of material k
643 : : // *****************************************************************************
644 : : {
645 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
646 : : std::array< std::array< tk::real, 3 >, 3 > gk;
647 : :
648 [ - + ]: 17794361 : if (solidx[k] > 0) {
649 : : // deformation gradient for solids
650 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
651 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
652 : 0 : gk[i][j] = state[deformIdx(nmat,solidx[k],i,j)];
653 : : }
654 : : }
655 : : else {
656 : : // empty vector for fluids
657 : 17794361 : gk = {{}};
658 : : }
659 : :
660 : 17794361 : return gk;
661 : : }
662 : :
663 : : std::array< std::array< tk::real, 3 >, 3 >
664 : 4160000 : getCauchyStress(
665 : : std::size_t nmat,
666 : : std::size_t k,
667 : : std::size_t ncomp,
668 : : const std::vector< tk::real >& state )
669 : : // *****************************************************************************
670 : : // Get the elastic Cauchy stress tensor for a material at given location
671 : : //! \param[in] nmat Number of materials in this PDE system
672 : : //! \param[in] k Material id whose deformation gradient is required
673 : : //! \param[in] ncomp Number of components in the PDE system
674 : : //! \param[in] state State vector at location
675 : : //! \return Elastic Cauchy stress tensor (alpha * \sigma_ij) of material k
676 : : // *****************************************************************************
677 : : {
678 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
679 : : std::array< std::array< tk::real, 3 >, 3 >
680 : 4160000 : asigk{{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
681 : :
682 : : // elastic Cauchy stress for solids
683 [ - + ]: 4160000 : if (solidx[k] > 0) {
684 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
685 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
686 : 0 : asigk[i][j] = state[ncomp +
687 : 0 : stressIdx(nmat, solidx[k], stressCmp[i][j])];
688 : : }
689 : : }
690 : :
691 : 4160000 : return asigk;
692 : : }
693 : :
694 : : bool
695 : 4034550 : haveSolid(
696 : : std::size_t nmat,
697 : : const std::vector< std::size_t >& solidx )
698 : : // *****************************************************************************
699 : : // Check whether we have solid materials in our problem
700 : : //! \param[in] nmat Number of materials in this PDE system
701 : : //! \param[in] solidx Material index indicator
702 : : //! \return true if we have at least one solid, false otherwise.
703 : : // *****************************************************************************
704 : : {
705 : : bool haveSolid = false;
706 [ + + ]: 12104400 : for (std::size_t k=0; k<nmat; ++k)
707 [ - + ]: 8069850 : if (solidx[k] > 0) haveSolid = true;
708 : :
709 : 4034550 : return haveSolid;
710 : : }
711 : :
712 : 9832293 : std::size_t numSolids(
713 : : std::size_t nmat,
714 : : const std::vector< std::size_t >& solidx )
715 : : // *****************************************************************************
716 : : // Count total number of solid materials in the problem
717 : : //! \param[in] nmat Number of materials in this PDE system
718 : : //! \param[in] solidx Material index indicator
719 : : //! \return Total number of solid materials in the problem
720 : : // *****************************************************************************
721 : : {
722 : : // count number of solid materials
723 : : std::size_t nsld(0);
724 [ + + ]: 32730909 : for (std::size_t k=0; k<nmat; ++k)
725 [ - + ]: 22898616 : if (solidx[k] > 0) ++nsld;
726 : :
727 : 9832293 : return nsld;
728 : : }
729 : :
730 : : } //inciter::
|
