Merge pull request #54 from gridley/rodWorthMatcher

Rod worth matcher

problems/rodMatchSerpentWorth/2d_rodded_lattice.geo

@@ -0,0 +1,84 @@
+Geometry.CopyMeshingMethod = 1;
+R = 78.58;
+H = 2*R + 6.75;
+num_segments = 14;
+pitch = R / num_segments;
+x = .237952211;
+fuel_rad = x * pitch;
+graph_rad = pitch;
+rodRad = 2.0; // CR radius
+lc = 1;
+lx = 0.4;
+ly = 2.5;
+
+Point(1) = {0, 0, 0, lc};
+Point(2) = {fuel_rad, 0, 0, lc};
+Point(3) = {graph_rad, 0, 0, lc};
+Point(4) = {0, H, 0, lc};
+Point(5) = {fuel_rad, H, 0, lc};
+Point(6) = {graph_rad, H, 0, lc};
+// fuel top
+Line(1) = {4, 5};
+// graph top
+Line(2) = {5, 6};
+// graph right edge
+Line(3) = {6, 3};
+// graph bottom
+Line(4) = {3, 2};
+// fuel bottom
+Line(5) = {1, 2};
+// fuel-graph interface
+Line(6) = {2, 5};
+// fuel left edge
+Line(7) = {4, 1};
+
+// Fuel
+Line Loop(8) = {1, -6, -5, -7};
+Plane Surface(9) = {8};
+
+// Moderator
+Line Loop(10) = {2, 3, 4, 6};
+Plane Surface(11) = {10};
+
+// Structured
+Transfinite Line{1, 5} = fuel_rad/lx;
+Transfinite Line{2, 4} = (graph_rad - fuel_rad)/lx;
+Transfinite Line{3, 6, 7} = H/ly;
+Transfinite Surface{9};
+Transfinite Surface{11};
+Recombine Surface{9};
+Recombine Surface{11};
+
+fuel_surfaces[] = {9};
+moder_surfaces[] = {11};
+fuel_tops[] = {1};
+fuel_bottoms[] = {5};
+moder_bottoms[] = {4};
+moder_tops[] = {2};
+For xindex In {1:num_segments-1}
+    new_f_surface = Translate {xindex*pitch, 0, 0} {
+      Duplicata { Surface{9}; }
+    };
+    fuel_surfaces += new_f_surface;
+    new_m_surface = Translate {xindex*pitch, 0, 0} {
+      Duplicata { Surface{11}; }
+    };
+    moder_surfaces += new_m_surface;
+    fuel_tops += {13 + (xindex - 1) * 8};
+    moder_tops += {17 + (xindex - 1) * 8};
+    fuel_bottoms += {15 + (xindex - 1) * 8};
+    moder_bottoms += {19 + (xindex - 1) * 8};
+    EndFor // xindex
+
+
+Physical Surface("cRod") =       { fuel_surfaces[0]            };
+Physical Line("cRod_top") =      { fuel_tops[0]                };
+Physical Line("cRod_bot") =      { fuel_bottoms[0]                        };
+// END NEW STUFF
+Physical Surface("fuel") =       { fuel_surfaces[{1:num_segments}]             };
+Physical Surface("moder") =      { moder_surfaces[]            };
+Physical Line("fuel_tops") =     { fuel_tops[{1:num_segments}]                 };
+Physical Line("moder_tops") =    { moder_tops[]                };
+Physical Line("fuel_bottoms") =  { fuel_bottoms[{1:num_segments}]              };
+Physical Line("moder_bottoms") = { moder_bottoms[]             };
+Physical Line("outer_wall") =    { 18 + 8 * (num_segments - 2) };

problems/rodMatchSerpentWorth/README.md

