Added VPP output options to flow channel

modules/thermal_hydraulics/doc/content/source/components/FlowChannel1Phase.md

@@ -24,6 +24,11 @@ the following parameters:
 - [!param](/Components/FlowChannel1Phase/initial_T)
 - [!param](/Components/FlowChannel1Phase/initial_vel)
 
+This component offers options to output quantities via vector post-processors:
+
+- [!param](/Components/FlowChannel1Phase/vpp_vars): Creates an [ElementValueSampler.md] with a vector for each of the listed variables (see below for a list of variables).
+- [!param](/Components/FlowChannel1Phase/create_flux_vpp): Creates a [NumericalFlux3EqnInternalValues.md], which creates a vector for each numerical flux component (mass, momentum, energy) at the internal sides.
+
 !syntax parameters /Components/FlowChannel1Phase
 
 ## Mesh id=mesh

modules/thermal_hydraulics/doc/content/source/vectorpostprocessors/NumericalFlux3EqnInternalValues.md

@@ -0,0 +1,20 @@
+# NumericalFlux3EqnInternalValues
+
+This vector post-processor computes the mass, momentum, and energy numerical fluxes for
+a [FlowChannel1Phase.md] at the internal sides. The following vectors are created:
+
+- `mass_flux`
+- `momentum_flux`
+- `energy_flux`
+
+!alert note
+For a given pair of adjacent elements, the numerical flux can be different on each
+side, due to non-conservative terms. In particular, if the cross-sectional area
+differs between the elements, a non-conservative flux appears. The approach taken
+here is to perform an arithmetic average of the fluxes on each side.
+
+!syntax parameters /VectorPostprocessors/NumericalFlux3EqnInternalValues
+
+!syntax inputs /VectorPostprocessors/NumericalFlux3EqnInternalValues
+
+!syntax children /VectorPostprocessors/NumericalFlux3EqnInternalValues

modules/thermal_hydraulics/include/components/Component1D.h

@@ -73,6 +73,7 @@ class Component1D : public GeneratedMeshComponent
 
 protected:
   virtual bool usingSecondOrderMesh() const override;
+  std::string sortBy() const;
 
   /// Map of end type to a list of connections
   std::map<Component1DConnection::EEndType, std::vector<Connection>> _connections;

modules/thermal_hydraulics/include/components/FlowChannel1Phase.h

@@ -21,10 +21,12 @@ class FlowChannel1Phase : public FlowChannel1PhaseBase
 
   FlowChannel1Phase(const InputParameters & params);
 
+  virtual void addMooseObjects() override;
   virtual const THM::FlowModelID & getFlowModelID() const override { return THM::FM_SINGLE_PHASE; }
   virtual std::vector<std::string> ICParameters() const override;
 
 protected:
   virtual void checkFluidProperties() const override;
   virtual std::string flowModelClassName() const override;
+  void addNumericalFluxVectorPostprocessor();
 };

modules/thermal_hydraulics/include/vectorpostprocessors/NumericalFlux3EqnInternalValues.h

@@ -0,0 +1,57 @@
+//* 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 "InternalSideVectorPostprocessor.h"
+#include "SamplerBase.h"
+
+class ADNumericalFlux3EqnBase;
+
+/**
+ * Computes internal fluxes for FlowChannel1Phase.
+ */
+class NumericalFlux3EqnInternalValues : public InternalSideVectorPostprocessor,
+                                        protected SamplerBase
+{
+public:
+  static InputParameters validParams();
+
+  NumericalFlux3EqnInternalValues(const InputParameters & parameters);
+
+  virtual void initialize() override;
+  virtual void execute() override;
+  virtual void finalize() override;
+
+  using SamplerBase::threadJoin;
+  virtual void threadJoin(const UserObject & y) override;
+
+protected:
+  /// Area in current element
+  const ADVariableValue & _A1;
+  /// Area in neighbor element
+  const ADVariableValue & _A2;
+
+  /// Reconstructed rho*A values in current element
+  const ADMaterialProperty<Real> & _rhoA1;
+  /// Reconstructed rho*u*A values in current element
+  const ADMaterialProperty<Real> & _rhouA1;
+  /// Reconstructed rho*E*A values in current element
+  const ADMaterialProperty<Real> & _rhoEA1;
+
+  /// Reconstructed rho*A values in neighbor element
+  const ADMaterialProperty<Real> & _rhoA2;
+  /// Reconstructed rho*u*A values in neighbor element
+  const ADMaterialProperty<Real> & _rhouA2;
+  /// Reconstructed rho*E*A values in neighbor element
+  const ADMaterialProperty<Real> & _rhoEA2;
+
+  /// Numerical flux user object
+  const ADNumericalFlux3EqnBase & _numerical_flux;
+};

