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