Merge pull request #249 from smpark7/fix-eigenvalue-scaling

Fix DNP source eigenvalue scaling and InScatter Jacobian formulation

include/actions/PrecursorAction.h

@@ -126,4 +126,6 @@ class PrecursorAction : public VariableNotAMooseObjectAction
 
   /// optional object name suffix
   std::string _object_suffix;
+
+  bool _is_loopapp;
 };

include/kernels/DelayedNeutronSource.h

@@ -24,7 +24,6 @@ class DelayedNeutronSource : public Kernel, public ScalarTransportBase
   unsigned int _num_precursor_groups;
   unsigned int _temp_id;
   const VariableValue & _temp;
-  Real _eigenvalue_scaling;
   std::vector<const VariableValue *> _pre_concs;
   std::vector<unsigned int> _pre_ids;
 };

include/kernels/PrecursorSource.h

@@ -26,4 +26,5 @@ class PrecursorSource : public Kernel, public ScalarTransportBase
   std::vector<const VariableValue *> _group_fluxes;
   std::vector<unsigned int> _flux_ids;
   Real _prec_scale;
+  Real _eigenvalue_scaling;
 };

problems/2021-cnrs-benchmark/phase-2/transient.i

@@ -20,6 +20,9 @@ dt = 0.00625        # Timestep size = 1 / freq / 200 = 0.00625 s
   sss2_input = true
   account_delayed = true
   integrate_p_by_parts = true
+  eigenvalue_scaling = 0.9927821802
+  ## Use the eigenvalue scaling factor below if running on a 40x40 mesh
+  # eigenvalue_scaling = 0.9926551482
 []
 
 [Mesh]
@@ -60,9 +63,6 @@ dt = 0.00625        # Timestep size = 1 / freq / 200 = 0.00625 s
   vacuum_boundaries = 'bottom left right top'
   create_temperature_var = false
   init_nts_from_file = true
-  eigenvalue_scaling = 0.9927821802
-  ## Use the eigenvalue scaling factor below if running on a 40x40 mesh
-  #  eigenvalue_scaling = 0.9926551482
 []
 
 [Precursors]

src/actions/NtAction.C

@@ -12,19 +12,12 @@
 #include "libmesh/enum_to_string.h"
 
 registerMooseAction("MoltresApp", NtAction, "add_kernel");
-
 registerMooseAction("MoltresApp", NtAction, "add_bc");
-
 registerMooseAction("MoltresApp", NtAction, "add_variable");
-
 registerMooseAction("MoltresApp", NtAction, "add_ic");
-
 registerMooseAction("MoltresApp", NtAction, "add_aux_variable");
-
 registerMooseAction("MoltresApp", NtAction, "add_aux_kernel");
-
 registerMooseAction("MoltresApp", NtAction, "check_copy_nodal_vars");
-
 registerMooseAction("MoltresApp", NtAction, "copy_nodal_vars");
 
 InputParameters
@@ -73,10 +66,9 @@ NtAction::validParams()
   params.addParam<std::vector<SubdomainName>>("pre_blocks", "The blocks the precursors live on.");
   params.addParam<Real>("eigenvalue_scaling",
                         1.0,
-                        "Artificial scaling factor for the fission "
-                        "source. Primarily introduced to make "
-                        "super/sub-critical systems exactly critical "
-                        "for the CNRS benchmark.");
+                        "Artificial scaling factor for the fission source. Primarily for "
+                        "introducing artificial reactivity to make super/subcritical systems "
+                        "exactly critical or to simulate reactivity insertions/withdrawals.");
   return params;
 }
 
@@ -262,7 +254,6 @@ NtAction::act()
         std::vector<std::string> include = {"temperature", "pre_concs"};
         params.applySpecificParameters(parameters(), include);
         params.set<unsigned int>("num_precursor_groups") = _num_precursor_groups;
-        params.set<Real>("eigenvalue_scaling") = getParam<Real>("eigenvalue_scaling");
 
         std::string kernel_name = "DelayedNeutronSource_" + var_name;
         _problem->addKernel("DelayedNeutronSource", kernel_name, params);

src/actions/PrecursorAction.C

@@ -52,6 +52,11 @@ PrecursorAction::validParams()
   params.addParam<bool>("nt_exp_form",
                         false,
                         "Whether concentrations should be in an expotential/logarithmic format.");
