Line data    Source code

       1             : /********************************************************************/
       2             : /*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
       3             : /*                             Cardinal                             */
       4             : /*                                                                  */
       5             : /*                  (c) 2021 UChicago Argonne, LLC                  */
       6             : /*                        ALL RIGHTS RESERVED                       */
       7             : /*                                                                  */
       8             : /*                 Prepared by UChicago Argonne, LLC                */
       9             : /*               Under Contract No. DE-AC02-06CH11357               */
      10             : /*                With the U. S. Department of Energy               */
      11             : /*                                                                  */
      12             : /*             Prepared by Battelle Energy Alliance, LLC            */
      13             : /*               Under Contract No. DE-AC07-05ID14517               */
      14             : /*                With the U. S. Department of Energy               */
      15             : /*                                                                  */
      16             : /*                 See LICENSE for full restrictions                */
      17             : /********************************************************************/
      18             : 
      19             : #ifdef ENABLE_OPENMC_COUPLING
      20             : 
      21             : #include "OpenMCVolumeCalculation.h"
      22             : #include "OpenMCCellAverageProblem.h"
      23             : #include "UserErrorChecking.h"
      24             : 
      25             : registerMooseObject("CardinalApp", OpenMCVolumeCalculation);
      26             : 
      27             : InputParameters
      28         258 : OpenMCVolumeCalculation::validParams()
      29             : {
      30         258 :   InputParameters params = GeneralUserObject::validParams();
      31         258 :   params += OpenMCBase::validParams();
      32         516 :   params.addParam<Point>("lower_left", "Lower left of the bounding box inside of which to "
      33             :     "compute volumes. If not specified, this will default to the lower left of the [Mesh] "
      34             :     "(which will NOT capture any OpenMC geometry that lies outside the [Mesh] extents.");
      35         516 :   params.addParam<Point>("upper_right", "Upper right of the bounding box inside of which to "
      36             :     "compute volumes. If not specified, this will default to the upper right of the [Mesh] "
      37             :     "(which will NOT capture any OpenMC geometry that lies outside the [Mesh] extents.");
      38         516 :   params.addRequiredRangeCheckedParam<unsigned int>("n_samples", "n_samples > 0",
      39             :     "Number of samples to use for the stochastic volume calculation");
      40         516 :   params.addRangeCheckedParam<Real>("trigger_threshold", "trigger_threshold > 0",
      41             :     "Trigger threshold to decide when to terminate the volume calculation");
      42             : 
      43         516 :   MooseEnum trigger("rel_err none", "none");
      44         516 :   params.addParam<MooseEnum>("trigger", getTallyTriggerEnum(),
      45             :     "Type of trigger to apply to decide when to terminate the volume calculation");
      46         258 :   params.addClassDescription("Stochastic volume calculation of the OpenMC cells which map to MOOSE");
      47         258 :   return params;
      48         258 : }
      49             : 
      50         129 : OpenMCVolumeCalculation::OpenMCVolumeCalculation(const InputParameters & parameters)
      51             :   : GeneralUserObject(parameters),
      52             :     OpenMCBase(this, parameters),
      53         129 :     _n_samples(getParam<unsigned int>("n_samples")),
      54         516 :     _trigger(getParam<MooseEnum>("trigger"))
      55             : {
      56         129 :   if (_trigger == "none")
      57         226 :     checkUnusedParam(parameters, "trigger_threshold", "not specifying any triggers");
      58          16 :   else if (_trigger == "rel_err")
      59             :   {
      60          32 :     checkRequiredParam(parameters, "trigger_threshold", "using a trigger");
      61          32 :     _trigger_threshold = getParam<Real>("trigger_threshold");
      62             :   }
      63             :   else
      64           0 :     mooseError("Unhandled trigger enum in OpenMCCellVolumeCalculation!");
      65             : 
      66         129 :   _scaling = _openmc_problem->scaling();
      67             : 
      68         258 :   if (isParamValid("lower_left"))
      69          20 :     _lower_left = getParam<Point>("lower_left");
      70         258 :   if (isParamValid("upper_right"))
      71          16 :     _upper_right = getParam<Point>("upper_right");
      72         129 : }
      73             : openmc::Position
      74         254 : OpenMCVolumeCalculation::position(const Point & pt) const
      75             : {
      76         254 :   openmc::Position p {pt(0), pt(1), pt(2)};
      77         254 :   return p;
      78             : }
      79             : 
      80             : void
      81           2 : OpenMCVolumeCalculation::resetVolumeCalculation()
      82             : {
      83           2 :   auto idx = openmc::model::volume_calcs.begin() + _calc_index;
      84             :   openmc::model::volume_calcs.erase(idx);
      85           2 : }
      86             : 
      87             : void
      88         129 : OpenMCVolumeCalculation::initializeVolumeCalculation()
      89             : {
      90             : 
      91         129 :   BoundingBox box = MeshTools::create_bounding_box(_openmc_problem->getMooseMesh().getMesh());
      92             : 
      93         258 :   if (!isParamValid("lower_left"))
      94         119 :     _lower_left = box.min();
      95         258 :   if (!isParamValid("upper_right"))
      96         121 :     _upper_right = box.max();
      97         129 :   if (_lower_left >= _upper_right)
      98           2 :     paramError("upper_right",
      99             :                "The 'upper_right' (",
     100             :                _upper_right(0),
     101             :                ", ",
     102             :                _upper_right(1),
     103             :                ", ",
     104             :                _upper_right(2),
     105             :                ") "
     106             :                "must be greater than the 'lower_left' (",
     107             :                _lower_left(0),
     108             :                ", ",
     109             :                _lower_left(1),
     110             :                ", ",
     111             :                _lower_left(2),
     112             :                ")!");
     113             : 
     114         127 :   _volume_calc.reset(new openmc::VolumeCalculation());
     115         127 :   _volume_calc->domain_type_ = openmc::VolumeCalculation::TallyDomain::CELL;
     116         127 :   _volume_calc->lower_left_ = position(_lower_left * _scaling);
     117         127 :   _volume_calc->upper_right_ = position(_upper_right * _scaling);
     118         127 :   _volume_calc->n_samples_ = _n_samples;
     119             : 
     120         127 :   if (_trigger == "rel_err")
     121             :   {
     122          16 :     _volume_calc->threshold_ = _trigger_threshold;
     123          16 :     _volume_calc->trigger_type_ = openmc::TriggerMetric::relative_error;
     124             :   }
     125             : 
     126         127 :   auto cell_to_elem = _openmc_problem->cellToElem();
     127             : 
     128             :   std::set<int> ids;
     129             :   _index_to_calc_index.clear();
     130             : 
     131             :   int i = 0;
     132        2239 :   for (const auto & c : cell_to_elem)
     133             :   {
     134        2112 :     auto index = c.first.first;
     135        2112 :     auto id = openmc::model::cells[index]->id_;
     136             : 
     137        2112 :     ids.insert(id);
     138        2112 :     if (!_index_to_calc_index.count(c.first.first))
     139        1526 :       _index_to_calc_index[c.first.first] = i++;
     140             :   }
     141             : 
     142         127 :   std::vector<int> domain_ids(ids.begin(), ids.end());
     143         127 :   _volume_calc->domain_ids_ = domain_ids;
     144         127 :   openmc::model::volume_calcs.push_back(*_volume_calc);
     145         127 :   _calc_index = openmc::model::volume_calcs.size();
     146         127 : }
     147             : 
     148             : void
     149         127 : OpenMCVolumeCalculation::computeVolumes()
     150             : {
     151         127 :   _console << "Running stochastic volume calculation... ";
     152         127 :   _results = _volume_calc->execute();
     153         127 : }
     154             : 
     155             : void
     156        2054 : OpenMCVolumeCalculation::cellVolume(const unsigned int & index, Real & volume, Real & std_dev) const
     157             : {
     158        2054 :   auto calc_index = _index_to_calc_index.at(index);
     159        2054 :   auto n_instances = openmc::model::cells[index]->n_instances_;
     160        2054 :   if (n_instances > 1)
     161         632 :     mooseDoOnce(mooseWarning(
     162             :         "OpenMC's stochastic volume calculation cannot individually measure volumes of cell "
     163             :         "INSTANCES. We assume that no cell instances are clipped by other cells (e.g. they are all "
     164             :         "identical) so that the volume of an individual instance is equal to the cell volume "
     165             :         "divided by number of instances. For most cases, this is correct - but if you have any "
     166             :         "instances which only partially exist in the geometry, this will give INCORRECT volumes."));
     167             : 
     168             :   // means add
     169        2052 :   volume = _results[calc_index].volume[0] / (_scaling * _scaling * _scaling) / n_instances;
     170             : 
     171             :   // but standard deviations only add as the variances
     172        2052 :   std_dev = _results[calc_index].volume[1] / (_scaling * _scaling * _scaling) / std::sqrt(n_instances);
     173        2052 : }
     174             : 
     175             : #endif