modules/thermal_hydraulics/src/components/Component1D.C

@@ -164,3 +164,25 @@ Component1D::getConnections(Component1DConnection::EEndType end_type) const
   else
     mooseError(name(), ": Invalid end type (", end_type, ").");
 }
+
+std::string
+Component1D::sortBy() const
+{
+  // choose the dominant direction
+  std::string dominant_direction = "x";
+  const Real x_abs = std::abs(_dir(0));
+  const Real y_abs = std::abs(_dir(1));
+  const Real z_abs = std::abs(_dir(2));
+  Real max_value = x_abs;
+  if (y_abs > max_value)
+  {
+    dominant_direction = "y";
+    max_value = y_abs;
+  }
+  if (z_abs > max_value)
+  {
+    dominant_direction = "z";
+    max_value = z_abs;
+  }
+  return dominant_direction;
+}

modules/thermal_hydraulics/src/components/FlowChannel1Phase.C

@@ -10,6 +10,7 @@
 #include "FlowChannel1Phase.h"
 #include "FlowModelSinglePhase.h"
 #include "SinglePhaseFluidProperties.h"
+#include "THMNames.h"
 
 registerMooseObject("ThermalHydraulicsApp", FlowChannel1Phase);
 
@@ -27,6 +28,10 @@ FlowChannel1Phase::validParams()
       "scaling_factor_1phase",
       sf_1phase,
       "Scaling factors for each single phase variable (rhoA, rhouA, rhoEA)");
+  params.addParam<bool>(
+      "create_flux_vpp",
+      false,
+      "If true, create a VectorPostprocessor with the the mass, momentum, and energy side fluxes");
 
   params.addParamNamesToGroup("scaling_factor_1phase", "Numerical scheme");
   params.addClassDescription("1-phase 1D flow channel");
@@ -57,3 +62,25 @@ FlowChannel1Phase::ICParameters() const
 {
   return {"initial_p", "initial_T", "initial_vel"};
 }
+
+void
+FlowChannel1Phase::addMooseObjects()
+{
+  FlowChannel1PhaseBase::addMooseObjects();
+
+  if (getParam<bool>("create_flux_vpp"))
+    addNumericalFluxVectorPostprocessor();
+}
+
+void
+FlowChannel1Phase::addNumericalFluxVectorPostprocessor()
+{
+  const std::string class_name = "NumericalFlux3EqnInternalValues";
+  InputParameters params = _factory.getValidParams(class_name);
+  params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
+  params.set<UserObjectName>("numerical_flux") = _numerical_flux_name;
+  params.set<std::vector<VariableName>>("A_linear") = {THM::AREA_LINEAR};
+  params.set<MooseEnum>("sort_by") = sortBy();
+  params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
+  getTHMProblem().addVectorPostprocessor(class_name, name() + "_flux_vpp", params);
+}

modules/thermal_hydraulics/src/components/FlowChannelBase.C

@@ -98,6 +98,8 @@ FlowChannelBase::validParams()
                         "If true, when there are multiple heat transfer components connected to "
                         "this flow channel, use their index for naming related quantities; "
                         "otherwise, use the name of the heat transfer component.");
