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             : #include "TallyBase.h"
      25             : #include "UserErrorChecking.h"
      26             : 
      27             : registerMooseObject("CardinalApp", FDTallyGradAux);
      28             : 
      29             : InputParameters
      30          26 : FDTallyGradAux::validParams()
      31             : {
      32          26 :   auto params = OpenMCVectorAuxKernel::validParams();
      33          26 :   params.addClassDescription(
      34             :       "An auxkernel which approximates tally gradients at element centroids using "
      35             :       "forward finite differences.");
      36          52 :   params.addRequiredParam<MooseEnum>(
      37             :       "score",
      38          52 :       getSingleTallyScoreEnum(),
      39             :       "The tally score this auxkernel should approximate the gradient of.");
      40          52 :   params.addParam<unsigned int>("ext_filter_bin",
      41          52 :                                 0,
      42             :                                 "The filter bin for the case where any filters are added to this "
      43             :                                 "tally with [Filters] (bin indices start at 0). This parameter "
      44             :                                 "should be specified if you wish to extract the relative error "
      45             :                                 "of a different non-spatial tally bin.");
      46          52 :   params.addParam<std::string>(
      47             :       "tally",
      48             :       "The name of the tally to fetch the score variable from. Only required if "
      49             :       "your problem contains multiple tallies which accumulate the same score.");
      50             : 
      51          52 :   params.addRelationshipManager("ElementSideNeighborLayers",
      52             :                                 Moose::RelationshipManagerType::ALGEBRAIC |
      53             :                                     Moose::RelationshipManagerType::GEOMETRIC,
      54          36 :                                 [](const InputParameters &, InputParameters & rm_params)
      55          36 :                                 { rm_params.set<unsigned short>("layers") = 2; });
      56             : 
      57          26 :   return params;
      58           0 : }
      59             : 
      60          16 : FDTallyGradAux::FDTallyGradAux(const InputParameters & parameters)
      61             :   : OpenMCVectorAuxKernel(parameters),
      62          16 :     _bin_index(getParam<unsigned int>("ext_filter_bin")),
      63          16 :     _sum_y_y_t(RealEigenMatrix::Zero(3, 3)),
      64          32 :     _sum_y_du_dy(RealEigenVector::Zero(3))
      65             : {
      66          16 :   if (_var.feType() != FEType(libMesh::CONSTANT, libMesh::MONOMIAL_VEC))
      67           2 :     paramError("variable",
      68             :                "FDTallyGradAux only supports CONSTANT MONOMIAL_VEC shape functions. Please "
      69             :                "ensure that 'variable' is of type MONOMIAL_VEC and order CONSTANT.");
      70             : 
      71          14 :   auto score = getScore("score");
      72          12 :   auto tally_name = tallyByScore(score, "tally");
      73             :   ;
      74             : 
      75          12 :   auto score_vars = _openmc_problem->getTallyScoreVariables(score, tally_name, _tid);
      76          12 :   auto score_bins = _openmc_problem->getTallyScoreVariableValues(score, tally_name, _tid);
      77             :   auto neighbor_score_bins =
      78          12 :       _openmc_problem->getTallyScoreNeighborVariableValues(score, tally_name, _tid);
      79             : 
      80          12 :   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          10 :   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          10 :   _tally_val = score_bins[_bin_index];
      95          10 :   _tally_neighbor_val = neighbor_score_bins[_bin_index];
      96          20 : }
      97             : 
      98             : void
      99        1536 : FDTallyGradAux::compute()
     100             : {
     101        1536 :   auto elem_c = _current_elem->true_centroid();
     102       10752 :   for (unsigned int side = 0; side < _current_elem->n_sides(); side++)
     103             :   {
     104        9216 :     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        9216 :     if (!neighbor)
     110         576 :       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        8640 :     if (!hasBlocks(neighbor->subdomain_id()))
     118           0 :       continue;
     119             : 
     120        8640 :     _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        8640 :     auto y_prime = neighbor->true_centroid() - elem_c;
     125        8640 :     RealEigenVector y_prime_eig(3);
     126        8640 :     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        8640 :     _sum_y_du_dy +=
     131        8640 :         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        8640 :     _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        3072 :   RealEigenVector val = _sum_y_y_t.fullPivLu().solve(_sum_y_du_dy);
     139             :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     140        1536 :   _var.setDofValue(val(0), 0);
     141        1536 :   _var.setDofValue(val(1), 1);
     142        1536 :   _var.setDofValue(val(2), 2);
     143             : 
     144        1536 :   _sum_y_y_t = RealEigenMatrix::Zero(3, 3);
     145        1536 :   _sum_y_du_dy = RealEigenVector::Zero(3);
     146        1536 : }
     147             : 
     148             : #endif