@@ -0,0 +1,8 @@
+The solve here is only done in eigenvalue mode. The objective is to change the absorbtion of the control rod until
+the same reactivity worth experienced at the MSRE is obtained in Moltres.
+
+In this top directory, the radius of the 2D case found the match Serpent's unrodded reactivity was found to be 78.58cm.
+Serpent found k=1.045 +- 0.0010. This simulation yields k= 1.045177e+00.
+
+In the adjustAbsorb directory, the control rod's absorbtion factor was found to match the MSRE results. To reproduce results,
+the .msh file generated by gmsh must be copied into that directory.

problems/rodMatchSerpentWorth/adjustAbsorb/README.md

@@ -0,0 +1,6 @@
+Serpent found that with one rod inserted, k fell to k=1.017+-0.001. Rod insertion means that the rod has a position
+23.622 cm above the bottom of the core. We know this from MSRE Design and Operations Report part 1 where the control rod
+system was described.
+
+The multiplicative factor to achieve the ORNL-specified rod worth of 2.8% \Delta k/k was found to be 14.22. This was found with
+a few secant method iterations on paper while I was reading stuff on sciencedirect.

problems/rodMatchSerpentWorth/adjustAbsorb/rodded.i

@@ -0,0 +1,231 @@
+flow_velocity=21.7 # cm/s. See MSRE-properties.ods
+nt_scale=1e13
+ini_temp=922
+diri_temp=922
+
+[GlobalParams]
+  num_groups = 4
+  num_precursor_groups = 6
+  use_exp_form = false
+  group_fluxes = 'group1 group2 group3 group4'
+  temperature = temp
+  sss2_input = true
+  pre_concs = 'pre1 pre2 pre3 pre4 pre5 pre6'
+  account_delayed = false
+[]
+
+[Mesh]
+  file = '2d_rodded_lattice.msh'
+[../]
+
+[Problem]
+  coord_type = RZ
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [./temp]
+    initial_condition = ${ini_temp}
+    scaling = 1e-4
+  [../]
+[]
+
+[Precursors]
+  [./pres]
+    var_name_base = pre
+    block = 'fuel'
+    outlet_boundaries = 'fuel_tops'
+    u_def = 0
+    v_def = ${flow_velocity}
+    w_def = 0
+    nt_exp_form = false
+    family = MONOMIAL
+    order = CONSTANT
+    # jac_test = true
+  [../]
+[]
+
+[Nt]
+  var_name_base = group
+  vacuum_boundaries = 'fuel_bottoms fuel_tops moder_bottoms moder_tops outer_wall cRod_top cRod_bot'
+  create_temperature_var = false
+  eigen = true
+[]
+
+[Kernels]
+  [./temp_source_fuel]
+    type = FissionHeatSource
+    variable = temp
+    nt_scale=${nt_scale}
+    block = 'fuel'
+    power = 7.5e6
+    tot_fissions = tot_fissions
+  [../]
+  [./temp_source_mod]
+    type = GammaHeatSource
+    variable = temp
+    gamma = .0144 # Cammi .0144
+    block = 'moder'
+    average_fission_heat = 'tot_fissions'
+  [../]
+  [./temp_diffusion]
+    type = MatDiffusion
+    D_name = 'k'
+    variable = temp
+  [../]
+  [./temp_advection_fuel]
+    type = ConservativeTemperatureAdvection
+    velocity = '0 ${flow_velocity} 0'
+    variable = temp
+    block = 'fuel'
+  [../]
+[]
+
+[BCs]
+  [./temp_diri_cg]
+    boundary = 'moder_bottoms fuel_bottoms outer_wall'
+    type = DirichletBC
+    variable = temp
+    value = ${diri_temp}
+  [../]
+  [./temp_advection_outlet]
+    boundary = 'fuel_tops'
+    type = TemperatureOutflowBC
+    variable = temp
+    velocity = '0 ${flow_velocity} 0'
+  [../]
+[]
+
+[Materials]
+  [./fuel]
+    type = GenericMoltresMaterial
+    property_tables_root = '../../../tutorial/step01_groupConstants/MSREProperties/msre_gentry_4gfuel_'
+    interp_type = 'spline'
+    block = 'fuel'
+    prop_names = 'k cp'
+    prop_values = '.0553 1967' # Robertson MSRE technical report @ 922 K
+  [../]
+  [./rho_fuel]
+    type = DerivativeParsedMaterial
+    f_name = rho
+    function = '2.146e-3 * exp(-1.8 * 1.18e-4 * (temp - 922))'
+    args = 'temp'
+    derivative_order = 1
+    block = 'fuel'
+  [../]
+  [./moder]
+    type = GenericMoltresMaterial
+    property_tables_root = '../../../tutorial/step01_groupConstants/MSREProperties/msre_gentry_4gmoder_'
+    interp_type = 'spline'
+    prop_names = 'k cp'
+    prop_values = '.312 1760' # Cammi 2011 at 908 K
+    block = 'moder'
+  [../]
+  [./rho_moder]
+    type = DerivativeParsedMaterial
+    f_name = rho
+    function = '1.86e-3 * exp(-1.8 * 1.0e-5 * (temp - 922))'
+    args = 'temp'
+    derivative_order = 1
+    block = 'moder'
+  [../]
+  [./cRod]
+    type = RoddedMaterial
+    property_tables_root = '../../../tutorial/step01_groupConstants/MSREProperties/msre_gentry_4gmoder_'
+    interp_type = 'spline'
+    prop_names = 'k cp'
+    prop_values = '.312 1760' # Cammi 2011 at 908 K
+    block = 'cRod'
+    rodDimension = 'y'
+    rodPosition = 23.622
+    absorb_factor = 14.22 # how much more absorbing than usual in absorbing region?
+  [../]
+  [./rho_crod]
+    type = DerivativeParsedMaterial
+    f_name = rho
+    function = '1.86e-3 * exp(-1.8 * 1.0e-5 * (temp - 922))'
+    args = 'temp'
+    derivative_order = 1
+    block = 'cRod'
+  [../]
+[]
+
+[Executioner]
+  type = InversePowerMethod
+  max_power_iterations = 50
+  xdiff = 'group1diff'
+
+  bx_norm = 'bnorm'
+  k0 = 1.5
+  pfactor = 1e-2
+  l_max_its = 100
+
+  # solve_type = 'PJFNK'
+  solve_type = 'NEWTON'
+  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
+  petsc_options_iname = '-pc_type -sub_pc_type -pc_factor_shift_type -pc_factor_shift_amount'
+  petsc_options_value = 'asm lu NONZERO 1e-10'
+[]
+
+[Preconditioning]
+  [./SMP]
+    type = SMP
+    full = true
+  [../]
+[]
+
+[Postprocessors]
+  [./bnorm]
+    type = ElmIntegTotFissNtsPostprocessor
+    execute_on = linear
+  [../]
+  [./tot_fissions]
+    type = ElmIntegTotFissPostprocessor
+    execute_on = linear
+  [../]
+  [./group1norm]
+    type = ElementIntegralVariablePostprocessor
+    variable = group1
+    execute_on = linear
+  [../]
+  [./group1max]
+    type = NodalMaxValue
+    variable = group1
+    execute_on = timestep_end
+  [../]
+  [./group1diff]
+    type = ElementL2Diff
+    variable = group1
+    execute_on = 'linear timestep_end'
+    use_displaced_mesh = false
+  [../]
+  [./group2norm]
+    type = ElementIntegralVariablePostprocessor
+    variable = group2
+    execute_on = linear
+  [../]
+  [./group2max]
+    type = NodalMaxValue
+    variable = group2
+    execute_on = timestep_end
+  [../]
+  [./group2diff]
+    type = ElementL2Diff
+    variable = group2
+    execute_on = 'linear timestep_end'
+    use_displaced_mesh = false
+  [../]
+[]
+
+[Outputs]
+  print_perf_log = true
+  print_linear_residuals = true
+  csv = true
+  exodus = true
+  execute_on = 'final'
+  file_base = 'out'
+[]
+
+[Debug]
+  show_var_residual_norms = true
+[]