+  params.addParam<std::vector<VariableName>>(
+      "vpp_vars", {}, "Variables to add in an ElementValueSampler");
 
   params.setDocString(
       "orientation",
@@ -321,6 +323,18 @@ FlowChannelBase::addMooseObjects()
     }
   }
 
+  const auto & vpp_vars = getParam<std::vector<VariableName>>("vpp_vars");
+  if (vpp_vars.size() > 0)
+  {
+    const std::string class_name = "ElementValueSampler";
+    InputParameters params = _factory.getValidParams(class_name);
+    params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
+    params.set<std::vector<VariableName>>("variable") = vpp_vars;
+    params.set<MooseEnum>("sort_by") = sortBy();
+    params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
+    getTHMProblem().addVectorPostprocessor(class_name, name() + "_vars_vpp", params);
+  }
+
   _flow_model->addMooseObjects();
 
   for (const auto & closures : _closures_objects)

modules/thermal_hydraulics/src/vectorpostprocessors/NumericalFlux3EqnInternalValues.C

@@ -0,0 +1,90 @@
+//* 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
+
+#include "NumericalFlux3EqnInternalValues.h"
+#include "THMIndicesVACE.h"
+#include "ADNumericalFlux3EqnBase.h"
+
+registerMooseObject("ThermalHydraulicsApp", NumericalFlux3EqnInternalValues);
+
+InputParameters
+NumericalFlux3EqnInternalValues::validParams()
+{
+  InputParameters params = InternalSideVectorPostprocessor::validParams();
+  params += SamplerBase::validParams();
+  params.addRequiredCoupledVar("A_linear", "Cross-sectional area, linear");
+  params.addRequiredParam<UserObjectName>("numerical_flux", "Name of numerical flux user object");
+  params.addClassDescription("Computes internal fluxes for FlowChannel1Phase.");
+  return params;
+}
+
+NumericalFlux3EqnInternalValues::NumericalFlux3EqnInternalValues(const InputParameters & parameters)
+  : InternalSideVectorPostprocessor(parameters),
+    SamplerBase(parameters, this, _communicator),
+    _A1(adCoupledValue("A_linear")),
+    _A2(adCoupledNeighborValue("A_linear")),
+    _rhoA1(getADMaterialProperty<Real>("rhoA")),
+    _rhouA1(getADMaterialProperty<Real>("rhouA")),
+    _rhoEA1(getADMaterialProperty<Real>("rhoEA")),
+    _rhoA2(getNeighborADMaterialProperty<Real>("rhoA")),
+    _rhouA2(getNeighborADMaterialProperty<Real>("rhouA")),
+    _rhoEA2(getNeighborADMaterialProperty<Real>("rhoEA")),
+    _numerical_flux(getUserObject<ADNumericalFlux3EqnBase>("numerical_flux"))
+{
+  std::vector<std::string> var_names(THMVACE1D::N_FLUX_OUTPUTS);
+  var_names[THMVACE1D::MASS] = "mass_flux";
+  var_names[THMVACE1D::MOMENTUM] = "momentum_flux";
+  var_names[THMVACE1D::ENERGY] = "energy_flux";
+
+  SamplerBase::setupVariables(var_names);
+}
+
+void
+NumericalFlux3EqnInternalValues::initialize()
+{
+  SamplerBase::initialize();
+}
+
+void
+NumericalFlux3EqnInternalValues::execute()
+{
+  // Assume we are in 1D, and internal sides have only a single quadrature point
+  const unsigned int _qp = 0;
+
+  std::vector<ADReal> U1 = {_rhoA1[_qp], _rhouA1[_qp], _rhoEA1[_qp], _A1[_qp]};
+  std::vector<ADReal> U2 = {_rhoA2[_qp], _rhouA2[_qp], _rhoEA2[_qp], _A2[_qp]};
+
+  const Real nLR_dot_d = _current_side * 2 - 1.0;
+
+  const std::vector<ADReal> & flux_elem_ad =
+      _numerical_flux.getFlux(_current_side, _current_elem->id(), true, U1, U2, nLR_dot_d);
+  const std::vector<ADReal> & flux_neig_ad =
+      _numerical_flux.getFlux(_current_side, _current_elem->id(), false, U1, U2, nLR_dot_d);
+
+  // Convert vector to non-AD
+  std::vector<Real> flux(flux_elem_ad.size());
+  for (const auto i : index_range(flux))
+    flux[i] = MetaPhysicL::raw_value(0.5 * (flux_elem_ad[i] + flux_neig_ad[i]));
+
+  SamplerBase::addSample(_q_point[_qp], _current_elem->id(), flux);
+}
+
+void
+NumericalFlux3EqnInternalValues::finalize()
+{
+  SamplerBase::finalize();
+}
+
+void
+NumericalFlux3EqnInternalValues::threadJoin(const UserObject & y)
+{
+  const auto & vpp = static_cast<const NumericalFlux3EqnInternalValues &>(y);
+
+  SamplerBase::threadJoin(vpp);
+}

