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 "FDTallyGradAux.h"
      22             : 
      23             : #include "CardinalEnums.h"
      24             : 
      25             : registerMooseObject("CardinalApp", FDTallyGradAux);
      26             : 
      27             : InputParameters
      28          10 : FDTallyGradAux::validParams()
      29             : {
      30          10 :   auto params = OpenMCVectorAuxKernel::validParams();
      31          10 :   params.addClassDescription(
      32             :       "An auxkernel which approximates tally gradients at element centroids using "
      33             :       "forward finite differences.");
      34          20 :   params.addRequiredParam<MooseEnum>(
      35             :       "score",
      36          20 :       getSingleTallyScoreEnum(),
      37             :       "The tally score this auxkernel should approximate the gradient of.");
      38          30 :   params.addParam<unsigned int>("ext_filter_bin",
      39          20 :                                 0,
      40             :                                 "The filter bin for the case where any filters are added to this "
      41             :                                 "tally with [Filters] (bin indices start at 0). This parameter "
      42             :                                 "should be specified if you wish to extract the relative error "
      43             :                                 "of a different non-spatial tally bin.");
      44             : 
      45          20 :   params.addRelationshipManager("ElementSideNeighborLayers",
      46             :                                 Moose::RelationshipManagerType::ALGEBRAIC |
      47             :                                     Moose::RelationshipManagerType::GEOMETRIC,
      48          12 :                                 [](const InputParameters &, InputParameters & rm_params)
      49          12 :                                 { rm_params.set<unsigned short>("layers") = 2; });
      50             : 
      51          10 :   return params;
      52           0 : }
      53             : 
      54           8 : FDTallyGradAux::FDTallyGradAux(const InputParameters & parameters)
      55             :   : OpenMCVectorAuxKernel(parameters),
      56           8 :     _bin_index(getParam<unsigned int>("ext_filter_bin")),
      57           8 :     _sum_y_y_t(RealEigenMatrix::Zero(3, 3)),
      58          16 :     _sum_y_du_dy(RealEigenVector::Zero(3))
      59             : {
      60           8 :   if (_var.feType() != FEType(libMesh::CONSTANT, libMesh::MONOMIAL_VEC))
      61           2 :     paramError("variable",
      62             :                "FDTallyGradAux only supports CONSTANT MONOMIAL_VEC shape functions. Please "
      63             :                "ensure that 'variable' is of type MONOMIAL_VEC and order CONSTANT.");
      64             : 
      65          12 :   std::string score = getParam<MooseEnum>("score");
      66             :   std::replace(score.begin(), score.end(), '_', '-');
      67             : 
      68             :   // Error check and fetch the tally score.
      69           6 :   if (!_openmc_problem->hasScore(score))
      70           2 :     paramError("score",
      71           2 :                "The problem does not contain any score named " +
      72           2 :                    std::string(getParam<MooseEnum>("score")) +
      73             :                    "! Please "
      74             :                    "ensure that one of your [Tallies] is scoring the requested score.");
      75             : 
      76           4 :   auto score_vars = _openmc_problem->getTallyScoreVariables(score, _tid);
      77           4 :   auto score_bins = _openmc_problem->getTallyScoreVariableValues(score, _tid);
      78           4 :   auto neighbor_score_bins = _openmc_problem->getTallyScoreNeighborVariableValues(score, _tid);
      79             : 
      80           4 :   if (_bin_index >= score_bins.size())
      81           2 :     paramError("ext_filter_bin",
      82             :                "The external filter bin provided is invalid for the number of "
      83           2 :                "external filter bins (" +
      84           2 :                    std::to_string(score_bins.size()) +
      85             :                    ") "
      86           2 :                    "applied to " +
      87           2 :                    std::string(getParam<MooseEnum>("score")) + "!");
      88             : 
      89           2 :   if (score_vars[_bin_index]->feType() != FEType(libMesh::CONSTANT, libMesh::MONOMIAL))
      90           0 :     paramError(
      91             :         "score",
      92             :         "FDTallyGradAux only supports CONSTANT MONOMIAL shape functions for tally variables.");
      93             : 
      94           2 :   _tally_val = score_bins[_bin_index];
      95           2 :   _tally_neighbor_val = neighbor_score_bins[_bin_index];
      96           2 : }
      97             : 
      98             : void
      99         512 : FDTallyGradAux::compute()
     100             : {
     101         512 :   auto elem_c = _current_elem->true_centroid();
     102        3584 :   for (unsigned int side = 0; side < _current_elem->n_sides(); side++)
     103             :   {
     104        3072 :     const Elem * neighbor = _current_elem->neighbor_ptr(side);
     105             : 
     106             :     // If the neighbor is null, the current side is a boundary and we can skip
     107             :     // it. Geometric ghosting ensures that element neighbors are on the
     108             :     // processors that need them during compute().
     109        3072 :     if (!neighbor)
     110         192 :       continue;
     111             : 
     112           0 :     if (!neighbor->active())
     113           0 :       continue;
     114             : 
     115             :     // Can skip the gradient computation for this side if the neighbor isn't on
     116             :     // this block.
     117        2880 :     if (!hasBlocks(neighbor->subdomain_id()))
     118           0 :       continue;
     119             : 
     120        2880 :     _openmc_problem->reinitNeighbor(_current_elem, side, _tid);
     121             : 
     122             :     // Fetch the vector pointing from the current element's centroid to the
     123             :     // neighbor's centroid (y').
     124        2880 :     auto y_prime = neighbor->true_centroid() - elem_c;
     125        2880 :     RealEigenVector y_prime_eig(3);
     126        2880 :     y_prime_eig << y_prime(0), y_prime(1), y_prime(2);
     127             : 
     128             :     // Compute du/dy along the direction pointing towards the neighbor's centroid.
     129             :     // Add to the b vector.
     130        2880 :     _sum_y_du_dy +=
     131        2880 :         y_prime_eig * ((*_tally_neighbor_val)[0] - (*_tally_val)[0]) / y_prime.norm_sq();
     132             :     // Compute the outer product between y' and y'.T.
     133             :     // Add to the A matrix.
     134        2880 :     _sum_y_y_t += y_prime_eig * y_prime_eig.transpose();
     135             :   }
     136             : 
     137             :   // Solve for an estimate of the tally gradient at the current element's centroid.
     138        1024 :   RealEigenVector val = _sum_y_y_t.fullPivLu().solve(_sum_y_du_dy);
     139             :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     140         512 :   _var.setDofValue(val(0), 0);
     141         512 :   _var.setDofValue(val(1), 1);
     142         512 :   _var.setDofValue(val(2), 2);
     143             : 
     144         512 :   _sum_y_y_t = RealEigenMatrix::Zero(3, 3);
     145         512 :   _sum_y_du_dy = RealEigenVector::Zero(3);
     146         512 : }
     147             : 
     148             : #endif