+  params.addParam<Real>("eigenvalue_scaling",
+                        1.0,
+                        "Artificial scaling factor for the fission source. Primarily for "
+                        "introducing artificial reactivity to make super/subcritical systems "
+                        "exactly critical or to simulate reactivity insertions/withdrawals.");
   params.addParam<bool>("jac_test",
                         false,
                         "Whether we're testing the Jacobian and should use some "
@@ -72,7 +77,7 @@ PrecursorAction::validParams()
                                               "from the base list of blocks. Replaces the 'block'"
                                               "parameter when initializing kernels.");
   params.addParam<MultiAppName>("multi_app", "Multiapp name for looping precursors.");
-  params.addParam<bool>("is_loopapp", "if circulating precursors, whether this is loop app");
+  params.addParam<bool>("is_loopapp", false, "if circulating precursors, whether this is loop app");
   params.addParam<bool>("eigen", false, "whether neutronics is in eigenvalue calculation mode");
   params.addParam<NonlinearVariableName>("outlet_vel",
                                          "Name of the velocity variable normal to the "
@@ -87,7 +92,8 @@ PrecursorAction::PrecursorAction(const InputParameters & params)
     _num_precursor_groups(getParam<unsigned int>("num_precursor_groups")),
     _var_name_base(getParam<std::string>("var_name_base")),
     _num_groups(getParam<unsigned int>("num_groups")),
-    _object_suffix(getParam<std::string>("object_suffix"))
+    _object_suffix(getParam<std::string>("object_suffix")),
+    _is_loopapp(getParam<bool>("is_loopapp"))
 {
   if (getParam<bool>("loop_precursors"))
   {
@@ -170,7 +176,7 @@ PrecursorAction::act()
 
     // transfers
     else if (_current_task == "add_transfer" && getParam<bool>("loop_precursors") &&
-        !getParam<bool>("is_loopapp"))
+        (!_is_loopapp))
     {
       // Set up MultiAppTransfer to simulate precursor looped flow into and
       // out of the reactor core
@@ -181,7 +187,7 @@ PrecursorAction::act()
   // Add outflow rate postprocessor for Navier-Stokes velocities in the main
   // app if precursors are looped
   if (_current_task == "add_postprocessor" && getParam<bool>("loop_precursors") &&
-      isParamValid("uvel") && !getParam<bool>("is_loopapp"))
+      isParamValid("uvel") && (!_is_loopapp))
     addCoolantOutflowPostprocessor();
 }
 
@@ -210,6 +216,7 @@ PrecursorAction::addPrecursorSource(const unsigned & op, const std::string & var
     std::vector<std::string> include = {"temperature", "group_fluxes"};
     params.applySpecificParameters(parameters(), include);
     params.set<bool>("use_exp_form") = getParam<bool>("nt_exp_form");
+    params.set<Real>("eigenvalue_scaling") = getParam<Real>("eigenvalue_scaling");
 
     std::string kernel_name = "PrecursorSource_" + var_name + "_" + _object_suffix;
     _problem->addKernel("PrecursorSource", kernel_name, params);

src/kernels/DelayedNeutronSource.C

@@ -13,10 +13,6 @@ DelayedNeutronSource::validParams()
                                             "concentrations. These MUST be listed by increasing "
                                             "group number.");
   params.addRequiredParam<unsigned int>("group_number","neutron energy group number for chi_d");
-  params.addParam<Real>("eigenvalue_scaling", 1.0, "Artificial scaling factor for the fission "
-                                                   "source. Primarily introduced to make "
-                                                   "super/sub-critical systems exactly critical "
-                                                   "for the CNRS benchmark.");
   return params;
 }
 
@@ -25,12 +21,11 @@ DelayedNeutronSource::DelayedNeutronSource(const InputParameters & parameters)
     ScalarTransportBase(parameters),
     _decay_constant(getMaterialProperty<std::vector<Real>>("decay_constant")),
     _d_decay_constant_d_temp(getMaterialProperty<std::vector<Real>>("d_decay_constant_d_temp")),
-    _group(getParam<unsigned int>("group_number")),
+    _group(getParam<unsigned int>("group_number") - 1),
     _chi_d(getMaterialProperty<std::vector<Real>>("chi_d")),
     _num_precursor_groups(getParam<unsigned int>("num_precursor_groups")),
     _temp_id(coupled("temperature")),
