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 "ADJunctionParallelChannels1PhaseUserObject.h"
      11             : #include "SinglePhaseFluidProperties.h"
      12             : #include "THMIndicesVACE.h"
      13             : #include "VolumeJunction1Phase.h"
      14             : #include "NumericalFlux3EqnBase.h"
      15             : #include "Numerics.h"
      16             : #include "LoggingInterface.h"
      17             : #include "metaphysicl/parallel_numberarray.h"
      18             : #include "metaphysicl/parallel_dualnumber.h"
      19             : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
      20             : #include "libmesh/parallel_algebra.h"
      21             : 
      22             : registerMooseObject("ThermalHydraulicsApp", ADJunctionParallelChannels1PhaseUserObject);
      23             : 
      24             : InputParameters
      25         651 : ADJunctionParallelChannels1PhaseUserObject::validParams()
      26             : {
      27         651 :   InputParameters params = ADVolumeJunction1PhaseUserObject::validParams();
      28        1302 :   params.addRequiredParam<std::string>("component_name", "Name of the associated component");
      29        1302 :   params.addRequiredParam<RealVectorValue>("dir_c0", "Direction of the first connection");
      30             : 
      31         651 :   params.addClassDescription("Computes and caches flux and residual vectors for a 1-phase junction "
      32             :                              "that connects flow channels that are parallel");
      33             : 
      34         651 :   return params;
      35           0 : }
      36             : 
      37         354 : ADJunctionParallelChannels1PhaseUserObject::ADJunctionParallelChannels1PhaseUserObject(
      38         354 :     const InputParameters & params)
      39             :   : ADVolumeJunction1PhaseUserObject(params),
      40             : 
      41         354 :     _p(getADMaterialProperty<Real>("p")),
      42             : 
      43         708 :     _dir_c0(getParam<RealVectorValue>("dir_c0")),
      44             : 
      45         354 :     _stored_pA(_n_connections),
      46         354 :     _areas(_n_connections),
      47         354 :     _is_inlet(_n_connections),
      48             : 
      49        1062 :     _component_name(getParam<std::string>("component_name"))
      50             : {
      51         354 : }
      52             : 
      53             : void
      54       17288 : ADJunctionParallelChannels1PhaseUserObject::initialize()
      55             : {
      56       17288 :   ADVolumeJunction1PhaseUserObject::initialize();
      57             : 
      58             :   _c_in.clear();
      59             :   _c_out.clear();
      60       17288 : }
      61             : 
      62             : void
      63       21110 : ADJunctionParallelChannels1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
      64             : {
      65       21110 :   ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(c);
      66             : 
      67       21110 :   const Real din = _normal[c];
      68       21110 :   const Point di = _dir[0];
      69             :   const Point ni = di * din;
      70             : 
      71             :   const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
      72             :   const auto & rhouV = _cached_junction_var_values[VolumeJunction1Phase::RHOUV_INDEX];
      73             :   const auto & rhovV = _cached_junction_var_values[VolumeJunction1Phase::RHOVV_INDEX];
      74             :   const auto & rhowV = _cached_junction_var_values[VolumeJunction1Phase::RHOWV_INDEX];
      75             :   const auto & rhoEV = _cached_junction_var_values[VolumeJunction1Phase::RHOEV_INDEX];
      76             : 
      77       21110 :   const ADReal vJ = THM::v_from_rhoA_A(rhoV, _volume);
      78             :   const ADRealVectorValue rhouV_vec(rhouV, rhovV, rhowV);
      79       21110 :   const ADReal rhouV2 = rhouV_vec * rhouV_vec;
      80       21110 :   const ADReal eJ = rhoEV / rhoV - 0.5 * rhouV2 / (rhoV * rhoV);
      81       21110 :   const ADReal pJ = _fp.p_from_v_e(vJ, eJ);
      82             : 
      83       42220 :   _residual[VolumeJunction1Phase::RHOUV_INDEX] -= pJ * ni(0) * _A[0];
      84       42220 :   _residual[VolumeJunction1Phase::RHOVV_INDEX] -= pJ * ni(1) * _A[0];
      85       42220 :   _residual[VolumeJunction1Phase::RHOWV_INDEX] -= pJ * ni(2) * _A[0];
      86             : 
      87       21110 :   if (c == 0)
      88             :   {
      89        9641 :     if (MooseUtils::absoluteFuzzyEqual(_rhouA[0], 0))
      90           0 :       _d_flow = _dir[0];
      91             :     else
      92             :       // FIXME: _d_flow should be again ADRealVectorValue when we have parallel comm on
      93             :       // ADRealVectorValue
      94        9641 :       _d_flow =
      95        9641 :           _dir[0] * MetaPhysicL::raw_value(_rhouA[0]) / std::abs(MetaPhysicL::raw_value(_rhouA[0]));
      96             :   }
      97             : 
      98       21110 :   if (!THM::areParallelVectors(_dir_c0, _dir[0]))
      99           2 :     mooseError(_component_name,
     100             :                ": Connected flow channels are not parallel, use VolumeJunction1Phase "
     101             :                "component instead.");
     102             : 
     103       21108 :   _areas[c] = _A[0];
     104       21108 :   _is_inlet[c] = THM::isOutlet(_rhouA[0], _normal[c]);
     105       42216 :   _stored_pA[c] = _p[0] * _A[0];
     106       21108 : }
     107             : 
     108             : void
     109        3576 : ADJunctionParallelChannels1PhaseUserObject::threadJoin(const UserObject & uo)
     110             : {
     111        3576 :   ADVolumeJunction1PhaseUserObject::threadJoin(uo);
     112             : 
     113             :   const auto & jpc_uo = static_cast<const ADJunctionParallelChannels1PhaseUserObject &>(uo);
     114             : 
     115             :   // Store the data computed/retrieved in the other threads
     116        5970 :   for (unsigned int i = 0; i < jpc_uo._connection_indices.size(); i++)
     117             :   {
     118        2394 :     const unsigned int c = jpc_uo._connection_indices[i];
     119             : 
     120        2394 :     _areas[c] = jpc_uo._areas[c];
     121        2394 :     _is_inlet[c] = jpc_uo._is_inlet[c];
     122        2394 :     _stored_pA[c] = jpc_uo._stored_pA[c];
     123        2394 :     if (c == 0)
     124         892 :       _d_flow = jpc_uo._d_flow;
     125             :   }
     126        3576 : }
     127             : 
     128             : void
     129       13709 : ADJunctionParallelChannels1PhaseUserObject::finalize()
     130             : {
     131       44091 :   for (unsigned int i = 0; i < _n_connections; i++)
     132             :   {
     133       30382 :     processor_id_type owner_proc = _processor_ids[i];
     134       30382 :     comm().broadcast(_areas[i], owner_proc, true);
     135       30382 :     comm().broadcast(_stored_pA[i], owner_proc, true);
     136             :     // because std::vector<bool> is very special
     137       30382 :     bool b = _is_inlet[i];
     138       30382 :     comm().broadcast(b, owner_proc, true);
     139       30382 :     _is_inlet[i] = b;
     140             :   }
     141       13709 :   comm().broadcast(_d_flow, _processor_ids[0], true);
     142             : 
     143       13709 :   ADReal Ain_total = 0;
     144       13709 :   ADReal pAin_total = 0;
     145       13709 :   ADReal pAout_total = 0;
     146       13709 :   ADReal Aout_total = 0;
     147       44091 :   for (unsigned int c = 0; c < _n_connections; c++)
     148             :   {
     149       30382 :     if (_is_inlet[c] == true)
     150             :     {
     151       15191 :       Ain_total += _areas[c];
     152       15191 :       pAin_total += _stored_pA[c];
     153       15191 :       _c_in.push_back(c);
     154             :     }
     155             :     else
     156             :     {
     157       15191 :       Aout_total += _areas[c];
     158       15191 :       pAout_total += _stored_pA[c];
     159       15191 :       _c_out.push_back(c);
     160             :     }
     161             :   }
     162             : 
     163       13709 :   ADReal p_wall = 0;
     164       13709 :   ADReal A_wall = 0;
     165             :   RealVectorValue d_wall;
     166             : 
     167       13709 :   if (Aout_total > Ain_total)
     168             :   {
     169        5314 :     if (Aout_total > 1e-15)
     170             :     {
     171        5314 :       p_wall = pAout_total / Aout_total;
     172        5314 :       _c_wall = _c_out;
     173             :     }
     174             :     else
     175             :     {
     176           0 :       p_wall = 0;
     177             :     }
     178        5314 :     A_wall = Aout_total - Ain_total;
     179        5314 :     d_wall = _d_flow;
     180             :   }
     181             :   else
     182             :   {
     183        8395 :     if (Ain_total > 1e-15)
     184             :     {
     185        8395 :       p_wall = pAin_total / Ain_total;
     186        8395 :       _c_wall = _c_in;
     187             :     }
     188             :     else
     189             :     {
     190           0 :       p_wall = 0;
     191             :     }
     192        8395 :     A_wall = Ain_total - Aout_total;
     193        8395 :     d_wall = -_d_flow;
     194             :   }
     195             : 
     196       82254 :   for (unsigned int i = 0; i < _n_scalar_eq; i++)
     197       68545 :     comm().sum(_residual[i]);
     198             : 
     199       13709 :   _residual[VolumeJunction1Phase::RHOUV_INDEX] -= d_wall(0) * p_wall * A_wall;
     200       13709 :   _residual[VolumeJunction1Phase::RHOVV_INDEX] -= d_wall(1) * p_wall * A_wall;
     201       13709 :   _residual[VolumeJunction1Phase::RHOWV_INDEX] -= d_wall(2) * p_wall * A_wall;
     202       13709 : }