LCOV - code coverage report
Current view: top level - src/meshgenerators - SCMQuadAssemblyMeshGenerator.C (source / functions) Hit Total Coverage
Test: idaholab/moose subchannel: #33359 (04c914) with base 7b3324 Lines: 320 333 96.1 %
Date: 2026-07-16 14:35:07 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

          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 "SCMQuadAssemblyMeshGenerator.h"
      11             : #include "SubChannelMesh.h"
      12             : 
      13             : #include "libmesh/edge_edge2.h"
      14             : 
      15             : #include <algorithm>
      16             : #include <cmath>
      17             : #include <limits>
      18             : #include <memory>
      19             : #include <numeric>
      20             : 
      21             : registerMooseObject("SubChannelApp", SCMQuadAssemblyMeshGenerator);
      22             : registerMooseObjectRenamed("SubChannelApp",
      23             :                            SCMQuadSubChannelMeshGenerator,
      24             :                            "06/30/2027 24:00",
      25             :                            SCMQuadAssemblyMeshGenerator);
      26             : registerMooseObjectRenamed("SubChannelApp",
      27             :                            QuadSubChannelMeshGenerator,
      28             :                            "06/30/2027 24:00",
      29             :                            SCMQuadAssemblyMeshGenerator);
      30             : registerMooseObjectRenamed("SubChannelApp",
      31             :                            SCMQuadPinMeshGenerator,
      32             :                            "06/30/2027 24:00",
      33             :                            SCMQuadAssemblyMeshGenerator);
      34             : registerMooseObjectRenamed("SubChannelApp",
      35             :                            QuadPinMeshGenerator,
      36             :                            "06/30/2027 24:00",
      37             :                            SCMQuadAssemblyMeshGenerator);
      38             : 
      39             : InputParameters
      40         428 : SCMQuadAssemblyMeshGenerator::validParams()
      41             : {
      42         428 :   InputParameters params = MeshGenerator::validParams();
      43         428 :   params.addClassDescription("Creates one mesh containing both 1D subchannels and 1D pins in a "
      44             :                              "square lattice arrangement");
      45             : 
      46         856 :   params.addRequiredParam<Real>("pitch", "Pitch [m]");
      47         856 :   params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
      48         856 :   params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
      49         856 :   params.addRequiredParam<Real>("heated_length", "Heated length [m]");
      50         856 :   params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
      51             : 
      52         856 :   params.addParam<std::vector<Real>>(
      53             :       "spacer_z", {}, "Axial location of spacers/vanes/mixing vanes [m]");
      54         856 :   params.addParam<std::vector<Real>>(
      55             :       "spacer_k", {}, "K-loss coefficient of spacers/vanes/mixing vanes [-]");
      56             : 
      57         856 :   params.addParam<std::vector<Real>>("z_blockage",
      58         856 :                                      std::vector<Real>({0.0, 0.0}),
      59             :                                      "Axial location of blockage (inlet, outlet) [m]");
      60         856 :   params.addParam<std::vector<unsigned int>>("index_blockage",
      61         856 :                                              std::vector<unsigned int>({0}),
      62             :                                              "Index of subchannels affected by blockage");
      63         856 :   params.addParam<std::vector<Real>>("reduction_blockage",
      64         856 :                                      std::vector<Real>({1.0}),
      65             :                                      "Area reduction of subchannels affected by blockage");
      66         856 :   params.addParam<std::vector<Real>>(
      67         856 :       "k_blockage", std::vector<Real>({0.0}), "Form loss coefficient of blocked subchannels");
      68             : 
      69         856 :   params.addParam<Real>("Kij", 0.5, "Lateral form loss coefficient [-]");
      70         856 :   params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
      71         856 :   params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
      72         856 :   params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");
      73             : 
      74         856 :   params.addRequiredParam<Real>(
      75             :       "side_gap",
      76             :       "The side gap, not to be confused with the gap between pins; this refers to the gap next "
      77             :       "to the duct or else the distance between the subchannel centroid and the duct wall. "
      78             :       "distance(edge pin center, duct wall) = pitch / 2 + side_gap [m]");
      79             : 
      80         856 :   params.addParam<unsigned int>("subchannel_block_id", 0, "Subchannel block id");
      81         856 :   params.addParam<unsigned int>("pin_block_id", 1, "Fuel Pin block id");
      82             : 
      83         428 :   return params;
      84           0 : }
      85             : 
      86         215 : SCMQuadAssemblyMeshGenerator::SCMQuadAssemblyMeshGenerator(const InputParameters & params)
      87             :   : MeshGenerator(params),
      88         215 :     _unheated_length_entry(getParam<Real>("unheated_length_entry")),
      89         430 :     _heated_length(getParam<Real>("heated_length")),
      90         430 :     _unheated_length_exit(getParam<Real>("unheated_length_exit")),
      91         430 :     _spacer_z(getParam<std::vector<Real>>("spacer_z")),
      92         430 :     _spacer_k(getParam<std::vector<Real>>("spacer_k")),
      93         430 :     _z_blockage(getParam<std::vector<Real>>("z_blockage")),
      94         430 :     _index_blockage(getParam<std::vector<unsigned int>>("index_blockage")),
      95         430 :     _reduction_blockage(getParam<std::vector<Real>>("reduction_blockage")),
      96         430 :     _k_blockage(getParam<std::vector<Real>>("k_blockage")),
      97         430 :     _kij(getParam<Real>("Kij")),
      98         430 :     _pitch(getParam<Real>("pitch")),
      99         430 :     _pin_diameter(getParam<Real>("pin_diameter")),
     100         430 :     _n_cells(getParam<unsigned int>("n_cells")),
     101         430 :     _nx(getParam<unsigned int>("nx")),
     102         430 :     _ny(getParam<unsigned int>("ny")),
     103         215 :     _n_channels(0),
     104         215 :     _n_gaps(0),
     105         215 :     _n_pins(0),
     106         430 :     _side_gap(getParam<Real>("side_gap")),
     107         430 :     _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
     108         645 :     _pin_block_id(getParam<unsigned int>("pin_block_id"))
     109             : {
     110         215 :   const Real total_length = _unheated_length_entry + _heated_length + _unheated_length_exit;
     111             : 
     112         215 :   if (_n_cells == 0)
     113           0 :     paramError("n_cells", "The number of axial cells must be greater than zero");
     114             : 
     115         215 :   if (total_length <= 0.0)
     116           0 :     mooseError("Total bundle length must be greater than zero");
     117             : 
     118         215 :   if (_nx == 0 || _ny == 0)
     119           0 :     mooseError("The number of subchannels must be greater than zero in each direction");
     120             : 
     121         215 :   if (_nx < 2 && _ny < 2)
     122           2 :     mooseError("The number of subchannels cannot be less than 2 in both directions. "
     123             :                "Smallest assembly allowed is either 2X1 or 1X2.");
     124             : 
     125         213 :   _n_channels = _nx * _ny;
     126         213 :   _n_gaps = (_nx - 1) * _ny + (_ny - 1) * _nx;
     127         213 :   _n_pins = (_nx - 1) * (_ny - 1);
     128             : 
     129         213 :   if (_spacer_z.size() != _spacer_k.size())
     130           0 :     mooseError("Size of vector spacer_z should equal size of spacer_k");
     131             : 
     132         615 :   for (const auto spacer_z : _spacer_z)
     133         402 :     if (spacer_z < 0.0 || spacer_z > total_length)
     134           0 :       paramError("spacer_z", "Spacer locations must be between zero and total bundle length");
     135             : 
     136         213 :   if (_z_blockage.size() != 2)
     137           0 :     paramError("z_blockage", "Size of vector z_blockage must be 2");
     138             : 
     139         213 :   if (_z_blockage.front() > _z_blockage.back())
     140           0 :     paramError("z_blockage", "z_blockage inlet location must not exceed outlet location");
     141             : 
     142         213 :   if (!_index_blockage.empty() &&
     143         213 :       *std::max_element(_index_blockage.begin(), _index_blockage.end()) > (_n_channels - 1))
     144           0 :     paramError("index_blockage", "Blocked subchannel index exceeds valid subchannel range");
     145             : 
     146         213 :   if (!_reduction_blockage.empty() &&
     147         213 :       *std::max_element(_reduction_blockage.begin(), _reduction_blockage.end()) > 1.0)
     148           0 :     paramError("reduction_blockage", "Area reduction of blocked subchannels cannot exceed 1");
     149             : 
     150         213 :   if ((_index_blockage.size() > _n_channels) || (_reduction_blockage.size() > _n_channels) ||
     151             :       (_k_blockage.size() > _n_channels))
     152           0 :     mooseError("Sizes of blockage-related vectors cannot exceed total number of subchannels");
     153             : 
     154         213 :   if ((_index_blockage.size() != _reduction_blockage.size()) ||
     155         426 :       (_index_blockage.size() != _k_blockage.size()) ||
     156             :       (_reduction_blockage.size() != _k_blockage.size()))
     157           0 :     mooseError("index_blockage, reduction_blockage, and k_blockage must have equal size");
     158             : 
     159         213 :   SubChannelMesh::generateZGrid(
     160         213 :       _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);
     161             : 
     162         213 :   initializeChannelData();
     163         213 : }
     164             : 
     165             : void
     166         213 : SCMQuadAssemblyMeshGenerator::initializeChannelData()
     167             : {
     168             :   // Defining the total length from 3 axial sections
     169         213 :   const Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
     170             : 
     171             :   // Defining the position of the spacer grid in the numerical solution array
     172             :   std::vector<int> spacer_cell;
     173         615 :   for (const auto & elem : _spacer_z)
     174         402 :     spacer_cell.emplace_back(std::round(elem * _n_cells / L));
     175             : 
     176             :   // Defining the arrays for axial resistances
     177         213 :   std::vector<Real> kgrid(_n_cells + 1, 0.0);
     178         213 :   _k_grid.resize(_n_channels, std::vector<Real>(_n_cells + 1));
     179             : 
     180             :   // Summing the spacer resistance to the 1D grid resistance array
     181         615 :   for (unsigned int index = 0; index < spacer_cell.size(); index++)
     182         402 :     kgrid[spacer_cell[index]] += _spacer_k[index];
     183             : 
     184             :   // Creating the 2D grid resistance array
     185        3837 :   for (unsigned int i = 0; i < _n_channels; i++)
     186        3624 :     _k_grid[i] = kgrid;
     187             : 
     188             :   // Add blockage resistance to the 2D grid resistance array
     189         213 :   const Real dz = L / _n_cells;
     190        4410 :   for (unsigned int i = 0; i < _n_cells + 1; i++)
     191        4197 :     if ((dz * i >= _z_blockage.front() && dz * i <= _z_blockage.back()))
     192             :     {
     193             :       unsigned int index = 0;
     194         426 :       for (const auto & i_ch : _index_blockage)
     195             :       {
     196         213 :         _k_grid[i_ch][i] += _k_blockage[index];
     197         213 :         index++;
     198             :       }
     199             :     }
     200             : 
     201             :   // Defining the size of the maps
     202         213 :   _gap_to_chan_map.resize(_n_gaps);
     203         213 :   _gap_to_pin_map.resize(_n_gaps);
     204         213 :   _gapnodes.resize(_n_gaps);
     205         213 :   _chan_to_gap_map.resize(_n_channels);
     206         213 :   _chan_to_pin_map.resize(_n_channels);
     207         213 :   _pin_to_chan_map.resize(_n_pins);
     208         213 :   _sign_id_crossflow_map.resize(_n_channels);
     209         213 :   _gij_map.resize(_n_cells + 1);
     210         213 :   _subchannel_position.resize(_n_channels);
     211             : 
     212        3837 :   for (unsigned int i = 0; i < _n_channels; i++)
     213             :   {
     214        3624 :     _subchannel_position[i].reserve(3);
     215       14496 :     for (unsigned int j = 0; j < 3; j++)
     216       10872 :       _subchannel_position.at(i).push_back(0.0);
     217             :   }
     218             : 
     219        4410 :   for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     220        4197 :     _gij_map[iz].reserve(_n_gaps);
     221             : 
     222             :   // Defining the signs for positive and negative flows
     223         213 :   const Real positive_flow = 1.0;
     224         213 :   const Real negative_flow = -1.0;
     225             : 
     226             :   // Defining the subchannel types
     227         213 :   _subch_type.resize(_n_channels);
     228        1005 :   for (unsigned int iy = 0; iy < _ny; iy++)
     229        4416 :     for (unsigned int ix = 0; ix < _nx; ix++)
     230             :     {
     231        3624 :       const unsigned int i_ch = _nx * iy + ix;
     232        3411 :       const bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
     233        6825 :                              (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
     234        3624 :       const bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);
     235             : 
     236             :       // Two channels side by side (a 1x2/2x1 grid) only occurs in verification cases; treat both
     237             :       // as CENTER. The smallest physical assembly is 2x2 subchannels around a single pin.
     238        3624 :       if (_n_channels == 2)
     239          18 :         _subch_type[i_ch] = EChannelType::CENTER;
     240        3606 :       else if (_n_channels == 4)
     241         104 :         _subch_type[i_ch] = EChannelType::CORNER;
     242        3502 :       else if (is_corner)
     243         712 :         _subch_type[i_ch] = EChannelType::CORNER;
     244        2790 :       else if (is_edge)
     245        1594 :         _subch_type[i_ch] = EChannelType::EDGE;
     246             :       else
     247        1196 :         _subch_type[i_ch] = EChannelType::CENTER;
     248             :     }
     249             : 
     250             :   // Index the east-west gaps.
     251         213 :   unsigned int i_gap = 0;
     252        1005 :   for (unsigned int iy = 0; iy < _ny; iy++)
     253        3624 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     254             :     {
     255        2832 :       const unsigned int i_ch = _nx * iy + ix;
     256        2832 :       const unsigned int j_ch = _nx * iy + (ix + 1);
     257        2832 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     258        2832 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     259        2832 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     260        2832 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     261        2832 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     262             : 
     263             :       // Make a gap size map.
     264        2832 :       if (iy == 0 || iy == _ny - 1)
     265        1264 :         _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
     266             :       else
     267        1568 :         _gij_map[0].push_back(_pitch - _pin_diameter);
     268             : 
     269        2832 :       ++i_gap;
     270             :     }
     271             : 
     272             :   // Index the north-south gaps.
     273         792 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     274        3355 :     for (unsigned int ix = 0; ix < _nx; ix++)
     275             :     {
     276        2776 :       const unsigned int i_ch = _nx * iy + ix;
     277        2776 :       const unsigned int j_ch = _nx * (iy + 1) + ix;
     278        2776 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     279        2776 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     280        2776 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     281        2776 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     282        2776 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     283             : 
     284             :       // Make a gap size map.
     285        2776 :       if (ix == 0 || ix == _nx - 1)
     286        1155 :         _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
     287             :       else
     288        1621 :         _gij_map[0].push_back(_pitch - _pin_diameter);
     289             : 
     290        2776 :       ++i_gap;
     291             :     }
     292             : 
     293        4197 :   for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
     294        3984 :     _gij_map[iz] = _gij_map[0];
     295             : 
     296             :   // Make pin to channel map.
     297         792 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     298        2776 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     299             :     {
     300        2197 :       const unsigned int i_pin = (_nx - 1) * iy + ix;
     301        2197 :       const unsigned int i_chan_1 = _nx * iy + ix;
     302        2197 :       const unsigned int i_chan_2 = _nx * (iy + 1) + ix;
     303        2197 :       const unsigned int i_chan_3 = _nx * (iy + 1) + (ix + 1);
     304        2197 :       const unsigned int i_chan_4 = _nx * iy + (ix + 1);
     305             : 
     306        2197 :       _pin_to_chan_map[i_pin].push_back(i_chan_1);
     307        2197 :       _pin_to_chan_map[i_pin].push_back(i_chan_2);
     308        2197 :       _pin_to_chan_map[i_pin].push_back(i_chan_3);
     309        2197 :       _pin_to_chan_map[i_pin].push_back(i_chan_4);
     310             :     }
     311             : 
     312             :   // Set the subchannel positions so that the center of the assembly is the zero point.
     313        1005 :   for (unsigned int iy = 0; iy < _ny; iy++)
     314        4416 :     for (unsigned int ix = 0; ix < _nx; ix++)
     315             :     {
     316        3624 :       const unsigned int i_ch = _nx * iy + ix;
     317        3624 :       const Real offset_x = (_nx - 1) * _pitch / 2.0;
     318        3624 :       const Real offset_y = (_ny - 1) * _pitch / 2.0;
     319        3624 :       _subchannel_position[i_ch][0] = _pitch * ix - offset_x;
     320        3624 :       _subchannel_position[i_ch][1] = _pitch * iy - offset_y;
     321             :     }
     322             : 
     323         213 :   if (_n_pins > 0)
     324             :   {
     325             :     // Make channel to pin map.
     326         984 :     for (unsigned int iy = 0; iy < _ny; iy++) // row
     327        4386 :       for (unsigned int ix = 0; ix < _nx; ix++)
     328             :       {
     329        3606 :         const unsigned int i_ch = _nx * iy + ix;
     330             : 
     331             :         // Corners contact 1/4 of one pin.
     332        3606 :         if (iy == 0 && ix == 0)
     333         204 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     334        3402 :         else if (iy == _ny - 1 && ix == 0)
     335         204 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     336        3198 :         else if (iy == 0 && ix == _nx - 1)
     337         204 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     338        2994 :         else if (iy == _ny - 1 && ix == _nx - 1)
     339         204 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     340             :         // Sides contact 1/4 of two pins.
     341        2790 :         else if (iy == 0)
     342             :         {
     343         425 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     344         425 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     345             :         }
     346        2365 :         else if (iy == _ny - 1)
     347             :         {
     348         425 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     349         425 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     350             :         }
     351        1940 :         else if (ix == 0)
     352             :         {
     353         372 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     354         372 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     355             :         }
     356        1568 :         else if (ix == _nx - 1)
     357             :         {
     358         372 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     359         372 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     360             :         }
     361             :         // Interior channels contact 1/4 of four pins.
     362             :         else
     363             :         {
     364        1196 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     365        1196 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     366        1196 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     367        1196 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     368             :         }
     369             :       }
     370             : 
     371             :     // Make gap to pin map.
     372        5803 :     for (unsigned int ig = 0; ig < _n_gaps; ig++)
     373             :     {
     374        5599 :       const auto i_ch = _gap_to_chan_map[ig].first;
     375        5599 :       const auto j_ch = _gap_to_chan_map[ig].second;
     376        5599 :       const auto & i_pins = _chan_to_pin_map[i_ch];
     377        5599 :       const auto & j_pins = _chan_to_pin_map[j_ch];
     378             : 
     379             :       // Initialize with default values.
     380             :       _gap_to_pin_map[ig] = {10000, 10000};
     381             : 
     382       20765 :       for (unsigned int i : i_pins)
     383       59918 :         for (unsigned int j : j_pins)
     384       44752 :           if (i == j)
     385             :           {
     386        8788 :             if (_gap_to_pin_map[ig].first == 10000)
     387             :             {
     388        5599 :               _gap_to_pin_map[ig].first = i;
     389        5599 :               _gap_to_pin_map[ig].second = i;
     390             :             }
     391             :             else
     392        3189 :               _gap_to_pin_map[ig].second = i;
     393             :           }
     394             :     }
     395             :   }
     396             : 
     397             :   // Reduce reserved memory in the channel-to-gap map.
     398        3837 :   for (auto & gap : _chan_to_gap_map)
     399             :     gap.shrink_to_fit();
     400             : 
     401             :   // Reduce reserved memory in the channel-to-pin map.
     402        3837 :   for (auto & pin : _chan_to_pin_map)
     403             :     pin.shrink_to_fit();
     404             : 
     405             :   // Reduce reserved memory in the pin-to-channel map.
     406        2410 :   for (auto & pin : _pin_to_chan_map)
     407             :     pin.shrink_to_fit();
     408         213 : }
     409             : 
     410             : void
     411         213 : SCMQuadAssemblyMeshGenerator::buildSubchannelMesh(MeshBase & mesh_base,
     412             :                                                   BoundaryInfo & boundary_info)
     413             : {
     414         213 :   mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * _ny * _nx);
     415         213 :   mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * _ny * _nx);
     416             : 
     417         213 :   _nodes.resize(_nx * _ny);
     418             : 
     419         213 :   const Real offset_x = (_nx - 1) * _pitch / 2.0;
     420         213 :   const Real offset_y = (_ny - 1) * _pitch / 2.0;
     421             : 
     422             :   // Add the points in the shape of a rectilinear grid. The grid is regular on the xy-plane with a
     423             :   // spacing of `pitch` between points. The grid along z is irregular to account for pin spacers.
     424             :   // Store pointers in the _nodes array so we can keep track of which points are in which channels.
     425         213 :   dof_id_type node_id = mesh_base.n_nodes();
     426        1005 :   for (unsigned int iy = 0; iy < _ny; iy++)
     427        4416 :     for (unsigned int ix = 0; ix < _nx; ix++)
     428             :     {
     429        3624 :       const unsigned int i_ch = _nx * iy + ix;
     430        3624 :       _nodes[i_ch].reserve(_n_cells + 1);
     431             : 
     432       70077 :       for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     433       66453 :         _nodes[i_ch].push_back(mesh_base.add_point(
     434      132906 :             Point(_pitch * ix - offset_x, _pitch * iy - offset_y, _z_grid[iz]), node_id++));
     435             :     }
     436             : 
     437             :   // Add the elements which in this case are 2-node edges that link each subchannel's nodes
     438             :   // vertically.
     439         213 :   dof_id_type elem_id = mesh_base.n_elem();
     440        1005 :   for (unsigned int iy = 0; iy < _ny; iy++)
     441        4416 :     for (unsigned int ix = 0; ix < _nx; ix++)
     442       66453 :       for (unsigned int iz = 0; iz < _n_cells; iz++)
     443             :       {
     444       62829 :         Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
     445       62829 :         elem->subdomain_id() = _subchannel_block_id;
     446       62829 :         elem->set_id(elem_id++);
     447             : 
     448       62829 :         const unsigned int i_ch = _nx * iy + ix;
     449       62829 :         elem->set_node(0, _nodes[i_ch][iz]);
     450       62829 :         elem->set_node(1, _nodes[i_ch][iz + 1]);
     451             : 
     452       62829 :         if (iz == 0)
     453        3624 :           boundary_info.add_side(elem, 0, 0);
     454       62829 :         if (iz == _n_cells - 1)
     455        3624 :           boundary_info.add_side(elem, 1, 1);
     456             :       }
     457             : 
     458         213 :   mesh_base.subdomain_name(_subchannel_block_id) = "subchannel";
     459         213 : }
     460             : 
     461             : void
     462         204 : SCMQuadAssemblyMeshGenerator::buildPinMesh(MeshBase & mesh_base)
     463             : {
     464         204 :   mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * (_ny - 1) * (_nx - 1));
     465         204 :   mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * (_ny - 1) * (_nx - 1));
     466             : 
     467         204 :   _pin_nodes.resize((_nx - 1) * (_ny - 1));
     468             : 
     469             :   // Add the points in the shape of a rectilinear grid. The grid is regular on the xy-plane with a
     470             :   // spacing of `pitch` between points. The grid along z is also regular. Store pointers in the
     471             :   // _pin_nodes array so we can keep track of which points are in which pins.
     472         204 :   const Real offset_x = (_nx - 2) * _pitch / 2.0;
     473         204 :   const Real offset_y = (_ny - 2) * _pitch / 2.0;
     474             : 
     475         204 :   dof_id_type node_id = mesh_base.n_nodes();
     476         780 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     477        2773 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     478             :     {
     479        2197 :       const unsigned int i_pin = (_nx - 1) * iy + ix;
     480        2197 :       _pin_nodes[i_pin].reserve(_n_cells + 1);
     481             : 
     482       42135 :       for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     483       39938 :         _pin_nodes[i_pin].push_back(mesh_base.add_point(
     484       79876 :             Point(_pitch * ix - offset_x, _pitch * iy - offset_y, _z_grid[iz]), node_id++));
     485             :     }
     486             : 
     487             :   // Add the elements which in this case are 2-node edges that link each pin's nodes vertically.
     488         204 :   dof_id_type elem_id = mesh_base.n_elem();
     489         780 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     490        2773 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     491       39938 :       for (unsigned int iz = 0; iz < _n_cells; iz++)
     492             :       {
     493       37741 :         Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
     494       37741 :         elem->subdomain_id() = _pin_block_id;
     495       37741 :         elem->set_id(elem_id++);
     496             : 
     497       37741 :         const unsigned int i_pin = (_nx - 1) * iy + ix;
     498       37741 :         elem->set_node(0, _pin_nodes[i_pin][iz]);
     499       37741 :         elem->set_node(1, _pin_nodes[i_pin][iz + 1]);
     500             :       }
     501             : 
     502         204 :   mesh_base.subdomain_name(_pin_block_id) = "fuel_pins";
     503         204 : }
     504             : 
     505             : void
     506         213 : SCMQuadAssemblyMeshGenerator::transferMetadata(QuadSubChannelMesh & sch_mesh)
     507             : {
     508             :   // Move the metadata into QuadSubChannelMesh.
     509         213 :   sch_mesh._unheated_length_entry = _unheated_length_entry;
     510         213 :   sch_mesh._heated_length = _heated_length;
     511         213 :   sch_mesh._unheated_length_exit = _unheated_length_exit;
     512         213 :   sch_mesh._z_grid = _z_grid;
     513         213 :   sch_mesh._k_grid = _k_grid;
     514         213 :   sch_mesh._spacer_z = _spacer_z;
     515         213 :   sch_mesh._spacer_k = _spacer_k;
     516         213 :   sch_mesh._z_blockage = _z_blockage;
     517         213 :   sch_mesh._index_blockage = _index_blockage;
     518         213 :   sch_mesh._reduction_blockage = _reduction_blockage;
     519         213 :   sch_mesh._kij = _kij;
     520         213 :   sch_mesh._pitch = _pitch;
     521         213 :   sch_mesh._pin_diameter = _pin_diameter;
     522         213 :   sch_mesh._n_cells = _n_cells;
     523             : 
     524         213 :   sch_mesh._nx = _nx;
     525         213 :   sch_mesh._ny = _ny;
     526         213 :   sch_mesh._n_channels = _n_channels;
     527         213 :   sch_mesh._n_gaps = _n_gaps;
     528         213 :   sch_mesh._n_pins = _n_pins;
     529         213 :   sch_mesh._side_gap = _side_gap;
     530             : 
     531         213 :   sch_mesh._nodes = _nodes;
     532         213 :   sch_mesh._pin_nodes = _pin_nodes;
     533         213 :   sch_mesh._gapnodes = _gapnodes;
     534         213 :   sch_mesh._gap_to_chan_map = _gap_to_chan_map;
     535         213 :   sch_mesh._gap_to_pin_map = _gap_to_pin_map;
     536         213 :   sch_mesh._chan_to_gap_map = _chan_to_gap_map;
     537         213 :   sch_mesh._chan_to_pin_map = _chan_to_pin_map;
     538         213 :   sch_mesh._pin_to_chan_map = _pin_to_chan_map;
     539         213 :   sch_mesh._sign_id_crossflow_map = _sign_id_crossflow_map;
     540         213 :   sch_mesh._gij_map = _gij_map;
     541         213 :   sch_mesh._subchannel_position = _subchannel_position;
     542         213 :   sch_mesh._subch_type = _subch_type;
     543             : 
     544         213 :   sch_mesh._duct_mesh_exist = false;
     545         213 :   sch_mesh._pin_mesh_exist = (_n_pins > 0);
     546         213 : }
     547             : 
     548             : std::unique_ptr<MeshBase>
     549         213 : SCMQuadAssemblyMeshGenerator::generate()
     550             : {
     551         213 :   auto mesh_base = buildMeshBaseObject();
     552             :   BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
     553             : 
     554         213 :   mesh_base->set_spatial_dimension(3);
     555             : 
     556         213 :   buildSubchannelMesh(*mesh_base, boundary_info);
     557             : 
     558         213 :   if (_n_pins > 0)
     559         204 :     buildPinMesh(*mesh_base);
     560             : 
     561         213 :   boundary_info.sideset_name(0) = "inlet";
     562         213 :   boundary_info.sideset_name(1) = "outlet";
     563         213 :   boundary_info.nodeset_name(0) = "inlet";
     564         213 :   boundary_info.nodeset_name(1) = "outlet";
     565             : 
     566         213 :   mesh_base->prepare_for_use();
     567             : 
     568         213 :   auto & sch_mesh = static_cast<QuadSubChannelMesh &>(*_mesh);
     569         213 :   transferMetadata(sch_mesh);
     570         213 :   sch_mesh.computeAssemblyHydraulicParameters();
     571             : 
     572         213 :   return mesh_base;
     573           0 : }

Generated by: LCOV version 1.14