FlowModel1PhaseUtils.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "SinglePhaseFluidProperties.h"
#include "THMIndicesVACE.h"
#include "MooseVariable.h"

#include "libmesh/elem.h"

namespace FlowModel1PhaseUtils
{

template <bool is_ad>
std::vector<GenericReal<is_ad>>
computePrimitiveSolutionVector(const std::vector<GenericReal<is_ad>> & U,
                               const SinglePhaseFluidProperties & fp)
{
  const auto & rhoA = U[THMVACE1D::RHOA];
  const auto & rhouA = U[THMVACE1D::RHOUA];
  const auto & rhoEA = U[THMVACE1D::RHOEA];
  const auto & A = U[THMVACE1D::AREA];
  const auto n_passives = U.size() - THMVACE1D::N_FLUX_INPUTS;

  const auto rho = rhoA / A;
  const auto vel = rhouA / rhoA;
  const auto v = 1.0 / rho;
  const auto e = rhoEA / rhoA - 0.5 * vel * vel;
  const auto p = fp.p_from_v_e(v, e);
  const auto T = fp.T_from_v_e(v, e);

  std::vector<GenericReal<is_ad>> W(THMVACE1D::N_PRIM_VARS + n_passives);
  W[THMVACE1D::PRESSURE] = fp.p_from_v_e(v, e);
  W[THMVACE1D::VELOCITY] = vel;
  W[THMVACE1D::TEMPERATURE] = fp.T_from_v_e(v, e);
  for (const auto i : make_range(n_passives))
    W[THMVACE1D::N_PRIM_VARS + i] = U[THMVACE1D::N_FLUX_INPUTS + i] / A;

  return W;
}

template <bool is_ad>
std::vector<GenericReal<is_ad>>
computeConservativeSolutionVector(const std::vector<GenericReal<is_ad>> & W,
                                  const GenericReal<is_ad> & A,
                                  const SinglePhaseFluidProperties & fp)
{
  const auto & p = W[THMVACE1D::PRESSURE];
  const auto & T = W[THMVACE1D::TEMPERATURE];
  const auto & vel = W[THMVACE1D::VELOCITY];
  const auto n_passives = W.size() - THMVACE1D::N_PRIM_VARS;

  const ADReal rho = fp.rho_from_p_T(p, T);
  const ADReal e = fp.e_from_p_rho(p, rho);
  const ADReal E = e + 0.5 * vel * vel;

  std::vector<GenericReal<is_ad>> U(THMVACE1D::N_FLUX_INPUTS + n_passives);
  U[THMVACE1D::RHOA] = rho * A;
  U[THMVACE1D::RHOUA] = U[THMVACE1D::RHOA] * vel;
  U[THMVACE1D::RHOEA] = U[THMVACE1D::RHOA] * E;
  U[THMVACE1D::AREA] = A;
  for (const auto i : make_range(n_passives))
    U[THMVACE1D::N_FLUX_INPUTS + i] = W[THMVACE1D::N_PRIM_VARS + i] * A;

  return U;
}

template <bool is_ad>
std::vector<GenericReal<is_ad>>
computeFluxFromPrimitive(const std::vector<GenericReal<is_ad>> & W,
                         const GenericReal<is_ad> & A,
                         const SinglePhaseFluidProperties & fp)
{
  const auto & p = W[THMVACE1D::PRESSURE];
  const auto & T = W[THMVACE1D::TEMPERATURE];
  const auto & vel = W[THMVACE1D::VELOCITY];
  const auto n_passives = W.size() - THMVACE1D::N_PRIM_VARS;

  const auto rho = fp.rho_from_p_T(p, T);
  const auto e = fp.e_from_p_rho(p, rho);
  const auto E = e + 0.5 * vel * vel;

  std::vector<ADReal> F(THMVACE1D::N_FLUX_OUTPUTS + n_passives, 0.0);
  F[THMVACE1D::MASS] = rho * vel * A;
  F[THMVACE1D::MOMENTUM] = (rho * vel * vel + p) * A;
  F[THMVACE1D::ENERGY] = vel * (rho * E + p) * A;
  for (const auto i : make_range(n_passives))
    F[THMVACE1D::N_FLUX_OUTPUTS + i] = vel * W[THMVACE1D::N_PRIM_VARS + i] * A;

  return F;
}

template <bool is_ad>
std::vector<GenericReal<is_ad>>
getElementalSolutionVector(const Elem * elem,
                           const std::vector<MooseVariable *> & U_vars,
                           bool is_implicit)
{
  mooseAssert(elem, "The supplied element is a nullptr.");

  std::vector<GenericReal<is_ad>> U(U_vars.size(), 0.0);

  if (is_implicit)
  {
    for (const auto i : make_range(U_vars.size()))
    {
      mooseAssert(U_vars[i], "The supplied variable is a nullptr.");
      U[i] = U_vars[i]->getElementalValue(elem);

      if (i != THMVACE1D::AREA)
      {
        std::vector<dof_id_type> dof_indices;
        U_vars[i]->dofMap().dof_indices(elem, dof_indices, U_vars[i]->number());
        Moose::derivInsert(U[i].derivatives(), dof_indices[0], 1.0);
      }
    }
  }
  else
  {
    for (const auto i : make_range(U_vars.size()))
      U[i] = U_vars[i]->getElementalValueOld(elem);
  }

  return U;
}
}