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 "OpenMCBase.h"
      22             : #include "openmc/eigenvalue.h"
      23             : #include "openmc/math_functions.h"
      24             : #include "libmesh/elem.h"
      25             : 
      26             : InputParameters
      27       13911 : OpenMCBase::validParams()
      28             : {
      29       13911 :   InputParameters params = emptyInputParameters();
      30       13911 :   return params;
      31             : }
      32             : 
      33        6326 : OpenMCBase::OpenMCBase(const ParallelParamObject * moose_object, const InputParameters & parameters)
      34        6326 :   : _openmc_problem(&getOpenMCProblem(*moose_object))
      35             : {
      36        6324 : }
      37             : 
      38             : Real
      39          96 : OpenMCBase::stdev(const double & mean, const double & sum_sq, unsigned int realizations) const
      40             : {
      41             :   return realizations > 1
      42         192 :              ? std::sqrt(std::max(0.0, (sum_sq / realizations - mean * mean) / (realizations - 1)))
      43          96 :              : 0.0;
      44             : }
      45             : 
      46             : Real
      47        1204 : OpenMCBase::kMean(const eigenvalue::EigenvalueEnum estimator) const
      48             : {
      49        1204 :   int n = openmc::simulation::n_realizations;
      50             :   const auto & gt = openmc::simulation::global_tallies;
      51             : 
      52        1204 :   switch (estimator)
      53             :   {
      54             :     case eigenvalue::collision:
      55          56 :       return gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n;
      56             :     case eigenvalue::absorption:
      57          48 :       return gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n;
      58             :     case eigenvalue::tracklength:
      59          48 :       return gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n;
      60        1052 :     case eigenvalue::combined:
      61             :     {
      62        1052 :       if (n <= 3)
      63           0 :         mooseError("Cannot compute combined k-effective estimate with fewer than 4 realizations!\n"
      64             :                    "Please change the estimator type to either 'collision', 'tracklength', or "
      65             :                    "'absorption'.");
      66             : 
      67             :       double k_eff[2];
      68        1052 :       openmc::openmc_get_keff(k_eff);
      69        1052 :       return k_eff[0];
      70             :     }
      71           0 :     default:
      72           0 :       mooseError("Internal error: Unhandled EigenvalueEnum!");
      73             :   }
      74             : }
      75             : 
      76             : Real
      77         668 : OpenMCBase::kStandardDeviation(const eigenvalue::EigenvalueEnum estimator) const
      78             : {
      79             :   const auto & gt = openmc::simulation::global_tallies;
      80         668 :   int n = openmc::simulation::n_realizations;
      81             : 
      82             :   double t_n1 = 1.0;
      83         668 :   if (openmc::settings::confidence_intervals)
      84             :   {
      85             :     double alpha = 1.0 - openmc::CONFIDENCE_LEVEL;
      86           0 :     t_n1 = openmc::t_percentile(1.0 - alpha / 2.0, n - 1);
      87             :   }
      88             : 
      89         668 :   switch (estimator)
      90             :   {
      91             :     case eigenvalue::collision:
      92             :     {
      93          32 :       double mean = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n;
      94          32 :       double sum_sq = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM_SQ);
      95          32 :       return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations);
      96             :     }
      97             :     case eigenvalue::absorption:
      98             :     {
      99          32 :       double mean = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n;
     100          32 :       double sum_sq = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM_SQ);
     101          32 :       return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations);
     102             :     }
     103             :     case eigenvalue::tracklength:
     104             :     {
     105          32 :       double mean = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n;
     106          32 :       double sum_sq = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM_SQ);
     107          32 :       return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations);
     108             :     }
     109         572 :     case eigenvalue::combined:
     110             :     {
     111         572 :       if (n <= 3)
     112           0 :         mooseError("Cannot compute combined k-effective standard deviation with fewer than 4 "
     113             :                    "realizations!");
     114             : 
     115             :       double k_eff[2];
     116         572 :       openmc::openmc_get_keff(k_eff);
     117         572 :       return k_eff[1];
     118             :     }
     119           0 :     default:
     120           0 :       mooseError("Internal error: Unhandled StandardDeviationEnum!");
     121             :   }
     122             : }
     123             : 
     124             : OpenMCCellAverageProblem &
     125        6414 : getOpenMCProblem(const ParallelParamObject & moose_object)
     126             : {
     127             :   OpenMCCellAverageProblem * const openmc_problem =
     128        6414 :       dynamic_cast<OpenMCCellAverageProblem *>(&moose_object.getMooseApp().feProblem());
     129        6414 :   if (!openmc_problem)
     130           2 :     moose_object.mooseError(
     131             :         "This object can only be used with problems of type 'OpenMCCellAverageProblem'!");
     132        6412 :   return *openmc_problem;
     133             : }
     134             : 
     135             : #endif