Steady State Results do not seem right.

[Jeremyb8707bigbrain]

Hey everyone I figured out how to run a Steady and the results arent exactly what I am looking for nor expected.
I believe this is me absolutely missing something that is crucial with scaling but who knows.

My results showed numbers around 1.1e-06 for all variables, with some spatial variance but not enough that made any sense to me. I know I have VacuumConcBCs, I would like to measure leakage effects with differing graphite thicknesses but it doesnt seem either that my mesh supports criticality or there is some scaling thing I am not familar with. Input file and picture of mesh below.

In OpenMC this same mesh with all reflective boundaries got a k_inf of around 1.65, so I was wondering if there was a better boundary option perhaps or something I could do to preserve leakage effects but not kill my neutrons instantly?

initial_temp = 973.0

[GlobalParams]
    num_groups = 4
    num_precursor_groups = 6
    group_fluxes = 'fast semi-fast epithermal thermal'
    account_delayed = false
    temperature = ${initial_temp}
[]

[Variables]
    [fast]
        order = FIRST
        family = LAGRANGE
        initial_condition = 1
    []
    [semi-fast]
        order = FIRST
        family = LAGRANGE
        initial_condition = 1
    []
    [epithermal]
        order = FIRST
        family = LAGRANGE
        initial_condition = 1
    []
    [thermal]
        order = FIRST
        family = LAGRANGE
        initial_condition = 1
    []
[]

[Mesh]
    file = '/workspace/moltensaltsim/Meshes/8cm-60cm.e'
[]

[Materials]
    [fuel_salt]
        type = MoltresJsonMaterial
        base_file = '/workspace/moltensaltsim/Material_Data/Moltres_JSONs/cluster_8cm.json'
        material_key = 'Molten_Salt'
        block = 'fuel_salt'
        interp_type = 'spline'
    []
    [graphite]
        type = MoltresJsonMaterial
        base_file = '/workspace/moltensaltsim/Material_Data/Moltres_JSONs/cluster_8cm.json'
        material_key = 'Graphite'
        block = 'graphite'
        interp_type = 'spline'
    []
[]

[BCs]
    [fast_vacuum]
        type = VacuumConcBC
        boundary = '10000 15008 15007'
        variable = fast
    []
    [semi-fast_vacuum]
        type = VacuumConcBC
        boundary = '10000 15008 15007'
        variable = semi-fast
    []
    [epithermal_vacuum]
        type = VacuumConcBC
        boundary = '10000 15008 15007'
        variable = epithermal
    []
    [thermal_vacuum]
        type = VacuumConcBC
        boundary = '10000 15008 15007'
        variable = thermal
    []
[]

[Kernels]
    [sigmar_fast]
        type = SigmaR
        variable = fast
        group_number = 1
    []
    [diffusion_fast]
        type = GroupDiffusion
        variable = fast
        group_number = 1
    []
    [scatter_fast]
        type = InScatter
        variable = fast
        group_number = 1
    []
    [fission_fast]
        type = CoupledFissionKernel
        variable = fast
        group_number = 1
        block = 'fuel_salt'
    []
    [sigmar_semi-fast]
        type = SigmaR
        variable = semi-fast
        group_number = 2
    []
    [diffusion_semi-fast]
        type = GroupDiffusion
        variable = semi-fast
        group_number = 2
    []
    [scatter_semi-fast]
        type = InScatter
        variable = semi-fast
        group_number = 2
    []
    [fission_semi-fast]
        type = CoupledFissionKernel
        variable = semi-fast
        group_number = 2
        block = 'fuel_salt'
    []
    [sigmar_epithermal]
        type = SigmaR
        variable = epithermal
        group_number = 3
    []
    [diffusion_epithermal]
        type = GroupDiffusion
        variable = epithermal
        group_number = 3
    []
    [scatter_epithermal]
        type = InScatter
        variable = epithermal
        group_number = 3
    []
    [fission_epithermal]
        type = CoupledFissionKernel
        variable = epithermal
        group_number = 3
        block = 'fuel_salt'
    []
    [sigmar_thermal]
        type = SigmaR
        variable = thermal
        group_number = 4
    []
    [diffusion_thermal]
        type = GroupDiffusion
        variable = thermal
        group_number = 4
    []
    [scatter_thermal]
        type = InScatter
        variable = thermal
        group_number = 4
    []
    [fission_thermal]
        type = CoupledFissionKernel
        variable = thermal
        group_number = 4
        block = 'fuel_salt'
    []
[]

[Executioner]
    type = Steady
    solve_type = NEWTON
