Picard iterations set up correctly?

[megan-crocker]

Hi all,

I am essentially trying to set up a tightly coupled neutronics-thermal simulation using Picard iterations, to iterate back and forth between OpenMC and MOOSE. Currently I have the scripts set up as below. I want the workflow to essentially perform an OpenMC solve, then send the heat source to MOOSE, send temp back to OpenMC and solve to convergence. I want to solve to steady state so I haven't included the time derivative (could this be the problem? the tutorials tend to use this).

I am getting an annoying issue where the solve is 'converging' within the first timestep. I have tried reducing the convergence tolerance and find that I just get oscillations and no actual convergence. This is my first time using this capability and I am just wondering if I have set everything up correctly? I have removed the large functions and blocks to make it easier to read.

Thank you :)

MOOSE.i

  [Mesh]
    [file]
      type = FileMeshGenerator
      file = mesh_tagged_m.msh
    []
  
  [Variables]
    [temp] 
      initial_condition = 673.15
    []
  []
  
  [Problem]
    kernel_coverage_check = false
  []
  
  [Functions]
    [ss316_k]
      type = PiecewiseLinear
      y = '14.25...'
      x = '290...'
    []
  
    [ss316_cp]
      type = PiecewiseLinear
      y =  '478.6...'
      x = '290...'
    []
  
    [ss316_cte]
      type = PiecewiseLinear
      y =  '1.5125505779752358e-05...'
      x = '290...'
    []
  
    [tungsten_k]
      type = PiecewiseLinear
      y =  '173.178...'
      x = '293.15...'
    []
  
    [tungsten_cp]
      type = PiecewiseLinear
      y =  '129.09...'
      x = '293.15...
    []
  
    [tungsten_cte]
      type = PiecewiseLinear
      y =  '4.484830611478225e-06...'
      x = '293.15...'
    []
  
    [li4_si04_k]
      type = PiecewiseLinear
      y =  '0.776...'
      x = '290...'
    []
  
    [li4_si04_cp]
      type = PiecewiseLinear
      y =  '916.99...'
      x = '290...'
    []
  
    [be12_ti_cp]
      type = PiecewiseLinear
      y =  '1433.845...'
      x = '323.15...'
    []
  []
  
  [AuxVariables]
    [heat_source]
      family = MONOMIAL
      order = CONSTANT
    []
  
    [k_output]
      family = MONOMIAL
      order = CONSTANT
    []
  
    [cp_output]
      family = MONOMIAL
      order = CONSTANT
    []
  []
  
  [AuxKernels]
    [k_output_aux]
      type = MaterialRealAux
      variable = k_output
      property = thermal_conductivity
      execute_on = timestep_end
    []
  
    [cp_output_aux]
      type = MaterialRealAux
      variable = cp_output
      property = specific_heat
      execute_on = timestep_end
    []
  []
  
  [Kernels]
    [hc]
      type = HeatConduction
      variable = temp
    []
  
    [heat]
      type = CoupledForce
      variable = temp
      v = heat_source
    []
  []
  
  [BCs]
    [first_wall_heat_flux]
      type = NeumannBC
      variable = temp
      boundary = 'first_wall_flux_surface'
      value = 0.5e6
    []
  
    [coolant_temp]
      type = DirichletBC
      boundary = 'end_surfaces fw_structure_coolant breeder_coolants structure_breeder_coolant_interface'
      value = 673.15
      variable = temp
    []
  
    [convective_bc_coolants]
      type = ConvectiveHeatFluxBC
      variable = temp
      boundary = 'end_surfaces fw_structure_coolant breeder_coolants structure_breeder_coolant_interface'
      T_infinity = 673.15
      heat_transfer_coefficient = 6000
    []
  []
  
  [Materials]
    [ss316_heat]
      type = HeatConductionMaterial
      temp = temp
      thermal_conductivity_temperature_function = ss316_k
      specific_heat_temperature_function = ss316_cp
      block = 'structure_front structure_rear'
    []
  
    [ss316_density]
      type = GenericConstantMaterial
      prop_names = 'density'
      prop_values = 7750
      block = 'structure_front structure_rear'
    []
  
    [tungsten_heat]
      type = HeatConductionMaterial
      temp = temp
      thermal_conductivity_temperature_function = tungsten_k
      specific_heat_temperature_function = tungsten_cp
      block = 'first_wall'
    []
  
    [tungsten_density]
      type = GenericConstantMaterial
      prop_names = 'density'
      prop_values = 19200
      block = 'first_wall'
    []
    
    [li4si04_heat]
      type = HeatConductionMaterial
      temp = temp
      thermal_conductivity_temperature_function = li4_si04_k
      specific_heat_temperature_function = li4_si04_cp
      block = 'Breeder_00_0 Breeder_01_0...'
    []
  
    [li4si04_density]
      type = GenericConstantMaterial
      prop_names = 'density'
      prop_values = 1584
      block = 'Breeder_00_0 Breeder_01_0...'
    []
  
    [be12ti_heat]
      type = HeatConductionMaterial
      temp = temp
      thermal_conductivity = 45.0
      specific_heat_temperature_function = be12_ti_cp
      block = 'Multiplier_0 Multiplier_1 Multiplier_2 Multiplier_3 Multiplier_4 Multiplier_5'
    []
  
    [be12ti_density]
      type = GenericConstantMaterial
      prop_names = 'density'
      prop_values = 2250
      block = 'Multiplier_0 Multiplier_1 Multiplier_2 Multiplier_3 Multiplier_4 Multiplier_5'
    []
  
    [he_heat]
      type = HeatConductionMaterial
      temp = temp
      thermal_conductivity = 250
      specific_heat = 20
      block = 'coolant_01...'
    []
  
    [he_density]
      type = GenericConstantMaterial
      prop_names = 'density'
      prop_values = 12.0
      block = 'coolant_01...'
    []
  []
  
  [Preconditioning]
    [SMP]
      type = SMP
      full = true
    []
  []
  
  [Executioner]
    type = Transient
    solve_type = 'NEWTON'
  
    end_time = 1
  
    multiapp_fixed_point_convergence = fp_conv
  
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = 'hypre boomeramg'
  
    verbose = true
  
    l_max_its = 1000
    nl_max_its = 200
  
    nl_abs_tol = 1e-10
    fixed_point_rel_tol = 1e-6
    fixed_point_abs_tol = 1e-10
  []
  
  [Outputs]
    exodus = true
    csv = true
    file_base = 'cosim_moose_thermal'
  []
  
  [MultiApps]
    [openmc]
      type = TransientMultiApp
      input_files = '7_openmc_cosim.i'
      execute_on = timestep_begin
    []
  []
  
  [Transfers]
    [heat_source_from_openmc]
      type = MultiAppGeneralFieldShapeEvaluationTransfer
      from_multi_app = openmc
      variable = heat_source
      source_variable = heating
      from_postprocessors_to_be_preserved = heating_total
      to_postprocessors_to_be_preserved = source_integral
    []
  
    [temp_to_openmc]
      type = MultiAppGeneralFieldShapeEvaluationTransfer
      to_multi_app = openmc
      variable = temp
      source_variable = temp
      execute_on = 'initial timestep_begin'
    []
  []
  
  [Convergence]
    [fp_conv]
      type = DefaultMultiAppFixedPointConvergence
      fixed_point_max_its = 100
      custom_pp = source_integral
      direct_pp_value = false
      custom_rel_tol = 0.00001
      verbose = true
    []
  []
  
  [Postprocessors]
    [source_integral] # openmc heating
      type = ElementIntegralVariablePostprocessor
      variable = heat_source
      execute_on = transfer
    []
  
    [first_wall_area]
      type = AreaPostprocessor
      boundary = 'first_wall_flux_surface'
      execute_on = timestep_end
    []
  
    [first_wall_flux_power]
      type = ParsedPostprocessor
      expression = '0.5e6*first_wall_area'
      execute_on = timestep_end
      pp_names = 'first_wall_area'
    []
  
    [heating_first_wall_thermal]
      type = ElementIntegralVariablePostprocessor
      variable = heat_source
      block = 'first_wall'
      execute_on = timestep_end
    []
  
    [heating_structure_thermal]
      type = ElementIntegralVariablePostprocessor
      variable = heat_source
      block = 'structure_front structure_rear'
      execute_on = timestep_end
    []
  
    [heating_all_multiplier_thermal]
      type = ElementIntegralVariablePostprocessor
      variable = heat_source
      block = 'Multiplier_0 Multiplier_1 Multiplier_2 Multiplier_3 Multiplier_4 Multiplier_5'
      execute_on = timestep_end
    []
  
    [heating_all_breeder_thermal]
      type = ElementIntegralVariablePostprocessor
      variable = heat_source
      block = 'Breeder_00_0 ...'
      execute_on = timestep_end
    []
  []

