``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170``` ```// ***************************************************************************** /*! \file src/PDE/CompFlow/Problem/VorticalFlow.cpp \copyright 2012-2015 J. Bakosi, 2016-2018 Los Alamos National Security, LLC., 2019-2020 Triad National Security, LLC. All rights reserved. See the LICENSE file for details. \brief Problem configuration for the compressible flow equations \details This file defines a Problem policy class for the compressible flow equations, defined in PDE/CompFlow/CompFlow.h. See PDE/CompFlow/Problem.h for general requirements on Problem policy classes for CompFlow. */ // ***************************************************************************** #include "VorticalFlow.hpp" #include "Inciter/InputDeck/InputDeck.hpp" #include "FieldOutput.hpp" namespace inciter { extern ctr::InputDeck g_inputdeck; } // ::inciter using inciter::CompFlowProblemVorticalFlow; tk::SolutionFn::result_type CompFlowProblemVorticalFlow::solution( ncomp_t system, [[maybe_unused]] ncomp_t ncomp, tk::real x, tk::real y, tk::real z, tk::real ) // ***************************************************************************** //! Evaluate analytical solution at (x,y,z,t) for all components //! \param[in] system Equation system index, i.e., which compressible //! flow equation system we operate on among the systems of PDEs //! \param[in] ncomp Number of scalar components in this PDE system //! \param[in] x X coordinate where to evaluate the solution //! \param[in] y Y coordinate where to evaluate the solution //! \param[in] z Z coordinate where to evaluate the solution //! \return Values of all components evaluated at (x) //! \note The function signature must follow tk::SolutionFn // ***************************************************************************** { Assert( ncomp == ncomp, "Number of scalar components must be " + std::to_string(ncomp) ); using tag::param; using tag::compflow; // manufactured solution parameters const auto& a = g_inputdeck.get< param, compflow, tag::alpha >()[ system ]; const auto& b = g_inputdeck.get< param, compflow, tag::beta >()[ system ]; const auto& p0 = g_inputdeck.get< param, compflow, tag::p0 >()[ system ]; // ratio of specific heats tk::real g = g_inputdeck.get< param, compflow, tag::gamma >()[ system ][0]; // velocity const tk::real ru = a*x - b*y; const tk::real rv = b*x + a*y; const tk::real rw = -2.0*a*z; // total specific energy const tk::real rE = (ru*ru+rv*rv+rw*rw)/2.0 + (p0-2.0*a*a*z*z)/(g-1.0); return {{ 1.0, ru, rv, rw, rE }}; } std::vector< std::string > CompFlowProblemVorticalFlow::fieldNames( ncomp_t ) const // ***************************************************************************** // Return field names to be output to file //! \return Vector of strings labelling fields output in file // ***************************************************************************** { auto n = CompFlowFieldNames(); n.push_back( "density_analytical" ); n.push_back( "x-velocity_analytical" ); n.push_back( "y-velocity_analytical" ); n.push_back( "z-velocity_analytical" ); n.push_back( "specific_total_energy_analytical" ); n.push_back( "pressure_analytical" ); return n; } std::vector< std::vector< tk::real > > CompFlowProblemVorticalFlow::fieldOutput( ncomp_t system, ncomp_t, ncomp_t offset, std::size_t nunk, std::size_t rdof, tk::real, tk::real, const std::vector< tk::real >&, const std::array< std::vector< tk::real >, 3 >& coord, const tk::Fields& U ) const // ***************************************************************************** // Return field output going to file //! \param[in] system Equation system index, i.e., which compressible //! flow equation system we operate on among the systems of PDEs //! \param[in] offset System offset specifying the position of the system of //! PDEs among other systems //! \param[in] nunk Number of unknowns to extract //! \param[in] coord Mesh node coordinates //! \param[in] U Solution vector at recent time step //! \return Vector of vectors to be output to file // ***************************************************************************** { // manufactured solution parameters const auto& a = g_inputdeck.get< tag::param, tag::compflow, tag::alpha >()[system]; const auto& b = g_inputdeck.get< tag::param, tag::compflow, tag::beta >()[system]; const auto& p0 = g_inputdeck.get< tag::param, tag::compflow, tag::p0 >()[system]; // ratio of specific heats tk::real g = g_inputdeck.get< tag::param, tag::compflow, tag::gamma >()[system][0]; auto out = CompFlowFieldOutput( system, offset, nunk, rdof, U ); const auto r = U.extract( 0*rdof, offset );<--- Variable 'r' is assigned a value that is never used. const auto ru = U.extract( 1*rdof, offset ); const auto rv = U.extract( 2*rdof, offset ); const auto rw = U.extract( 3*rdof, offset ); const auto re = U.extract( 4*rdof, offset ); // mesh node coordinates const auto& x = coord[0]; const auto& y = coord[1]; const auto& z = coord[2]; out.push_back( std::vector< tk::real >( nunk, 1.0 ) ); std::vector< tk::real > u = ru; for (std::size_t i=0; i v = rv; for (std::size_t i=0; i w = rw; for (std::size_t i=0; i E = re; for (std::size_t i=0; i P( nunk, 0.0 ); for (std::size_t i=0; i CompFlowProblemVorticalFlow::names( ncomp_t ) const // ***************************************************************************** // Return names of integral variables to be output to diagnostics file //! \return Vector of strings labelling integral variables output // ***************************************************************************** { return { "r", "ru", "rv", "rw", "re" }; } ```