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 : :
21 : : namespace inciter {
22 : :
23 : : extern ctr::InputDeck g_inputdeck;
24 : :
25 : 209 : void initializeMaterialEoS( std::vector< EOS >& mat_blk )
26 : : // *****************************************************************************
27 : : // Initialize the material block with configured EOS
28 : : //! \param[in,out] mat_blk Material block that gets initialized
29 : : // *****************************************************************************
30 : : {
31 : : // EoS initialization
32 : 209 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
33 : 209 : const auto& matprop = g_inputdeck.get< tag::material >();
34 : 209 : const auto& matidxmap = g_inputdeck.get< tag::matidxmap >();
35 [ + + ]: 778 : for (std::size_t k=0; k<nmat; ++k) {
36 : 569 : auto mateos = matprop[matidxmap.get< tag::eosidx >()[k]].get<tag::eos>();
37 [ + - ]: 569 : mat_blk.emplace_back(mateos, EqType::multimat, k);
38 : : }
39 : 209 : }
40 : :
41 : : bool
42 : 7048 : cleanTraceMultiMat(
43 : : tk::real t,
44 : : std::size_t nelem,
45 : : const std::vector< EOS >& mat_blk,
46 : : const tk::Fields& geoElem,
47 : : std::size_t nmat,
48 : : tk::Fields& U,
49 : : tk::Fields& P )
50 : : // *****************************************************************************
51 : : // Clean up the state of trace materials for multi-material PDE system
52 : : //! \param[in] t Physical time
53 : : //! \param[in] nelem Number of elements
54 : : //! \param[in] mat_blk EOS material block
55 : : //! \param[in] geoElem Element geometry array
56 : : //! \param[in] nmat Number of materials in this PDE system
57 : : //! \param[in/out] U High-order solution vector which gets modified
58 : : //! \param[in/out] P High-order vector of primitives which gets modified
59 : : //! \return Boolean indicating if an unphysical material state was found
60 : : // *****************************************************************************
61 : : {
62 : 7048 : const auto ndof = g_inputdeck.get< tag::ndof >();
63 : 7048 : const auto rdof = g_inputdeck.get< tag::rdof >();
64 : 7048 : std::size_t ncomp = U.nprop()/rdof;
65 : 7048 : auto al_eps = 1.0e-02;
66 : 7048 : auto neg_density = false;
67 : :
68 [ + - ]: 7048 : std::vector< tk::real > ugp(ncomp, 0.0);
69 : :
70 [ + + ]: 2575548 : for (std::size_t e=0; e<nelem; ++e)
71 : : {
72 : : // find material in largest quantity, and determine if pressure
73 : : // relaxation is needed. If it is, determine materials that need
74 : : // relaxation, and the total volume of these materials.
75 [ + - ]: 5137000 : std::vector< int > relaxInd(nmat, 0);
76 : 2568500 : auto almax(0.0), relaxVol(0.0);
77 : 2568500 : std::size_t kmax = 0;
78 [ + + ]: 8850960 : for (std::size_t k=0; k<nmat; ++k)
79 : : {
80 [ + - ]: 6282460 : auto al = U(e, volfracDofIdx(nmat, k, rdof, 0));
81 [ + + ]: 6282460 : if (al > almax)
82 : : {
83 : 3937154 : almax = al;
84 : 3937154 : kmax = k;
85 : : }
86 [ + + ]: 2345306 : else if (al < al_eps)
87 : : {
88 : 2316852 : relaxInd[k] = 1;
89 : 2316852 : relaxVol += al;
90 : : }
91 : : }
92 : 2568500 : relaxInd[kmax] = 1;
93 : 2568500 : relaxVol += almax;
94 : :
95 : : // get conserved quantities
96 [ + - ]: 5137000 : std::vector< tk::real > B(rdof, 0.0);
97 : 2568500 : B[0] = 1.0;
98 [ + - ]: 2568500 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
99 : :
100 [ + - ]: 2568500 : auto u = P(e, velocityDofIdx(nmat, 0, rdof, 0));
101 [ + - ]: 2568500 : auto v = P(e, velocityDofIdx(nmat, 1, rdof, 0));
102 [ + - ]: 2568500 : auto w = P(e, velocityDofIdx(nmat, 2, rdof, 0));
103 [ + - ]: 2568500 : auto pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0))/almax;
104 [ + - ]: 2568500 : auto gmax = getDeformGrad(nmat, kmax, ugp);
105 : 2568500 : auto tmax = mat_blk[kmax].compute< EOS::temperature >(
106 [ + - ][ + - ]: 2568500 : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
107 [ + - ]: 2568500 : U(e, energyDofIdx(nmat, kmax, rdof, 0)), almax, gmax );
108 : :
109 : 2568500 : tk::real p_target(0.0), d_al(0.0), d_arE(0.0);
110 : : //// get equilibrium pressure
111 : : //std::vector< tk::real > kmat(nmat, 0.0);
112 : : //tk::real ratio(0.0);
113 : : //for (std::size_t k=0; k<nmat; ++k)
114 : : //{
115 : : // auto arhok = U(e, densityDofIdx(nmat,k,rdof,0));
116 : : // auto alk = U(e, volfracDofIdx(nmat,k,rdof,0));
117 : : // auto apk = P(e, pressureDofIdx(nmat,k,rdof,0)3);
118 : : // auto ak = mat_blk[k].compute< EOS::soundspeed >(arhok, apk, alk, k);
119 : : // kmat[k] = arhok * ak * ak;
120 : :
121 : : // p_target += alk * apk / kmat[k];
122 : : // ratio += alk * alk / kmat[k];
123 : : //}
124 : : //p_target /= ratio;
125 : : //p_target = std::max(p_target, 1e-14);
126 : 2568500 : p_target = std::max(pmax, 1e-14);
127 : :
128 : : // 1. Correct minority materials and store volume/energy changes
129 [ + + ]: 8850960 : for (std::size_t k=0; k<nmat; ++k)
130 : : {
131 [ + - ]: 6282460 : auto alk = U(e, volfracDofIdx(nmat, k, rdof, 0));
132 [ + - ]: 6282460 : auto pk = P(e, pressureDofIdx(nmat, k, rdof, 0)) / alk;
133 : : // for positive volume fractions
134 [ + + ]: 6282460 : if (matExists(alk))
135 : : {
136 : : // check if volume fraction is lesser than threshold (al_eps) and
137 : : // if the material (effective) pressure is negative. If either of
138 : : // these conditions is true, perform pressure relaxation.
139 [ + + ]: 5459093 : if ((alk < al_eps) ||
140 [ + + ][ + + ]: 8081857 : (pk < mat_blk[k].compute< EOS::min_eff_pressure >(1e-12,
141 [ + - ][ + - ]: 2622764 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk))
142 : : /*&& (std::fabs((pk-pmax)/pmax) > 1e-08)*/)
143 : : {
144 : : //auto gk = gamma< tag::multimat >(0, k);
145 : :
146 : 213567 : tk::real alk_new(0.0);
147 : : //// volume change based on polytropic expansion/isentropic compression
148 : : //if (pk > p_target)
149 : : //{
150 : : // alk_new = std::pow((pk/p_target), (1.0/gk)) * alk;
151 : : //}
152 : : //else
153 : : //{
154 : : // auto arhok = U(e, densityDofIdx(nmat, k, rdof, 0));
155 : : // auto ck = eos_soundspeed< tag::multimat >(0, arhok, alk*pk,
156 : : // alk, k);
157 : : // auto kk = arhok * ck * ck;
158 : : // alk_new = alk - (alk*alk/kk) * (p_target-pk);
159 : : //}
160 : 213567 : alk_new = alk;
161 : :
162 : : // determine target relaxation pressure
163 : 213567 : auto prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
164 [ + - ][ + - ]: 213567 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk_new);
165 : 213567 : prelax = std::max(prelax, p_target);
166 : :
167 : : // energy change
168 [ + - ]: 213567 : auto rhomat = U(e, densityDofIdx(nmat, k, rdof, 0))
169 : 213567 : / alk_new;
170 : 213567 : std::array< std::array< tk::real, 3 >, 3 > gmat {{
171 : : {{1, 0, 0}},
172 : : {{0, 1, 0}},
173 : : {{0, 0, 1}} }};
174 [ + - ]: 213567 : auto rhoEmat = mat_blk[k].compute< EOS::totalenergy >(rhomat, u, v, w,
175 : : prelax, gmat);
176 : :
177 : : // volume-fraction and total energy flux into majority material
178 : 213567 : d_al += (alk - alk_new);
179 [ + - ]: 213567 : d_arE += (U(e, energyDofIdx(nmat, k, rdof, 0))
180 : 213567 : - alk_new * rhoEmat);
181 : :
182 : : // update state of trace material
183 [ + - ]: 213567 : U(e, volfracDofIdx(nmat, k, rdof, 0)) = alk_new;
184 [ + - ]: 213567 : U(e, energyDofIdx(nmat, k, rdof, 0)) = alk_new*rhoEmat;
185 [ + - ]: 213567 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk_new*prelax;
186 [ + - ]: 213567 : resetSolidTensors(nmat, k, e, alk_new*prelax, U, P);
187 : : }
188 : : }
189 : : // check for unbounded volume fractions
190 [ + + ][ - + ]: 3446131 : else if (alk < 0.0 || !std::isfinite(alk))
[ + + ]
191 : : {
192 [ + - ]: 20 : auto rhok = mat_blk[k].compute< EOS::density >(p_target,
193 : 10 : std::max(1e-8,tmax));
194 [ + - ]: 10 : if (std::isfinite(alk)) d_al += (alk - 1e-14);
195 : : // update state of trace material
196 [ + - ]: 10 : U(e, volfracDofIdx(nmat, k, rdof, 0)) = 1e-14;
197 [ + - ]: 10 : U(e, densityDofIdx(nmat, k, rdof, 0)) = 1e-14 * rhok;
198 [ + - ]: 10 : U(e, energyDofIdx(nmat, k, rdof, 0)) = 1e-14
199 [ + - ]: 10 : * mat_blk[k].compute< EOS::totalenergy >(rhok, u, v, w, p_target,
200 : : gmax);
201 [ + - ]: 10 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = 1e-14 *
202 : : p_target;
203 [ + - ]: 10 : resetSolidTensors(nmat, k, e, 1e-14*p_target, U, P);
204 [ + + ]: 40 : for (std::size_t i=1; i<rdof; ++i) {
205 [ + - ]: 30 : U(e, volfracDofIdx(nmat, k, rdof, i)) = 0.0;
206 [ + - ]: 30 : U(e, densityDofIdx(nmat, k, rdof, i)) = 0.0;
207 [ + - ]: 30 : U(e, energyDofIdx(nmat, k, rdof, i)) = 0.0;
208 [ + - ]: 30 : P(e, pressureDofIdx(nmat, k, rdof, i)) = 0.0;
209 : : }
210 : : }
211 : : else {
212 : : // determine target relaxation pressure
213 : 3446121 : auto prelax = mat_blk[k].compute< EOS::min_eff_pressure >(1e-10,
214 [ + - ][ + - ]: 3446121 : U(e, densityDofIdx(nmat, k, rdof, 0)), alk);
215 : 3446121 : prelax = std::max(prelax, p_target);
216 [ + - ]: 3446121 : auto rhok = U(e, densityDofIdx(nmat, k, rdof, 0)) / alk;
217 [ + - ]: 3446121 : auto gk = getDeformGrad(nmat, k, ugp);
218 : : // update state of trace material
219 [ + - ]: 3446121 : U(e, energyDofIdx(nmat, k, rdof, 0)) = alk
220 [ + - ]: 3446121 : * mat_blk[k].compute< EOS::totalenergy >( rhok, u, v, w, prelax, gk );
221 [ + - ]: 3446121 : P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk *
222 : : prelax;
223 [ + - ]: 3446121 : resetSolidTensors(nmat, k, e, alk*prelax, U, P);
224 [ + + ]: 6963918 : for (std::size_t i=1; i<rdof; ++i) {
225 [ + - ]: 3517797 : U(e, energyDofIdx(nmat, k, rdof, i)) = 0.0;
226 [ + - ]: 3517797 : P(e, pressureDofIdx(nmat, k, rdof, i)) = 0.0;
227 : : }
228 : : }
229 : : }
230 : :
231 : : // 2. Based on volume change in majority material, compute energy change
232 : : //auto gmax = gamma< tag::multimat >(0, kmax);
233 : : //auto pmax_new = pmax * std::pow(almax/(almax+d_al), gmax);
234 : : //auto rhomax_new = U(e, densityDofIdx(nmat, kmax, rdof, 0))
235 : : // / (almax+d_al);
236 : : //auto rhoEmax_new = eos_totalenergy< tag::multimat >(0, rhomax_new, u,
237 : : // v, w, pmax_new, kmax);
238 : : //auto d_arEmax_new = (almax+d_al) * rhoEmax_new
239 : : // - U(e, energyDofIdx(nmat, kmax, rdof, 0));
240 : :
241 [ + - ]: 2568500 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)) += d_al;
242 : : //U(e, energyDofIdx(nmat, kmax, rdof, 0)) += d_arEmax_new;
243 : :
244 : : // 2. Flux energy change into majority material
245 [ + - ]: 2568500 : U(e, energyDofIdx(nmat, kmax, rdof, 0)) += d_arE;
246 [ + - ]: 2568500 : P(e, pressureDofIdx(nmat, kmax, rdof, 0)) =
247 : 5137000 : mat_blk[kmax].compute< EOS::pressure >(
248 [ + - ][ + - ]: 2568500 : U(e, densityDofIdx(nmat, kmax, rdof, 0)), u, v, w,
249 [ + - ]: 2568500 : U(e, energyDofIdx(nmat, kmax, rdof, 0)),
250 [ + - ]: 2568500 : U(e, volfracDofIdx(nmat, kmax, rdof, 0)), kmax, gmax );
251 : :
252 : : // enforce unit sum of volume fractions
253 : 2568500 : auto alsum = 0.0;
254 [ + + ]: 8850960 : for (std::size_t k=0; k<nmat; ++k)
255 [ + - ]: 6282460 : alsum += U(e, volfracDofIdx(nmat, k, rdof, 0));
256 : :
257 [ + + ]: 8850960 : for (std::size_t k=0; k<nmat; ++k) {
258 [ + - ]: 6282460 : U(e, volfracDofIdx(nmat, k, rdof, 0)) /= alsum;
259 [ + - ]: 6282460 : U(e, densityDofIdx(nmat, k, rdof, 0)) /= alsum;
260 [ + - ]: 6282460 : U(e, energyDofIdx(nmat, k, rdof, 0)) /= alsum;
261 [ + - ]: 6282460 : P(e, pressureDofIdx(nmat, k, rdof, 0)) /= alsum;
262 : : }
263 : :
264 : : //// bulk quantities
265 : : //auto rhoEb(0.0), rhob(0.0), volb(0.0);
266 : : //for (std::size_t k=0; k<nmat; ++k)
267 : : //{
268 : : // if (relaxInd[k] > 0.0)
269 : : // {
270 : : // rhoEb += U(e, energyDofIdx(nmat,k,rdof,0));
271 : : // volb += U(e, volfracDofIdx(nmat,k,rdof,0));
272 : : // rhob += U(e, densityDofIdx(nmat,k,rdof,0));
273 : : // }
274 : : //}
275 : :
276 : : //// 2. find mixture-pressure
277 : : //tk::real pmix(0.0), den(0.0);
278 : : //pmix = rhoEb - 0.5*rhob*(u*u+v*v+w*w);
279 : : //for (std::size_t k=0; k<nmat; ++k)
280 : : //{
281 : : // auto gk = gamma< tag::multimat >(0, k);
282 : : // auto Pck = pstiff< tag::multimat >(0, k);
283 : :
284 : : // pmix -= U(e, volfracDofIdx(nmat,k,rdof,0)) * gk * Pck *
285 : : // relaxInd[k] / (gk-1.0);
286 : : // den += U(e, volfracDofIdx(nmat,k,rdof,0)) * relaxInd[k]
287 : : // / (gk-1.0);
288 : : //}
289 : : //pmix /= den;
290 : :
291 : : //// 3. correct energies
292 : : //for (std::size_t k=0; k<nmat; ++k)
293 : : //{
294 : : // if (relaxInd[k] > 0.0)
295 : : // {
296 : : // auto alk_new = U(e, volfracDofIdx(nmat,k,rdof,0));
297 : : // U(e, energyDofIdx(nmat,k,rdof,0)) = alk_new *
298 : : // eos_totalenergy< tag::multimat >(0, rhomat[k], u, v, w, pmix,
299 : : // k);
300 : : // P(e, pressureDofIdx(nmat, k, rdof, 0)) = alk_new * pmix;
301 : : // }
302 : : //}
303 : :
304 [ + - ]: 2568500 : pmax = P(e, pressureDofIdx(nmat, kmax, rdof, 0)) /
305 [ + - ]: 2568500 : U(e, volfracDofIdx(nmat, kmax, rdof, 0));
306 : :
307 : : // check for unphysical state
308 [ + + ]: 8850960 : for (std::size_t k=0; k<nmat; ++k)
309 : : {
310 [ + - ]: 6282460 : auto alpha = U(e, volfracDofIdx(nmat, k, rdof, 0));
311 [ + - ]: 6282460 : auto arho = U(e, densityDofIdx(nmat, k, rdof, 0));
312 [ + - ]: 6282460 : auto apr = P(e, pressureDofIdx(nmat, k, rdof, 0));
313 : :
314 : : // lambda for screen outputs
315 : 0 : auto screenout = [&]()
316 : : {
317 : 0 : std::cout << "Physical time: " << t << std::endl;
318 : 0 : std::cout << "Element centroid: " << geoElem(e,1) << ", "
319 : 0 : << geoElem(e,2) << ", " << geoElem(e,3) << std::endl;
320 : 0 : std::cout << "Material-id: " << k << std::endl;
321 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
322 : 0 : std::cout << "Partial density: " << arho << std::endl;
323 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
324 : 0 : std::cout << "Major pressure: " << pmax << " (mat " << kmax << ")"
325 : 0 : << std::endl;
326 : 0 : std::cout << "Major temperature: " << tmax << " (mat " << kmax << ")"
327 : 0 : << std::endl;
328 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
329 : 0 : << std::endl;
330 : 6282460 : };
331 : :
332 [ - + ]: 6282460 : if (arho < 0.0)
333 : : {
334 : 0 : neg_density = true;
335 [ - - ]: 0 : screenout();
336 : : }
337 : : }
338 : : }
339 : 14096 : return neg_density;
340 : : }
341 : :
342 : : tk::real
343 : 335 : timeStepSizeMultiMat(
344 : : const std::vector< EOS >& mat_blk,
345 : : const std::vector< int >& esuf,
346 : : const tk::Fields& geoFace,
347 : : const tk::Fields& geoElem,
348 : : const std::size_t nelem,
349 : : std::size_t nmat,
350 : : const tk::Fields& U,
351 : : const tk::Fields& P )
352 : : // *****************************************************************************
353 : : // Time step restriction for multi material cell-centered schemes
354 : : //! \param[in] mat_blk EOS material block
355 : : //! \param[in] esuf Elements surrounding elements array
356 : : //! \param[in] geoFace Face geometry array
357 : : //! \param[in] geoElem Element geometry array
358 : : //! \param[in] nelem Number of elements
359 : : //! \param[in] nmat Number of materials in this PDE system
360 : : //! \param[in] U High-order solution vector
361 : : //! \param[in] P High-order vector of primitives
362 : : //! \return Maximum allowable time step based on cfl criterion
363 : : // *****************************************************************************
364 : : {
365 : 335 : const auto ndof = g_inputdeck.get< tag::ndof >();
366 : 335 : const auto rdof = g_inputdeck.get< tag::rdof >();
367 : 335 : std::size_t ncomp = U.nprop()/rdof;
368 : 335 : std::size_t nprim = P.nprop()/rdof;
369 : :
370 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
371 [ + - ]: 670 : std::vector< tk::real > delt(U.nunk(), 0.0);
372 [ + - ][ + - ]: 670 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
373 : :
374 : : // compute maximum characteristic speed at all internal element faces
375 [ + + ]: 526965 : for (std::size_t f=0; f<esuf.size()/2; ++f)
376 : : {
377 : 526630 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
378 : 526630 : auto er = esuf[2*f+1];
379 : :
380 : : std::array< tk::real, 3 > fn;
381 [ + - ]: 526630 : fn[0] = geoFace(f,1);
382 [ + - ]: 526630 : fn[1] = geoFace(f,2);
383 [ + - ]: 526630 : fn[2] = geoFace(f,3);
384 : :
385 : : // left element
386 : :
387 : : // Compute the basis function for the left element
388 [ + - ]: 526630 : std::vector< tk::real > B_l(rdof, 0.0);
389 : 526630 : B_l[0] = 1.0;
390 : :
391 : : // get conserved quantities
392 [ + - ]: 526630 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
393 : : // get primitive quantities
394 [ + - ]: 526630 : pgp = eval_state(nprim, rdof, ndof, el, P, B_l);
395 : :
396 : : // advection velocity
397 : 526630 : u = pgp[velocityIdx(nmat, 0)];
398 : 526630 : v = pgp[velocityIdx(nmat, 1)];
399 : 526630 : w = pgp[velocityIdx(nmat, 2)];
400 : :
401 [ + - ][ + - ]: 526630 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ + - ]
402 : :
403 : : // acoustic speed
404 : 526630 : a = 0.0;
405 [ + + ]: 1579890 : for (std::size_t k=0; k<nmat; ++k)
406 : : {
407 [ + + ]: 1053260 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
408 [ + - ]: 546463 : auto gk = getDeformGrad(nmat, k, ugp);
409 [ + - ]: 546463 : gk = tk::rotateTensor(gk, fn);
410 [ + - ]: 1639389 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
411 : 546463 : ugp[densityIdx(nmat, k)],
412 : 1092926 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
413 : : }
414 : : }
415 : :
416 [ + - ]: 526630 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
417 : :
418 : : // right element
419 [ + + ]: 526630 : if (er > -1) {
420 : 406010 : std::size_t eR = static_cast< std::size_t >( er );
421 : :
422 : : // Compute the basis function for the right element
423 [ + - ]: 406010 : std::vector< tk::real > B_r(rdof, 0.0);
424 : 406010 : B_r[0] = 1.0;
425 : :
426 : : // get conserved quantities
427 [ + - ]: 406010 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
428 : : // get primitive quantities
429 [ + - ]: 406010 : pgp = eval_state( nprim, rdof, ndof, eR, P, B_r);
430 : :
431 : : // advection velocity
432 : 406010 : u = pgp[velocityIdx(nmat, 0)];
433 : 406010 : v = pgp[velocityIdx(nmat, 1)];
434 : 406010 : w = pgp[velocityIdx(nmat, 2)];
435 : :
436 [ + - ][ + - ]: 406010 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ + - ]
437 : :
438 : : // acoustic speed
439 : 406010 : a = 0.0;
440 [ + + ]: 1218030 : for (std::size_t k=0; k<nmat; ++k)
441 : : {
442 [ + + ]: 812020 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
443 [ + - ]: 421988 : auto gk = getDeformGrad(nmat, k, ugp);
444 [ + - ]: 421988 : gk = tk::rotateTensor(gk, fn);
445 [ + - ]: 1265964 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
446 : 421988 : ugp[densityIdx(nmat, k)],
447 : 843976 : pgp[pressureIdx(nmat, k)], ugp[volfracIdx(nmat, k)], k, gk ) );
448 : : }
449 : : }
450 : :
451 [ + - ]: 406010 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
452 : :
453 : 406010 : delt[eR] += std::max( dSV_l, dSV_r );
454 : : } else {
455 : 120620 : dSV_r = dSV_l;
456 : : }
457 : :
458 : 526630 : delt[el] += std::max( dSV_l, dSV_r );
459 : : }
460 : :
461 : 335 : tk::real mindt = std::numeric_limits< tk::real >::max();
462 : :
463 : : // compute allowable dt
464 [ + + ]: 229115 : for (std::size_t e=0; e<nelem; ++e)
465 : : {
466 [ + - ]: 228780 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
467 : : }
468 : :
469 : 670 : return mindt;
470 : : }
471 : :
472 : : tk::real
473 : 225 : timeStepSizeMultiMatFV(
474 : : const std::vector< EOS >& mat_blk,
475 : : const tk::Fields& geoElem,
476 : : std::size_t nelem,
477 : : std::size_t nmat,
478 : : const tk::Fields& U,
479 : : const tk::Fields& P,
480 : : std::vector< tk::real >& local_dte )
481 : : // *****************************************************************************
482 : : // Time step restriction for multi material cell-centered FV scheme
483 : : //! \param[in] mat_blk Material EOS block
484 : : //! \param[in] geoElem Element geometry array
485 : : //! \param[in] nelem Number of elements
486 : : //! \param[in] nmat Number of materials in this PDE system
487 : : //! \param[in] U High-order solution vector
488 : : //! \param[in] P High-order vector of primitives
489 : : //! \param[in,out] local_dte Time step size for each element (for local
490 : : //! time stepping)
491 : : //! \return Maximum allowable time step based on cfl criterion
492 : : // *****************************************************************************
493 : : {
494 : 225 : const auto ndof = g_inputdeck.get< tag::ndof >();
495 : 225 : const auto rdof = g_inputdeck.get< tag::rdof >();
496 : 225 : std::size_t ncomp = U.nprop()/rdof;
497 : 225 : std::size_t nprim = P.nprop()/rdof;
498 : :
499 [ + - ][ + - ]: 450 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
500 : 225 : tk::real mindt = std::numeric_limits< tk::real >::max();
501 : :
502 : : // loop over all elements
503 [ + + ]: 84325 : for (std::size_t e=0; e<nelem; ++e)
504 : : {
505 : : // basis function at centroid
506 [ + - ]: 84100 : std::vector< tk::real > B(rdof, 0.0);
507 : 84100 : B[0] = 1.0;
508 : :
509 : : // get conserved quantities
510 [ + - ]: 84100 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
511 : : // get primitive quantities
512 [ + - ]: 84100 : pgp = eval_state(nprim, rdof, ndof, e, P, B);
513 : :
514 : : // magnitude of advection velocity
515 : 84100 : auto u = pgp[velocityIdx(nmat, 0)];
516 : 84100 : auto v = pgp[velocityIdx(nmat, 1)];
517 : 84100 : auto w = pgp[velocityIdx(nmat, 2)];
518 : 84100 : auto vmag = std::sqrt(tk::dot({u, v, w}, {u, v, w}));
519 : :
520 : : // acoustic speed
521 : 84100 : tk::real a = 0.0;
522 [ + + ]: 252300 : for (std::size_t k=0; k<nmat; ++k)
523 : : {
524 [ + + ]: 168200 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
525 [ + - ]: 308808 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
526 : 102936 : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
527 : 205872 : ugp[volfracIdx(nmat, k)], k ) );
528 : : }
529 : : }
530 : :
531 : : // characteristic wave speed
532 : 84100 : auto v_char = vmag + a;
533 : :
534 : : // characteristic length (radius of insphere)
535 [ + - ][ + - ]: 84100 : auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
536 : 84100 : /std::sqrt(24.0);
537 : :
538 : : // element dt
539 : 84100 : local_dte[e] = dx/v_char;
540 : 84100 : mindt = std::min(mindt, local_dte[e]);
541 : : }
542 : :
543 : 450 : return mindt;
544 : : }
545 : :
546 : : void
547 : 3659698 : resetSolidTensors(
548 : : std::size_t nmat,
549 : : std::size_t k,
550 : : std::size_t e,
551 : : tk::real apr_target,
552 : : tk::Fields& U,
553 : : tk::Fields& P )
554 : : // *****************************************************************************
555 : : // Reset the solid tensors
556 : : //! \param[in] nmat Number of materials in this PDE system
557 : : //! \param[in] k Material id whose deformation gradient is required
558 : : //! \param[in] e Id of element whose solution is to be limited
559 : : //! \param[in] apr_target Target partial pressure (alpha_k * p_k)
560 : : //! \param[in/out] U High-order solution vector which gets modified
561 : : //! \param[in/out] P High-order vector of primitives which gets modified
562 : : // *****************************************************************************
563 : : {
564 : 3659698 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
565 : 3659698 : const auto rdof = g_inputdeck.get< tag::rdof >();
566 : :
567 [ + + ]: 3659698 : if (solidx[k] > 0) {
568 [ + + ]: 212256 : for (std::size_t i=0; i<3; ++i) {
569 [ + + ]: 636768 : for (std::size_t j=0; j<3; ++j) {
570 : : // deformation gradient reset
571 [ + + ]: 477576 : if (i==j) U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = 1.0;
572 : 318384 : else U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, 0)) = 0.0;
573 : :
574 : : // Cauchy-stress reset
575 [ + + ]: 636768 : if (i==j) P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, 0))
576 : 318384 : = -apr_target;
577 : 318384 : else P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, 0)) = 0.0;
578 : :
579 : : // high-order reset
580 [ + + ]: 1910304 : for (std::size_t l=1; l<rdof; ++l) {
581 : 1432728 : U(e, deformDofIdx(nmat, solidx[k], i, j, rdof, l)) = 0.0;
582 : 1432728 : P(e, stressDofIdx(nmat, solidx[k], stressCmp[i][j], rdof, l)) = 0.0;
583 : : }
584 : : }
585 : : }
586 : : }
587 : 3659698 : }
588 : :
589 : : std::array< std::array< tk::real, 3 >, 3 >
590 : 19430506 : getDeformGrad(
591 : : std::size_t nmat,
592 : : std::size_t k,
593 : : const std::vector< tk::real >& state )
594 : : // *****************************************************************************
595 : : // Get the inverse deformation gradient tensor for a material at given location
596 : : //! \param[in] nmat Number of materials in this PDE system
597 : : //! \param[in] k Material id whose deformation gradient is required
598 : : //! \param[in] state State vector at location
599 : : //! \return Inverse deformation gradient tensor (g_k) of material k
600 : : // *****************************************************************************
601 : : {
602 : 19430506 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
603 : : std::array< std::array< tk::real, 3 >, 3 > gk;
604 : :
605 [ + + ]: 19430506 : if (solidx[k] > 0) {
606 : : // deformation gradient for solids
607 [ + + ]: 4872120 : for (std::size_t i=0; i<3; ++i) {
608 [ + + ]: 14616360 : for (std::size_t j=0; j<3; ++j)
609 : 10962270 : gk[i][j] = state[deformIdx(nmat,solidx[k],i,j)];
610 : : }
611 : : }
612 : : else {
613 : : // empty vector for fluids
614 : 18212476 : gk = {{}};
615 : : }
616 : :
617 : 19430506 : return gk;
618 : : }
619 : :
620 : : std::array< std::array< tk::real, 3 >, 3 >
621 : 675360 : getCauchyStress(
622 : : std::size_t nmat,
623 : : std::size_t k,
624 : : std::size_t ncomp,
625 : : const std::vector< tk::real >& state )
626 : : // *****************************************************************************
627 : : // Get the elastic Cauchy stress tensor for a material at given location
628 : : //! \param[in] nmat Number of materials in this PDE system
629 : : //! \param[in] k Material id whose deformation gradient is required
630 : : //! \param[in] ncomp Number of components in the PDE system
631 : : //! \param[in] state State vector at location
632 : : //! \return Elastic Cauchy stress tensor (alpha * \sigma_ij) of material k
633 : : // *****************************************************************************
634 : : {
635 : 675360 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
636 : : std::array< std::array< tk::real, 3 >, 3 >
637 : 675360 : asigk{{ {{0,0,0}}, {{0,0,0}}, {{0,0,0}} }};
638 : :
639 : : // elastic Cauchy stress for solids
640 [ + + ]: 675360 : if (solidx[k] > 0) {
641 [ + + ]: 1843200 : for (std::size_t i=0; i<3; ++i) {
642 [ + + ]: 5529600 : for (std::size_t j=0; j<3; ++j)
643 : 4147200 : asigk[i][j] = state[ncomp +
644 : 4147200 : stressIdx(nmat, solidx[k], stressCmp[i][j])];
645 : : }
646 : : }
647 : :
648 : 675360 : return asigk;
649 : : }
650 : :
651 : : bool
652 : 7374171 : haveSolid(
653 : : std::size_t nmat,
654 : : const std::vector< std::size_t >& solidx )
655 : : // *****************************************************************************
656 : : // Check whether we have solid materials in our problem
657 : : //! \param[in] nmat Number of materials in this PDE system
658 : : //! \param[in] solidx Material index indicator
659 : : //! \return true if we have at least one solid, false otherwise.
660 : : // *****************************************************************************
661 : : {
662 : 7374171 : bool haveSolid = false;
663 [ + + ]: 22122528 : for (std::size_t k=0; k<nmat; ++k)
664 [ + + ]: 14748357 : if (solidx[k] > 0) haveSolid = true;
665 : :
666 : 7374171 : return haveSolid;
667 : : }
668 : :
669 : 2098350 : std::size_t numSolids(
670 : : std::size_t nmat,
671 : : const std::vector< std::size_t >& solidx )
672 : : // *****************************************************************************
673 : : // Count total number of solid materials in the problem
674 : : //! \param[in] nmat Number of materials in this PDE system
675 : : //! \param[in] solidx Material index indicator
676 : : //! \return Total number of solid materials in the problem
677 : : // *****************************************************************************
678 : : {
679 : : // count number of solid materials
680 : 2098350 : std::size_t nsld(0);
681 [ + + ]: 6967200 : for (std::size_t k=0; k<nmat; ++k)
682 [ + + ]: 4868850 : if (solidx[k] > 0) ++nsld;
683 : :
684 : 2098350 : return nsld;
685 : : }
686 : :
687 : : } //inciter::
|
