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 "MultiControlDrumFunction.h"
      11             : 
      12             : registerMooseObject("ReactorApp", MultiControlDrumFunction);
      13             : 
      14             : InputParameters
      15          54 : MultiControlDrumFunction::validParams()
      16             : {
      17          54 :   InputParameters params = Function::validParams();
      18          54 :   params.addClassDescription(
      19             :       "A function that returns an absorber fraction for multiple control drums application.");
      20         108 :   params.addRequiredParam<MeshGeneratorName>(
      21             :       "mesh_generator",
      22             :       "Name of the mesh generator to be used to retrieve control drums information.");
      23         108 :   params.addRequiredParam<std::vector<Real>>(
      24             :       "angular_speeds", "Vector of angular rotation speeds of the control drum.");
      25         108 :   params.addRequiredParam<std::vector<Real>>(
      26             :       "start_angles",
      27             :       "Vector of initial angular positions of the beginning of the absorber components.");
      28         108 :   params.addRequiredParam<std::vector<Real>>(
      29             :       "angle_ranges", "Vector of angular ranges of the beginning of the absorber components.");
      30         108 :   params.addParam<Real>(
      31         108 :       "rotation_start_time", 0.0, "The time point at which the control drums start rotating.");
      32         108 :   params.addParam<Real>("rotation_end_time",
      33         108 :                         std::numeric_limits<Real>::max(),
      34             :                         "The time point at which the control drums stop rotating.");
      35         108 :   params.addParam<bool>(
      36         108 :       "use_control_drum_id", true, "Whether extra element id user_control_drum_id is used.");
      37         108 :   params.addParamNamesToGroup("use_control_drum_id", "Advanced");
      38             : 
      39          54 :   return params;
      40           0 : }
      41             : 
      42          32 : MultiControlDrumFunction::MultiControlDrumFunction(const InputParameters & parameters)
      43             :   : Function(parameters),
      44             :     MeshMetaDataInterface(this),
      45             :     ElementIDInterface(this),
      46           0 :     _mesh_generator(getParam<MeshGeneratorName>("mesh_generator")),
      47          64 :     _angular_speeds(getParam<std::vector<Real>>("angular_speeds")),
      48          64 :     _start_angles(getParam<std::vector<Real>>("start_angles")),
      49          64 :     _angle_ranges(getParam<std::vector<Real>>("angle_ranges")),
      50          64 :     _rotation_start_time(getParam<Real>("rotation_start_time")),
      51          64 :     _rotation_end_time(getParam<Real>("rotation_end_time")),
      52          64 :     _use_control_drum_id(getParam<bool>("use_control_drum_id")),
      53          32 :     _control_drum_id(_use_control_drum_id ? getElementIDByName("control_drum_id")
      54             :                                           : DofObject::invalid_id),
      55          32 :     _control_drum_positions(
      56          32 :         getMeshProperty<std::vector<Point>>("control_drum_positions", _mesh_generator)),
      57          32 :     _control_drums_azimuthal_meta(getMeshProperty<std::vector<std::vector<Real>>>(
      58          32 :         "control_drums_azimuthal_meta", _mesh_generator))
      59             : {
      60          32 :   if (_angular_speeds.size() != _start_angles.size())
      61           2 :     paramError("start_angles",
      62             :                "the size of 'start_angles' must be identical to the size of 'angular_speeds'");
      63          30 :   if (_angular_speeds.size() != _angle_ranges.size())
      64           2 :     paramError("angle_ranges",
      65             :                "the size of 'angle_ranges' must be identical to the size of 'angular_speeds'");
      66          28 :   if (_angular_speeds.size() != _control_drum_positions.size())
      67           2 :     paramError("angular_speeds",
      68             :                "the size of this parameter must be identical to the control drum number recorded "
      69             :                "in the MeshMetaData.");
      70          26 :   if (_rotation_end_time <= _rotation_start_time)
      71           2 :     paramError("rotation_end_time", "this parameter must be larger than rotation_start_time.");
      72          24 : }
      73             : 
      74             : Real
      75        6912 : MultiControlDrumFunction::value(Real t, const Point & p) const
      76             : {
      77             :   // Calculate the effective rotation time
      78        6912 :   const Real t_eff = t < _rotation_start_time
      79        6912 :                          ? 0.0
      80        6912 :                          : (t > _rotation_end_time ? (_rotation_end_time - _rotation_start_time)
      81             :                                                    : t - _rotation_start_time);
      82             :   // Find the closest control drum for a given Point p; only needed if control drum id is not used.
      83             :   std::vector<std::pair<Real, unsigned int>> meshcontrol_drum_dist_vec;
      84        6912 :   if (!_use_control_drum_id)
      85             :   {
      86       24192 :     for (unsigned int i = 0; i < _control_drum_positions.size(); i++)
      87             :     {
      88       20736 :       Real cd_dist = std::sqrt(Utility::pow<2>(p(0) - _control_drum_positions[i](0)) +
      89       20736 :                                Utility::pow<2>(p(1) - _control_drum_positions[i](1)));
      90       20736 :       meshcontrol_drum_dist_vec.push_back(std::make_pair(cd_dist, i));
      91             :     }
      92        3456 :     std::sort(meshcontrol_drum_dist_vec.begin(), meshcontrol_drum_dist_vec.end());
      93             :   }
      94        6912 :   const unsigned int cd_id = _use_control_drum_id
      95        6912 :                                  ? (_control_drum_id > 0 ? _control_drum_id - 1 : 0)
      96        6912 :                                  : meshcontrol_drum_dist_vec.front().second;
      97        6912 :   const auto & cd_pos = _control_drum_positions[cd_id];
      98             : 
      99        6912 :   Real dynamic_start_angle = (_angular_speeds[cd_id] * t_eff + _start_angles[cd_id]) / 180.0 * M_PI;
     100        6912 :   Real dynamic_end_angle = dynamic_start_angle + _angle_ranges[cd_id] / 180.0 * M_PI;
     101             : 
     102        6912 :   dynamic_start_angle = atan2(std::sin(dynamic_start_angle),
     103        6912 :                               std::cos(dynamic_start_angle)) /
     104        6912 :                         M_PI * 180.0; // quick way to move to -M_PI to M_PI
     105        6912 :   dynamic_end_angle = atan2(std::sin(dynamic_end_angle),
     106        6912 :                             std::cos(dynamic_end_angle)) /
     107        6912 :                       M_PI * 180.0; // quick way to move to -M_PI to M_PI
     108             : 
     109             :   // Locate the first mesh azimuthal angle that is greater than dynamic_start_angle
     110        6912 :   auto start_bound = std::upper_bound(_control_drums_azimuthal_meta[cd_id].begin(),
     111        6912 :                                       _control_drums_azimuthal_meta[cd_id].end(),
     112             :                                       dynamic_start_angle);
     113             :   // Locate the first mesh azimuthal angle that is greater than dynamic_end_angle
     114        6912 :   auto end_bound = std::upper_bound(_control_drums_azimuthal_meta[cd_id].begin(),
     115        6912 :                                     _control_drums_azimuthal_meta[cd_id].end(),
     116             :                                     dynamic_end_angle);
     117             : 
     118             :   // Locate the elements that contain the start/end angles.
     119             :   // This part seems long; but it only solves one simple issue -> transition from M_PI to -M_PI
     120             : 
     121             :   // The two azimuthal ends of the element that the start angle intercepts;
     122             :   // and the two azimuthal ends of the element that the end angle intercepts;
     123             :   Real start_low, start_high, end_low, end_high;
     124             :   // This means that the dynamic_start_angle is greater than the lowest mesh azimuthal angle and
     125             :   // lower than the greatest mesh azimuthal angle. Namely, dynamic_start_angle does not cause the
     126             :   // "transition from M_PI to -M_PI" issue.
     127        6912 :   if (start_bound != _control_drums_azimuthal_meta[cd_id].begin() &&
     128             :       start_bound != _control_drums_azimuthal_meta[cd_id].end())
     129             :   {
     130        6816 :     start_low = *(start_bound - 1);
     131        6816 :     start_high = *start_bound;
     132             :   }
     133             :   // On the contrary, this means the dynamic_start_angle intercepts an element across the
     134             :   // M_PI(-M_PI) azimuthal angle. Namely, start_high is actually lower than start_low because of
     135             :   // this transition.
     136             :   else
     137             :   {
     138          96 :     start_high = _control_drums_azimuthal_meta[cd_id].front();
     139          96 :     start_low = _control_drums_azimuthal_meta[cd_id].back();
     140             :   }
     141             : 
     142             :   // This means that dynamic_end_angle is greater than the lowest mesh azimuthal angle and lower
     143             :   // than the greatest mesh azimuthal angle. Namely, dynamic_end_angle does not cause the
     144             :   // "transition from M_PI to -M_PI" issue.
     145        6912 :   if (end_bound != _control_drums_azimuthal_meta[cd_id].begin() &&
     146             :       end_bound != _control_drums_azimuthal_meta[cd_id].end())
     147             :   {
     148        6528 :     end_low = *(end_bound - 1);
     149        6528 :     end_high = *end_bound;
     150             :   }
     151             :   // On the contrary, this means the dynamic_end_angle intercepts an element across the
     152             :   // M_PI(-M_PI) azimuthal angle. Namely, end_high is actually lower than end_low because of
     153             :   // this transition.
     154             :   else
     155             :   {
     156         384 :     end_high = _control_drums_azimuthal_meta[cd_id].front();
     157         384 :     end_low = _control_drums_azimuthal_meta[cd_id].back();
     158             :   }
     159             : 
     160             :   // azimuthal_p is the azimuthal angle of the point whose functon value needs to be calculated.
     161        6912 :   Real azimuthal_p = atan2(p(1) - cd_pos(1), p(0) - cd_pos(0)) / M_PI * 180;
     162             : 
     163             :   // If there is not periodical transition from M_PI to -M_PI, we always have:
     164             :   // start_low < start_high < end_low < end_high
     165             :   // In the presence of "transition from M_PI to -M_PI", the relation above is compromised.
     166             :   // Based on how the relation is compromised, we can tell where the transition occurs.
     167             : 
     168             :   // Most trivial scenario -> no transition from M_PI to -M_PI is involved (it happens to a pure
     169             :   // reflector element)
     170        6912 :   if (end_high >= start_low)
     171             :   {
     172             :     // the point is located in an element that is purely absorber.
     173        5952 :     if (azimuthal_p >= start_high && azimuthal_p <= end_low)
     174             :       return 100.0;
     175             :     // the point is located in an element that is purely reflector.
     176        5034 :     else if (azimuthal_p <= start_low || azimuthal_p >= end_high)
     177             :       return 0.0;
     178             :     // the point is located in a mixed element intercepted by dynamic_start_angle so that volume
     179             :     // fraction is calculated.
     180         744 :     else if (azimuthal_p < start_high && azimuthal_p > start_low)
     181             :     {
     182         372 :       Real start_interval = start_high - start_low;
     183         372 :       Real stab_interval = start_high - dynamic_start_angle;
     184         372 :       return stab_interval / start_interval * 100.0;
     185             :     }
     186             :     // the point is located in a mixed element intercepted by dynamic_end_angle so that volume
     187             :     // fraction is calculated.
     188             :     else
     189             :     {
     190         372 :       Real end_interval = end_high - end_low;
     191         372 :       Real endab_interval = dynamic_end_angle - end_low;
     192         372 :       return endab_interval / end_interval * 100.0;
     193             :     }
     194             :   }
     195             :   // The element intercepted by dynamic_end_angle has the "transition from M_PI to -M_PI" issue.
     196             :   // This is still relatively simple because only one of the mixed element is affected.
     197         960 :   else if (end_low >= end_high)
     198             :   {
     199             :     // (Not affected) the point is located in an element that is purely absorber.
     200         384 :     if (azimuthal_p >= start_high && azimuthal_p <= end_low)
     201             :       return 100.0;
     202             :     // (Affected) the point is located in an element that is purely reflector.
     203         286 :     else if (azimuthal_p <= start_low && azimuthal_p >= end_high)
     204             :       return 0.0;
     205             :     // (Not affected) the point is located in a mixed element intercepted by dynamic_start_angle so
     206             :     // that volume fraction is calculated.
     207          48 :     else if (azimuthal_p < start_high && azimuthal_p > start_low)
     208             :     {
     209          24 :       Real start_interval = start_high - start_low;
     210          24 :       Real stab_interval = start_high - dynamic_start_angle;
     211          24 :       return stab_interval / start_interval * 100.0;
     212             :     }
     213             :     // (Affected) the point is located in a mixed element intercepted by dynamic_end_angle so that
     214             :     // volume fraction is calculated. As this element is also intercepted by azimuthal_angle = M_PI
     215             :     // (-M_PI), an extra 2 * M_PI (360 degrees) shift is used to correct this transition effect.
     216             :     else
     217             :     {
     218          24 :       Real end_interval = end_high - end_low + 360.0;
     219          24 :       Real endab_interval = (dynamic_end_angle - end_low >= 0.0)
     220          24 :                                 ? (dynamic_end_angle - end_low)
     221             :                                 : (dynamic_end_angle - end_low + 360.0);
     222          24 :       return endab_interval / end_interval * 100.0;
     223             :     }
     224             :   }
     225             :   // The "transition from M_PI to -M_PI" happens to an pure absorber element
     226         576 :   else if (start_high >= end_low)
     227             :   {
     228             :     // (Affected) the point is located in an element that is purely absorber.
     229         480 :     if (azimuthal_p >= start_high || azimuthal_p <= end_low)
     230             :       return 100.0;
     231             :     // (Affected) the point is located in an element that is purely reflector.
     232         372 :     else if (azimuthal_p >= end_high && azimuthal_p <= start_low)
     233             :       return 0.0;
     234             :     // (Not affected) the point is located in a mixed element intercepted by dynamic_start_angle so
     235             :     // that volume fraction is calculated.
     236          60 :     else if (azimuthal_p < start_high && azimuthal_p > start_low)
     237             :     {
     238          30 :       Real start_interval = start_high - start_low;
     239          30 :       Real stab_interval = start_high - dynamic_start_angle;
     240          30 :       return stab_interval / start_interval * 100.0;
     241             :     }
     242             :     // (Not affected) the point is located in a mixed element intercepted by dynamic_end_angle so
     243             :     // that volume fraction is calculated.
     244             :     else
     245             :     {
     246          30 :       Real end_interval = end_high - end_low;
     247          30 :       Real endab_interval = dynamic_end_angle - end_low;
     248          30 :       return endab_interval / end_interval * 100.0;
     249             :     }
     250             :   }
     251             :   // The element intercepted by dynamic_start_angle has the "transition from M_PI to -M_PI" issue.
     252             :   // This is still relatively simple because only one of the mixed element is affected.
     253             :   else
     254             :   {
     255             :     // (Not affected) the point is located in an element that is purely absorber.
     256          96 :     if (azimuthal_p >= start_high && azimuthal_p <= end_low)
     257             :       return 100.0;
     258             :     // (Affected) the point is located in an element that is purely reflector.
     259          84 :     else if (azimuthal_p >= end_high && azimuthal_p <= start_low)
     260             :       return 0.0;
     261             :     // (Affected) the point is located in a mixed element intercepted by dynamic_start_angle so that
     262             :     // volume fraction is calculated. As this element is also intercepted by azimuthal_angle = M_PI
     263             :     // (-M_PI), an extra 2 * M_PI (360 degrees) shift is used to correct this transition effect.
     264          12 :     else if (azimuthal_p > start_low || azimuthal_p < start_high)
     265             :     {
     266           6 :       Real start_interval = start_high - start_low + 360.0;
     267           6 :       Real stab_interval = (start_high - dynamic_start_angle >= 0.0)
     268           6 :                                ? (start_high - dynamic_start_angle)
     269             :                                : (start_high - dynamic_start_angle + 360.0);
     270           6 :       return stab_interval / start_interval * 100.0;
     271             :     }
     272             :     // (Not affected) the point is located in a mixed element intercepted by dynamic_end_angle so
     273             :     // that volume fraction is calculated.
     274             :     else
     275             :     {
     276           6 :       Real end_interval = end_high - end_low;
     277           6 :       Real endab_interval = dynamic_end_angle - end_low;
     278           6 :       return endab_interval / end_interval * 100.0;
     279             :     }
     280             :   }
     281        6912 : }