openmc.i

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh_tagged_m.msh
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  scaling = 100.0
  initial_properties = xml
  volume_calculation = vol
  verbose = true
  source_strength = 8.9e16
  lowest_cell_level = 1
  temperature_blocks = 'structure_front structure_rear first_wall... '

  [Tallies]
    [issf]
      type = MeshTally
      mesh_template = mesh_tagged_m.msh
      score = 'heating H3_production'
      output = unrelaxed_tally_std_dev
    []
  []
[]

[UserObjects]
  [vol]
    type = OpenMCVolumeCalculation
    n_samples = 1000000
  []

  [moab]
    type = MoabSkinner
    verbose = true

    temperature = temp
    n_temperature_bins = 50
    temperature_min = 0
    temperature_max = 2500
    build_graveyard = true
    output_skins = true
    material_names = 'SS316LN_Front SS316LN_Rear Tungsten...'
  []
[]

[Postprocessors]
  [heating_total]
    type = ElementIntegralVariablePostprocessor
    variable = heating
  []

  [tritium_production_total]
    type = ElementIntegralVariablePostprocessor
    variable = H3_production
  []

  [heating_first_wall_neutronics]
    type = ElementIntegralVariablePostprocessor
    variable = heating
    block = 'first_wall'
  []

  [heating_structure_neutronics]
    type = ElementIntegralVariablePostprocessor
    variable = heating
    block = 'structure_front structure_rear'
  []

  [heating_all_multiplier_neutronics]
    type = ElementIntegralVariablePostprocessor
    variable = heating
    block = 'Multiplier_0 Multiplier_1 Multiplier_2 Multiplier_3 Multiplier_4 Multiplier_5'
  []

  [heating_all_breeder_neutronics]
    type = ElementIntegralVariablePostprocessor
    variable = heating
    block = 'Breeder_00_0 Breeder_01_0...'
  []

  [tritium_error]
    type = TallyRelativeError
    tally_score = H3_production
    value_type = average
  []

  [heating_error]
    type = TallyRelativeError
    tally_score = heating
    value_type = average
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
[]