[]

[Outputs]
    exodus = true
[]

[Jeremyb8707bigbrain]

I know real MSRs use reflector rings and I was just wondering if I had to model those also just for this slice of the core or would it be better to increase core size, or play with the boundaries?

[Jeremyb8707bigbrain]

@smpark7

The Nt and Eigenvalue mode works great now. My only issue is the k_eff is unphysical (keff of 60). Is it possible that using the NtAction it reverses the order of which the groups are assumed? For example:

initial_temp = 973.0

[GlobalParams]
    num_groups = 4
    num_precursor_groups = 6
    account_delayed = false
    temperature = ${initial_temp}
    sss2_input = false
    use_exp_form = false
[]

[Mesh]
    file = '/workspace/moltensaltsim/Meshes/8cm-60cm.e'
[]

[Materials]
    [fuel_salt]
        type = MoltresJsonMaterial
        base_file = '/workspace/moltensaltsim/Material_Data/Moltres_JSONs/cluster_8cm.json'
        material_key = 'Molten_Salt'
        block = 'fuel_salt'
        interp_type = 'spline'
    []
    [graphite]
        type = MoltresJsonMaterial
        base_file = '/workspace/moltensaltsim/Material_Data/Moltres_JSONs/cluster_8cm.json'
        material_key = 'Graphite'
        block = 'graphite'
        interp_type = 'spline'
    []
[]

[Nt]
    var_name_base = g_
    vacuum_boundaries = '10000 15008 15007'
    pre_blocks = '0'
    create_temperature_var = false
    eigen = true
    fission_blocks = 'fuel_salt'
[]

[Executioner]
    type = Eigenvalue
    solve_type = NEWTON
[]

[Outputs]
    exodus = true
[]

I mean for g_1 to be the fast group, but is it possible that it thinks its the thermal group? I only ask this because my results look like this
group1:

group2:

group3 and group4 are basically 0.
If so that may break the Group Constant ordering and thermal neutrons are using fast neutron diffusion coefficients. Any help here would be appreciated.

[Jeremyb8707bigbrain]

Also realized that my REMXS values for my graphite are really small (e-03 to e-05), I did use pure C12 but I was curious if I missed anything else.

I know we calculate REMXS through the total cross sections and the nu-scatter matrix, would these values be different if we did absorption and nu-scatter, and we added absorption and the out scatter values together?

[smpark7]

If you used OpenMC, sss2_input should be true. This would have transposed your group-to-group scattering matrix

[smpark7]

Your graphite REMXS seems typical. The neutron group numbering 1-4 is consistent with your original fast-thermal ordering. Your scattering matrix was transposed which is why your neutrons were not down-scattering after being born in group 1 and 2.

I recommend using interp_type=linear or monotone_cubic. Spline regularly exhibits unphysical oscillations.

For your executioner settings, I recommend these from $MOLTRES/tutorial/eigenvalue/nts-action.i:

[Executioner]
  type = Eigenvalue
  initial_eigenvalue = 1
  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'

  automatic_scaling = true
  compute_scaling_once = false
  resid_vs_jac_scaling_param = 0.1

  line_search = none
[]

If it still gives weird eigenvalues, you can try including:

[Problem]
  type = EigenProblem
  bx_norm = bnorm
[]

[Postprocessors]
  [bnorm]
    type = ElmIntegTotFissNtsPostprocessor
    block = 'fuel_salt'
    execute_on = linear
  []
[]

[Jeremyb8707bigbrain]

Ok this is working a lot better! I have figured out my vacuum boundaries are probably way too aggressive because its giving an eigenvalue of around 0.08. To be fair I do have all outside boundaries being vacuumed.

If I added another hexagon ring of just pure graphite im wondering if that will reflect some of the neutrons, and then let others leak. Will have to experiment! Thank you very much!

[Jeremyb8707bigbrain]

I added extra graphite around the core to act as a reflector andddddddd the eigenvalue went negative.
Iteration 1 eigenvalue = -0.0040497075

This only occured when I added the extra reflectors, which in theory should increase k_eff.... however it went negative.
g1

g2

g3

g4

My eigenvalue stuff is exactly like you said, which worked great when there was no reflector, even though is was 0.08.
Im curious the cause this time, my reflective OpenMC calculations gave it around 1.65 so im wondering the difference here.

[smpark7]

Try removing all vacuum BCs and see if you get close to OpenMC. All BCs are reflecting by default

[smpark7]

You also didn't set the temperature in the material blocks.

Edit: Ignore this comment. I didn't see the temperature parameter in GlobalParams.