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