Problems with time-step convergence

[TLWise]

I am attempting to model the following coupled equations:

dN1/dt = -1e-18 * N1

dN2/dt = 1e-18 * N1 - 1e-7 * N2

The solution iterations of the problem do not converge for the Newton method timesteps. Would you happen to have a recommendation for how I should formulate the timestepper for this type of problem?

[GiudGiud]

Hello

These coefficients are tiny.
Try automatic scaling first to see if it does the trick.

And if not then it gets complicated. The right way to do this is to use matrix exponentiation techniques iirc.

[GiudGiud]

ICs are setting the value at time 0.

I don't see how they would not. Can you make a minimal non working example that shows them failing?

In doing so you'll probably identify the problem

[TLWise]

I will share the code here. Thanks!
Simple_238U_Transmutation_Decay.i.txt

[GiudGiud]

Which one is set with N1=1e24 ?

And what are these units of time? Seems like big dts to me

[TLWise]

~one year in seconds for the dt. U238 IC is the number of atoms U238 (the value > 1e24) for less than one kg of U238.

[GiudGiud]

I ll try to run this tomorrow.

A few things:

• the simple predictor is likely not doing anything good here. You need to create a custom predictor, or better yet, a custom preconditioner since it's known how to solve the Bateman equations
• the first order monomial variables are not a good choice at this point. This is going to be constant as there is no spatial variation, so constant monomial suffices

[GiudGiud]

ok this took me a while but I think this should work.

• I fixed the lambda coefficients. I got them from google AI actually so maybe they are wrong now haha
• I re-wrote the equation in a logarithm form. in that form the terms are all much more reasonable numerically. I used AD-parsed materials to write the reaction terms. It's not the only option to do this. Note that we no longer support a 0-valued initial guess

Please let me know if you made progress on a custom predictor or a custom preconditioner

# Uranium decay chain by transmutation and atomic decay

# mpiexec -n 4 ../../crab-opt -i Simple_U238_Transmutation_Decay.i --dump

#************************* Begin Decay Constants ************************************
lambda_U238  = 4.9e-18  # [1/sec] ***** U238 *****
lambda_Th234 = 3.33e-7 # [1/sec] ***** Th234 *****
lambda_Pa234 = 2.8e-5 # [1/sec] ***** Pa234 *****
#lambda_234U = 8.97e0  # [1/sec] ***** 234U *****
#************************* End Decay Constants ************************************

# t = 3.15e17 # seconds


