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