Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "SinglePhaseFluidProperties.h"
      13             : #include "THMIndicesVACE.h"
      14             : #include "MooseVariable.h"
      15             : 
      16             : #include "libmesh/elem.h"
      17             : 
      18             : namespace FlowModel1PhaseUtils
      19             : {
      20             : 
      21             : /**
      22             :  * Computes the primitive solution vector from the conservative solution vector
      23             :  *
      24             :  * @param[in] U   Conservative solution vector
      25             :  * @param[in] fp  Fluid properties object
      26             :  */
      27             : template <bool is_ad>
      28             : std::vector<GenericReal<is_ad>>
      29    18301372 : computePrimitiveSolutionVector(const std::vector<GenericReal<is_ad>> & U,
      30             :                                const SinglePhaseFluidProperties & fp)
      31             : {
      32             :   const auto & rhoA = U[THMVACE1D::RHOA];
      33             :   const auto & rhouA = U[THMVACE1D::RHOUA];
      34             :   const auto & rhoEA = U[THMVACE1D::RHOEA];
      35             :   const auto & A = U[THMVACE1D::AREA];
      36             : 
      37             :   const auto rho = rhoA / A;
      38             :   const auto vel = rhouA / rhoA;
      39    18301372 :   const auto v = 1.0 / rho;
      40    18301372 :   const auto e = rhoEA / rhoA - 0.5 * vel * vel;
      41    18301372 :   const auto p = fp.p_from_v_e(v, e);
      42    18301372 :   const auto T = fp.T_from_v_e(v, e);
      43             : 
      44    18301372 :   std::vector<GenericReal<is_ad>> W(THMVACE1D::N_PRIM_VARS);
      45    18301372 :   W[THMVACE1D::PRESSURE] = fp.p_from_v_e(v, e);
      46    18301372 :   W[THMVACE1D::VELOCITY] = vel;
      47    18301372 :   W[THMVACE1D::TEMPERATURE] = fp.T_from_v_e(v, e);
      48             : 
      49    18301372 :   return W;
      50             : }
      51             : 
      52             : /**
      53             :  * Computes the conservative solution vector from the primitive solution vector
      54             :  *
      55             :  * @param[in] W   Primitive solution vector
      56             :  * @param[in] A   Cross-sectional area
      57             :  * @param[in] fp  Fluid properties object
      58             :  */
      59             : template <bool is_ad>
      60             : std::vector<GenericReal<is_ad>>
      61     4719466 : computeConservativeSolutionVector(const std::vector<GenericReal<is_ad>> & W,
      62             :                                   const GenericReal<is_ad> & A,
      63             :                                   const SinglePhaseFluidProperties & fp)
      64             : {
      65             :   const auto & p = W[THMVACE1D::PRESSURE];
      66             :   const auto & T = W[THMVACE1D::TEMPERATURE];
      67             :   const auto & vel = W[THMVACE1D::VELOCITY];
      68             : 
      69     4719466 :   const ADReal rho = fp.rho_from_p_T(p, T);
      70     4719466 :   const ADReal e = fp.e_from_p_rho(p, rho);
      71     4719466 :   const ADReal E = e + 0.5 * vel * vel;
      72             : 
      73     4719466 :   std::vector<GenericReal<is_ad>> U(THMVACE1D::N_FLUX_INPUTS);
      74     4719466 :   U[THMVACE1D::RHOA] = rho * A;
      75     4719466 :   U[THMVACE1D::RHOUA] = U[THMVACE1D::RHOA] * vel;
      76     4719466 :   U[THMVACE1D::RHOEA] = U[THMVACE1D::RHOA] * E;
      77     4719466 :   U[THMVACE1D::AREA] = A;
      78             : 
      79     4719466 :   return U;
      80             : }
      81             : 
      82             : /**
      83             :  * Computes the numerical flux vector from the primitive solution vector
      84             :  *
      85             :  * @param[in] W   Primitive solution vector
      86             :  * @param[in] A   Cross-sectional area
      87             :  * @param[in] fp  Fluid properties object
      88             :  */
      89             : template <bool is_ad>
      90             : std::vector<GenericReal<is_ad>>
      91           2 : computeFluxFromPrimitive(const std::vector<GenericReal<is_ad>> & W,
      92             :                          const GenericReal<is_ad> & A,
      93             :                          const SinglePhaseFluidProperties & fp)
      94             : {
      95             :   const auto & p = W[THMVACE1D::PRESSURE];
      96             :   const auto & T = W[THMVACE1D::TEMPERATURE];
      97             :   const auto & vel = W[THMVACE1D::VELOCITY];
      98             : 
      99           2 :   const auto rho = fp.rho_from_p_T(p, T);
     100           2 :   const auto e = fp.e_from_p_rho(p, rho);
     101           2 :   const auto E = e + 0.5 * vel * vel;
     102             : 
     103           2 :   std::vector<ADReal> F(THMVACE1D::N_FLUX_OUTPUTS, 0.0);
     104           2 :   F[THMVACE1D::MASS] = rho * vel * A;
     105           2 :   F[THMVACE1D::MOMENTUM] = (rho * vel * vel + p) * A;
     106           2 :   F[THMVACE1D::ENERGY] = vel * (rho * E + p) * A;
     107             : 
     108           2 :   return F;
     109             : }
     110             : 
     111             : /**
     112             :  * Gets the elemental conservative solution vector
     113             :  *
     114             :  * @param[in] elem         Element
     115             :  * @param[in] U_vars       Vector of conservative variable pointers
     116             :  * @param[in] is_implicit  Is implicit?
     117             :  */
     118             : template <bool is_ad>
     119             : std::vector<GenericReal<is_ad>>
     120    13647016 : getElementalSolutionVector(const Elem * elem,
     121             :                            const std::vector<MooseVariable *> & U_vars,
     122             :                            bool is_implicit)
     123             : {
     124             :   mooseAssert(elem, "The supplied element is a nullptr.");
     125             : 
     126    13647016 :   std::vector<GenericReal<is_ad>> U(THMVACE1D::N_FLUX_INPUTS, 0.0);
     127             : 
     128    13647016 :   if (is_implicit)
     129             :   {
     130    68235080 :     for (unsigned int i = 0; i < THMVACE1D::N_FLUX_INPUTS; i++)
     131             :     {
     132             :       mooseAssert(U_vars[i], "The supplied variable is a nullptr.");
     133    54588064 :       U[i] = U_vars[i]->getElementalValue(elem);
     134             :     }
     135             : 
     136             :     std::vector<dof_id_type> dof_indices;
     137             : 
     138    13647016 :     const std::vector<unsigned int> ind = {THMVACE1D::RHOA, THMVACE1D::RHOUA, THMVACE1D::RHOEA};
     139    54588064 :     for (unsigned int j = 0; j < ind.size(); j++)
     140             :     {
     141    40941048 :       const auto i = ind[j];
     142    40941048 :       U_vars[i]->dofMap().dof_indices(elem, dof_indices, U_vars[i]->number());
     143    40941048 :       Moose::derivInsert(U[i].derivatives(), dof_indices[0], 1.0);
     144             :     }
     145             :   }
     146             :   else
     147             :   {
     148           0 :     for (unsigned int i = 0; i < THMVACE1D::N_FLUX_INPUTS; i++)
     149           0 :       U[i] = U_vars[i]->getElementalValueOld(elem);
     150             :   }
     151             : 
     152    13647016 :   return U;
     153             : }
     154             : }