modules/thermal_hydraulics/test/tests/components/flow_channel_1phase/gold/output_vpps_out_pipe_flux_vpp_0002.csv

@@ -0,0 +1,5 @@
+energy_flux,id,mass_flux,momentum_flux,x,y,z
+42923519.842431,0,142.59242584674,109018.77462793,0.2,0,1.2246467991474e-17
+42926546.468124,1,142.60899170696,108990.26418095,0.4,0,2.4492935982947e-17
+42929591.848237,2,142.62563642346,108961.75041735,0.6,0,3.6739403974421e-17
+42932655.998424,3,142.64236005822,108933.23332136,0.8,0,4.8985871965894e-17

modules/thermal_hydraulics/test/tests/components/flow_channel_1phase/gold/output_vpps_out_pipe_vars_vpp_0002.csv

@@ -0,0 +1,6 @@
+T,id,p,x,y,z
+291.32385895141,0,90222.60983281,0.1,0,6.1232339957368e-18
+291.30184724105,1,90182.891036046,0.3,0,1.836970198721e-17
+291.27977673955,2,90143.131282025,0.5,0,3.0616169978684e-17
+291.25764992243,3,90103.333225525,0.7,0,4.2862637970157e-17
+291.23546671384,4,90063.496755943,0.9,0,5.5109105961631e-17

modules/thermal_hydraulics/test/tests/components/flow_channel_1phase/output_vpps.i

@@ -0,0 +1,75 @@
+[GlobalParams]
+  gravity_vector = '0 0 0'
+  scaling_factor_1phase = '1 1 1e-5'
+[]
+
+[FluidProperties]
+  [fp]
+    type = IdealGasFluidProperties
+  []
+[]
+
+[Closures]
+  [simple_closures]
+    type = Closures1PhaseSimple
+  []
+[]
+
+[Components]
+  [inlet]
+    type = InletStagnationPressureTemperature1Phase
+    input = 'pipe:in'
+    p0 = 1e5
+    T0 = 300
+  []
+  [pipe]
+    type = FlowChannel1Phase
+    position = '0 0 0'
+    orientation = '1 0 0'
+    length = 1.0
+    n_elems = 5
+    A = 1.0
+    initial_T = 300
+    initial_p = 1e5
+    initial_vel = 0
+    fp = fp
+    closures = simple_closures
+    f = 0
+    vpp_vars = 'p T'
+    create_flux_vpp = true
+  []
+  [outlet]
+    type = Outlet1Phase
+    input = 'pipe:out'
+    p = 0.9e5
+  []
+[]
+
+[Preconditioning]
+  [pc]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  scheme = bdf2
+
+  start_time = 0
+  dt = 1
+  num_steps = 1
+
+  solve_type = NEWTON
+  nl_rel_tol = 1e-8
+  nl_abs_tol = 1e-8
+  nl_max_its = 10
+
+  l_tol = 1e-3
+  l_max_its = 10
+[]
+
+[Outputs]
+  csv = true
+  execute_on = 'FINAL'
+[]