Skip to content

Commit

Permalink
Add FDTallyGradAux.
Browse files Browse the repository at this point in the history 
Ref. neams-th-coe#1031
  • Loading branch information
@nuclearkevin
nuclearkevin committed Jan 26, 2025
1 parent 5f39b93 commit 98651e8
Show file tree
Hide file tree
Showing 2 changed files with 198 additions and 0 deletions.
63 changes: 63 additions & 0 deletions include/auxkernels/FDTallyGradAux.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
/********************************************************************/
/* SOFTWARE COPYRIGHT NOTIFICATION */
/* Cardinal */
/* */
/* (c) 2021 UChicago Argonne, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by UChicago Argonne, LLC */
/* Under Contract No. DE-AC02-06CH11357 */
/* With the U. S. Department of Energy */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See LICENSE for full restrictions */
/********************************************************************/

#pragma once

#include "OpenMCAuxKernel.h"

/**
* A class which approximates gradients of constant monomial tallies using forward finite
* differences. The gradient computation is based on:
* K. N. Stolte and P. V. Tsvetkov (2023), Annals of Nuclear Energy, 182, 109617.
* https://doi.org/10.1016/j.anucene.2022.109617
*/
class FDTallyGradAux : public OpenMCVectorAuxKernel
{
public:
static InputParameters validParams();

FDTallyGradAux(const InputParameters & parameters);

/// We handle computing and storing the variable value manually.
virtual void compute() override;

protected:
/// Need to override computeValue() to avoid creating a pure-virtual class.
virtual RealVectorValue computeValue() override { return RealVectorValue(0.0, 0.0, 0.0); }

/// The external filter bin index for the score.
const unsigned int _bin_index;

/// The element's tally value.
const VariableValue * _tally_val;

/// The neighboring element's tally value.
const VariableValue * _tally_neighbor_val;

/**
* The sum of outer products of y' with itself, where y' = x_j - x_i.
* x_j is the centroid of the neighboring element j of element i.
* x_i is the centroid of element i.
*/
RealEigenMatrix _sum_y_y_t;

/**
* The sum of the finite difference approximation of u multiplied by y'.
*/
RealEigenVector _sum_y_du_dy;
};
135 changes: 135 additions & 0 deletions src/auxkernels/FDTallyGradAux.C
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
/********************************************************************/
/* SOFTWARE COPYRIGHT NOTIFICATION */
/* Cardinal */
/* */
/* (c) 2021 UChicago Argonne, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by UChicago Argonne, LLC */
/* Under Contract No. DE-AC02-06CH11357 */
/* With the U. S. Department of Energy */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See LICENSE for full restrictions */
/********************************************************************/

#ifdef ENABLE_OPENMC_COUPLING

#include "FDTallyGradAux.h"

#include "CardinalEnums.h"

registerMooseObject("CardinalApp", FDTallyGradAux);

InputParameters
FDTallyGradAux::validParams()
{
auto params = OpenMCVectorAuxKernel::validParams();
params.addClassDescription(
"An auxkernel which approximates tally gradients at element centroids using "
"forward finite differences.");
params.addRequiredParam<MooseEnum>("score",
getSingleTallyScoreEnum(), "The tally score this auxkernel should approximate the gradient of.");
params.addParam<unsigned int>("ext_filter_bin", 0, "The non-spatial filter bin for the tally score.");

params.addRelationshipManager("ElementSideNeighborLayers",
Moose::RelationshipManagerType::ALGEBRAIC |
Moose::RelationshipManagerType::GEOMETRIC,
[](const InputParameters &, InputParameters & rm_params)
{ rm_params.set<unsigned short>("layers") = 2; });

return params;
}

FDTallyGradAux::FDTallyGradAux(const InputParameters & parameters)
: OpenMCVectorAuxKernel(parameters),
_bin_index(getParam<unsigned int>("ext_filter_bin")),
_sum_y_y_t(RealEigenMatrix::Zero(3, 3)),
_sum_y_du_dy(RealEigenVector::Zero(3))
{
if (_var.feType() != FEType(CONSTANT, MONOMIAL_VEC))
paramError(
"variable",
"FDTallyGradAux only supports CONSTANT MONOMIAL_VEC shape functions. Please "
"ensure that 'variable' is of type CONSTANT MONOMIAL_VEC.");

std::string score = getParam<MooseEnum>("score");
std::replace(score.begin(), score.end(), '_', '-');

// Error check and fetch the tally score.
const auto & all_scores = _openmc_problem->getTallyScores();
if (std::find(all_scores.begin(), all_scores.end(), score) == all_scores.end())
paramError("score", "The problem does not contain any score named " +
std::string(getParam<MooseEnum>("score")) +
"! Please "
"ensure that one of your [Tallies] is scoring the requested score.");

auto score_vars = getTallyScoreVariables(score);
auto score_bins = getTallyScoreVariableValues(score);
auto neighbor_score_bins = getTallyScoreNeighborVariableValues(score);

if (_bin_index >= score_bins.size())
paramError("ext_filter_bin",
"The external filter bin for the score " + std::string(getParam<MooseEnum>("score")) +
" is larger than the number of external filter bins!");

if (score_vars[_bin_index]->feType() != FEType(CONSTANT, MONOMIAL))
paramError("score", "FDTallyGradAux only supports CONSTANT MONOMIAL shape functions for tally variables.");

_tally_val = score_bins[_bin_index];
_tally_neighbor_val = neighbor_score_bins[_bin_index];
}

void
FDTallyGradAux::compute()
{
auto elem_c = _current_elem->true_centroid();
for (unsigned int side = 0; side < _current_elem->n_sides(); side++)
{
const Elem * neighbor = _current_elem->neighbor_ptr(side);

// If the neighbor is null, the current side is a boundary and we can skip
// it. Geometric ghosting ensures that element neighbors are on the
// processors that need them during compute().
if (!neighbor)
continue;

if (!neighbor->active())
continue;

// Can skip the gradient computation for this side if the neighbor isn't on
// this block.
if (!hasBlocks(neighbor->subdomain_id()))
continue;

_openmc_problem->reinitNeighbor(_current_elem, side, _tid);

// Fetch the vector pointing from the current element's centroid to the
// neighbor's centroid (y').
auto y_prime = neighbor->true_centroid() - elem_c;
RealEigenVector y_prime_eig(3);
y_prime_eig << y_prime(0), y_prime(1), y_prime(2);

// Compute du/dy along the direction pointing towards the neighbor's centroid.
// Add to the b vector.
_sum_y_du_dy += y_prime_eig * ((*_tally_neighbor_val)[0] - (*_tally_val)[0]) / y_prime.norm_sq();
// Compute the outer product between y' and y'.T.
// Add to the A matrix.
_sum_y_y_t += y_prime_eig * y_prime_eig.transpose();
}

// Solve for an estimate of the tally gradient at the current element's centroid.
RealEigenVector val = _sum_y_y_t.fullPivLu().solve(_sum_y_du_dy);
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
_var.setDofValue(val(0), 0);
_var.setDofValue(val(1), 1);
_var.setDofValue(val(2), 2);

_sum_y_y_t = RealEigenMatrix::Zero(3, 3);
_sum_y_du_dy = RealEigenVector::Zero(3);
}

#endif

0 comments on commit 98651e8

Please sign in to comment.