<commit_before>//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* 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

// Navier-Stokes includes
#include "ConservedVarMaterial.h"
#include "NS.h"

// FluidProperties includes
#include "SinglePhaseFluidProperties.h"

// MOOSE includes
#include "AuxiliarySystem.h"

registerADMooseObject("NavierStokesApp", ConservedVarMaterial);

namespace nms = NS;

defineADValidParams(
    ConservedVarMaterial,
    VarMaterialBase,
    params.addRequiredCoupledVar(nms::density, "density");
    params.addRequiredCoupledVar(nms::total_energy_density, "total fluid energy");
    params.addRequiredCoupledVar(nms::momentum_x, "rhou");
    params.addCoupledVar(nms::momentum_y, "rhov");
    params.addCoupledVar(nms::momentum_z, "rhow");
    params.addCoupledVar(
        nms::T_fluid,
        "optionally use for approximating second spatial derivatives for stabilization");
    params.addClassDescription("Provides access to variables for a conserved variable set "
      "of density, total fluid energy, and momentum"););

using MetaPhysicL::raw_value;

ConservedVarMaterial::ConservedVarMaterial(const InputParameters & params)
  : VarMaterialBase(params),
    _var_rho(adCoupledValue(nms::density)),
    _var_total_energy_density(adCoupledValue(nms::total_energy_density)),
    _var_rho_u(adCoupledValue(nms::momentum_x)),
    _var_rho_v(isCoupled(nms::momentum_y) ? adCoupledValue(nms::momentum_y) : _ad_zero),
    _var_rho_w(isCoupled(nms::momentum_z) ? adCoupledValue(nms::momentum_z) : _ad_zero),
    _var_grad_rho(adCoupledGradient(nms::density)),
    _var_grad_rho_et(adCoupledGradient(nms::total_energy_density)),
    _var_grad_rho_u(adCoupledGradient(nms::momentum_x)),
    _var_grad_rho_v(isCoupled(nms::momentum_y) ? adCoupledGradient(nms::momentum_y) : _ad_grad_zero),
    _var_grad_rho_w(isCoupled(nms::momentum_z) ? adCoupledGradient(nms::momentum_z) : _ad_grad_zero),
    _var_grad_grad_rho(adCoupledSecond(nms::density)),
    _var_grad_grad_rho_u(adCoupledSecond(nms::momentum_x)),
    _var_grad_grad_rho_v(isCoupled(nms::momentum_y) ? adCoupledSecond(nms::momentum_y) : _ad_second_zero),
    _var_grad_grad_rho_w(isCoupled(nms::momentum_z) ? adCoupledSecond(nms::momentum_z) : _ad_second_zero),
    _var_rho_dot(_is_transient ? adCoupledDot(nms::density) : _ad_zero),
    _var_rho_et_dot(_is_transient ? adCoupledDot(nms::total_energy_density) : _ad_zero),
    _var_rho_u_dot(_is_transient ? adCoupledDot(nms::momentum_x) : _ad_zero),
    _var_rho_v_dot(isCoupled(nms::momentum_y) && _is_transient ? adCoupledDot(nms::momentum_y)
                                                          : _ad_zero),
    _var_rho_w_dot(isCoupled(nms::momentum_z) && _is_transient ? adCoupledDot(nms::momentum_z)
                                                          : _ad_zero),
    _var_grad_grad_T_fluid(isCoupled(nms::T_fluid) ? adCoupledSecond(nms::T_fluid)
                                                   : _ad_second_zero)
{
  warnAuxiliaryVariables();
}

bool
ConservedVarMaterial::coupledAuxiliaryVariables() const
{
  auto & sys = _fe_problem.getAuxiliarySystem();
  return sys.hasVariable(nms::density) || sys.hasVariable(nms::total_energy_density) || sys.hasVariable(nms::momentum_x) ||
      sys.hasVariable(nms::momentum_y) || sys.hasVariable(nms::momentum_z);
}