[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1.0
    nx = 1
    ny = 1
#    bias_x = 1.0125
  []
#  [refine]
#    type = RefineBlockGenerator
#    input = basic_mesh
#    block = 'Fluid'
#    refinement = '0'
#  []
[]

# These variables represent the logarithms
[Variables]
  [U238]
    order = CONSTANT
    family = MONOMIAL
  []
  [Th234]
    order = CONSTANT
    family = MONOMIAL
  []
  [Pa234]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[ICs]
  [U238]
    type = ConstantIC
    value = '${fparse log(2.53e31)}'
    variable = U238
  []
  [Th234]
    type = ConstantIC
    value = -1
    variable = Th234
  []
  [Pa234]
    type = ConstantIC
    value = -1
    variable = Pa234
  []
[]

[Kernels]

#************************* U238 ********************************************
  [./chngrate_U238]
    type     = TimeDerivative
    variable = U238
  [../]
  [rr_tran_decay_U238]
    type = BodyForce
    variable = U238
    value = ${fparse -lambda_U238}
  []
#************************* Th234 ********************************************
  [./chngrate_Th234]
    type     = TimeDerivative
    variable = Th234
  [../]
  [rr_coupl_Th234_to_U238]
    type = ADMatBodyForce
    variable = Th234
    material_property = 'gain_Th34_from_U8'
  []
  [rr_tran_decay_Th234]
    type = BodyForce
    variable = Th234
    value = ${fparse -lambda_Th234}
  []
#************************* Pa234 ********************************************
  [./chngrate_Pa234]
    type     = TimeDerivative
    variable = Pa234
  [../]
  [rr_coupl_Pa234_to_Th234]
    type = ADMatBodyForce
    variable = Pa234
    material_property = 'gain_Pa34_from_Th34'
  []
  [rr_tran_decay_Pa234]
    type = BodyForce
    variable = Pa234
    value = ${fparse -lambda_Pa234}
  []
[]

[AuxVariables]

[]

[AuxKernels]

[]

[BCs]

[]

[Materials]
  [gain_Th34_from_U8]
    type = ADParsedMaterial
    expression = "lambda_U238 * exp(U238) / exp(Th234)"
    property_name = 'gain_Th34_from_U8'
    coupled_variables = 'U238 Th234'
    functor_symbols = "lambda_U238"
    functor_names = '${lambda_U238}'
    outputs = 'all'
  []
  [gain_Th34]
    type = ADParsedMaterial
    expression = "lambda_Th234 * exp(Th234) / exp(Pa234)"
    property_name = 'gain_Pa34_from_Th34'
    coupled_variables = 'Th234 Pa234'
    constant_names = lambda_Th234
    constant_expressions = '${lambda_Th234}'
    outputs = 'all'
  []
[]

[Functions]
  [conditional_function]
     type =  ParsedFunction
     expression = 't < 101536005'
  []
 [U238_loading_function]
     type =  ParsedFunction
     expression = 't < 101536005'
 []
[]

# [Controls]
#   [initial_step]
#     type = ConditionalFunctionEnableControl
#     conditional_function = conditional_function
#     disable_objects = 'Kernels::rr_tran_decay_U238 Kernels::rr_coupl_Th234_to_U238 Kernels::rr_tran_decay_Th234'
# #    enable_objects = 'Kernels::rr_tran_decay_U2380'
#     execute_on = 'INITIAL TIMESTEP_END'
#   []
# []

[Preconditioning]
  [SMP]
    type = SMP
    full = true
    # petsc_options_iname = '-pc_type -pc_factor_shift_type'
    # petsc_options_value = 'lu       NONZERO'
  []
[]

#**************************************** FEM Example Executioner *****************************************
[Executioner]
#  execute_on = 'initial timestep_begin'
  type = Transient
  solve_type = 'NEWTON'
#  solve_type = 'PJFNK'

  # nl_rel_tol = 1e-5
  nl_abs_tol = 1e-8
  nl_max_its = 1000
  l_max_its = 10000
  nl_forced_its = 2

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  line_search = 'none'

#  automatic_scaling = true
#  line_search = none
#  ignore_variables_for_autoscaling = ''
#  compute_scaling_once = false
#  scaling_group_variables = 'U238 Th234 Pa234'  # Automatic Scaling of variables for better initial timestep guess.
#  end_time = 3600 # 1 hour
#  dtmax = 300 # 5 minutes
#  dtmin = 1.0 # 1 sec
#  end_time = 2.592e5 # 30 days
#  dtmax = 3600 # 1 hour
#  dtmin = 60 # 1 minute
#  end_time = 3.15e7 # 1 year
#  dtmax = 86400 # 1 day
#  dtmin = 3600  # 1 hour
#  end_time = 3153600000 # 1000 years
#  dtmax = 3153600000 # 100 yrs
#  dtmin = 315360  # 0.01 yrs
#  scheme = 'bdf2'
# dt = 2e7
  [TimeStepper]
    type = IterationAdaptiveDT
    # dt = ${fparse 31536000 / 10}
    dt = 2e2
    optimal_iterations = 10
    growth_factor = 1.25
  []

  end_time = ${fparse 100 * 365.25 * 24 * 3600}
#  num_steps = 100000
#  dt = 3153600000
  # [Predictor]
  #   type = SimplePredictor
  #   scale = 1.0
  # []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [./U238]
    type = ElementAverageValue
    variable = U238
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./Th234]
    type = ElementAverageValue
    variable = Th234
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./Pa234]
    type = ElementAverageValue
    variable = Pa234
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

[Outputs]
  exodus = true
  [out]
    type = CSV
#    execute_on = FINAL
  []
[]

[AuxVariables]
  [U238_out]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = ParsedAux
      expression = 'exp(U238)'
      coupled_variables = 'U238'
    []
  []
  [Th234_out]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = ParsedAux
      expression = 'exp(Th234)'
      coupled_variables = 'Th234'
    []
  []
  [Pa234_out]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = ParsedAux
      expression = 'exp(Pa234)'
      coupled_variables = 'Pa234'
    []
  []
[]

[TLWise]

Thank you for taking time to understand what I provided and for applying your knowledge and efforts to finding a potential solution to my proposed problem. I look forward to applying what you developed for use in my other future efforts.