LCOV - code coverage report
Current view: top level - src/meshgenerators - SCMQuadAssemblyMeshGenerator.C (source / functions) Hit Total Coverage
Test: idaholab/moose subchannel: #33187 (5aa0b2) with base d7c4bd Lines: 320 333 96.1 %
Date: 2026-06-30 12:24:57 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         410 : SCMQuadAssemblyMeshGenerator::validParams()
      41             : {
      42         410 :   InputParameters params = MeshGenerator::validParams();
      43         410 :   params.addClassDescription("Creates one mesh containing both 1D subchannels and 1D pins in a "
      44             :                              "square lattice arrangement");
      45             : 
      46         820 :   params.addRequiredParam<Real>("pitch", "Pitch [m]");
      47         820 :   params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
      48         820 :   params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
      49         820 :   params.addRequiredParam<Real>("heated_length", "Heated length [m]");
      50         820 :   params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
      51             : 
      52         820 :   params.addParam<std::vector<Real>>(
      53             :       "spacer_z", {}, "Axial location of spacers/vanes/mixing vanes [m]");
      54         820 :   params.addParam<std::vector<Real>>(
      55             :       "spacer_k", {}, "K-loss coefficient of spacers/vanes/mixing vanes [-]");
      56             : 
      57         820 :   params.addParam<std::vector<Real>>("z_blockage",
      58         820 :                                      std::vector<Real>({0.0, 0.0}),
      59             :                                      "Axial location of blockage (inlet, outlet) [m]");
      60         820 :   params.addParam<std::vector<unsigned int>>("index_blockage",
      61         820 :                                              std::vector<unsigned int>({0}),
      62             :                                              "Index of subchannels affected by blockage");
      63         820 :   params.addParam<std::vector<Real>>("reduction_blockage",
      64         820 :                                      std::vector<Real>({1.0}),
      65             :                                      "Area reduction of subchannels affected by blockage");
      66         820 :   params.addParam<std::vector<Real>>(
      67         820 :       "k_blockage", std::vector<Real>({0.0}), "Form loss coefficient of blocked subchannels");
      68             : 
      69         820 :   params.addParam<Real>("Kij", 0.5, "Lateral form loss coefficient [-]");
      70         820 :   params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
      71         820 :   params.addRequiredParam<unsigned int>("nx", "Number of channels in the x direction [-]");
      72         820 :   params.addRequiredParam<unsigned int>("ny", "Number of channels in the y direction [-]");
      73             : 
      74         820 :   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         820 :   params.addParam<unsigned int>("subchannel_block_id", 0, "Subchannel block id");
      81         820 :   params.addParam<unsigned int>("pin_block_id", 1, "Fuel Pin block id");
      82             : 
      83         410 :   return params;
      84           0 : }
      85             : 
      86         206 : SCMQuadAssemblyMeshGenerator::SCMQuadAssemblyMeshGenerator(const InputParameters & params)
      87             :   : MeshGenerator(params),
      88         206 :     _unheated_length_entry(getParam<Real>("unheated_length_entry")),
      89         412 :     _heated_length(getParam<Real>("heated_length")),
      90         412 :     _unheated_length_exit(getParam<Real>("unheated_length_exit")),
      91         412 :     _spacer_z(getParam<std::vector<Real>>("spacer_z")),
      92         412 :     _spacer_k(getParam<std::vector<Real>>("spacer_k")),
      93         412 :     _z_blockage(getParam<std::vector<Real>>("z_blockage")),
      94         412 :     _index_blockage(getParam<std::vector<unsigned int>>("index_blockage")),
      95         412 :     _reduction_blockage(getParam<std::vector<Real>>("reduction_blockage")),
      96         412 :     _k_blockage(getParam<std::vector<Real>>("k_blockage")),
      97         412 :     _kij(getParam<Real>("Kij")),
      98         412 :     _pitch(getParam<Real>("pitch")),
      99         412 :     _pin_diameter(getParam<Real>("pin_diameter")),
     100         412 :     _n_cells(getParam<unsigned int>("n_cells")),
     101         412 :     _nx(getParam<unsigned int>("nx")),
     102         412 :     _ny(getParam<unsigned int>("ny")),
     103         206 :     _n_channels(0),
     104         206 :     _n_gaps(0),
     105         206 :     _n_pins(0),
     106         412 :     _side_gap(getParam<Real>("side_gap")),
     107         412 :     _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
     108         618 :     _pin_block_id(getParam<unsigned int>("pin_block_id"))
     109             : {
     110         206 :   const Real total_length = _unheated_length_entry + _heated_length + _unheated_length_exit;
     111             : 
     112         206 :   if (_n_cells == 0)
     113           0 :     paramError("n_cells", "The number of axial cells must be greater than zero");
     114             : 
     115         206 :   if (total_length <= 0.0)
     116           0 :     mooseError("Total bundle length must be greater than zero");
     117             : 
     118         206 :   if (_nx == 0 || _ny == 0)
     119           0 :     mooseError("The number of subchannels must be greater than zero in each direction");
     120             : 
     121         206 :   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         204 :   _n_channels = _nx * _ny;
     126         204 :   _n_gaps = (_nx - 1) * _ny + (_ny - 1) * _nx;
     127         204 :   _n_pins = (_nx - 1) * (_ny - 1);
     128             : 
     129         204 :   if (_spacer_z.size() != _spacer_k.size())
     130           0 :     mooseError("Size of vector spacer_z should equal size of spacer_k");
     131             : 
     132         597 :   for (const auto spacer_z : _spacer_z)
     133         393 :     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         204 :   if (_z_blockage.size() != 2)
     137           0 :     paramError("z_blockage", "Size of vector z_blockage must be 2");
     138             : 
     139         204 :   if (_z_blockage.front() > _z_blockage.back())
     140           0 :     paramError("z_blockage", "z_blockage inlet location must not exceed outlet location");
     141             : 
     142         204 :   if (!_index_blockage.empty() &&
     143         204 :       *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         204 :   if (!_reduction_blockage.empty() &&
     147         204 :       *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         204 :   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         204 :   if ((_index_blockage.size() != _reduction_blockage.size()) ||
     155         408 :       (_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         204 :   SubChannelMesh::generateZGrid(
     160         204 :       _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);
     161             : 
     162         204 :   initializeChannelData();
     163         204 : }
     164             : 
     165             : void
     166         204 : SCMQuadAssemblyMeshGenerator::initializeChannelData()
     167             : {
     168             :   // Defining the total length from 3 axial sections
     169         204 :   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         597 :   for (const auto & elem : _spacer_z)
     174         393 :     spacer_cell.emplace_back(std::round(elem * _n_cells / L));
     175             : 
     176             :   // Defining the arrays for axial resistances
     177         204 :   std::vector<Real> kgrid(_n_cells + 1, 0.0);
     178         204 :   _k_grid.resize(_n_channels, std::vector<Real>(_n_cells + 1));
     179             : 
     180             :   // Summing the spacer resistance to the 1D grid resistance array
     181         597 :   for (unsigned int index = 0; index < spacer_cell.size(); index++)
     182         393 :     kgrid[spacer_cell[index]] += _spacer_k[index];
     183             : 
     184             :   // Creating the 2D grid resistance array
     185        3747 :   for (unsigned int i = 0; i < _n_channels; i++)
     186        3543 :     _k_grid[i] = kgrid;
     187             : 
     188             :   // Add blockage resistance to the 2D grid resistance array
     189         204 :   const Real dz = L / _n_cells;
     190        4365 :   for (unsigned int i = 0; i < _n_cells + 1; i++)
     191        4161 :     if ((dz * i >= _z_blockage.front() && dz * i <= _z_blockage.back()))
     192             :     {
     193             :       unsigned int index = 0;
     194         408 :       for (const auto & i_ch : _index_blockage)
     195             :       {
     196         204 :         _k_grid[i_ch][i] += _k_blockage[index];
     197         204 :         index++;
     198             :       }
     199             :     }
     200             : 
     201             :   // Defining the size of the maps
     202         204 :   _gap_to_chan_map.resize(_n_gaps);
     203         204 :   _gap_to_pin_map.resize(_n_gaps);
     204         204 :   _gapnodes.resize(_n_gaps);
     205         204 :   _chan_to_gap_map.resize(_n_channels);
     206         204 :   _chan_to_pin_map.resize(_n_channels);
     207         204 :   _pin_to_chan_map.resize(_n_pins);
     208         204 :   _sign_id_crossflow_map.resize(_n_channels);
     209         204 :   _gij_map.resize(_n_cells + 1);
     210         204 :   _subchannel_position.resize(_n_channels);
     211             : 
     212        3747 :   for (unsigned int i = 0; i < _n_channels; i++)
     213             :   {
     214        3543 :     _subchannel_position[i].reserve(3);
     215       14172 :     for (unsigned int j = 0; j < 3; j++)
     216       10629 :       _subchannel_position.at(i).push_back(0.0);
     217             :   }
     218             : 
     219        4365 :   for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     220        4161 :     _gij_map[iz].reserve(_n_gaps);
     221             : 
     222             :   // Defining the signs for positive and negative flows
     223         204 :   const Real positive_flow = 1.0;
     224         204 :   const Real negative_flow = -1.0;
     225             : 
     226             :   // Defining the subchannel types
     227         204 :   _subch_type.resize(_n_channels);
     228         969 :   for (unsigned int iy = 0; iy < _ny; iy++)
     229        4308 :     for (unsigned int ix = 0; ix < _nx; ix++)
     230             :     {
     231        3543 :       const unsigned int i_ch = _nx * iy + ix;
     232        3339 :       const bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
     233        6681 :                              (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
     234        3543 :       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        3543 :       if (_n_channels == 2)
     239          18 :         _subch_type[i_ch] = EChannelType::CENTER;
     240        3525 :       else if (_n_channels == 4)
     241         104 :         _subch_type[i_ch] = EChannelType::CORNER;
     242        3421 :       else if (is_corner)
     243         676 :         _subch_type[i_ch] = EChannelType::CORNER;
     244        2745 :       else if (is_edge)
     245        1558 :         _subch_type[i_ch] = EChannelType::EDGE;
     246             :       else
     247        1187 :         _subch_type[i_ch] = EChannelType::CENTER;
     248             :     }
     249             : 
     250             :   // Index the east-west gaps.
     251         204 :   unsigned int i_gap = 0;
     252         969 :   for (unsigned int iy = 0; iy < _ny; iy++)
     253        3543 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     254             :     {
     255        2778 :       const unsigned int i_ch = _nx * iy + ix;
     256        2778 :       const unsigned int j_ch = _nx * iy + (ix + 1);
     257        2778 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     258        2778 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     259        2778 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     260        2778 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     261        2778 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     262             : 
     263             :       // Make a gap size map.
     264        2778 :       if (iy == 0 || iy == _ny - 1)
     265        1228 :         _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
     266             :       else
     267        1550 :         _gij_map[0].push_back(_pitch - _pin_diameter);
     268             : 
     269        2778 :       ++i_gap;
     270             :     }
     271             : 
     272             :   // Index the north-south gaps.
     273         765 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     274        3283 :     for (unsigned int ix = 0; ix < _nx; ix++)
     275             :     {
     276        2722 :       const unsigned int i_ch = _nx * iy + ix;
     277        2722 :       const unsigned int j_ch = _nx * (iy + 1) + ix;
     278        2722 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     279        2722 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     280        2722 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     281        2722 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     282        2722 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     283             : 
     284             :       // Make a gap size map.
     285        2722 :       if (ix == 0 || ix == _nx - 1)
     286        1119 :         _gij_map[0].push_back((_pitch - _pin_diameter) / 2.0 + _side_gap);
     287             :       else
     288        1603 :         _gij_map[0].push_back(_pitch - _pin_diameter);
     289             : 
     290        2722 :       ++i_gap;
     291             :     }
     292             : 
     293        4161 :   for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
     294        3957 :     _gij_map[iz] = _gij_map[0];
     295             : 
     296             :   // Make pin to channel map.
     297         765 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     298        2722 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     299             :     {
     300        2161 :       const unsigned int i_pin = (_nx - 1) * iy + ix;
     301        2161 :       const unsigned int i_chan_1 = _nx * iy + ix;
     302        2161 :       const unsigned int i_chan_2 = _nx * (iy + 1) + ix;
     303        2161 :       const unsigned int i_chan_3 = _nx * (iy + 1) + (ix + 1);
     304        2161 :       const unsigned int i_chan_4 = _nx * iy + (ix + 1);
     305             : 
     306        2161 :       _pin_to_chan_map[i_pin].push_back(i_chan_1);
     307        2161 :       _pin_to_chan_map[i_pin].push_back(i_chan_2);
     308        2161 :       _pin_to_chan_map[i_pin].push_back(i_chan_3);
     309        2161 :       _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         969 :   for (unsigned int iy = 0; iy < _ny; iy++)
     314        4308 :     for (unsigned int ix = 0; ix < _nx; ix++)
     315             :     {
     316        3543 :       const unsigned int i_ch = _nx * iy + ix;
     317        3543 :       const Real offset_x = (_nx - 1) * _pitch / 2.0;
     318        3543 :       const Real offset_y = (_ny - 1) * _pitch / 2.0;
     319        3543 :       _subchannel_position[i_ch][0] = _pitch * ix - offset_x;
     320        3543 :       _subchannel_position[i_ch][1] = _pitch * iy - offset_y;
     321             :     }
     322             : 
     323         204 :   if (_n_pins > 0)
     324             :   {
     325             :     // Make channel to pin map.
     326         948 :     for (unsigned int iy = 0; iy < _ny; iy++) // row
     327        4278 :       for (unsigned int ix = 0; ix < _nx; ix++)
     328             :       {
     329        3525 :         const unsigned int i_ch = _nx * iy + ix;
     330             : 
     331             :         // Corners contact 1/4 of one pin.
     332        3525 :         if (iy == 0 && ix == 0)
     333         195 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     334        3330 :         else if (iy == _ny - 1 && ix == 0)
     335         195 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     336        3135 :         else if (iy == 0 && ix == _nx - 1)
     337         195 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     338        2940 :         else if (iy == _ny - 1 && ix == _nx - 1)
     339         195 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     340             :         // Sides contact 1/4 of two pins.
     341        2745 :         else if (iy == 0)
     342             :         {
     343         416 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     344         416 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     345             :         }
     346        2329 :         else if (iy == _ny - 1)
     347             :         {
     348         416 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     349         416 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     350             :         }
     351        1913 :         else if (ix == 0)
     352             :         {
     353         363 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     354         363 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     355             :         }
     356        1550 :         else if (ix == _nx - 1)
     357             :         {
     358         363 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     359         363 :           _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        1187 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     365        1187 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     366        1187 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     367        1187 :           _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     368             :         }
     369             :       }
     370             : 
     371             :     // Make gap to pin map.
     372        5686 :     for (unsigned int ig = 0; ig < _n_gaps; ig++)
     373             :     {
     374        5491 :       const auto i_ch = _gap_to_chan_map[ig].first;
     375        5491 :       const auto j_ch = _gap_to_chan_map[ig].second;
     376        5491 :       const auto & i_pins = _chan_to_pin_map[i_ch];
     377        5491 :       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       20441 :       for (unsigned int i : i_pins)
     383       59270 :         for (unsigned int j : j_pins)
     384       44320 :           if (i == j)
     385             :           {
     386        8644 :             if (_gap_to_pin_map[ig].first == 10000)
     387             :             {
     388        5491 :               _gap_to_pin_map[ig].first = i;
     389        5491 :               _gap_to_pin_map[ig].second = i;
     390             :             }
     391             :             else
     392        3153 :               _gap_to_pin_map[ig].second = i;
     393             :           }
     394             :     }
     395             :   }
     396             : 
     397             :   // Reduce reserved memory in the channel-to-gap map.
     398        3747 :   for (auto & gap : _chan_to_gap_map)
     399             :     gap.shrink_to_fit();
     400             : 
     401             :   // Reduce reserved memory in the channel-to-pin map.
     402        3747 :   for (auto & pin : _chan_to_pin_map)
     403             :     pin.shrink_to_fit();
     404             : 
     405             :   // Reduce reserved memory in the pin-to-channel map.
     406        2365 :   for (auto & pin : _pin_to_chan_map)
     407             :     pin.shrink_to_fit();
     408         204 : }
     409             : 
     410             : void
     411         204 : SCMQuadAssemblyMeshGenerator::buildSubchannelMesh(MeshBase & mesh_base,
     412             :                                                   BoundaryInfo & boundary_info)
     413             : {
     414         204 :   mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * _ny * _nx);
     415         204 :   mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * _ny * _nx);
     416             : 
     417         204 :   _nodes.resize(_nx * _ny);
     418             : 
     419         204 :   const Real offset_x = (_nx - 1) * _pitch / 2.0;
     420         204 :   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         204 :   dof_id_type node_id = mesh_base.n_nodes();
     426         969 :   for (unsigned int iy = 0; iy < _ny; iy++)
     427        4308 :     for (unsigned int ix = 0; ix < _nx; ix++)
     428             :     {
     429        3543 :       const unsigned int i_ch = _nx * iy + ix;
     430        3543 :       _nodes[i_ch].reserve(_n_cells + 1);
     431             : 
     432       69672 :       for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     433       66129 :         _nodes[i_ch].push_back(mesh_base.add_point(
     434      132258 :             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         204 :   dof_id_type elem_id = mesh_base.n_elem();
     440         969 :   for (unsigned int iy = 0; iy < _ny; iy++)
     441        4308 :     for (unsigned int ix = 0; ix < _nx; ix++)
     442       66129 :       for (unsigned int iz = 0; iz < _n_cells; iz++)
     443             :       {
     444       62586 :         Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
     445       62586 :         elem->subdomain_id() = _subchannel_block_id;
     446       62586 :         elem->set_id(elem_id++);
     447             : 
     448       62586 :         const unsigned int i_ch = _nx * iy + ix;
     449       62586 :         elem->set_node(0, _nodes[i_ch][iz]);
     450       62586 :         elem->set_node(1, _nodes[i_ch][iz + 1]);
     451             : 
     452       62586 :         if (iz == 0)
     453        3543 :           boundary_info.add_side(elem, 0, 0);
     454       62586 :         if (iz == _n_cells - 1)
     455        3543 :           boundary_info.add_side(elem, 1, 1);
     456             :       }
     457             : 
     458         204 :   mesh_base.subdomain_name(_subchannel_block_id) = "subchannel";
     459         204 : }
     460             : 
     461             : void
     462         195 : SCMQuadAssemblyMeshGenerator::buildPinMesh(MeshBase & mesh_base)
     463             : {
     464         195 :   mesh_base.reserve_elem(mesh_base.n_elem() + _n_cells * (_ny - 1) * (_nx - 1));
     465         195 :   mesh_base.reserve_nodes(mesh_base.n_nodes() + (_n_cells + 1) * (_ny - 1) * (_nx - 1));
     466             : 
     467         195 :   _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         195 :   const Real offset_x = (_nx - 2) * _pitch / 2.0;
     473         195 :   const Real offset_y = (_ny - 2) * _pitch / 2.0;
     474             : 
     475         195 :   dof_id_type node_id = mesh_base.n_nodes();
     476         753 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     477        2719 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     478             :     {
     479        2161 :       const unsigned int i_pin = (_nx - 1) * iy + ix;
     480        2161 :       _pin_nodes[i_pin].reserve(_n_cells + 1);
     481             : 
     482       41955 :       for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     483       39794 :         _pin_nodes[i_pin].push_back(mesh_base.add_point(
     484       79588 :             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         195 :   dof_id_type elem_id = mesh_base.n_elem();
     489         753 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     490        2719 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     491       39794 :       for (unsigned int iz = 0; iz < _n_cells; iz++)
     492             :       {
     493       37633 :         Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
     494       37633 :         elem->subdomain_id() = _pin_block_id;
     495       37633 :         elem->set_id(elem_id++);
     496             : 
     497       37633 :         const unsigned int i_pin = (_nx - 1) * iy + ix;
     498       37633 :         elem->set_node(0, _pin_nodes[i_pin][iz]);
     499       37633 :         elem->set_node(1, _pin_nodes[i_pin][iz + 1]);
     500             :       }
     501             : 
     502         195 :   mesh_base.subdomain_name(_pin_block_id) = "fuel_pins";
     503         195 : }
     504             : 
     505             : void
     506         204 : SCMQuadAssemblyMeshGenerator::transferMetadata(QuadSubChannelMesh & sch_mesh)
     507             : {
     508             :   // Move the metadata into QuadSubChannelMesh.
     509         204 :   sch_mesh._unheated_length_entry = _unheated_length_entry;
     510         204 :   sch_mesh._heated_length = _heated_length;
     511         204 :   sch_mesh._unheated_length_exit = _unheated_length_exit;
     512         204 :   sch_mesh._z_grid = _z_grid;
     513         204 :   sch_mesh._k_grid = _k_grid;
     514         204 :   sch_mesh._spacer_z = _spacer_z;
     515         204 :   sch_mesh._spacer_k = _spacer_k;
     516         204 :   sch_mesh._z_blockage = _z_blockage;
     517         204 :   sch_mesh._index_blockage = _index_blockage;
     518         204 :   sch_mesh._reduction_blockage = _reduction_blockage;
     519         204 :   sch_mesh._kij = _kij;
     520         204 :   sch_mesh._pitch = _pitch;
     521         204 :   sch_mesh._pin_diameter = _pin_diameter;
     522         204 :   sch_mesh._n_cells = _n_cells;
     523             : 
     524         204 :   sch_mesh._nx = _nx;
     525         204 :   sch_mesh._ny = _ny;
     526         204 :   sch_mesh._n_channels = _n_channels;
     527         204 :   sch_mesh._n_gaps = _n_gaps;
     528         204 :   sch_mesh._n_pins = _n_pins;
     529         204 :   sch_mesh._side_gap = _side_gap;
     530             : 
     531         204 :   sch_mesh._nodes = _nodes;
     532         204 :   sch_mesh._pin_nodes = _pin_nodes;
     533         204 :   sch_mesh._gapnodes = _gapnodes;
     534         204 :   sch_mesh._gap_to_chan_map = _gap_to_chan_map;
     535         204 :   sch_mesh._gap_to_pin_map = _gap_to_pin_map;
     536         204 :   sch_mesh._chan_to_gap_map = _chan_to_gap_map;
     537         204 :   sch_mesh._chan_to_pin_map = _chan_to_pin_map;
     538         204 :   sch_mesh._pin_to_chan_map = _pin_to_chan_map;
     539         204 :   sch_mesh._sign_id_crossflow_map = _sign_id_crossflow_map;
     540         204 :   sch_mesh._gij_map = _gij_map;
     541         204 :   sch_mesh._subchannel_position = _subchannel_position;
     542         204 :   sch_mesh._subch_type = _subch_type;
     543             : 
     544         204 :   sch_mesh._duct_mesh_exist = false;
     545         204 :   sch_mesh._pin_mesh_exist = (_n_pins > 0);
     546         204 : }
     547             : 
     548             : std::unique_ptr<MeshBase>
     549         204 : SCMQuadAssemblyMeshGenerator::generate()
     550             : {
     551         204 :   auto mesh_base = buildMeshBaseObject();
     552             :   BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
     553             : 
     554         204 :   mesh_base->set_spatial_dimension(3);
     555             : 
     556         204 :   buildSubchannelMesh(*mesh_base, boundary_info);
     557             : 
     558         204 :   if (_n_pins > 0)
     559         195 :     buildPinMesh(*mesh_base);
     560             : 
     561         204 :   boundary_info.sideset_name(0) = "inlet";
     562         204 :   boundary_info.sideset_name(1) = "outlet";
     563         204 :   boundary_info.nodeset_name(0) = "inlet";
     564         204 :   boundary_info.nodeset_name(1) = "outlet";
     565             : 
     566         204 :   mesh_base->prepare_for_use();
     567             : 
     568         204 :   auto & sch_mesh = static_cast<QuadSubChannelMesh &>(*_mesh);
     569         204 :   transferMetadata(sch_mesh);
     570         204 :   sch_mesh.computeAssemblyHydraulicParameters();
     571             : 
     572         204 :   return mesh_base;
     573           0 : }

Generated by: LCOV version 1.14