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             : #include "ADGateValve1PhaseUserObject.h"
      11             : #include "THMIndicesVACE.h"
      12             : #include "ADNumericalFlux3EqnBase.h"
      13             : #include "Function.h"
      14             : #include "Numerics.h"
      15             : #include "metaphysicl/parallel_numberarray.h"
      16             : #include "metaphysicl/parallel_dualnumber.h"
      17             : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
      18             : #include "libmesh/parallel_algebra.h"
      19             : 
      20             : registerMooseObject("ThermalHydraulicsApp", ADGateValve1PhaseUserObject);
      21             : 
      22             : const std::vector<std::pair<std::string, unsigned int>>
      23             :     ADGateValve1PhaseUserObject::_varname_eq_index_pairs{
      24             :         std::pair<std::string, unsigned int>("rhoA", THMVACE1D::MASS),
      25             :         std::pair<std::string, unsigned int>("rhouA", THMVACE1D::MOMENTUM),
      26             :         std::pair<std::string, unsigned int>("rhoEA", THMVACE1D::ENERGY)};
      27             : 
      28             : InputParameters
      29          78 : ADGateValve1PhaseUserObject::validParams()
      30             : {
      31          78 :   InputParameters params = ADFlowJunctionUserObject::validParams();
      32             : 
      33         156 :   params.addRequiredParam<Real>("open_area_fraction",
      34             :                                 "Fraction of possible flow area that is open");
      35         156 :   params.addParam<Real>("open_area_fraction_min", 1e-10, "Minimum open area fraction");
      36             : 
      37         156 :   params.addRequiredCoupledVar("A", "Cross-sectional area of connected flow channels");
      38         156 :   params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
      39         156 :   params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
      40         156 :   params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
      41             : 
      42         156 :   params.addRequiredParam<UserObjectName>("numerical_flux", "Numerical flux user object name");
      43             : 
      44         156 :   params.addRequiredParam<std::string>("component_name", "Name of the associated component");
      45             : 
      46         156 :   params.addClassDescription("Gate valve user object for 1-phase flow");
      47             : 
      48         156 :   params.declareControllable("open_area_fraction");
      49             : 
      50          78 :   return params;
      51           0 : }
      52             : 
      53          42 : ADGateValve1PhaseUserObject::ADGateValve1PhaseUserObject(const InputParameters & params)
      54             :   : ADFlowJunctionUserObject(params),
      55             : 
      56          42 :     _f_open(getParam<Real>("open_area_fraction")),
      57          84 :     _f_open_min(getParam<Real>("open_area_fraction_min")),
      58             : 
      59          42 :     _A(adCoupledValue("A")),
      60          42 :     _rhoA(adCoupledValue("rhoA")),
      61          42 :     _rhouA(adCoupledValue("rhouA")),
      62          42 :     _rhoEA(adCoupledValue("rhoEA")),
      63             : 
      64          42 :     _rhoA_jvar(coupled("rhoA")),
      65          42 :     _rhouA_jvar(coupled("rhouA")),
      66          42 :     _rhoEA_jvar(coupled("rhoEA")),
      67             : 
      68          84 :     _p(getADMaterialProperty<Real>("p")),
      69             : 
      70          42 :     _numerical_flux(getUserObject<ADNumericalFlux3EqnBase>("numerical_flux")),
      71             : 
      72          84 :     _component_name(getParam<std::string>("component_name")),
      73             : 
      74          84 :     _dir(getMaterialProperty<RealVectorValue>("direction")),
      75             : 
      76          42 :     _solutions(_n_connections),
      77          42 :     _fluxes(_n_connections, std::vector<ADReal>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
      78          42 :     _dof_indices(_n_connections, std::vector<dof_id_type>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
      79             : 
      80          42 :     _stored_p(_n_connections),
      81             : 
      82          42 :     _elem_ids(_n_connections),
      83          42 :     _local_side_ids(_n_connections),
      84             : 
      85          42 :     _areas(_n_connections),
      86          84 :     _directions(_n_connections)
      87             : {
      88          42 : }
      89             : 
      90             : void
      91        4766 : ADGateValve1PhaseUserObject::initialize()
      92             : {
      93        4766 :   _connection_indices.clear();
      94             : 
      95        4766 :   std::vector<ADReal> zero(THMVACE1D::N_FLUX_INPUTS, ADReal(0.));
      96       14298 :   for (auto & s : _solutions)
      97        9532 :     s = zero;
      98       19064 :   for (unsigned int i = 0; i < THMVACE1D::N_FLUX_OUTPUTS; i++)
      99             :   {
     100       14298 :     _fluxes[0][i] = 0;
     101       14298 :     _fluxes[1][i] = 0;
     102             :   }
     103        4766 : }
     104             : 
     105             : void
     106        5740 : ADGateValve1PhaseUserObject::execute()
     107             : {
     108             :   // Get connection index
     109        5740 :   const unsigned int c = getBoundaryIDIndex();
     110        5740 :   _connection_indices.push_back(c);
     111             : 
     112             :   // Store DoF indices
     113       22960 :   for (const auto & varname_eq_index_pair : _varname_eq_index_pairs)
     114             :   {
     115       17220 :     MooseVariable * var = getVar(varname_eq_index_pair.first, 0);
     116             :     auto && dofs = var->dofIndices();
     117             :     mooseAssert(dofs.size() == 1, "There should be only one DoF index.");
     118       17220 :     _dof_indices[c][varname_eq_index_pair.second] = dofs[0];
     119             :   }
     120             : 
     121             :   // Store solution vector for connection
     122        5740 :   std::vector<ADReal> U(THMVACE1D::N_FLUX_INPUTS, ADReal(0.));
     123        5740 :   U[THMVACE1D::RHOA] = _rhoA[0];
     124        5740 :   U[THMVACE1D::RHOUA] = _rhouA[0];
     125        5740 :   U[THMVACE1D::RHOEA] = _rhoEA[0];
     126        5740 :   U[THMVACE1D::AREA] = _A[0];
     127        5740 :   _solutions[c] = U;
     128             : 
     129             :   // Store element ID and local side ID for connection
     130        5740 :   _elem_ids[c] = _current_elem->id();
     131        5740 :   _local_side_ids[c] = _current_side;
     132             : 
     133             :   // Store direction and area of channel at junction (used for error-checking)
     134        5740 :   _areas[c] = _A[0];
     135        5740 :   _directions[c] = _dir[0];
     136             : 
     137        5740 :   _stored_p[c] = _p[0];
     138        5740 : }
     139             : 
     140             : void
     141         482 : ADGateValve1PhaseUserObject::threadJoin(const UserObject & uo)
     142             : {
     143             :   const auto & junction_uo = static_cast<const ADGateValve1PhaseUserObject &>(uo);
     144             : 
     145             :   // Store the data computed/retrieved in the other threads
     146        1205 :   for (unsigned int i = 0; i < junction_uo._connection_indices.size(); i++)
     147             :   {
     148         723 :     const unsigned int c = junction_uo._connection_indices[i];
     149             : 
     150         723 :     _dof_indices[c] = junction_uo._dof_indices[c];
     151         723 :     _solutions[c] = junction_uo._solutions[c];
     152         723 :     _elem_ids[c] = junction_uo._elem_ids[c];
     153         723 :     _local_side_ids[c] = junction_uo._local_side_ids[c];
     154         723 :     _areas[c] = junction_uo._areas[c];
     155         723 :     _directions[c] = junction_uo._directions[c];
     156         723 :     _stored_p[c] = junction_uo._stored_p[c];
     157             :   }
     158         482 : }
     159             : 
     160             : void
     161        4284 : ADGateValve1PhaseUserObject::finalize()
     162             : {
     163             :   // Check direction compatibility
     164        4284 :   if (!THM::areParallelVectors(_directions[0], _directions[1]))
     165           2 :     mooseError(_component_name, ": The connected channels must be parallel at the junction.");
     166             : 
     167       12846 :   for (unsigned int i = 0; i < _n_connections; i++)
     168             :   {
     169        8564 :     processor_id_type owner_proc = _processor_ids[i];
     170        8564 :     comm().broadcast(_elem_ids[i], owner_proc, true);
     171        8564 :     comm().broadcast(_local_side_ids[i], owner_proc, true);
     172        8564 :     comm().broadcast(_solutions[i], owner_proc, true);
     173        8564 :     comm().broadcast(_directions[i], owner_proc, true);
     174        8564 :     comm().broadcast(_areas[i], owner_proc, true);
     175        8564 :     comm().broadcast(_stored_p[i], owner_proc, true);
     176             :   }
     177             : 
     178        4282 :   const Real & n1_dot_d1 = _normal[0];
     179        4282 :   const Real d1_dot_d2 = _directions[0] * _directions[1];
     180             : 
     181             :   const ADReal & A1 = _areas[0];
     182             :   const ADReal & A2 = _areas[1];
     183        4282 :   const ADReal A = std::min(A1, A2);
     184        4282 :   const ADReal A_flow = _f_open * A;
     185             : 
     186        4282 :   if (_f_open > _f_open_min)
     187             :   {
     188             :     // compute flow contribution
     189        2126 :     std::vector<ADReal> U_flow1 = _solutions[0];
     190        2126 :     std::vector<ADReal> U_flow2 = _solutions[1];
     191        2126 :     U_flow1[THMVACE1D::AREA] = A_flow;
     192        2126 :     U_flow2[THMVACE1D::AREA] = A_flow;
     193        8504 :     for (unsigned int i = 0; i < THMVACE1D::N_FLUX_OUTPUTS; i++)
     194             :     {
     195        6378 :       U_flow1[i] *= A_flow / A1;
     196        6378 :       U_flow2[i] *= A_flow / A2;
     197             :     }
     198        2126 :     U_flow2[THMVACE1D::RHOUA] *= d1_dot_d2;
     199             : 
     200        2126 :     _fluxes[0] = _numerical_flux.getFlux(
     201        2126 :         _local_side_ids[0], _elem_ids[0], true, U_flow1, U_flow2, n1_dot_d1);
     202             : 
     203        2126 :     _fluxes[1] = _numerical_flux.getFlux(
     204        2126 :         _local_side_ids[0], _elem_ids[0], false, U_flow1, U_flow2, n1_dot_d1);
     205        2126 :     _fluxes[1][THMVACE1D::RHOA] *= d1_dot_d2;
     206        2126 :     _fluxes[1][THMVACE1D::RHOEA] *= d1_dot_d2;
     207        2126 :   }
     208             : 
     209             :   // compute wall contribution
     210       12846 :   for (unsigned int c = 0; c < _n_connections; c++)
     211             :   {
     212        8564 :     const ADReal A_wall = _areas[c] - A_flow;
     213        8564 :     _fluxes[c][THMVACE1D::RHOUA] += _stored_p[c] * A_wall;
     214             :   }
     215        4282 : }
     216             : 
     217             : const std::vector<ADReal> &
     218       17076 : ADGateValve1PhaseUserObject::getFlux(const unsigned int & connection_index) const
     219             : {
     220             :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     221             : 
     222       17076 :   return _fluxes[connection_index];
     223             : }