void
ConservedVarMaterial::setNonlinearProperties()
{
  _rho[_qp] = _var_rho[_qp];
  _total_energy_density[_qp] = _var_total_energy_density[_qp];
  _momentum[_qp] = {_var_rho_u[_qp], _var_rho_v[_qp], _var_rho_w[_qp]};

  _grad_rho[_qp] = _var_grad_rho[_qp];
  _grad_rho_et[_qp] = _var_grad_rho_et[_qp];
  _grad_rho_u[_qp] = _var_grad_rho_u[_qp];
  _grad_rho_v[_qp] = _var_grad_rho_v[_qp];
  _grad_rho_w[_qp] = _var_grad_rho_w[_qp];

  _drho_dt[_qp] = _var_rho_dot[_qp];
  _drho_et_dt[_qp] = _var_rho_et_dot[_qp];
  _drho_u_dt[_qp] = _var_rho_u_dot[_qp];
  _drho_v_dt[_qp] = _var_rho_v_dot[_qp];
  _drho_w_dt[_qp] = _var_rho_w_dot[_qp];

  // TODO: Figure out how to approximate second spatial derivatives better
  //
  // The Quasi-linear (strong) residual part of the stabilization term has effectively a laplace
  // operator acting on T_fluid - so we need the second spatial derivatives of temperature. But if
  // the nonlinear variable is not T_fluid (e.g. energy instead) we get temperature via the fluid
  // properties system which doesn't support second derivatives. So with the current state of
  // things, we will not be able to get the full accuracy for the stabilization term - unless you
  // are using PrimitiveVarMaterial or MixedVarMaterial. However, when using linear
  // elements, the aux variables would have been giving us zero for the second derivatives anyway
  // when using regular cartesian grids.  So I would say that the loss of accuracy here is fine
  // for now. And the solution in the future will be to support some sort of approximation or
  // figure out how to get second derivs of fluid properties.
  _grad_grad_T_fluid[_qp] = _var_grad_grad_T_fluid[_qp];
}

void
ConservedVarMaterial::computeQpProperties()
{
  setNonlinearProperties();

  _velocity[_qp] = _momentum[_qp] / _rho[_qp];
  _speed[_qp] = computeSpeed();

  _v[_qp] = 1 / _rho[_qp];
  auto grad_vol = -1 / (_rho[_qp] * _rho[_qp]) * _grad_rho[_qp];
  _grad_vel_x[_qp] = _v[_qp] * _grad_rho_u[_qp] + _momentum[_qp](0) * grad_vol;
  _grad_vel_y[_qp] = _v[_qp] * _grad_rho_v[_qp] + _momentum[_qp](1) * grad_vol;
  _grad_vel_z[_qp] = _v[_qp] * _grad_rho_w[_qp] + _momentum[_qp](2) * grad_vol;

  _grad_grad_vel_x[_qp] = (_var_grad_grad_rho_u[_qp] - outer_product(_grad_rho[_qp], _grad_vel_x[_qp]) -
    outer_product(_grad_vel_x[_qp], _grad_rho[_qp]) - _velocity[_qp](0) * _var_grad_grad_rho[_qp]) / _rho[_qp];
  _grad_grad_vel_y[_qp] = (_var_grad_grad_rho_v[_qp] - outer_product(_grad_rho[_qp], _grad_vel_y[_qp]) -
    outer_product(_grad_vel_y[_qp], _grad_rho[_qp]) - _velocity[_qp](1) * _var_grad_grad_rho[_qp]) / _rho[_qp];
  _grad_grad_vel_z[_qp] = (_var_grad_grad_rho_w[_qp] - outer_product(_grad_rho[_qp], _grad_vel_z[_qp]) -
    outer_product(_grad_vel_z[_qp], _grad_rho[_qp]) - _velocity[_qp](2) * _var_grad_grad_rho[_qp]) / _rho[_qp];

  _specific_internal_energy[_qp] = _total_energy_density[_qp] / _rho[_qp] - (_velocity[_qp] * _velocity[_qp]) / 2;

  Real dTdvol = 0;
  Real dTde = 0;
  Real dpdvol = 0;
  Real dpde = 0;
  _T_fluid[_qp] = _fluid.T_from_v_e(_v[_qp], _specific_total_energy[_qp]);
  _pressure[_qp] = _fluid.p_from_v_e(_v[_qp], _specific_total_energy[_qp]);

  _specific_total_enthalpy[_qp] = (_total_energy_density[_qp] + _pressure[_qp]) / _rho[_qp];
  _total_enthalpy_density[_qp] = _rho[_qp] * _specific_total_enthalpy[_qp];

  auto grad_e = _grad_rho_et[_qp] * _v[_qp] + _total_energy_density[_qp] * grad_vol - _velocity[_qp](0) * _grad_vel_x[_qp] -
                _velocity[_qp](1) * _grad_vel_y[_qp] - _velocity[_qp](2) * _grad_vel_z[_qp];

  _grad_T_fluid[_qp] = grad_vol * dTdvol + grad_e * dTde;
  _grad_pressure[_qp] = grad_vol * dpdvol + grad_e * dpde;

  auto dvelx_dt = (_drho_u_dt[_qp] - _velocity[_qp](0) * _drho_dt[_qp]) / _rho[_qp];
  auto dvely_dt = (_drho_v_dt[_qp] - _velocity[_qp](1) * _drho_dt[_qp]) / _rho[_qp];
  auto dvelz_dt = (_drho_w_dt[_qp] - _velocity[_qp](2) * _drho_dt[_qp]) / _rho[_qp];
  auto de_dt = _drho_et_dt[_qp] / _rho[_qp] - _total_energy_density[_qp] / (_rho[_qp] * _rho[_qp]) * _drho_dt[_qp] -
    (_velocity[_qp](0) * dvelx_dt + _velocity[_qp](1) * dvely_dt + _velocity[_qp](2) * dvelz_dt);
  auto dvol_dt = -1 / (_rho[_qp] * _rho[_qp]) * _drho_dt[_qp];
  _dT_dt[_qp] = dTde * de_dt + dTdvol * dvol_dt;
}<commit_msg>Compute dpdvol and dpde in ConservedVarMaterial<commit_after>//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* 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