[Outputs]
  exodus = true
  csv = true
  file_base = 'cosim_openmc'
[]

[GiudGiud]

Hello

I believe the time steps are used to perform fixed point iterations (instead of nesting fixed point iterations in a time step) in Cardinal.
Do you have the log? I'd be interested to see what ' DefaultMultiAppFixedPointConvergence' reports, and what values source_integral takes throughout the simulation

[megan-crocker]

Hi, thank you for the quick reply. Here is the log for a simulation I performed with a convergence of 1e-5, it did about 24 iterations before I terminated it because I wasn't getting any convergence.
out_thermal_cosim.txt

[megan-crocker]

I did suspect that the issue may be related to the fact I've used a steady state solver to nest picard iterations within one timestep - is it better then to do a transient simulation and execute the multi-app for each step of the parent solve?

[GiudGiud]

Moose should be able to do either.
Cardinal I m not sure. They did report some lag issues and I don't think it was all addressed
@nuclearkevin thoughts?

[GiudGiud]

ok so source_integral, the convergence critierion, is computed from a quantity that varies stochastically no?
It's the integral value of 'heating' after its projection onto the MOOSE mesh before renormalization.

if so then its variation from the stochasticity of the OpenMC simulation would prevent convergence.

Can you instead base the convergence on a normalized quantity? And adapt the desired tolerance to the stochastic variation observed as well.

I think people probably just use the k_eff. This is not the best choice. I think the relative change in the fission source might be better

[nuclearkevin]

Moose should be able to do either. Cardinal I m not sure. They did report some lag issues and I don't think it was all addressed @nuclearkevin thoughts?

The issue you're thinking about only occurs when adaptivity is active (which we've partially fixed with the EXEC_POST_ADAPTIVITY) flag. Still need to extend that flag to the output system as well, but it hasn't been that high on my list of priorities at the moment.

[nuclearkevin]

@megan-crocker have you tried using constant relaxation on the mesh tally? Oscillations between a stochastic solver and a deterministic solver are expected (especially if the coupling between physics is weak) as the fields aren't statistically converged. Relaxation (which averages the response from the stochastic solver over Picard iterations) dampen these cycles and forces convergence. See the "Other Features" section here for more information.

[megan-crocker]

Okay, I will try this. Thank you! :)
I have currently set it up so that the iterations occur with a single steady state solve (one 'timestep') - is this the correct way to do it or is it preferable to iterate instead with evolving timesteps ? Or does it not make a difference?

[megan-crocker]

Unfortunately I am still getting oscillations even with a relaxation parameter of 0.3. The Picard residuals are reducing, though, where they weren't before. I don't fully understand this behaviour, but is a reduction of the Picard residual the usual way to measure convergence of a tightly coupled system? Perhaps I have just made an error using a stochastically varying quantity as a measure of convergence.

[GiudGiud]

but is a reduction of the Picard residual the usual way to measure convergence of a tightly coupled system?

There are many ways to measure convergence of the tightly coupled system. you picked a postprocessor on the norm of the power density from the MC calculation. I would not call this standard or usual though I don't have the Cardinal experience others may have. I suspect watching the convergence of the keff is a lot more usual.

[nuclearkevin]

Okay, I will try this. Thank you! :) I have currently set it up so that the iterations occur with a single steady state solve (one 'timestep') - is this the correct way to do it or is it preferable to iterate instead with evolving timesteps ? Or does it not make a difference?

When it comes to the way we often do Picard iterations in Cardinal, a Transient executioner in the main app is used with steady-state detection enabled (bypassing MOOSE's Picard system). Functionally, this should be the same as using the Picard system with convergence based on the Picard residual.

Unfortunately I am still getting oscillations even with a relaxation parameter of 0.3. The Picard residuals are reducing, though, where they weren't before. I don't fully understand this behaviour, but is a reduction of the Picard residual the usual way to measure convergence of a tightly coupled system? Perhaps I have just made an error using a stochastically varying quantity as a measure of convergence.

I'm not an expert on convergence of Picard iterations including a stochastic quantity, but from my understanding there will still be some residual statistical error. Ideally, you should measure convergence based on the temperature (which smooths out statistical fluctuations), though you'll find that convergence is quite slow even when measuring based on temperature (and so you'll need to pick an larger convergence threshold or select the right number Picard iterations). I'll defer to @meltawila when it comes to selecting the optimal Picard convergence parameters for a fusion case, as this is a little outside of my wheelhouse.