diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md new file mode 100644 index 0000000000..ce48accd20 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -0,0 +1,31 @@ +--- +name: Bug report +about: Report software bugs +title: "[Bug]: " +labels: Type:Bug +assignees: '' + +--- + +## Expected behavior +A clear and concise description of what you expected to happen. + +## Describe the bug +A clear and concise description of what the bug is. + +## To Reproduce +Steps to reproduce the behavior: +1. +2. +3. +4. + +## Screenshots/output +If applicable, add screenshots or program output to help explain your problem. + +## System Specifications: + - Version: + - Platform: + - Subsystem: + +## How can this issue be closed? diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md new file mode 100644 index 0000000000..b4eff87996 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/feature_request.md @@ -0,0 +1,16 @@ +--- +name: Feature request +about: Suggest an idea for this project +title: 'Feature: ' +labels: Status:1-New, Type:Feature +assignees: '' + +--- + +## Background and motivation + +## Description of idea + +## Implementation details + +## Potential snags diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000000..811d9c77a4 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,17 @@ +version: 2 +updates: + - package-ecosystem: "gitsubmodule" + directory: "/" + allow: + - dependency-name: "moose" + schedule: + interval: "monthly" + commit-message: + prefix: "auto_update_moose" + open-pull-requests-limit: 1 + labels: + - "Comp:Core" + - "Priority:1-Critical" + reviewers: + - "smpark7" + - "nglaser3" diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 0000000000..40ec999d2e --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,38 @@ +## Summary of changes + + + + +## Types of changes + + +- [ ] Bug fix (non-breaking change which fixes an issue) +- [ ] New feature (non-breaking change which adds functionality) +- [ ] Breaking change (fix or feature that would cause existing functionality to change) + +## Required for Merging +- [ ] I have read the **CONTRIBUTING** document. +- [ ] My code follows the code style of this project. +- [ ] My change requires a change to the documentation. +- [ ] I have updated the documentation accordingly. +- [ ] I have added tests to cover my changes. +- [ ] All new and existing tests passed. + - [ ] CI tests pass + - [ ] Local tests pass (including Serpent2 integration tests) + +## Associated Issues and PRs + + +- Issue: # + + +## Associated Developers + + +- Dev: @ + + +## Checklist for Reviewers + +Reviewers should use [this link](https://arfc.github.io/manual/guides/pull_requests) to get to the +Review Checklist before they begin their review. diff --git a/.gitmodules b/.gitmodules index a733c15e8c..e47e28d554 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,3 +4,4 @@ [submodule "moose"] path = moose url = ../../idaholab/moose + branch = devel \ No newline at end of file diff --git a/doc/content/bib/references.bib b/doc/content/bib/references.bib new file mode 100644 index 0000000000..1644211776 --- /dev/null +++ b/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}, +} diff --git a/doc/content/getting_started/theory.md b/doc/content/getting_started/theory.md index 154adff0e4..2a869e5011 100644 --- a/doc/content/getting_started/theory.md +++ b/doc/content/getting_started/theory.md @@ -16,7 +16,46 @@ the MOOSE finite element framework, enabling highly flexible and scalable reacto ## Multigroup Neutron Diffusion -More information to be added in future. +The neutron diffusion equation is an approximation to the Boltzmann transport equation, and is derived by taking the zeroth and first moment with respect to $\hat\Omega$, the direction of neutron travel. +A full derivation of the diffusion equation will not be presented here for conciseness. Importantly, when taking the first moment of the transport equation the angular flux is assumed to be linearly anisotropic. +The approximations reduce the phase space through the elimination of angular dependence, but the resulting equation also has reduced fidelity when compared to the transport equation, particularly in regions where the neutron flux has strong angular dependence. +A non-exhaustive list of regions in which the neutron flux strongly depends on angle is near material interfaces between materials with highly dissimilar neutronic properties, within strong absorbers, and within near-void regions. +As a deterministic method, the neutron diffusion method also requires discretization of the continuous energy dependence into energy groups consisting of non-overlapping, finite energy ranges across the entire energy spectrum. +This energy discretization creates a system of equations referred to as the multigroup neutron diffusion equations: + +!equation +\frac{1}{v_g} \frac{\partial\phi_g}{\partial t}-\nabla \cdot D_g \nabla \phi_g +\Sigma^R_{g} \phi_g =\sum^G_{g \neq g'} {\Sigma^s_{g'\rightarrow g} \phi_{g'}}+ \chi^p_{g} \sum^G_{g'=1} {\left(1- \beta \right) \nu_{g'} \Sigma^f_{g'}\phi_{g'} }+\chi^d_{g} \sum^I_i {\lambda_i C_i} + +The delayed neutron precursor distributions are governed by: + +!equation +\frac{\partial C_i}{\partial t}=\sum^G_{g'=1}{\beta_i \nu \Sigma^f_{g'}\phi_{g'}}-\lambda_i C_i-\vec{u} \cdot \nabla C_i + +Notably, there are two production terms of neutrons, the prompt fission source and delayed neutron precursor decay source. +The first term describes the neutrons immediately born from fission, and the second term describes the neutrons born from the radioactive decay of neutron-emitting radionuclides, commonly called delayed neutron precursors. +The multigroup neutron diffusion equations are generally impossible to solve analytically for realistic problems and are therefore typically solved with numerical methods. +Moltres utilizes the [Finite Element Method](https://mooseframework.inl.gov/help/faq/what_is_fem.html) (FEM) capabilities provided by the MOOSE framework. +In FEM, the spatial domain is discretized into finite mesh elements and the time dependence is modeled through discrete time steps and time integration methods. + +!table id=diff_table caption= Terms in multi-group neutron diffusion equation and their associated kernels +|Term in Multi-Group Diffusion Equation | Associated Kernel |Definition of Term| +| - | - | - | +| [!eq](\frac{1}{v_g} \frac{\partial\phi_g}{\partial t}) | [NtTimeDerivative](NtTimeDerivative.md) | Time rate of change of energy group g | +| [!eq](- \nabla \cdot D_g \nabla \phi_g) | [GroupDiffusion](GroupDiffusion.md) | Streaming term of energy group g | +| [!eq](\Sigma^R_{g} \phi_g) | [SigmaR](SigmaR.md) | Removal from energy group g | +| [!eq](\sum^G_{g \neq g'} {\Sigma^s_{g'\rightarrow g} \phi_{g'}}) | [InScatter](InScatter.md) | In-scattering into energy group g | +| [!eq](\chi^p_{g} \sum^G_{g'=1} {\left[1- \beta \right] \nu_{g'} \Sigma^f_{g'}\phi_{g'} } ) | [CoupledFissionKernel](CoupledFissionKernel.md) | Prompt fission neutron source | +| [!eq](\chi^d_{g} \sum^I_i {\lambda_i C_i}) | [DelayedNeutronSource](DelayedNeutronSource.md) | Delayed fission neutron source | + +!table id=prec_table caption= Terms in delayed neutron precursor equation and their associated kernels +| Term in Delayed Precursor Equation | Associated Kernel(s) | Definition of Term | +| - | - | - | +| [!eq](\frac{\partial C_i}{\partial t}) | [ScalarTransportTimeDerivative](ScalarTransportTimeDerivative.md) | Time rate of change of precursor population | +| [!eq](\sum^G_{g'=1}{\beta_i \nu \Sigma^f_{g'}\phi_{g'}}) | [PrecursorSource](PrecursorSource.md) | Production of precursor from fission +| [!eq](-\lambda_i C_i) | [PrecursorDecay](PrecursorDecay.md) | Loss of precusors due to radioactive decay | +| [!eq](-\vec{u} \cdot \nabla C_i) | [DGCoupledAdvection](DGCoupledAdvection.md), [DGFunctionConvection](DGFunctionConvection.md), [DGConvection](https://mooseframework.inl.gov/source/dgkernels/DGConvection.html) | Advection of the precursors | + +For the advective term in the delayed neutron precursor equation, $-\vec{u}\cdot \nabla C_i$, there are three kernels. DGCoupledAdvection is used when the velocity is a variable in the simulation, DGFunctionConvection is used when the velocity is a function, and DGConvection is for when the velocity is constant. ## Heat Transfer and Fluid Flow diff --git a/doc/content/source/actions/NtAction.md b/doc/content/source/actions/NtAction.md index c83824310f..0a47fd8c82 100644 --- a/doc/content/source/actions/NtAction.md +++ b/doc/content/source/actions/NtAction.md @@ -1,20 +1,42 @@ # NtAction -!alert construction title=Undocumented Action Class -The NtAction has not been documented. The content listed below should be used as a starting point for -documenting the class, which includes the typical automatic documentation associated with an Action; -however, what is contained is ultimately determined by what is necessary to make the documentation -clear for users. +!alert note +This action sets up ONLY the variables, kernels, and BCs for neutron diffusion, +NOT the variables, kernels, and BCs for precursor tracking. +To generate these, see [PrecursorAction](PrecursorAction.md). !syntax description /Nt/NtAction ## Overview -!! Replace these lines with information regarding the NtAction action. +```NtAction``` greatly reduces input file length and complexity by automatically +setting up variables, kernels, and BCs for the neutron diffusion equations. +When including only the required input parameters, this action constructs the +neutron group variables and their associated kernels: +[SigmaR](SigmaR.md), [GroupDiffusion](GroupDiffusion.md), +[InScatter](InScatter.md), and [CoupledFissionKernel](CoupledFissionKernel.md). +With ```account_delayed``` set to ```True```, +[DelayedNeutronSource](DelayedNeutronSource.md) is also constructed, +otherwise this kernel is not constructed and the simulation does not consider +the delayed neutron precursor effects on the neutron flux distributions. + +For more information regarding the use of ```NtAction``` please refer to the +tutorials located [here](tutorials.md), specifically the +Multiphysics Reactor +Simulations+ section. ## Example Input File Syntax -!! Describe and include an example of how to use the NtAction action. +An example input file without the ```NtAction```, showing only the portion +affected by the action: + +!listing tutorial/eigenvalue/nts.i + start=Variables end=Materials remove=Precursors + + +And then the same information created with the ```NtAction```: + +!listing tutorial/eigenvalue/nts-action.i + block=Nt !syntax description /Nt/NtAction diff --git a/doc/content/source/bcs/VacuumConcBC.md b/doc/content/source/bcs/VacuumConcBC.md index 8d20ba5d27..a2d13cd142 100644 --- a/doc/content/source/bcs/VacuumConcBC.md +++ b/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 diff --git a/include/bcs/VacuumConcBC.h b/include/bcs/VacuumConcBC.h index 7170ad1e44..0b92826f68 100644 --- a/include/bcs/VacuumConcBC.h +++ b/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; }; diff --git a/include/kernels/INSBoussinesqBodyForce.h b/include/kernels/INSBoussinesqBodyForce.h index 7d2f64b8e3..838d949f1b 100644 --- a/include/kernels/INSBoussinesqBodyForce.h +++ b/include/kernels/INSBoussinesqBodyForce.h @@ -19,7 +19,7 @@ class INSBoussinesqBodyForce : public Kernel protected: virtual Real computeQpResidual() override; virtual Real computeQpJacobian() override; - virtual Real computeQpOffDiagJacobian(unsigned jvar); + virtual Real computeQpOffDiagJacobian(unsigned jvar) override; // Parameters const VariableValue & _T; diff --git a/include/materials/SATauMaterial.h b/include/materials/SATauMaterial.h index 5d5cc328af..3e368b12f6 100644 --- a/include/materials/SATauMaterial.h +++ b/include/materials/SATauMaterial.h @@ -54,6 +54,150 @@ class SATauMaterialTempl : public T using T::_velocity; using T::_grad_velocity; using T::_tau; + using Coupleable::adCoupledValue; + using Coupleable::adCoupledGradient; + using Coupleable::coupledValue; + using Coupleable::adCoupledDot; }; typedef SATauMaterialTempl SATauMaterial; + +template +InputParameters +SATauMaterialTempl::validParams() +{ + InputParameters params = T::validParams(); + params.addClassDescription( + "This is the material class used to compute the stabilization parameter tau_viscosity " + "for the Spalart-Allmaras turbulent viscosity equation."); + params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); + params.addRequiredCoupledVar("wall_distance_var", "Wall distance aux variable name"); + params.addParam("use_ft2_term", false, "Whether to apply the f_t2 term in the equation"); + return params; +} + +template +SATauMaterialTempl::SATauMaterialTempl(const InputParameters & parameters) + : T(parameters), + _sigma(2.0/3.0), + _cb1(0.1355), + _cb2(0.622), + _kappa(0.41), + _cw1(_cb1 / _kappa / _kappa + (1 + _cb2) / _sigma), + _cw2(0.3), + _cw3(2.0), + _cv1(7.1), + _ct1(1.0), + _ct2(2.0), + _ct3(1.2), + _ct4(0.5), + _mu_tilde(adCoupledValue("mu_tilde")), + _grad_mu(adCoupledGradient("mu_tilde")), + _tau_visc(this->template declareADProperty("tau_viscosity")), + _wall_dist(coupledValue("wall_distance_var")), + _use_ft2_term(this->template getParam("use_ft2_term")), + _strain_mag(this->template declareADProperty("strain_mag")), + _convection_strong_residual(this->template + declareADProperty("convection_strong_residual")), + _destruction_strong_residual(this->template + declareADProperty("destruction_strong_residual")), + _diffusion_strong_residual(this->template + declareADProperty("diffusion_strong_residual")), + _source_strong_residual(this->template + declareADProperty("source_strong_residual")), + _time_strong_residual(this->template + declareADProperty("time_strong_residual")), + _visc_strong_residual(this->template + declareADProperty("visc_strong_residual")) +{ +} + +template +void +SATauMaterialTempl::subdomainSetup() +{ + T::subdomainSetup(); + + if (_has_transient) + _visc_dot = & adCoupledDot("mu_tilde"); + else + _visc_dot = nullptr; +} + +template +void +SATauMaterialTempl::computeQpProperties() +{ + T::computeQpProperties(); + + // Compute strain rate and vorticity magnitudes + _strain_mag[_qp] = 2.0 * Utility::pow<2>(_grad_velocity[_qp](0, 0)) + + 2.0 * Utility::pow<2>(_grad_velocity[_qp](1, 1)) + + 2.0 * Utility::pow<2>(_grad_velocity[_qp](2, 2)) + + Utility::pow<2>(_grad_velocity[_qp](0, 2) + _grad_velocity[_qp](2, 0)) + + Utility::pow<2>(_grad_velocity[_qp](0, 1) + _grad_velocity[_qp](1, 0)) + + Utility::pow<2>(_grad_velocity[_qp](1, 2) + _grad_velocity[_qp](2, 1)); + _strain_mag[_qp] = std::sqrt(_strain_mag[_qp] + 1e-16); + ADReal vorticity_mag = + Utility::pow<2>(_grad_velocity[_qp](0, 2) - _grad_velocity[_qp](2, 0)) + + Utility::pow<2>(_grad_velocity[_qp](0, 1) - _grad_velocity[_qp](1, 0)) + + Utility::pow<2>(_grad_velocity[_qp](1, 2) - _grad_velocity[_qp](2, 1)); + vorticity_mag = std::sqrt(vorticity_mag + 1e-16); + + // Compute relevant parameters for the SA equation + const Real d = std::max(_wall_dist[_qp], 1e-16); // Avoid potential division by zero + const ADReal chi = _mu_tilde[_qp] / _mu[_qp]; + const ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(_cv1)); + const ADReal fv2 = 1.0 - chi / (1. + chi * fv1); + const ADReal S_tilde = + vorticity_mag + _mu_tilde[_qp] * fv2 / (_kappa * _kappa * d * d * _rho[_qp]); + const ADReal S = S_tilde + 2 * std::min(0.0, _strain_mag[_qp] - vorticity_mag); + ADReal r; + if (S_tilde <= 0.0) // Avoid potential division by zero + r = 10.; + else + r = std::min(_mu_tilde[_qp] / (S_tilde * _kappa * _kappa * d * d * _rho[_qp]), 10.0); + const ADReal g = r + _cw2 * (Utility::pow<6>(r) - r); + const ADReal fw = g * std::pow((1. + Utility::pow<6>(_cw3)) / + (Utility::pow<6>(g) + Utility::pow<6>(_cw3)), + 1.0 / 6.0); + + // Compute strong forms of the SA equation + if (_use_ft2_term) // Whether to apply the f_t2 term in the SA equation + { + const ADReal ft2 = _ct3 * std::exp(-_ct4 * chi * chi); + _destruction_strong_residual[_qp] = + (_cw1 * fw - _cb1 * ft2 / _kappa / _kappa) * Utility::pow<2>(_mu_tilde[_qp] / d); + _source_strong_residual[_qp] = -(1 - ft2) * _rho[_qp] * _cb1 * S * _mu_tilde[_qp]; + } + else + { + _destruction_strong_residual[_qp] = _cw1 * fw * Utility::pow<2>(_mu_tilde[_qp] / d); + _source_strong_residual[_qp] = -_rho[_qp] * _cb1 * S * _mu_tilde[_qp]; + } + _convection_strong_residual[_qp] = _rho[_qp] * _velocity[_qp] * _grad_mu[_qp]; + _diffusion_strong_residual[_qp] = -1.0 / _sigma * _cb2 * (_grad_mu[_qp] * _grad_mu[_qp]); + if (_has_transient) + _time_strong_residual[_qp] = (*_visc_dot)[_qp] * _rho[_qp]; + _visc_strong_residual[_qp] = _has_transient ? _time_strong_residual[_qp] : 0.0; + _visc_strong_residual[_qp] += (_convection_strong_residual[_qp] + + _destruction_strong_residual[_qp] + + _diffusion_strong_residual[_qp] + + _source_strong_residual[_qp]); + + // Compute the tau stabilization parameter for mu_tilde SUPG stabilization + const auto nu_visc = (_mu[_qp] + _mu_tilde[_qp]) / _rho[_qp] / _sigma; + const auto transient_part = _has_transient ? 4.0 / (_dt * _dt) : 0.0; + const auto speed = NS::computeSpeed(_velocity[_qp]); + _tau_visc[_qp] = _alpha / std::sqrt(transient_part + + (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + + 9.0 * (4.0 * nu_visc / (_hmax * _hmax)) * + 4.0 * (nu_visc / (_hmax * _hmax))); + + // Replace the nu value in the tau stabilization parameter for INS SUPG stabilization + const auto nu = (_mu[_qp] + _mu_tilde[_qp] * fv1) / _rho[_qp]; + _tau[_qp] = _alpha / std::sqrt(transient_part + + (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + + 9.0 * (4.0 * nu / (_hmax * _hmax)) * + 4.0 * (nu / (_hmax * _hmax))); +} diff --git a/include/materials/SATauStabilized3Eqn.h b/include/materials/SATauStabilized3Eqn.h index 41a02d5172..8fb53afef6 100644 --- a/include/materials/SATauStabilized3Eqn.h +++ b/include/materials/SATauStabilized3Eqn.h @@ -1,13 +1,10 @@ #pragma once #include "SATauMaterial.h" +#include "INSADStabilized3Eqn.h" -class INSADMaterial; class INSADStabilized3Eqn; -#include "SATauMaterial.h" -#include "INSADStabilized3Eqn.h" - class SATauStabilized3Eqn : public SATauMaterialTempl { public: diff --git a/moose b/moose index b38351138b..cedfd055b0 160000 --- a/moose +++ b/moose @@ -1 +1 @@ -Subproject commit b38351138bf10bd150c21926f36a6a286e278a66 +Subproject commit cedfd055b01b8ccebc50e5d10a89847caa8e24ae diff --git a/src/actions/NtAction.C b/src/actions/NtAction.C index d9fd95f9f1..3ff027cca2 100644 --- a/src/actions/NtAction.C +++ b/src/actions/NtAction.C @@ -27,7 +27,8 @@ NtAction::validParams() params.addRequiredParam("num_precursor_groups", "specifies the total number of precursors to create"); - params.addRequiredParam("var_name_base", "specifies the base name of the variables"); + params.addRequiredParam("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("nt_ic_function", "An initial condition function for the neutrons."); - params.addParam>("vacuum_boundaries", - "The boundaries on which to apply vacuum boundaries."); + params.addParam>( + "vacuum_boundaries", + "The boundaries on which to apply vacuum boundaries."); + MooseEnum vacuum_bc_type("marshak mark milne", "marshak"); + params.addParam("vacuum_bc_type", vacuum_bc_type, + "Whether to apply Marshak, Mark, or Milne vacuum boundary conditions. Defaults to Marshak."); params.addParam( "create_temperature_var", true, "Whether to create the temperature variable."); params.addParam( @@ -147,6 +152,7 @@ NtAction::act() params.set("variable") = var_name; if (isParamValid("use_exp_form")) params.set("use_exp_form") = getParam("use_exp_form"); + params.set("vacuum_bc_type") = getParam("vacuum_bc_type"); std::string bc_name = "VacuumConcBC_" + var_name; _problem->addBoundaryCondition("VacuumConcBC", bc_name, params); } diff --git a/src/auxkernels/WallShearStressAux.C b/src/auxkernels/WallShearStressAux.C index 046948785b..43ba0c13b6 100644 --- a/src/auxkernels/WallShearStressAux.C +++ b/src/auxkernels/WallShearStressAux.C @@ -40,5 +40,5 @@ WallShearStressAux::computeValue() // Parallel component of prior velocity gradient grad_vel_wall -= (_normals[_qp] * grad_vel_wall) * _normals[_qp]; - return raw_value(_mu[_qp]) * std::norm(grad_vel_wall); + return raw_value(_mu[_qp]) * grad_vel_wall.norm(); } diff --git a/src/base/ScalarTransportBase.C b/src/base/ScalarTransportBase.C index 28a7d12801..57d1e418e8 100644 --- a/src/base/ScalarTransportBase.C +++ b/src/base/ScalarTransportBase.C @@ -7,7 +7,7 @@ ScalarTransportBase::validParams() { InputParameters params = emptyInputParameters(); params.addParam("use_exp_form", - true, + false, "Whether concentrations should be in an expotential/logarithmic format."); return params; } diff --git a/src/bcs/VacuumConcBC.C b/src/bcs/VacuumConcBC.C index cfb9a24bf6..71718bc413 100644 --- a/src/bcs/VacuumConcBC.C +++ b/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("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("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; } diff --git a/src/materials/SATauMaterial.C b/src/materials/SATauMaterial.C index 54eebf5972..30bf407569 100644 --- a/src/materials/SATauMaterial.C +++ b/src/materials/SATauMaterial.C @@ -3,144 +3,4 @@ registerMooseObject("MoltresApp", SATauMaterial); -template -InputParameters -SATauMaterialTempl::validParams() -{ - InputParameters params = T::validParams(); - params.addClassDescription( - "This is the material class used to compute the stabilization parameter tau_viscosity " - "for the Spalart-Allmaras turbulent viscosity equation."); - params.addRequiredCoupledVar("mu_tilde", "Spalart-Allmaras turbulence viscosity variable"); - params.addRequiredCoupledVar("wall_distance_var", "Wall distance aux variable name"); - params.addParam("use_ft2_term", false, "Whether to apply the f_t2 term in the equation"); - return params; -} - -template -SATauMaterialTempl::SATauMaterialTempl(const InputParameters & parameters) - : T(parameters), - _sigma(2.0/3.0), - _cb1(0.1355), - _cb2(0.622), - _kappa(0.41), - _cw1(_cb1 / _kappa / _kappa + (1 + _cb2) / _sigma), - _cw2(0.3), - _cw3(2.0), - _cv1(7.1), - _ct1(1.0), - _ct2(2.0), - _ct3(1.2), - _ct4(0.5), - _mu_tilde(this->template adCoupledValue("mu_tilde")), - _grad_mu(this->template adCoupledGradient("mu_tilde")), - _tau_visc(this->template declareADProperty("tau_viscosity")), - _wall_dist(this->template coupledValue("wall_distance_var")), - _use_ft2_term(this->template getParam("use_ft2_term")), - _strain_mag(this->template declareADProperty("strain_mag")), - _convection_strong_residual(this->template - declareADProperty("convection_strong_residual")), - _destruction_strong_residual(this->template - declareADProperty("destruction_strong_residual")), - _diffusion_strong_residual(this->template - declareADProperty("diffusion_strong_residual")), - _source_strong_residual(this->template - declareADProperty("source_strong_residual")), - _time_strong_residual(this->template - declareADProperty("time_strong_residual")), - _visc_strong_residual(this->template - declareADProperty("visc_strong_residual")) -{ -} - -template -void -SATauMaterialTempl::subdomainSetup() -{ - T::subdomainSetup(); - - if (_has_transient) - _visc_dot = & this->template adCoupledDot("mu_tilde"); - else - _visc_dot = nullptr; -} - -template -void -SATauMaterialTempl::computeQpProperties() -{ - T::computeQpProperties(); - - // Compute strain rate and vorticity magnitudes - _strain_mag[_qp] = 2.0 * Utility::pow<2>(_grad_velocity[_qp](0, 0)) + - 2.0 * Utility::pow<2>(_grad_velocity[_qp](1, 1)) + - 2.0 * Utility::pow<2>(_grad_velocity[_qp](2, 2)) + - Utility::pow<2>(_grad_velocity[_qp](0, 2) + _grad_velocity[_qp](2, 0)) + - Utility::pow<2>(_grad_velocity[_qp](0, 1) + _grad_velocity[_qp](1, 0)) + - Utility::pow<2>(_grad_velocity[_qp](1, 2) + _grad_velocity[_qp](2, 1)); - _strain_mag[_qp] = std::sqrt(_strain_mag[_qp] + 1e-16); - ADReal vorticity_mag = - Utility::pow<2>(_grad_velocity[_qp](0, 2) - _grad_velocity[_qp](2, 0)) + - Utility::pow<2>(_grad_velocity[_qp](0, 1) - _grad_velocity[_qp](1, 0)) + - Utility::pow<2>(_grad_velocity[_qp](1, 2) - _grad_velocity[_qp](2, 1)); - vorticity_mag = std::sqrt(vorticity_mag + 1e-16); - - // Compute relevant parameters for the SA equation - const Real d = std::max(_wall_dist[_qp], 1e-16); // Avoid potential division by zero - const ADReal chi = _mu_tilde[_qp] / _mu[_qp]; - const ADReal fv1 = Utility::pow<3>(chi) / (Utility::pow<3>(chi) + Utility::pow<3>(_cv1)); - const ADReal fv2 = 1.0 - chi / (1. + chi * fv1); - const ADReal S_tilde = - vorticity_mag + _mu_tilde[_qp] * fv2 / (_kappa * _kappa * d * d * _rho[_qp]); - const ADReal S = S_tilde + 2 * std::min(0.0, _strain_mag[_qp] - vorticity_mag); - ADReal r; - if (S_tilde <= 0.0) // Avoid potential division by zero - r = 10.; - else - r = std::min(_mu_tilde[_qp] / (S_tilde * _kappa * _kappa * d * d * _rho[_qp]), 10.0); - const ADReal g = r + _cw2 * (Utility::pow<6>(r) - r); - const ADReal fw = g * std::pow((1. + Utility::pow<6>(_cw3)) / - (Utility::pow<6>(g) + Utility::pow<6>(_cw3)), - 1.0 / 6.0); - - // Compute strong forms of the SA equation - if (_use_ft2_term) // Whether to apply the f_t2 term in the SA equation - { - const ADReal ft2 = _ct3 * std::exp(-_ct4 * chi * chi); - _destruction_strong_residual[_qp] = - (_cw1 * fw - _cb1 * ft2 / _kappa / _kappa) * Utility::pow<2>(_mu_tilde[_qp] / d); - _source_strong_residual[_qp] = -(1 - ft2) * _rho[_qp] * _cb1 * S * _mu_tilde[_qp]; - } - else - { - _destruction_strong_residual[_qp] = _cw1 * fw * Utility::pow<2>(_mu_tilde[_qp] / d); - _source_strong_residual[_qp] = -_rho[_qp] * _cb1 * S * _mu_tilde[_qp]; - } - _convection_strong_residual[_qp] = _rho[_qp] * _velocity[_qp] * _grad_mu[_qp]; - _diffusion_strong_residual[_qp] = -1.0 / _sigma * _cb2 * (_grad_mu[_qp] * _grad_mu[_qp]); - if (_has_transient) - _time_strong_residual[_qp] = (*_visc_dot)[_qp] * _rho[_qp]; - _visc_strong_residual[_qp] = _has_transient ? _time_strong_residual[_qp] : 0.0; - _visc_strong_residual[_qp] += (_convection_strong_residual[_qp] + - _destruction_strong_residual[_qp] + - _diffusion_strong_residual[_qp] + - _source_strong_residual[_qp]); - - // Compute the tau stabilization parameter for mu_tilde SUPG stabilization - const auto nu_visc = (_mu[_qp] + _mu_tilde[_qp]) / _rho[_qp] / _sigma; - const auto transient_part = _has_transient ? 4.0 / (_dt * _dt) : 0.0; - const auto speed = NS::computeSpeed(_velocity[_qp]); - _tau_visc[_qp] = _alpha / std::sqrt(transient_part + - (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + - 9.0 * (4.0 * nu_visc / (_hmax * _hmax)) * - 4.0 * (nu_visc / (_hmax * _hmax))); - - // Replace the nu value in the tau stabilization parameter for INS SUPG stabilization - const auto nu = (_mu[_qp] + _mu_tilde[_qp] * fv1) / _rho[_qp]; - _tau[_qp] = _alpha / std::sqrt(transient_part + - (2.0 * speed / _hmax) * (2.0 * speed / _hmax) + - 9.0 * (4.0 * nu / (_hmax * _hmax)) * - 4.0 * (nu / (_hmax * _hmax))); -} - template class SATauMaterialTempl; diff --git a/tests/bcs/gold/vacuum_bc_mark_out.e b/tests/bcs/gold/vacuum_bc_mark_out.e new file mode 100644 index 0000000000..0fe8a2fe3b Binary files /dev/null and b/tests/bcs/gold/vacuum_bc_mark_out.e differ diff --git a/tests/bcs/gold/vacuum_bc_milne_out.e b/tests/bcs/gold/vacuum_bc_milne_out.e new file mode 100644 index 0000000000..83e3ec48c1 Binary files /dev/null and b/tests/bcs/gold/vacuum_bc_milne_out.e differ diff --git a/tests/bcs/gold/vacuum_bc_out.e b/tests/bcs/gold/vacuum_bc_out.e new file mode 100644 index 0000000000..285895b3b8 Binary files /dev/null and b/tests/bcs/gold/vacuum_bc_out.e differ diff --git a/tests/bcs/tests b/tests/bcs/tests index f44f2d9ed4..e8b547e0e4 100644 --- a/tests/bcs/tests +++ b/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.' + [] [] diff --git a/tests/bcs/vacuum_bc.i b/tests/bcs/vacuum_bc.i new file mode 100644 index 0000000000..0d6f00de77 --- /dev/null +++ b/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 + [] +[] diff --git a/tests/bcs/xsdata.json b/tests/bcs/xsdata.json new file mode 100644 index 0000000000..f8705e06f3 --- /dev/null +++ b/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 + ] + } +} diff --git a/tests/pre/gold/pre_loop_ins_out.e b/tests/pre/gold/pre_loop_ins_out.e index b616343e92..2f2c478abb 100644 Binary files a/tests/pre/gold/pre_loop_ins_out.e and b/tests/pre/gold/pre_loop_ins_out.e differ diff --git a/tests/pre/gold/pre_loop_out.e b/tests/pre/gold/pre_loop_out.e index df40f9e886..6827324952 100644 Binary files a/tests/pre/gold/pre_loop_out.e and b/tests/pre/gold/pre_loop_out.e differ