// Navier-Stokes includes
#include "ConservedVarMaterial.h"
#include "NS.h"

// FluidProperties includes
#include "SinglePhaseFluidProperties.h"

// MOOSE includes
#include "AuxiliarySystem.h"

registerADMooseObject("NavierStokesApp", ConservedVarMaterial);

namespace nms = NS;

defineADValidParams(
    ConservedVarMaterial,
    VarMaterialBase,
    params.addRequiredCoupledVar(nms::density, "density");
    params.addRequiredCoupledVar(nms::total_energy_density, "total fluid energy");
    params.addRequiredCoupledVar(nms::momentum_x, "rhou");
    params.addCoupledVar(nms::momentum_y, "rhov");
    params.addCoupledVar(nms::momentum_z, "rhow");
    params.addCoupledVar(
        nms::T_fluid,
        "optionally use for approximating second spatial derivatives for stabilization");
    params.addClassDescription("Provides access to variables for a conserved variable set "
      "of density, total fluid energy, and momentum"););

using MetaPhysicL::raw_value;

ConservedVarMaterial::ConservedVarMaterial(const InputParameters & params)
  : VarMaterialBase(params),
    _var_rho(adCoupledValue(nms::density)),
    _var_total_energy_density(adCoupledValue(nms::total_energy_density)),
    _var_rho_u(adCoupledValue(nms::momentum_x)),
    _var_rho_v(isCoupled(nms::momentum_y) ? adCoupledValue(nms::momentum_y) : _ad_zero),
    _var_rho_w(isCoupled(nms::momentum_z) ? adCoupledValue(nms::momentum_z) : _ad_zero),
    _var_grad_rho(adCoupledGradient(nms::density)),
    _var_grad_rho_et(adCoupledGradient(nms::total_energy_density)),
    _var_grad_rho_u(adCoupledGradient(nms::momentum_x)),
    _var_grad_rho_v(isCoupled(nms::momentum_y) ? adCoupledGradient(nms::momentum_y) : _ad_grad_zero),
    _var_grad_rho_w(isCoupled(nms::momentum_z) ? adCoupledGradient(nms::momentum_z) : _ad_grad_zero),
    _var_grad_grad_rho(adCoupledSecond(nms::density)),
    _var_grad_grad_rho_u(adCoupledSecond(nms::momentum_x)),
    _var_grad_grad_rho_v(isCoupled(nms::momentum_y) ? adCoupledSecond(nms::momentum_y) : _ad_second_zero),
    _var_grad_grad_rho_w(isCoupled(nms::momentum_z) ? adCoupledSecond(nms::momentum_z) : _ad_second_zero),
    _var_rho_dot(_is_transient ? adCoupledDot(nms::density) : _ad_zero),
    _var_rho_et_dot(_is_transient ? adCoupledDot(nms::total_energy_density) : _ad_zero),
    _var_rho_u_dot(_is_transient ? adCoupledDot(nms::momentum_x) : _ad_zero),
    _var_rho_v_dot(isCoupled(nms::momentum_y) && _is_transient ? adCoupledDot(nms::momentum_y)
                                                          : _ad_zero),
    _var_rho_w_dot(isCoupled(nms::momentum_z) && _is_transient ? adCoupledDot(nms::momentum_z)
                                                          : _ad_zero),
    _var_grad_grad_T_fluid(isCoupled(nms::T_fluid) ? adCoupledSecond(nms::T_fluid)
                                                   : _ad_second_zero)
{
  warnAuxiliaryVariables();
}

