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     7045720 : 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     7045720 :   const auto n_passives = U.size() - THMVACE1D::N_FLUX_INPUTS;
      37             : 
      38             :   const auto rho = rhoA / A;
      39             :   const auto vel = rhouA / rhoA;
      40     7045720 :   const auto v = 1.0 / rho;
      41     7045720 :   const auto e = rhoEA / rhoA - 0.5 * vel * vel;
      42     7045720 :   const auto p = fp.p_from_v_e(v, e);
      43     7045720 :   const auto T = fp.T_from_v_e(v, e);
      44             : 
      45     7045720 :   std::vector<GenericReal<is_ad>> W(THMVACE1D::N_PRIM_VARS + n_passives);
      46     7045720 :   W[THMVACE1D::PRESSURE] = fp.p_from_v_e(v, e);
      47     7045720 :   W[THMVACE1D::VELOCITY] = vel;
      48     7045720 :   W[THMVACE1D::TEMPERATURE] = fp.T_from_v_e(v, e);
      49     7045720 :   for (const auto i : make_range(n_passives))
      50           0 :     W[THMVACE1D::N_PRIM_VARS + i] = U[THMVACE1D::N_FLUX_INPUTS + i] / A;
      51             : 
      52     7045720 :   return W;
      53           0 : }
      54             : 
      55             : /**
      56             :  * Computes the conservative solution vector from the primitive solution vector
      57             :  *
      58             :  * @param[in] W   Primitive solution vector
      59             :  * @param[in] A   Cross-sectional area
      60             :  * @param[in] fp  Fluid properties object
      61             :  */
      62             : template <bool is_ad>
      63             : std::vector<GenericReal<is_ad>>
      64     1786136 : computeConservativeSolutionVector(const std::vector<GenericReal<is_ad>> & W,
      65             :                                   const GenericReal<is_ad> & A,
      66             :                                   const SinglePhaseFluidProperties & fp)
      67             : {
      68             :   const auto & p = W[THMVACE1D::PRESSURE];
      69             :   const auto & T = W[THMVACE1D::TEMPERATURE];
      70             :   const auto & vel = W[THMVACE1D::VELOCITY];
      71     1786136 :   const auto n_passives = W.size() - THMVACE1D::N_PRIM_VARS;
      72             : 
      73     1786136 :   const ADReal rho = fp.rho_from_p_T(p, T);
      74     1786136 :   const ADReal e = fp.e_from_p_rho(p, rho);
      75     1786136 :   const ADReal E = e + 0.5 * vel * vel;
      76             : 
      77     1786136 :   std::vector<GenericReal<is_ad>> U(THMVACE1D::N_FLUX_INPUTS + n_passives);
      78     1786136 :   U[THMVACE1D::RHOA] = rho * A;
      79     1786136 :   U[THMVACE1D::RHOUA] = U[THMVACE1D::RHOA] * vel;
      80     1786136 :   U[THMVACE1D::RHOEA] = U[THMVACE1D::RHOA] * E;
      81     1786136 :   U[THMVACE1D::AREA] = A;
      82     1786136 :   for (const auto i : make_range(n_passives))
      83           0 :     U[THMVACE1D::N_FLUX_INPUTS + i] = W[THMVACE1D::N_PRIM_VARS + i] * A;
      84             : 
      85     1786136 :   return U;
      86           0 : }
      87             : 
      88             : /**
      89             :  * Computes the numerical flux vector from the primitive solution vector
      90             :  *
      91             :  * @param[in] W   Primitive solution vector
      92             :  * @param[in] A   Cross-sectional area
      93             :  * @param[in] fp  Fluid properties object
      94             :  */
      95             : template <bool is_ad>
      96             : std::vector<GenericReal<is_ad>>
      97           2 : computeFluxFromPrimitive(const std::vector<GenericReal<is_ad>> & W,
      98             :                          const GenericReal<is_ad> & A,
      99             :                          const SinglePhaseFluidProperties & fp)
     100             : {
     101             :   const auto & p = W[THMVACE1D::PRESSURE];
     102             :   const auto & T = W[THMVACE1D::TEMPERATURE];
     103             :   const auto & vel = W[THMVACE1D::VELOCITY];
     104           2 :   const auto n_passives = W.size() - THMVACE1D::N_PRIM_VARS;
     105             : 
     106           2 :   const auto rho = fp.rho_from_p_T(p, T);
     107           2 :   const auto e = fp.e_from_p_rho(p, rho);
     108           2 :   const auto E = e + 0.5 * vel * vel;
     109             : 
     110           2 :   std::vector<ADReal> F(THMVACE1D::N_FLUX_OUTPUTS + n_passives, 0.0);
     111           2 :   F[THMVACE1D::MASS] = rho * vel * A;
     112           2 :   F[THMVACE1D::MOMENTUM] = (rho * vel * vel + p) * A;
     113           2 :   F[THMVACE1D::ENERGY] = vel * (rho * E + p) * A;
     114           2 :   for (const auto i : make_range(n_passives))
     115           0 :     F[THMVACE1D::N_FLUX_OUTPUTS + i] = vel * W[THMVACE1D::N_PRIM_VARS + i] * A;
     116             : 
     117           2 :   return F;
     118           0 : }
     119             : 
     120             : /**
     121             :  * Gets the elemental conservative solution vector
     122             :  *
     123             :  * @param[in] elem         Element
     124             :  * @param[in] U_vars       Vector of conservative variable pointers
     125             :  * @param[in] is_implicit  Is implicit?
     126             :  */
     127             : template <bool is_ad>
     128             : std::vector<GenericReal<is_ad>>
     129     5269450 : getElementalSolutionVector(const Elem * elem,
     130             :                            const std::vector<MooseVariable *> & U_vars,
     131             :                            bool is_implicit)
     132             : {
     133             :   mooseAssert(elem, "The supplied element is a nullptr.");
     134             : 
     135     5269450 :   std::vector<GenericReal<is_ad>> U(U_vars.size(), 0.0);
     136             : 
     137     5269450 :   if (is_implicit)
     138             :   {
     139    26347250 :     for (const auto i : make_range(U_vars.size()))
     140             :     {
     141             :       mooseAssert(U_vars[i], "The supplied variable is a nullptr.");
     142    21077800 :       U[i] = U_vars[i]->getElementalValue(elem);
     143             : 
     144    21077800 :       if (i != THMVACE1D::AREA)
     145             :       {
     146             :         std::vector<dof_id_type> dof_indices;
     147    15808350 :         U_vars[i]->dofMap().dof_indices(elem, dof_indices, U_vars[i]->number());
     148    15808350 :         Moose::derivInsert(U[i].derivatives(), dof_indices[0], 1.0);
     149    15808350 :       }
     150             :     }
     151             :   }
     152             :   else
     153             :   {
     154           0 :     for (const auto i : make_range(U_vars.size()))
     155           0 :       U[i] = U_vars[i]->getElementalValueOld(elem);
     156             :   }
     157             : 
     158     5269450 :   return U;
     159           0 : }
     160             : }