Merge pull request #273 from smpark7/vacuum-bc

Add Mark and Milne vacuum BCs

doc/content/bib/references.bib

@@ -0,0 +1,36 @@
+
+@article{placzek_milnes_1947,
+	title = {Milne's {Problem} in {Transport} {Theory}},
+	volume = {72},
+	url = {https://link.aps.org/doi/10.1103/PhysRev.72.550},
+	doi = {10.1103/PhysRev.72.550},
+	abstract = {A modified derivation of the Wiener-Hopf solution of Milne's problem is given in a form suitable for application to problems in the theory of neutron diffusion.},
+	number = {7},
+	urldate = {2024-05-31},
+	journal = {Physical Review},
+	author = {Placzek, G. and Seidel, W.},
+	month = oct,
+	year = {1947},
+	note = {Publisher: American Physical Society},
+	pages = {550--555},
+	file = {APS Snapshot:C\:\\Users\\Sun Myung\\Zotero\\storage\\IJVWPMDS\\PhysRev.72.html:text/html;Full Text PDF:C\:\\Users\\Sun Myung\\Zotero\\storage\\TVUWZQ6J\\Placzek and Seidel - 1947 - Milne's Problem in Transport Theory.pdf:application/pdf},
+}
+
+@article{rulko_variational_1995,
+	title = {Variational {P1} {Approximations} of {General}-{Geometry} {Multigroup} {Transport} {Problems}},
+	volume = {121},
+	issn = {0029-5639},
+	url = {https://doi.org/10.13182/NSE121-393},
+	doi = {10.13182/NSE121-393},
+	abstract = {A variational approximation is developed for general-geometry multigroup transport problems with arbitrary anisotropic scattering. The variational principle is based on a functional that approximates a reaction rate in a subdomain of the system. In principle, approximations that result from this functional “optimally”determine such reaction rates. The functional contains an arbitrary parameter α and requires the approximate solutions of a forward and an adjoint transport problem. If the basis functions for the forward and adjoint solutions are chosen to be linear functions of the angular variable Ω, the functional yields the familiar multigroup P1 equations for all values of α. However, the boundary conditions that result from the functional depend on α. In particular, for problems with vacuum boundaries, one obtains the conventional mixed boundary condition, but with an extrapolation distance that depends continuously on α. The choice α = 0 yields a generalization of boundary conditions derived earlier by Federighi and Pomraning for a more limited class of problems. The choice α = 1 yields a generalization of boundary conditions derived previously by Davis for mono-energetic problems. Other boundary conditions are obtained by choosing different values of α. We discuss this indeterminancy of a in conjunction with numerical experiments.},
+	number = {3},
+	urldate = {2024-05-31},
+	journal = {Nuclear Science and Engineering},
+	author = {Rulko, Robert P. and Tomašević, Djordje and Larsen, Edward W.},
+	month = dec,
+	year = {1995},
+	note = {Publisher: Taylor \& Francis
+\_eprint: https://doi.org/10.13182/NSE121-393},
+	pages = {393--404},
+	file = {Full Text PDF:C\:\\Users\\Sun Myung\\Zotero\\storage\\8NHDX7DM\\Rulko et al. - 1995 - Variational P1 Approximations of General-Geometry .pdf:application/pdf},
+}

doc/content/source/bcs/VacuumConcBC.md

@@ -4,12 +4,25 @@
 
 ## Overview
 
-This object adds the $\frac{\phi}{4}-\frac{D_g}{2}\hat{n}\cdot\nabla\phi_g = 0$ vacuum boundary
+This object adds the $\phi_g+\alpha D_g\hat{n}\cdot\nabla\phi_g = 0$ vacuum boundary
 condition of the multigroup neutron diffusion equations. The weak form after applying integration
 by parts to the neutron diffusion term ([GroupDiffusion](/GroupDiffusion.md)) is:
 
 !equation
-\int_{\partial V}\psi D_g\nabla\phi_g\cdot\hat{n}\ dS = \int_{\partial V}\psi\frac{\phi_g}{2}\ dS
+-\int_{\partial V}\psi D_g\nabla\phi_g\cdot\hat{n}\ dS =
+\int_{\partial V}\psi\frac{\phi_g}{\alpha}\ dS
+
+The value of $\alpha$ varies depending on the BC type selected using the `vacuum_bc_type`
+parameter. The available options are:
+
+1. Marshak: $\alpha=2$
+2. Mark: $\alpha=\sqrt{3}$
+3. Milne: $\alpha=3\times 0.710446$
+
+The Marshak and Mark BCs are derived from vacuum boundary condition approximations using
+$P_1$ and $S_2$ transport methods, respectively. The Milne BC is derived from the exact analytical
+solution of the Milne problem [!citep](placzek_milnes_1947) and has been proven to be accurate for
+a wide range of diffusive problems [!citep](rulko_variational_1995).
 
 ## Example Input File Syntax
 

include/bcs/VacuumConcBC.h

@@ -4,10 +4,13 @@
 #include "ScalarTransportBase.h"
 
 /**
- * Implements a simple VacuumConc BC for neutron diffusion on the boundary.
- * VacuumConc BC is defined as \f$ D\frac{du}{dn}+\frac{u}{2} = 0\f$, where u is neutron flux.
- * Hence, \f$ D\frac{du}{dn}=-\frac{u}{2} \f$ and \f$ -\frac{u}{2} \f$ is substituted into
- * the Neumann BC term produced from integrating the diffusion operator by parts.
+ * Implements a simple VacuumConc BC for neutron diffusion on vacuum boundaries.
+ * VacuumConc BC is defined as \f$ D\frac{du}{dn}+\frac{u}{\alpha} = 0\f$, where u is neutron flux.
+ * Hence, \f$ D\frac{du}{dn}=-\frac{u}{\alpha} \f$ and \f$ -\frac{u}{\alpha} \f$ is substituted
+ * into the Neumann BC term produced from integrating the diffusion operator by parts.
+ *
+ * The three types of vacuum BCs available are Marshak ($\alpha=2$), Mark ($\alpha=\sqrt{3}$), and
+ * Milne ($\alpha=3\times 0.710446$). VacuumConcBC defaults to Marshak if `bc_type` is not set.
  */
 class VacuumConcBC : public IntegratedBC, public ScalarTransportBase
 {
@@ -22,9 +25,20 @@ class VacuumConcBC : public IntegratedBC, public ScalarTransportBase
 
 protected:
   virtual Real computeQpResidual() override;
-
   virtual Real computeQpJacobian() override;
 
-  /// Ratio of u to du/dn
+  enum BC_TYPE
+  {
+    MARSHAK,
+    MARK,
+    MILNE
+  };
+
+  // Ratio of u to du/dn
   Real _alpha;
+
+  // Milne vacuum boundary extrapolation coefficient
+  // Derived from the exact analytical solution to the Milne problem. See MooseDocs-based
+  // documentation for more information.
+  Real _milne_extrapolation_coefficient = 3 * 0.710446;
 };

src/actions/NtAction.C

@@ -27,7 +27,8 @@ NtAction::validParams()
 
   params.addRequiredParam<unsigned int>("num_precursor_groups",
                                         "specifies the total number of precursors to create");
-  params.addRequiredParam<std::string>("var_name_base", "specifies the base name of the variables");
+  params.addRequiredParam<std::string>("var_name_base",
+                                       "specifies the base name of the variables");
   params.addRequiredCoupledVar("temperature", "Name of temperature variable");
   params.addCoupledVar("pre_concs",
                        "All the variables that hold the precursor concentrations. "
@@ -42,8 +43,12 @@ NtAction::validParams()
                         "random initial conditions for the precursors.");
   params.addParam<FunctionName>("nt_ic_function",
                                 "An initial condition function for the neutrons.");
-  params.addParam<std::vector<BoundaryName>>("vacuum_boundaries",
-                                             "The boundaries on which to apply vacuum boundaries.");
+  params.addParam<std::vector<BoundaryName>>(
+      "vacuum_boundaries",
+      "The boundaries on which to apply vacuum boundaries.");
+  MooseEnum vacuum_bc_type("marshak mark milne", "marshak");
+  params.addParam<MooseEnum>("vacuum_bc_type", vacuum_bc_type,
+      "Whether to apply Marshak, Mark, or Milne vacuum boundary conditions. Defaults to Marshak.");
   params.addParam<bool>(
       "create_temperature_var", true, "Whether to create the temperature variable.");
   params.addParam<bool>(
@@ -147,6 +152,7 @@ NtAction::act()
         params.set<NonlinearVariableName>("variable") = var_name;
         if (isParamValid("use_exp_form"))
           params.set<bool>("use_exp_form") = getParam<bool>("use_exp_form");
+        params.set<MooseEnum>("vacuum_bc_type") = getParam<MooseEnum>("vacuum_bc_type");
         std::string bc_name = "VacuumConcBC_" + var_name;
         _problem->addBoundaryCondition("VacuumConcBC", bc_name, params);
       }

src/bcs/VacuumConcBC.C

@@ -7,22 +7,37 @@ VacuumConcBC::validParams()
 {
   InputParameters params = IntegratedBC::validParams();
   params += ScalarTransportBase::validParams();
+  MooseEnum vacuum_bc_type("marshak mark milne", "marshak");
+  params.addParam<MooseEnum>("vacuum_bc_type", vacuum_bc_type,
+      "Whether to apply Marshak, Mark, or Milne vacuum boundary conditions. Defaults to Marshak.");
   return params;
 }
 
 VacuumConcBC::VacuumConcBC(const InputParameters & parameters)
   : IntegratedBC(parameters), ScalarTransportBase(parameters)
 {
+  switch (getParam<MooseEnum>("vacuum_bc_type"))
+  {
+    case MARSHAK:
+      _alpha = 2.;
+      break;
+    case MARK:
+      _alpha = std::sqrt(3.);
+      break;
+    case MILNE:
+      _alpha = _milne_extrapolation_coefficient;
+      break;
+  }
 }
 
 Real
 VacuumConcBC::computeQpResidual()
 {
-  return _test[_i][_qp] * computeConcentration(_u, _qp) / 2.;
+  return _test[_i][_qp] * computeConcentration(_u, _qp) / _alpha;
 }
 
 Real
 VacuumConcBC::computeQpJacobian()
 {
-  return _test[_i][_qp] * computeConcentrationDerivative(_u, _phi, _j, _qp) / 2.;
+  return _test[_i][_qp] * computeConcentrationDerivative(_u, _phi, _j, _qp) / _alpha;
 }

tests/bcs/tests

@@ -1,7 +1,47 @@
 [Tests]
-  [./coupled_bc]
+  [coupled_bc]
     type = 'Exodiff'
     input = 'coupled_bc_test.i'
     exodiff = 'coupled_bc_test_out.e'
-  [../]
+    requirement = 'The system shall be able to impose PostprocessorCoupledInflowBC and
+                   CoupledOutflowBC for advection with flow from coupled velocity variables.'
+  []
+  [vacuum_bc]
+    type = 'Exodiff'
+    input = 'vacuum_bc.i'
+    exodiff = 'vacuum_bc_out.e'
+    requirement = 'The system shall impose Marshak vacuum boundary conditions by default if
+                   `vacuum_bc_type` is not set.'
+  []
+  [vacuum_bc_marshak]
+    type = 'Exodiff'
+    input = 'vacuum_bc.i'
+    # Left boundary value = 1 + 2 * 0.001 = 1.002
+    # Right boundary value = 2 * 0.001 = 0.002
+    cli_args = 'BCs/vacuum/vacuum_bc_type=marshak'
+    exodiff = 'vacuum_bc_out.e'
+    requirement = 'The system shall impose Marshak vacuum boundary conditions.'
+  []
+  [vacuum_bc_mark]
+    type = 'Exodiff'
+    input = 'vacuum_bc.i'
+    # Left boundary value = 1 + sqrt(3) * 0.001 = 1.0017320508
+    # Right boundary value = sqrt(3) * 0.001 = 0.0017320508
+    cli_args = 'BCs/vacuum/vacuum_bc_type=mark
+                BCs/dirichlet/value=1.0017320508
+                Outputs/out/file_base=vacuum_bc_mark_out'
+    exodiff = 'vacuum_bc_mark_out.e'
+    requirement = 'The system shall impose Mark vacuum boundary conditions.'
+  []
+  [vacuum_bc_milne]
+    type = 'Exodiff'
+    input = 'vacuum_bc.i'
+    # Left boundary value = 1 + 3 * 0.710446 * 0.001 = 1.002131338
+    # Right boundary value = 3 * 0.710446 * 0.001 = 0.002131338
+    cli_args = 'BCs/vacuum/vacuum_bc_type=milne
+                BCs/dirichlet/value=1.0017320508
+                Outputs/out/file_base=vacuum_bc_milne_out'
+    exodiff = 'vacuum_bc_milne_out.e'
+    requirement = 'The system shall impose Milne vacuum boundary conditions.'
+  []
 []

tests/bcs/vacuum_bc.i

@@ -0,0 +1,68 @@
+[GlobalParams]
+  use_exp_form = false
+[]
+
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 2
+  []
+[]
+
+[Variables]
+  [flux]
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = GroupDiffusion
+    variable = flux
+    group_number = 1
+  []
+[]
+
+[BCs]
+  [dirichlet]
+    type = DirichletBC
+    variable = flux
+    boundary = 'left'
+    value = 1.002
+  []
+  [vacuum]
+    type = VacuumConcBC
+    variable = flux
+    boundary = 'right'
+    vacuum_bc_type = marshak
+  []
+[]
+
+[Materials]
+  [mat]
+    type = MoltresJsonMaterial
+    base_file = 'xsdata.json'
+    material_key = 0
+    interp_type = none
+    num_groups = 1
+    num_precursor_groups = 1
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = PJFNK
+[]
+
+[Postprocessors]
+  [boundary_value]
+    type = SideAverageValue
+    variable = flux
+    boundary = 'right'
+  []
+[]
+
+[Outputs]
+  [out]
+    type = Exodus
+  []
+[]

tests/bcs/xsdata.json

@@ -0,0 +1,45 @@
+{
+    "0": {
+        "900": {
+            "BETA_EFF": [
+                0.0
+            ],
+            "CHI_D": [
+                0.0
+            ],
+            "CHI_P": [
+                0.0
+            ],
+            "CHI_T": [
+                0.0
+            ],
+            "DECAY_CONSTANT": [
+                0.0
+            ],
+            "DIFFCOEF": [
+                1e-3
+            ],
+            "FISSE": [
+                0.0
+            ],
+            "FISSXS": [
+                0.0
+            ],
+            "GTRANSFXS": [
+                0.0
+            ],
+            "NSF": [
+                0.0
+            ],
+            "RECIPVEL": [
+                0.0
+            ],
+            "REMXS": [
+                0.0
+            ]
+        },
+        "temp": [
+            900
+        ]
+    }
+}