bool
ConservedVarMaterial::coupledAuxiliaryVariables() const
{
  auto & sys = _fe_problem.getAuxiliarySystem();
  return sys.hasVariable(nms::density) || sys.hasVariable(nms::total_energy_density) || sys.hasVariable(nms::momentum_x) ||
      sys.hasVariable(nms::momentum_y) || sys.hasVariable(nms::momentum_z);
}

void
ConservedVarMaterial::setNonlinearProperties()
{
  _rho[_qp] = _var_rho[_qp];
  _total_energy_density[_qp] = _var_total_energy_density[_qp];
  _momentum[_qp] = {_var_rho_u[_qp], _var_rho_v[_qp], _var_rho_w[_qp]};

  _grad_rho[_qp] = _var_grad_rho[_qp];
  _grad_rho_et[_qp] = _var_grad_rho_et[_qp];
  _grad_rho_u[_qp] = _var_grad_rho_u[_qp];
  _grad_rho_v[_qp] = _var_grad_rho_v[_qp];
  _grad_rho_w[_qp] = _var_grad_rho_w[_qp];

  _drho_dt[_qp] = _var_rho_dot[_qp];
  _drho_et_dt[_qp] = _var_rho_et_dot[_qp];
  _drho_u_dt[_qp] = _var_rho_u_dot[_qp];
  _drho_v_dt[_qp] = _var_rho_v_dot[_qp];
  _drho_w_dt[_qp] = _var_rho_w_dot[_qp];

  // TODO: Figure out how to approximate second spatial derivatives better
  //
  // The Quasi-linear (strong) residual part of the stabilization term has effectively a laplace
  // operator acting on T_fluid - so we need the second spatial derivatives of temperature. But if
  // the nonlinear variable is not T_fluid (e.g. energy instead) we get temperature via the fluid
  // properties system which doesn't support second derivatives. So with the current state of
  // things, we will not be able to get the full accuracy for the stabilization term - unless you
  // are using PrimitiveVarMaterial or MixedVarMaterial. However, when using linear
  // elements, the aux variables would have been giving us zero for the second derivatives anyway
  // when using regular cartesian grids.  So I would say that the loss of accuracy here is fine
  // for now. And the solution in the future will be to support some sort of approximation or
  // figure out how to get second derivs of fluid properties.
  _grad_grad_T_fluid[_qp] = _var_grad_grad_T_fluid[_qp];
}

