Branch data Line data Source code
1 : : // ***************************************************************************** 2 : : /*! 3 : : \file src/PDE/CompFlow/Problem/FieldOutput.cpp 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 Field outputs for single-material equation solver 9 : : \details This file defines functions for field quantites to be output to 10 : : files for compressible single-material equations. 11 : : */ 12 : : // ***************************************************************************** 13 : : 14 : : #include "FieldOutput.hpp" 15 : : #include "ContainerUtil.hpp" 16 : : #include "History.hpp" 17 : : #include "Inciter/InputDeck/InputDeck.hpp" 18 : : #include "EoS/GetMatProp.hpp" 19 : : #include "ConfigureCompFlow.hpp" 20 : : 21 : : namespace inciter { 22 : : 23 : : extern ctr::InputDeck g_inputdeck; 24 : : 25 : 2654 : std::map< std::string, tk::GetVarFn > CompFlowOutVarFn() 26 : : // ***************************************************************************** 27 : : // Return a map that associates user-specified strings to functions 28 : : //! \return Map that associates user-specified strings to functions that compute 29 : : //! relevant quantities to be output to file 30 : : // ***************************************************************************** 31 : : { 32 : 2654 : std::map< std::string, tk::GetVarFn > OutFnMap; 33 : : 34 : : // Allowed strings for user-def field output vars 35 [ + - ][ + - ]: 2654 : OutFnMap["density"] = compflow::densityOutVar; [ + - ] 36 [ + - ][ + - ]: 2654 : OutFnMap["x-velocity"] = compflow::velocityOutVar<0>; [ + - ] 37 [ + - ][ + - ]: 2654 : OutFnMap["y-velocity"] = compflow::velocityOutVar<1>; [ + - ] 38 [ + - ][ + - ]: 2654 : OutFnMap["z-velocity"] = compflow::velocityOutVar<2>; [ + - ] 39 [ + - ][ + - ]: 2654 : OutFnMap["specific_total_energy"] = compflow::specificTotalEnergyOutVar; [ + - ] 40 [ + - ][ + - ]: 2654 : OutFnMap["volumetric_total_energy"] = compflow::volumetricTotalEnergyOutVar; [ + - ] 41 [ + - ][ + - ]: 2654 : OutFnMap["x-momentum"] = compflow::momentumOutVar<0>; [ + - ] 42 [ + - ][ + - ]: 2654 : OutFnMap["y-momentum"] = compflow::momentumOutVar<1>; [ + - ] 43 [ + - ][ + - ]: 2654 : OutFnMap["z-momentum"] = compflow::momentumOutVar<2>; [ + - ] 44 [ + - ][ + - ]: 2654 : OutFnMap["pressure"] = compflow::pressureOutVar; [ + - ] 45 : : 46 : 2654 : return OutFnMap; 47 : : } 48 : : 49 : 906 : std::vector< std::string > CompFlowSurfNames() 50 : : // ***************************************************************************** 51 : : // Return surface field names to be output to file 52 : : //! \note Every surface will output these fields. 53 : : //! \return Vector of strings labelling surface fields output in file 54 : : // ***************************************************************************** 55 : : { 56 : 906 : std::vector< std::string > n; 57 : : 58 [ + - ][ + - ]: 906 : n.push_back( "density_numerical" ); 59 [ + - ][ + - ]: 906 : n.push_back( "x-velocity_numerical" ); 60 [ + - ][ + - ]: 906 : n.push_back( "y-velocity_numerical" ); 61 [ + - ][ + - ]: 906 : n.push_back( "z-velocity_numerical" ); 62 [ + - ][ + - ]: 906 : n.push_back( "specific_total_energy_numerical" ); 63 [ + - ][ + - ]: 906 : n.push_back( "pressure_numerical" ); 64 : : 65 : 906 : return n; 66 : : } 67 : : 68 : : std::vector< std::vector< tk::real > > 69 : 906 : CompFlowSurfOutput( const std::vector< EOS >& mat_blk, 70 : : const std::map< int, std::vector< std::size_t > >& bnd, 71 : : const tk::Fields& U ) 72 : : // ***************************************************************************** 73 : : // Return surface field output going to file 74 : : //! \param[in] bnd Boundary node/elem lists mapped to side set ids 75 : : //! \param[in] U Solution vector at recent time step 76 : : //! \return Vector of vectors of solution along side sets to be output to file 77 : : // ***************************************************************************** 78 : : { 79 : 906 : std::vector< std::vector< tk::real > > out; 80 : : 81 : : // extract field output along side sets requested 82 [ + + ]: 934 : for (auto s : g_inputdeck.get< tag::field_output, tag::sideset >()) { 83 : : // get node list for side set requested 84 [ + - ]: 28 : auto b = bnd.find(static_cast<int>(s)); 85 [ - + ]: 28 : if (b == end(bnd)) continue; 86 : 28 : const auto& nodes = b->second; 87 [ + - ]: 56 : std::vector< tk::real > surfaceSol( nodes.size() ); 88 : 28 : auto i = out.size(); 89 [ + - ]: 28 : out.insert( end(out), 6, surfaceSol ); 90 : 28 : std::size_t j = 0; 91 [ + + ]: 3766 : for (auto n : nodes) { 92 [ + - ]: 3738 : const auto u = U.extract( n ); 93 [ - + ][ - - ]: 3738 : Assert( u.size() == 5, "Size mismatch" ); [ - - ][ - - ] 94 : 3738 : out[i+0][j] = u[0]; 95 : 3738 : out[i+1][j] = u[1]/u[0]; 96 : 3738 : out[i+2][j] = u[2]/u[0]; 97 : 3738 : out[i+3][j] = u[3]/u[0]; 98 : 3738 : out[i+4][j] = u[4]/u[0]; 99 [ + - ]: 7476 : out[i+5][j] = mat_blk[0].compute< EOS::pressure >( u[0], u[1]/u[0], 100 : 3738 : u[2]/u[0], u[3]/u[0], u[4] ); 101 : 3738 : ++j; 102 : : } 103 : : } 104 : : 105 : 906 : return out; 106 : : } 107 : : 108 : : std::vector< std::vector< tk::real > > 109 : 906 : CompFlowElemSurfOutput( 110 : : const std::vector< EOS >& mat_blk, 111 : : const std::map< int, std::vector< std::size_t > >& bface, 112 : : const std::vector< std::size_t >& triinpoel, 113 : : const tk::Fields& U ) 114 : : // ***************************************************************************** 115 : : // Return element surface field output (on triangle faces) going to file 116 : : //! \param[in] mat_blk Material EOS block 117 : : //! \param[in] bface Boundary-faces mapped to side set ids 118 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids 119 : : //! \param[in] U Solution vector at recent time step 120 : : //! \return Vector of vectors of solution on side set faces to be output to file 121 : : // ***************************************************************************** 122 : : { 123 : 906 : std::vector< std::vector< tk::real > > out; 124 : : 125 : : // extract field output along side sets requested 126 [ + + ]: 934 : for (auto s : g_inputdeck.get< tag::field_output, tag::sideset >()) { 127 : : // get face list for side set requested 128 [ + - ]: 28 : auto b = bface.find(static_cast<int>(s)); 129 [ - + ]: 28 : if (b == end(bface)) continue; 130 : 28 : const auto& faces = b->second; 131 [ + - ]: 56 : std::vector< tk::real > surfaceSol( faces.size() ); 132 : 28 : auto i = out.size(); 133 [ + - ]: 28 : out.insert( end(out), 6, surfaceSol ); 134 : 28 : std::size_t j = 0; 135 [ + + ]: 5460 : for (auto f : faces) { 136 : : // access solutions at nodes 137 [ + - ]: 10864 : auto UA = U.extract(triinpoel[f*3+0]); 138 [ + - ]: 10864 : auto UB = U.extract(triinpoel[f*3+1]); 139 [ + - ]: 5432 : auto UC = U.extract(triinpoel[f*3+2]); 140 : : 141 : 5432 : auto rho = (UA[0] + UB[0] + UC[0]) / 3.0; 142 : 5432 : auto u = (UA[1]/UA[0] + UB[1]/UB[0] + UC[1]/UC[0]) / 3.0; 143 : 5432 : auto v = (UA[2]/UA[0] + UB[2]/UB[0] + UC[2]/UC[0]) / 3.0; 144 : 5432 : auto w = (UA[3]/UA[0] + UB[3]/UB[0] + UC[3]/UC[0]) / 3.0; 145 : 5432 : auto E = (UA[4]/UA[0] + UB[4]/UB[0] + UC[4]/UC[0]) / 3.0; 146 : : 147 : 5432 : out[i+0][j] = rho; 148 : 5432 : out[i+1][j] = u; 149 : 5432 : out[i+2][j] = v; 150 : 5432 : out[i+3][j] = w; 151 : 5432 : out[i+4][j] = E; 152 [ + - ]: 5432 : out[i+5][j] = mat_blk[0].compute< EOS::pressure >( rho, u, v, w, rho*E ); 153 : 5432 : ++j; 154 : : } 155 : : } 156 : : 157 : 906 : return out; 158 : : } 159 : : 160 : 30 : std::vector< std::string > CompFlowHistNames() 161 : : // ***************************************************************************** 162 : : // Return time history field names to be output to file 163 : : //! \note Every time history point will output these fields. 164 : : //! \return Vector of strings labelling time history fields output in file 165 : : // ***************************************************************************** 166 : : { 167 : 30 : std::vector< std::string > n; 168 : : 169 [ + - ][ + - ]: 30 : n.push_back( "density" ); 170 [ + - ][ + - ]: 30 : n.push_back( "x-velocity" ); 171 [ + - ][ + - ]: 30 : n.push_back( "y-velocity" ); 172 [ + - ][ + - ]: 30 : n.push_back( "z-velocity" ); 173 [ + - ][ + - ]: 30 : n.push_back( "energy" ); 174 [ + - ][ + - ]: 30 : n.push_back( "pressure" ); 175 : : 176 : 30 : return n; 177 : : } 178 : : 179 : : std::vector< std::vector< tk::real > > 180 : 254 : CompFlowHistOutput( const std::vector< EOS >& mat_blk, 181 : : const std::vector< HistData >& h, 182 : : const std::vector< std::size_t >& inpoel, 183 : : const tk::Fields& U ) 184 : : // ***************************************************************************** 185 : : // Return time history field output evaluated at time history points 186 : : //! \param[in] h History point data 187 : : //! \param[in] inpoel Mesh element connectivity 188 : : //! \param[in] U Solution vector at recent time step 189 : : //! \return Vector of vectors of solution variables evaluated in all history 190 : : //! points. Inner vector: variables, outer vector: points. 191 : : // ***************************************************************************** 192 : : { 193 [ + - ]: 254 : std::vector< std::vector< tk::real > > out( h.size() ); 194 : : 195 : 254 : std::size_t j = 0; 196 [ + + ]: 332 : for (const auto& p : h) { 197 : 78 : auto e = p.get< tag::elem >(); // host element id 198 : 78 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point 199 [ + - ]: 78 : out[j].resize( 6, 0.0 ); 200 [ + + ]: 390 : for (std::size_t i=0; i<4; ++i) { 201 [ + - ]: 312 : const auto u = U.extract( inpoel[e*4+i] ); 202 [ - + ][ - - ]: 312 : Assert( u.size() == 5, "Size mismatch" ); [ - - ][ - - ] 203 : 312 : out[j][0] += n[i] * u[0]; 204 : 312 : out[j][1] += n[i] * u[1]/u[0]; 205 : 312 : out[j][2] += n[i] * u[2]/u[0]; 206 : 312 : out[j][3] += n[i] * u[3]/u[0]; 207 : 312 : out[j][4] += n[i] * u[4]/u[0]; 208 [ + - ]: 624 : out[j][5] += n[i] * mat_blk[0].compute< EOS::pressure >( u[0], u[1]/u[0], 209 : 624 : u[2]/u[0], u[3]/u[0], u[4] ); 210 : : } 211 : 78 : ++j; 212 : : } 213 : : 214 : 254 : return out; 215 : : } 216 : : 217 : : } //inciter::