Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiSpecies/MiscMultiSpeciesFns.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-species system functions
9 : : \details This file defines functions that required for multi-species
10 : : compressible fluid dynamics.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <iostream>
15 : :
16 : : #include "MiscMultiSpeciesFns.hpp"
17 : : #include "Inciter/InputDeck/InputDeck.hpp"
18 : : #include "Integrate/Basis.hpp"
19 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
20 : : #include "EoS/GetMatProp.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : :
26 : 12 : void initializeSpeciesEoS( std::vector< EOS >& mat_blk )
27 : : // *****************************************************************************
28 : : // Initialize the species block with configured EOS
29 : : //! \param[in,out] mat_blk Material block that gets initialized
30 : : // *****************************************************************************
31 : : {
32 : : // EoS initialization
33 : : // ---------------------------------------------------------------------------
34 : 12 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
35 : 12 : const auto& matprop = g_inputdeck.get< tag::material >();
36 : 12 : const auto& matidxmap = g_inputdeck.get< tag::matidxmap >();
37 : : // assume only one type of species
38 : 12 : auto mateos = matprop[matidxmap.get< tag::eosidx >()[0]].get<tag::eos>();
39 [ + + ]: 36 : for (std::size_t k=0; k<nspec; ++k) {
40 [ + - ]: 24 : mat_blk.emplace_back(mateos, EqType::multispecies, k);
41 : : }
42 : 12 : }
43 : :
44 : : tk::real
45 : 0 : timeStepSizeMultiSpecies(
46 : : const std::vector< EOS >& mat_blk,
47 : : const std::vector< int >& esuf,
48 : : const tk::Fields& geoFace,
49 : : const tk::Fields& geoElem,
50 : : const std::size_t nelem,
51 : : std::size_t nspec,
52 : : const tk::Fields& U,
53 : : const tk::Fields& /*P*/ )
54 : : // *****************************************************************************
55 : : // Time step restriction for multi species cell-centered schemes
56 : : //! \param[in] mat_blk EOS species block
57 : : //! \param[in] esuf Elements surrounding elements array
58 : : //! \param[in] geoFace Face geometry array
59 : : //! \param[in] geoElem Element geometry array
60 : : //! \param[in] nelem Number of elements
61 : : //! \param[in] nspec Number of speciess in this PDE system
62 : : //! \param[in] U High-order solution vector
63 : : // //! \param[in] P High-order vector of primitives
64 : : //! \return Maximum allowable time step based on cfl criterion
65 : : // *****************************************************************************
66 : : {
67 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
68 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
69 : 0 : std::size_t ncomp = U.nprop()/rdof;
70 : :
71 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
72 [ - - ]: 0 : std::vector< tk::real > delt(U.nunk(), 0.0);
73 [ - - ]: 0 : std::vector< tk::real > ugp(ncomp, 0.0);
74 : :
75 : : // compute maximum characteristic speed at all internal element faces
76 [ - - ]: 0 : for (std::size_t f=0; f<esuf.size()/2; ++f)
77 : : {
78 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
79 : 0 : auto er = esuf[2*f+1];
80 : :
81 : : std::array< tk::real, 3 > fn;
82 [ - - ]: 0 : fn[0] = geoFace(f,1);
83 [ - - ]: 0 : fn[1] = geoFace(f,2);
84 [ - - ]: 0 : fn[2] = geoFace(f,3);
85 : :
86 : : // left element
87 : :
88 : : // Compute the basis function for the left element
89 [ - - ]: 0 : std::vector< tk::real > B_l(rdof, 0.0);
90 : 0 : B_l[0] = 1.0;
91 : :
92 : : // get conserved quantities
93 [ - - ]: 0 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
94 : :
95 : : // mixture density
96 : 0 : tk::real rhob(0.0);
97 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
98 : 0 : rhob += ugp[multispecies::densityIdx(nspec, k)];
99 : :
100 : : // advection velocity
101 : 0 : u = ugp[multispecies::momentumIdx(nspec, 0)]/rhob;
102 : 0 : v = ugp[multispecies::momentumIdx(nspec, 1)]/rhob;
103 : 0 : w = ugp[multispecies::momentumIdx(nspec, 2)]/rhob;
104 : :
105 [ - - ][ - - ]: 0 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ - - ]
106 : :
107 : : // acoustic speed
108 [ - - ]: 0 : auto p = mat_blk[0].compute< EOS::pressure >( rhob, u, v, w,
109 : 0 : ugp[multispecies::energyIdx(nspec, 0)] );
110 [ - - ]: 0 : a = mat_blk[0].compute< EOS::soundspeed >( rhob, p );
111 : :
112 [ - - ]: 0 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
113 : :
114 : : // right element
115 [ - - ]: 0 : if (er > -1) {
116 : 0 : std::size_t eR = static_cast< std::size_t >( er );
117 : :
118 : : // Compute the basis function for the right element
119 [ - - ]: 0 : std::vector< tk::real > B_r(rdof, 0.0);
120 : 0 : B_r[0] = 1.0;
121 : :
122 : : // get conserved quantities
123 [ - - ]: 0 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
124 : :
125 : : // mixture density
126 : 0 : rhob = 0.0;
127 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
128 : 0 : rhob += ugp[multispecies::densityIdx(nspec, k)];
129 : :
130 : : // advection velocity
131 : 0 : u = ugp[multispecies::momentumIdx(nspec, 0)]/rhob;
132 : 0 : v = ugp[multispecies::momentumIdx(nspec, 1)]/rhob;
133 : 0 : w = ugp[multispecies::momentumIdx(nspec, 2)]/rhob;
134 : :
135 [ - - ][ - - ]: 0 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
[ - - ]
136 : :
137 : : // acoustic speed
138 [ - - ]: 0 : p = mat_blk[0].compute< EOS::pressure >( rhob, u, v, w,
139 : 0 : ugp[multispecies::energyIdx(nspec, 0)] );
140 [ - - ]: 0 : a = mat_blk[0].compute< EOS::soundspeed >( rhob, p );
141 : :
142 [ - - ]: 0 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
143 : :
144 : 0 : delt[eR] += std::max( dSV_l, dSV_r );
145 : : } else {
146 : 0 : dSV_r = dSV_l;
147 : : }
148 : :
149 : 0 : delt[el] += std::max( dSV_l, dSV_r );
150 : : }
151 : :
152 : 0 : tk::real mindt = std::numeric_limits< tk::real >::max();
153 : :
154 : : // compute allowable dt
155 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
156 : : {
157 [ - - ]: 0 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
158 : : }
159 : :
160 : 0 : return mindt;
161 : : }
162 : :
163 : : } //inciter::
|