void
ConservedVarMaterial::computeQpProperties()
{
  setNonlinearProperties();

  _velocity[_qp] = _momentum[_qp] / _rho[_qp];
  _speed[_qp] = computeSpeed();

  _v[_qp] = 1 / _rho[_qp];
  auto grad_vol = -1 / (_rho[_qp] * _rho[_qp]) * _grad_rho[_qp];
  _grad_vel_x[_qp] = _v[_qp] * _grad_rho_u[_qp] + _momentum[_qp](0) * grad_vol;
  _grad_vel_y[_qp] = _v[_qp] * _grad_rho_v[_qp] + _momentum[_qp](1) * grad_vol;
  _grad_vel_z[_qp] = _v[_qp] * _grad_rho_w[_qp] + _momentum[_qp](2) * grad_vol;

  _grad_grad_vel_x[_qp] = (_var_grad_grad_rho_u[_qp] - outer_product(_grad_rho[_qp], _grad_vel_x[_qp]) -
    outer_product(_grad_vel_x[_qp], _grad_rho[_qp]) - _velocity[_qp](0) * _var_grad_grad_rho[_qp]) / _rho[_qp];
  _grad_grad_vel_y[_qp] = (_var_grad_grad_rho_v[_qp] - outer_product(_grad_rho[_qp], _grad_vel_y[_qp]) -
    outer_product(_grad_vel_y[_qp], _grad_rho[_qp]) - _velocity[_qp](1) * _var_grad_grad_rho[_qp]) / _rho[_qp];
  _grad_grad_vel_z[_qp] = (_var_grad_grad_rho_w[_qp] - outer_product(_grad_rho[_qp], _grad_vel_z[_qp]) -
    outer_product(_grad_vel_z[_qp], _grad_rho[_qp]) - _velocity[_qp](2) * _var_grad_grad_rho[_qp]) / _rho[_qp];

  _specific_internal_energy[_qp] = _total_energy_density[_qp] / _rho[_qp] - (_velocity[_qp] * _velocity[_qp]) / 2;

  Real dTdvol = 0;
  Real dTde = 0;
  ADReal dpdvol = 0;
  ADReal dpde = 0;
  _T_fluid[_qp] = _fluid.T_from_v_e(_v[_qp], _specific_internal_energy[_qp]);
  _fluid.p_from_v_e(_v[_qp], _specific_internal_energy[_qp], _pressure[_qp], dpdvol, dpde);

  _specific_total_enthalpy[_qp] = (_total_energy_density[_qp] + _pressure[_qp]) / _rho[_qp];
  _total_enthalpy_density[_qp] = _rho[_qp] * _specific_total_enthalpy[_qp];

  auto grad_e = _grad_rho_et[_qp] * _v[_qp] + _total_energy_density[_qp] * grad_vol - _velocity[_qp](0) * _grad_vel_x[_qp] -
                _velocity[_qp](1) * _grad_vel_y[_qp] - _velocity[_qp](2) * _grad_vel_z[_qp];

  _grad_T_fluid[_qp] = grad_vol * dTdvol + grad_e * dTde;
  _grad_pressure[_qp] = grad_vol * dpdvol + grad_e * dpde;

  auto dvelx_dt = (_drho_u_dt[_qp] - _velocity[_qp](0) * _drho_dt[_qp]) / _rho[_qp];
  auto dvely_dt = (_drho_v_dt[_qp] - _velocity[_qp](1) * _drho_dt[_qp]) / _rho[_qp];
  auto dvelz_dt = (_drho_w_dt[_qp] - _velocity[_qp](2) * _drho_dt[_qp]) / _rho[_qp];
  auto de_dt = _drho_et_dt[_qp] / _rho[_qp] - _total_energy_density[_qp] / (_rho[_qp] * _rho[_qp]) * _drho_dt[_qp] -
    (_velocity[_qp](0) * dvelx_dt + _velocity[_qp](1) * dvely_dt + _velocity[_qp](2) * dvelz_dt);
  auto dvol_dt = -1 / (_rho[_qp] * _rho[_qp]) * _drho_dt[_qp];
  _dT_dt[_qp] = dTde * de_dt + dTdvol * dvol_dt;
}
{'fraction_non_alphanumeric': 0.12655324029822212, 'fraction_numerical': 0.0034410246606767347, 'mean_word_length': 5.029197080291971, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 6, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '675243', 'n_tokens_mistral': 6172, 'n_tokens_neox': 5947, 'n_words': 993}