-    _temp(coupledValue("temperature")),
-    _eigenvalue_scaling(getParam<Real>("eigenvalue_scaling"))
+    _temp(coupledValue("temperature"))
 {
   unsigned int n = coupledComponents("pre_concs");
   if (!(n == _num_precursor_groups))
@@ -53,10 +48,7 @@ DelayedNeutronSource::computeQpResidual()
   for (unsigned int i = 0; i < _num_precursor_groups; ++i)
     r += -_decay_constant[_qp][i] * computeConcentration((*_pre_concs[i]), _qp);
 
-  if ((_eigenvalue_scaling != 1.0))
-    r /= _eigenvalue_scaling;
-
-  return _chi_d[_qp][_group-1] * _test[_i][_qp] * r;
+  return _chi_d[_qp][_group] * _test[_i][_qp] * r;
 }
 
 Real
@@ -84,8 +76,5 @@ DelayedNeutronSource::computeQpOffDiagJacobian(unsigned int jvar)
       jac += -_test[_i][_qp] * computeConcentration((*_pre_concs[i]), _qp) *
              _d_decay_constant_d_temp[_qp][i] * _phi[_j][_qp];
 
-  if ((_eigenvalue_scaling != 1.0))
-    jac /= _eigenvalue_scaling;
-
-  return _chi_d[_qp][_group-1] * jac;
+  return _chi_d[_qp][_group] * jac;
 }

src/kernels/InScatter.C

@@ -73,7 +73,7 @@ InScatter::computeQpOffDiagJacobian(unsigned int jvar)
   Real jac = 0;
   for (unsigned int i = 0; i < _num_groups; ++i)
   {
-    if (jvar == _flux_ids[i] && jvar != _group)
+    if (jvar == _flux_ids[i])
     {
       if (_sss2_input)
         jac += -_test[_i][_qp] * _gtransfxs[_qp][i * _num_groups + _group] *

src/kernels/PrecursorSource.C

@@ -16,6 +16,11 @@ PrecursorSource::validParams()
   params.addCoupledVar(
       "temperature", 800, "The temperature used to interpolate material properties.");
   params.addParam<Real>("prec_scale", 1, "The factor by which the neutron fluxes are scaled.");
+  params.addParam<Real>("eigenvalue_scaling",
+                        1.0,
+                        "Artificial scaling factor for the fission source. Primarily for "
+                        "introducing artificial reactivity to make super/subcritical systems "
+                        "exactly critical or to simulate reactivity insertions/withdrawals.");
   return params;
 }
 
@@ -30,7 +35,8 @@ PrecursorSource::PrecursorSource(const InputParameters & parameters)
     _precursor_group(getParam<unsigned int>("precursor_group_number") - 1),
     _temp(coupledValue("temperature")),
     _temp_id(coupled("temperature")),
-    _prec_scale(getParam<Real>("prec_scale"))
+    _prec_scale(getParam<Real>("prec_scale")),
+    _eigenvalue_scaling(getParam<Real>("eigenvalue_scaling"))
 {
   _group_fluxes.resize(_num_groups);
   _flux_ids.resize(_num_groups);
@@ -51,6 +57,9 @@ PrecursorSource::computeQpResidual()
          computeConcentration((*_group_fluxes[i]), _qp) * _prec_scale;
   }
 
+  if ((_eigenvalue_scaling != 1.0))
+    r /= _eigenvalue_scaling;
+
   return r;
 }
 
@@ -69,7 +78,7 @@ PrecursorSource::computeQpOffDiagJacobian(unsigned int jvar)
     {
       jac = -_test[_i][_qp] * _beta_eff[_qp][_precursor_group] * _nsf[_qp][i] *
             computeConcentrationDerivative((*_group_fluxes[i]), _phi, _j, _qp) * _prec_scale;
-      return jac;
+      break;
     }
 
   if (jvar == _temp_id)
@@ -80,8 +89,10 @@ PrecursorSource::computeQpOffDiagJacobian(unsigned int jvar)
                   computeConcentration((*_group_fluxes[i]), _qp) * _prec_scale +
               _d_beta_eff_d_temp[_qp][_precursor_group] * _phi[_j][_qp] * _nsf[_qp][i] *
                   computeConcentration((*_group_fluxes[i]), _qp) * _prec_scale);
-    return jac;
   }
 
-  return 0;
+  if ((_eigenvalue_scaling != 1.0))
+    jac /= _eigenvalue_scaling;
+
+  return jac;
 }