diff --git a/problems/rodMatchSerpentWorth/2d_rodded_lattice.geo b/problems/rodMatchSerpentWorth/2d_rodded_lattice.geo new file mode 100644 index 0000000000..667eadf97b --- /dev/null +++ b/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) }; diff --git a/problems/rodMatchSerpentWorth/README.md b/problems/rodMatchSerpentWorth/README.md new file mode 100644 index 0000000000..7d509fbcb9 --- /dev/null +++ b/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. diff --git a/problems/rodMatchSerpentWorth/adjustAbsorb/README.md b/problems/rodMatchSerpentWorth/adjustAbsorb/README.md new file mode 100644 index 0000000000..7262dfe61b --- /dev/null +++ b/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. diff --git a/problems/rodMatchSerpentWorth/adjustAbsorb/rodded.i b/problems/rodMatchSerpentWorth/adjustAbsorb/rodded.i new file mode 100644 index 0000000000..0661e519c7 --- /dev/null +++ b/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 +[] diff --git a/problems/rodMatchSerpentWorth/rodded.i b/problems/rodMatchSerpentWorth/rodded.i new file mode 100644 index 0000000000..24eaaa39d3 --- /dev/null +++ b/problems/rodMatchSerpentWorth/rodded.i @@ -0,0 +1,232 @@ +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 = 200.0 + absorb_factor = 1e6 # 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 +[] diff --git a/python/controlRodWorthMatcher.py b/python/controlRodWorthMatcher.py new file mode 100644 index 0000000000..d86caaa8bc --- /dev/null +++ b/python/controlRodWorthMatcher.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python3 +# This script matches control rod worths of the Moltres RoddedMaterial to +# results from a more detailed transport code. +# +# usage: +# controlRodWorthMatcher.py +# +# The file called rodWorthFile can have any number of rows. The last column +# is assumed to be the control rod worth IN RELATIVE REACTIVITY ie change in (k-1)/k +# as compared to the state with no control rods inserted. +# All preceding columns are assumed to belong to control rods 1, 2, 3, all +# the way up to n. This script will find the moltres absorbtion factor for +# each rod that most closely matches reactivities. +# +# NOTE: +# This script is unfinished. I'm putting it into the python directory so +# that others have an idea of the methodology I had in mind to determine +# the rod multiplicative factor for absorbtion for systems with many rods. +# +# The general idea would be to form a residual term that is squared reactivity +# differences between detailed transport (e.g. Serpent) and Moltres' results. +# A Jacobian is formed from perturbing each rod multiplicative factor. Then +# a Newton iteration is done to drive the residual to zero. Repeat until results +# are within uncertainty of the Monte Carlo neutronics. Note that forming the Jacobian +# for this residual should probably happen in parallel on a computing cluster. +# +import numpy +import subprocess +import sys + +inputFile = sys.argv[1] +rodWorthFile = sys.argv[2] + +if len(sys.argv) > 2: + raise Exception("Too many command line args.") + +with open(inputFile) as infile: + infileText = infile.read() + +with open(rodWorthFile) as rWFile: + for line in rWFile: + sline = line.split() + rodPos, worth = ( float(sline[0]), float(sline[1]) ) diff --git a/tutorial/step01_groupConstants/extractSerpent2GCs.py b/python/extractSerpent2GCs.py similarity index 100% rename from tutorial/step01_groupConstants/extractSerpent2GCs.py rename to python/extractSerpent2GCs.py diff --git a/tutorial/step01_groupConstants/README.md b/tutorial/step01_groupConstants/README.md index 6e632e5c45..e0caa70cc8 100644 --- a/tutorial/step01_groupConstants/README.md +++ b/tutorial/step01_groupConstants/README.md @@ -29,7 +29,10 @@ get parsed by PyNE in the extractSerpent2GCs script. The command to run in order to generate the moltres-compatible group constants is: -```./extractSerpent2GCs.py MSREProperties msre_gentry_4g tempMapping.txt universeMapping.txt fuel moder``` +```$MOLTRES/python/extractSerpent2GCs.py MSREProperties msre_gentry_4g tempMapping.txt universeMapping.txt fuel moder``` + +Where $MOLTRES is an environment variable leading to the install location of Moltres. An alternative would be to just add +the $MOLTRES/python directory to your path. The syntax requires the arbitrary directory name you'd like to create, then the arbitrary file base name that moltres will look at, then a file that maps branch names to temperatures, then a file that maps universe numbers from serpent to material names,