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