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 : :
72 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
73 : 0 : std::vector< tk::real > delt(U.nunk(), 0.0);
74 [ - - ][ - - ]: 0 : std::vector< tk::real > ugp(ncomp, 0.0);
75 : :
76 : : // compute maximum characteristic speed at all internal element faces
77 [ - - ]: 0 : for (std::size_t f=0; f<esuf.size()/2; ++f)
78 : : {
79 [ - - ]: 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
80 : 0 : auto er = esuf[2*f+1];
81 : :
82 : : std::array< tk::real, 3 > fn;
83 : : fn[0] = geoFace(f,1);
84 : : fn[1] = geoFace(f,2);
85 : : fn[2] = geoFace(f,3);
86 : :
87 : : // left element
88 : :
89 : : // Compute the basis function for the left element
90 [ - - ][ - - ]: 0 : std::vector< tk::real > B_l(rdof, 0.0);
91 : 0 : B_l[0] = 1.0;
92 : :
93 : : // get conserved quantities
94 [ - - ]: 0 : ugp = eval_state(ncomp, rdof, ndof, el, U, B_l);
95 : :
96 : : // initialize mixture
97 [ - - ]: 0 : Mixture mix(nspec, ugp, mat_blk);
98 : :
99 : : // mixture density
100 : 0 : auto rhob = mix.get_mix_density();
101 : :
102 : : // advection velocity
103 [ - - ]: 0 : u = ugp[multispecies::momentumIdx(nspec, 0)]/rhob;
104 : 0 : v = ugp[multispecies::momentumIdx(nspec, 1)]/rhob;
105 [ - - ]: 0 : w = ugp[multispecies::momentumIdx(nspec, 2)]/rhob;
106 : :
107 [ - - ]: 0 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
108 : :
109 : : // acoustic speed
110 [ - - ]: 0 : auto p = mix.pressure( rhob, u, v, w,
111 : : ugp[multispecies::energyIdx(nspec, 0)], mat_blk);
112 [ - - ]: 0 : a = mix.frozen_soundspeed( rhob, p, mat_blk );
113 : :
114 : 0 : dSV_l = geoFace(f,0) * (std::fabs(vn) + a);
115 : :
116 : : // right element
117 [ - - ]: 0 : if (er > -1) {
118 : 0 : std::size_t eR = static_cast< std::size_t >( er );
119 : :
120 : : // Compute the basis function for the right element
121 [ - - ][ - - ]: 0 : std::vector< tk::real > B_r(rdof, 0.0);
122 : 0 : B_r[0] = 1.0;
123 : :
124 : : // get conserved quantities
125 [ - - ]: 0 : ugp = eval_state( ncomp, rdof, ndof, eR, U, B_r);
126 : :
127 : : // initialize mixture
128 [ - - ]: 0 : Mixture mixr(nspec, ugp, mat_blk);
129 : :
130 : : // mixture density
131 : 0 : rhob = mixr.get_mix_density();
132 : :
133 : : // advection velocity
134 [ - - ]: 0 : u = ugp[multispecies::momentumIdx(nspec, 0)]/rhob;
135 : 0 : v = ugp[multispecies::momentumIdx(nspec, 1)]/rhob;
136 [ - - ]: 0 : w = ugp[multispecies::momentumIdx(nspec, 2)]/rhob;
137 : :
138 : 0 : vn = u*geoFace(f,1) + v*geoFace(f,2) + w*geoFace(f,3);
139 : :
140 : : // acoustic speed
141 [ - - ]: 0 : p = mixr.pressure( rhob, u, v, w,
142 : : ugp[multispecies::energyIdx(nspec, 0)], mat_blk );
143 [ - - ]: 0 : a = mixr.frozen_soundspeed( rhob, p, mat_blk );
144 : :
145 [ - - ]: 0 : dSV_r = geoFace(f,0) * (std::fabs(vn) + a);
146 : :
147 [ - - ]: 0 : delt[eR] += std::max( dSV_l, dSV_r );
148 : : } else {
149 : 0 : dSV_r = dSV_l;
150 : : }
151 : :
152 [ - - ]: 0 : delt[el] += std::max( dSV_l, dSV_r );
153 : : }
154 : :
155 : 0 : tk::real mindt = std::numeric_limits< tk::real >::max();
156 : :
157 : : // compute allowable dt
158 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
159 : : {
160 [ - - ]: 0 : mindt = std::min( mindt, geoElem(e,0)/delt[e] );
161 : : }
162 : :
163 [ - - ]: 0 : return mindt;
164 : : }
165 : :
166 : : } //inciter::
|
