Martensite grain and phase nucleation in austenite matrix

[sfarkya]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello,

I am new to this area and am trying to implement nucleation of martensite grains within an austenite matrix.

Model Setup

I start with a fully austenitic phase, with austenite grains defined as OPs gr0–gr4.

Based on the transformation temperatures (Ms, Mf), I define a probability of nucleation.

I reserved:

gr8 for nucleation

gr5–gr7 for remapping the newly formed nucleus into martensite grains

I use bulk free energies for FCC (austenite) and BCC (martensite) from TDB files to drive the η evolution.

For grain nucleation and grain growth, I am using the Grain Growth Model from MOOSE.

My current goal is to stabilize martensite grains after nucleation.

Issue

Nucleation appears to occur as expected. However, I am having trouble with the GrainTracker:

The GrainTracker does not seem to correctly map nucleated grains (from gr8) into the reserved OPs (gr5–gr8).

I do observe grain allocation to gr5–gr8, but I am unclear how this assignment is happening.

It appears that multiple OPs may be representing the same grain.

Questions

How can I ensure that after nucleation in gr8, the grains are correctly remapped into gr5–gr8?

What is the recommended approach to stabilizing a newly nucleated martensite grain?

Are there specific considerations when using the Grain Growth Model and GrainTracker together for phase transformation problems?

Any insights or guidance on using the Grain Growth Model or GrainTracker for this type of phase transformation problem would be greatly appreciated.

Thank you, and apologies if my terminology is not precise — I am still learning and trying to understand the framework better.

Attached my current script and execution:
main_out.csv

`# ============================================================
# Phase evolution only (eta): Austenite (FCC_A1) <-> BCC proxy
# Tracks FCC/BCC free energies as fields (placeholders)
# You will fill the expressions for G_fcc, G_bcc, and f_eta.
# ============================================================

# =========================
# Mesh
# =========================
[Mesh]
  type = GeneratedMesh
  dim = 2

  # smaller mesh for testing (9 x 9 µm^2 with dx = 0.05 µm -> 180 x 180)
  nx = 180
  ny = 180
  xmax = 9.0
  ymax = 9.0
  elem_type = QUAD4
  uniform_refine = 0
[]

[GlobalParams]
# Parameters used by several kernels that are defined globally to simplify input file
  op_num = 9 # Number of order parameters used
  var_name_base = gr # Base name of grains
  penalty = 1
[]

# =========================
# Variables (ONLY eta evolves)
# =========================
[Variables]

  [./PolycrystalVariables]
  [../]

  # Phase indicator: eta = 0 (austenite), eta = 1 (BCC/martensite proxy)
  [./eta]
    order = FIRST
    family = LAGRANGE
  [../]
[]

# =========================
# Initial conditions
# =========================
[ICs]
  [PolycrystalICs]
    [PolycrystalColoringIC]
      polycrystal_ic_uo = voronoi
    []
  [] 

  [./eta_ic]
    type = ConstantIC
    variable = eta
    value = 0.0
  [../]
  
    [./bnds]
    type = BndsCalcIC
    variable = bnds
  [../]

  # # --- tiny noise on reserved BCC OPs to break symmetry ---
  # [./gr5_noise]
  #   type = PolycrystalRandomIC
  #   op_index = 5
  #   seed = 5005
  #   variable = 'gr5'
  #   random_type = 'discrete'
  # []

  # [./gr6_noise]
  #   type = PolycrystalRandomIC
  #   op_index = 6
  #   seed = 5006
  #   variable = 'gr6'
  #   random_type = 'discrete'
  # []

  # [./gr7_noise]
  #   type = PolycrystalRandomIC
  #   op_index = 7
  #   seed = 5007
  #   variable = 'gr7'
  #   random_type = 'discrete'
  # []

[]
[UserObjects]
  [voronoi]
    type = PolycrystalVoronoi
    grain_num = 5 # Number of grains
    rand_seed = 10
    int_width = 0.3 # this is the width of the interface between the grains # ball park value # not sure if this is correct
    op_num = 5
  []
  [grain_tracker]
    type = GrainTracker
    variable = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
  []
  
  [./inserter]
    # step 1 
    # The inserter runs at the end of each time step to add nucleation events
    # that happend during the timestep (if it converged) to the list of nuclei
    type = DiscreteNucleationInserter
    hold_time = 0.00 # 1 # step 6 
    # Apply both a DiscreteNucleationForce and a Reaction kernel to a reserved order 
    # parameter and set the hold_time in the inserter to 0. This will cause the 
    # modification of the reserved OP during only a single timestep, 
    # which will be enough for the GrainTracker to pick up and remap the new nucleus.
    probability = P
    radius = 1 # 0.15 # 2 (works)
  [../]

  [./map]
    # step 3
    # The map UO runs at the beginning of a timestep and generates a per-element/qp
    # map of nucleus locations. The map is only regenerated if the mesh changed or
    # the list of nuclei was modified.
    # The map converts the nucleation points into finite area objects with a given radius.
    type = DiscreteNucleationMap
    periodic = eta
    inserter = inserter
    int_width = 0.3 # 0.1
    # int_width - Nucleus interface width for smooth nuclei
  [../]


[] 

# =========================
# AuxVariables (tracked fields)
# =========================
[AuxVariables]
  # Fixed composition for "phase evolution only"
  # [./c]
  #   order = FIRST
  #   family = LAGRANGE
  # []
  [./Xc]
    order = FIRST
    family = LAGRANGE
  []
  [./Yc_fcc]
    order = FIRST
    family = LAGRANGE
  []
  [./Yc_bcc]
    order = FIRST
    family = LAGRANGE
  []
  # grain tracking variables
    [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
  [var_indices]
    order = CONSTANT
    family = MONOMIAL
  []
  [./bnds]
    order = FIRST
    family = LAGRANGE
  [../]

  # Temperature schedule (from AuxKernel below)
  [./T]
    type = MooseVariable
  []

  [./Ffcc_total]
    order = FIRST
    family = MONOMIAL 
  []
  [./Fbcc_total]
    order = FIRST
    family = MONOMIAL
  []

  [./Ftotal_eta_aux]
    order = FIRST
    family = MONOMIAL
  []

  [./Ftotal_gr_aux]
    order = FIRST
    family = MONOMIAL
  []

  [./Ftotal_loc_aux]
    order = FIRST
    family = MONOMIAL
  []

  [./F_phase_grain_aux]
    order = FIRST
    family = MONOMIAL
  []

  [./Sfcc]
    order = FIRST
    family = LAGRANGE
  []

  [./Sbcc]
    order = FIRST
    family = LAGRANGE
  []

  [./L_fcc_aux]
    order = FIRST
    family = MONOMIAL
  []

  [./L_bcc_aux]
    order = FIRST
    family = MONOMIAL
  []

  # Optional: track local free energy density and d(f)/d(eta)
  [./f_eta_aux]
    order = FIRST
    family = LAGRANGE
  []
  [./dfdeta_aux]
    order = FIRST
    family = LAGRANGE
  []

[]



# =========================
# AuxKernels
# =========================
[AuxKernels]
  # ---- Fixed carbon (phase evolution only) ----
  # [./c_aux]
  #   type = ConstantAux
  #   variable = c
  #   value = 0.0050251 # for calphad x = 0.005
  #   execute_on = 'initial'
  # []
  [./Xc]
  type = ConstantAux
  variable = Xc
  value = 0.005
  execute_on = 'initial'
  []

  [./Yc_fcc_aux]
  type = ParsedAux
  variable = Yc_fcc
  coupled_variables = 'Xc'
  # FCC: 1 interstitial sublattice site per metal site
  expression = 'Xc/(1 - Xc)'
  execute_on = 'initial timestep_end'
[]

[./Yc_bcc_aux]
  type = ParsedAux
  variable = Yc_bcc
  coupled_variables = 'Xc'
  # BCC_A2 in most Fe-C tdbs: 3 interstitial sites per metal site
  expression = 'Xc/(3*(1 - Xc))'
  execute_on = 'initial timestep_end'
[]
  

  # ---- Quench schedule: linear cooling then hold ----
  # NOTE: Make sure units match your thermodynamics tables/expressions
  # expression = max(Tq, T0 - rate*t)
  # [./T_aux]
  #   type = ParsedAux
  #   variable = T
  #   use_xyzt = true
  #   expression = 'max(Tq, T0 - rate*t)'
  #   constant_names = 'T0 rate Tq'
  #   constant_expressions = '1300 20 300'
  #   execute_on = 'initial timestep_end'
  # []

  #[./T_aux]
  #  type = ParsedAux
  #  variable = T
  #  use_xyzt = true

  #  # Linear cool from T0 at rate until reaching Tq, then hold
  #  # expression = 'if(t < tq, T0 - rate*t, Tq)'
  #  expression = 'if(t < t_start, T0, Tinf + (T0 - Tinf)*exp(-(t - t_start)/tau))'
  #  constant_names = 'T0 Tinf t_start tau'

  #  # t_start ≈ 12.6 s
  #  # T(21.8 s) ≈ 560 C is matched with tau ≈ 13.19 s (for T0=900C, Tinf=223C)
  #  constant_expressions = '1173.15 496.15 12.6 13.19'

  #  # constant_names = 'T0 rate Tq tq'
  #  # # --- choose ONE of the following sets ---

  #  # # Quench to 223 C:
  #  # constant_expressions = '1173.15 20.0 516.15 32.85'

  #  # # Quench to 243 C:
  #  # # constant_expressions = '1173.15 20.0 516.15 32.85'

  #  # # Quench to 263 C:
  #  # # constant_expressions = '1173.15 20.0 536.15 31.85'

  #  execute_on = 'initial timestep_end'
  #[]

  [./T_aux]
    type = ParsedAux
    variable = T
    use_xyzt = true

    # Hold at T0 until t_start, then T = Tinf + (T0 - Tinf)*exp(-k*(t - t_start))
    # Matches Python plot_temperature_profiles_various_k.py: 850 C hold until 12.6 s, then exponential to 100 C, k = 0.12
    expression = 'if(t < t_start, T0, Tinf + (T0 - Tinf)*exp(-k*(t - t_start)))'

    constant_names = 'T0 Tinf t_start k'
    # T0 = 850 C, Tinf = 100 C, t_start = 12.6 s, k = 0.12 (1/s)
    constant_expressions = '1123.15 373.15 12.6 0.12'

    execute_on = 'initial timestep_end'
  []

  [./Ffcc_total]
  type = MaterialRealAux
  variable = Ffcc_total
  property = Ffcc
  execute_on = timestep_end
  [../]
  [./Fbcc_total]
  type = MaterialRealAux
  variable = Fbcc_total
  property = Fbcc
  execute_on = timestep_end
  [../]

  [./Ftotal_eta_aux]
  type = MaterialRealAux
  variable = Ftotal_eta_aux
  property = Ftotal_eta
  execute_on = timestep_end
  [../]

  [./Ftotal_gr_aux]
  type = MaterialRealAux
  variable = Ftotal_gr_aux
  property = Ftotal_gr
  execute_on = timestep_end
  [../]

  [./Ftotal_loc_aux]
  type = MaterialRealAux
  variable = Ftotal_loc_aux
  property = Ftotal_loc
  execute_on = timestep_end
  [../]

  [./F_phase_grain_aux]
  type = MaterialRealAux
  variable = F_phase_grain_aux
  property = F_phase_grain
  execute_on = timestep_end
  [../]

  [./Sfcc_aux]
    type = ParsedAux
    variable = Sfcc
    coupled_variables = 'gr0 gr1 gr2 gr3 gr4'
    expression = 'gr0+gr1+gr2+gr3+gr4'
    execute_on = 'initial timestep_end'
  []

  [./Sbcc_aux]
    type = ParsedAux
    variable = Sbcc
    coupled_variables = 'gr5 gr6 gr7 gr8'
    expression = 'gr5+gr6+gr7+gr8'
    execute_on = 'initial timestep_end'
  []

  [./L_fcc_aux]
  type = MaterialRealAux
  variable = L_fcc_aux
  property = L_fcc
  execute_on = timestep_end
  [../]

  [./L_bcc_aux]
  type = MaterialRealAux
  variable = L_bcc_aux
  property = L_bcc
  execute_on = timestep_end
  [../]

  # AuxKernel block, defining the equations used to calculate the auxvars
  [bnds_aux]
    # AuxKernel that calculates the GB term
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  [unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    flood_counter = grain_tracker
    field_display = UNIQUE_REGION
    execute_on = 'initial timestep_end'
  []
  [var_indices]
    type = FeatureFloodCountAux
    variable = var_indices
    flood_counter = grain_tracker
    field_display = VARIABLE_COLORING
    execute_on = 'initial timestep_end'
  []
[]

# =========================
# Boundary Conditions
# =========================
[BCs]
  [Periodic]
    [All]
      auto_direction = 'x y'
    []
  []
[]

# =========================
# Kernels (eta evolution only)
# =========================
[Kernels]
  [./detadt]
    type = TimeDerivative
    variable = eta
  [../]

  [./eta_interface]
    type = ACInterface
    variable = eta
    kappa_name = kappa_eta
    mob_name = L_eta
  [../]

  [./eta_bulk]
    type = AllenCahn
    variable = eta
    f_name = Ftotal_eta
    mob_name = L_eta
  [../]
  

  [dgr0dt]
    type = TimeDerivative
    variable = gr0
  []
  [ACInterface0]
    type = ACInterface
    variable = gr0
    coupled_variables = 'T gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    kappa_name = kappa_c
    mob_name = L_fcc
  []
  [ACBulk0]
    type = AllenCahn
    f_name = Ftotal_gr
    variable = gr0
    coupled_variables =  'eta gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [./penalty0]
    type = SwitchingFunctionPenalty
    variable = gr0
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]

  [dgr1dt]
    type = TimeDerivative
    variable = gr1
  []
  [ACInterface1]
    type = ACInterface
    variable = gr1
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [ACBulk1]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr1
    # v = 'gr0 gr2 gr3 gr4 gr5 gr6 gr7'
    coupled_variables ='eta gr0 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [./penalty1]
    type = SwitchingFunctionPenalty
    variable = gr1
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  # ------
  [dgr2dt]
    type = TimeDerivative
    variable = gr2
  []
  [ACInterface2]
    type = ACInterface
    variable = gr2
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr3 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [ACBulk2]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr2
    # v =  'gr0 gr1 gr3 gr4 gr5 gr6 gr7'
    coupled_variables =  'eta gr0 gr1 gr3 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [./penalty2]
    type = SwitchingFunctionPenalty
    variable = gr2
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  [dgr3dt]
    type = TimeDerivative
    variable = gr3
  []
  [ACInterface3]
    type = ACInterface
    variable = gr3
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [ACBulk3]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr3
    # v = 'gr0 gr1 gr2 gr4 gr5 gr6 gr7'
    coupled_variables = 'eta gr0 gr1 gr2 gr4 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [./penalty3]
    type = SwitchingFunctionPenalty
    variable = gr3
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  [dgr4dt]
    type = TimeDerivative
    variable = gr4
  []
  [ACInterface4]
    type = ACInterface
    variable = gr4
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr3 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [ACBulk4]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr4
    # v = 'gr0 gr1 gr2 gr3 gr5 gr6 gr7'
    coupled_variables = 'eta gr0 gr1 gr2 gr3 gr5 gr6 gr7 gr8'
    mob_name = L_fcc
  []
  [./penalty4]
    type = SwitchingFunctionPenalty
    variable = gr4
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  # --------------------------
  [dgr5dt]
    type = TimeDerivative
    variable = gr5
  []
  [ACInterface5]
    type = ACInterface
    variable = gr5
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr3 gr4 gr6 gr7 gr8'
    mob_name = L_bcc
  []
  [ACBulk5]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr5
    # v = 'gr0 gr1 gr2 gr3 gr4 gr6 gr7'
    coupled_variables = 'eta gr0 gr1 gr2 gr3 gr4 gr6 gr7 gr8'
    mob_name = L_bcc
  []
  [./penalty5]
    type = SwitchingFunctionPenalty
    variable = gr5
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  [dgr6dt]
    type = TimeDerivative
    variable = gr6
  []
  [ACInterface6]
    type = ACInterface
    variable = gr6
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr3 gr4 gr5 gr7 gr8'
    mob_name = L_bcc
  []
  [ACBulk6]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr6
    # v = 'gr0 gr1 gr2 gr3 gr4 gr5 gr7'
    coupled_variables = 'eta gr0 gr1 gr2 gr3 gr4 gr5 gr7 gr8'
    mob_name = L_bcc
  []
  [./penalty6]
    type = SwitchingFunctionPenalty
    variable = gr6
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  [dgr7dt]
    type = TimeDerivative
    variable = gr7
  []
  [ACInterface7]
    type = ACInterface
    variable = gr7
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr8'
    mob_name = L_bcc
  []
  [ACBulk7]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr7
    # v = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6'
    coupled_variables = 'eta gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr8'
    mob_name = L_bcc
  []
  [./penalty7]
    type = SwitchingFunctionPenalty
    variable = gr7
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  [dgr8dt]
    type = TimeDerivative
    variable = gr8
  []
  [ACInterface8]
    type = ACInterface
    variable = gr8
    kappa_name = kappa_c
    coupled_variables = 'T gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7'
    mob_name = L_bcc
  []
  [ACBulk8]
    type = AllenCahn
    f_name = Ftotal_gr
    # type = ACGrGrPoly
    variable = gr8
    # v = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6'
    coupled_variables = 'eta gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7'
    mob_name = L_bcc
  []
  [./penalty8]
    type = SwitchingFunctionPenalty
    variable = gr8
    etas = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    h_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
  [../]
  
  #option 1: use reaction kernel to nucleate gr8
  
  [rxn_gr8]
    type = Reaction
    variable = gr8
    rate = 1
  []
  [./nucleation_gr8]
    type = DiscreteNucleationForce
    map = map
    variable = gr8
  [../]
[]

# =========================
# Materials
# =========================
[Materials]

  # ---- Phase-field parameters for eta ----
  [./eta_params]
    type = GenericConstantMaterial
    prop_names  = 'L_eta kappa_eta' # convergence issues with L_eta=1
    prop_values = '1e-4 1.0' # 1e-3 very fast 
  [../]

  # https://mooseframework.inl.gov/source/materials/GBEvolution.html
  [./consts]
    type = GenericConstantMaterial
    prop_names = 'gamma_asymm'
    prop_values = '1.5' 
  []
  
  # https://github.com/idaholab/moose/blob/next/modules/phase_field/src/materials/GBEvolutionBase.C#L14
  # [./mobility] 
  #   type = GenericFunctionMaterial
  #   prop_names = 'L'
  #   # prop_values = '(4.0/3.0)*(mgb/wGB)'
  #   prop_values = '(4.0/3.0)*(((2.5e-6*1e-9)/(6.24150974e18*1e-36)*exp(-0.23/(8.617333262e-5*450)))/14.0)'
  # []
  # ---------------------------
  # mobility: L = (4/3) * (mgb / wGB)
  # (FORMULA UNCHANGED)
  #
  # IMPORTANT: wGB must be an actual length, not a bare "14".
  # We keep your "14" as number of interface points and convert it to µm using dx_um.
  # Keep your mgb expression exactly as-is (do not touch your conversions),
  # but divide by a physical wGB in the same length system (µm).
  # wGB_um = wGB_pts * dx_um
  # ---------------------------
  
  # [./mobility]
  #   type = GenericFunctionMaterial
  #   prop_names = 'L'
  #   #constant_names = 'x_um wGB_pts'
  #   #constant_expressions = '0.05 14' # multipied by 1e-3 
  #   prop_values = '4e-4*(4.0/3.0)*(((2.5e-6*1e-9)/(6.24150974e18*1e-36)*exp(-0.23/(8.617333262e-5*1100)))/(14*0.05))'
  # [../]


  # [./kappa_c]
  #   type = GenericFunctionMaterial
  #   prop_names = 'kappa_c'
  #   prop_values = '3*(3/4)*(0.708*6.24150974e18*1e-18)*14'
  # []
  # ---------------------------
  # kappa_c = (3/4) * gamma * wGB
  # (FORMULA UNCHANGED)
  #
  # Target units: J*um^2/mol (consistent with F in J/mol and x in µm)
  #
  # gamma input is in J/m^2 -> convert to J/um^2 by *1e-12
  # then convert J/um -> J*um^2/mol using molar volume:
  #   (J/um) * (1e18*Vm_m3_mol) = J*um^2/mol
  # ---------------------------
  [./kappa_c]
    type = GenericFunctionMaterial
    prop_names = 'kappa_c'  # multiply by 3
    prop_values = '(3.0/4.0)*(0.708*1e-12)*(14*0.05)*(1e18*7.10e-6)'

    # constant_names = 'gamma_J_m2 Vm_m3_mol dx_um wGB_pts'
    # constant_expressions = '0.708 7.10e-6 0.05 14'
    # kappa = (3/4)*gamma*wGB
    # gamma(J/m^2)->gamma(J/um^2): *1e-12
    # wGB(um): wGB_pts*dx_um
    # convert J/um -> J*um^2/mol: *(1e18*Vm_m3_mol)
    # prop_values = '(3.0/4.0)*(gamma_J_m2*1e-12)*(wGB_pts*dx_um)*(1e18*Vm_m3_mol)'
  [../]

  # mu = 6 * gamma / wGB
  # (FORMULA UNCHANGED)
  #
  # Target units: J/mol
  #
  # gamma(J/m^2)->gamma(J/um^2): *1e-12
  # gamma/wGB => J/um^3
  # convert J/um^3 -> J/mol using molar volume:
  #   (J/um^3) * (1e18*Vm_m3_mol) = J/mol
  # ---------------------------
  [./mu]
    type = GenericFunctionMaterial
    prop_names = 'mu'
    # constant_names = 'gamma_J_m2 Vm_m3_mol dx_um wGB_pts'
    # constant_expressions = '0.708 7.10e-6 0.05 14'
    prop_values = '(6.0*(0.708*1e-12)/(14*0.05))*(1e18*7.10e-6)'
    # prop_values = '(6.0*(gamma_J_m2*1e-12)/(wGB_pts*dx_um))*(1e18*Vm_m3_mol)'

    # mu = 6*gamma/wGB
    # gamma(J/m^2)->gamma(J/um^2): *1e-12
    # /wGB(um): /(wGB_pts*dx_um) => J/um^3
    # J/um^3 -> J/mol: *(1e18*Vm_m3_mol)
  [../]

  # [./mu]
  #   type = GenericFunctionMaterial
  #   prop_names = 'mu'
  #   prop_values = '(6*(0.708*6.24150974e18*1e-18))/14'
  # #   # prop_values = '6*(GBenergy/wGB)'
  # # []

  # ===================================
  # Eta switching and barrier function 
  # ===================================
  [./switching]
    type = SwitchingFunctionMaterial
    block = 0
    eta = eta
    h_order = SIMPLE
    outputs = exodus
  [../]

  [./barrier]
    type = BarrierFunctionMaterial
    block = 0
    eta = eta
    g_order = SIMPLE
    outputs = exodus
  [../]
 
  # =============================
  # Switching function Multiplhase 
  # =============================

  [./hsum]
    type = ParsedMaterial
    expression = h0+h1+h2+h3+h4+h5+h6+h7+h8
    property_name = hsum
    material_property_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
    # coupled_variables = 'c'
    outputs = exodus
  [../]
  [./switching0]
    type = SwitchingFunctionMaterial
    function_name = h0
    eta = gr0
    h_order = SIMPLE
  [../]
  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = gr1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = gr2
    h_order = SIMPLE
  [../] 
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = gr3
    h_order = SIMPLE
  [../]
  [./switching4]
    type = SwitchingFunctionMaterial
    function_name = h4
    eta = gr4
    h_order = SIMPLE
  [../]
  [./switching5]
    type = SwitchingFunctionMaterial
    function_name = h5
    eta = gr5
    h_order = SIMPLE
  [../]
  [./switching6]
    type = SwitchingFunctionMaterial
    function_name = h6
    eta = gr6
    h_order = SIMPLE
  [../]
  [./switching7]
    type = SwitchingFunctionMaterial
    function_name = h7
    eta = gr7
    h_order = SIMPLE
  [../]
  [./switching8]
    type = SwitchingFunctionMaterial
    function_name = h8
    eta = gr8
    h_order = SIMPLE
  [../]
# =============================================================
# Mobility
  [./mobility_fcc]
    type = DerivativeParsedMaterial
    property_name = L_fcc
    coupled_variables = 'T'
    expression = '4e-4*(4.0/3.0)*(((2.5e-6*1e-9)/(6.24150974e18*1e-36)*exp(-0.23/(8.617333262e-5*T)))/(14*0.05))'
  [../]
    
  [./mobility_bcc]
    type = DerivativeParsedMaterial
    property_name = L_bcc
    coupled_variables = 'T'
    expression = '1e-3*(4.0/3.0)*(((2.5e-6*1e-9)/(6.24150974e18*1e-36)*exp(-0.23/(8.617333262e-5*T)))/(14*0.05))'
  [../]
# ============================================================
# Interface bulk energy
# ============================================================
  [./total_free_energy]
  type = DerivativeParsedMaterial
  block = 0
  property_name = Ftotal_loc
  material_property_names = 'mu gamma_asymm h'
  coupled_variables = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8 eta'
  expression = '
    Sfcc := gr0^2+gr1^2+gr2^2+gr3^2+gr4^2;
    Sbcc := gr5^2+gr6^2+gr7^2+gr8^2;
    Wfcc := (gr0^4/4-gr0^2/2)+(gr1^4/4-gr1^2/2)+(gr2^4/4-gr2^2/2)+(gr3^4/4-gr3^2/2)+(gr4^4/4-gr4^2/2);
    Pfcc := gr0^2*gr1^2+gr0^2*gr2^2+gr0^2*gr3^2+gr0^2*gr4^2
          + gr1^2*gr2^2+gr1^2*gr3^2+gr1^2*gr4^2
          + gr2^2*gr3^2+gr2^2*gr4^2
          + gr3^2*gr4^2;
    Wbcc := (gr5^4/4-gr5^2/2)+(gr6^4/4-gr6^2/2)+(gr7^4/4-gr7^2/2)+(gr8^4/4-gr8^2/2);
    Pbcc := gr5^2*gr6^2+gr5^2*gr7^2+gr5^2*gr8^2
          + gr6^2*gr7^2+gr6^2*gr8^2
          + gr7^2*gr8^2;
    Ffcc := (1-h)*mu*(Wfcc + gamma_asymm*Pfcc + 0.25);
    Fbcc := h*mu*(Wbcc + gamma_asymm*Pbcc + 0.25);
    Ffcc + Fbcc 
  '
#  Ffcc + Fbcc + Fcross
# Fcross := mu*gamma_asymm*(Sfcc*Sbcc);
# # choose ONE of the following two lines:
# # (A) always-on cross-family penalty (your current intent)
# Fcross := mu*gamma_asymm*(Sfcc*Sbcc);
# # (B) gated cross-family penalty (only active in eta transition zone)
# # Fcross := mu*gamma_asymm*h*(1-h)*(Sfcc*Sbcc);
# Ffcc + Fbcc + Fcross

  derivative_order = 2
  enable_jit = true
  outputs = exodus
  [../]

  # # Combined local free energy for all grains (gr0-gr8) with phase switching
  # [./total_free_energy]
  #   type = DerivativeParsedMaterial
  #   block = 0
  #   property_name = Ftotal_loc
  #   material_property_names = 'mu gamma_asymm h'
  #   coupled_variables = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8 eta'
  #   # expression = '(1-h)*mu*(((gr0^4/4)-(gr0^2/2))+((gr1^4/4)-(gr1^2/2))+((gr2^4/4)-(gr2^2/2))+((gr3^4/4)-(gr3^2/2))+((gr4^4/4)-(gr4^2/2))+((gr5^4/4)-(gr5^2/2))+((gr6^4/4)-(gr6^2/2))+((gr7^4/4)-(gr7^2/2))+gamma_asymm*(gr0^2*gr1^2+gr0^2*gr2^2+gr0^2*gr3^2+gr0^2*gr4^2+gr0^2*gr5^2+gr0^2*gr6^2+gr0^2*gr7^2+gr1^2*gr2^2+gr1^2*gr3^2+gr1^2*gr4^2+gr1^2*gr5^2+gr1^2*gr6^2+gr1^2*gr7^2+gr2^2*gr3^2+gr2^2*gr4^2+gr2^2*gr5^2+gr2^2*gr6^2+gr2^2*gr7^2+gr3^2*gr4^2+gr3^2*gr5^2+gr3^2*gr6^2+gr3^2*gr7^2+gr4^2*gr5^2+gr4^2*gr6^2+gr4^2*gr7^2+gr5^2*gr6^2+gr5^2*gr7^2+gr6^2*gr7^2)+0.25) + h*mu*(((gr8^4/4)-(gr8^2/2))+gamma_asymm*gr8^2*(gr0^2+gr1^2+gr2^2+gr3^2+gr4^2+gr5^2+gr6^2+gr7^2)+0.25)'
  #   expression = '(1-h)*mu*(
  #                  ((gr0^4/4)-(gr0^2/2))+((gr1^4/4)-(gr1^2/2))+((gr2^4/4)-(gr2^2/2))+((gr3^4/4)-(gr3^2/2))+((gr4^4/4)-(gr4^2/2))
  #                  + gamma_asymm*(
  #                    gr0^2*gr1^2+gr0^2*gr2^2+gr0^2*gr3^2+gr0^2*gr4^2
  #                   +gr1^2*gr2^2+gr1^2*gr3^2+gr1^2*gr4^2
  #                   +gr2^2*gr3^2+gr2^2*gr4^2
  #                   +gr3^2*gr4^2
  #                  )
  #                  + 0.25
  #                )
  #                + h*mu*(
  #                  ((gr5^4/4)-(gr5^2/2))+((gr6^4/4)-(gr6^2/2))+((gr7^4/4)-(gr7^2/2))+((gr8^4/4)-(gr8^2/2))
  #                  + gamma_asymm*(
  #                    gr5^2*gr6^2+gr5^2*gr7^2+gr5^2*gr8^2
  #                   +gr6^2*gr7^2+gr6^2*gr8^2
  #                   +gr7^2*gr8^2
  #                  )
  #                  + 0.25
  #                )
  #                + mu*gamma_asymm*(
  #                  (gr0^2+gr1^2+gr2^2+gr3^2+gr4^2)
  #                 *(gr5^2+gr6^2+gr7^2+gr8^2)
  #                )'
  #   
  #   derivative_order = 2
  #   outputs = exodus
  # [../]

  # much slower - could be made fast if written in C.
  # [./total_free_energy] 
  #   type = DerivativeParsedMaterial
  #   property_name = Ftotal_loc
  #   # Use material properties directly (no duplicate mu expression here)
  #   material_property_names = 'mu gamma_asymm'
  #   coupled_variables = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
  #   expression = 'mu*(((gr0^4/4)-(gr0^2/2))+((gr1^4/4)-(gr1^2/2))+((gr2^4/4)-(gr2^2/2))+((gr3^4/4)-(gr3^2/2))+((gr4^4/4)-(gr4^2/2))+((gr5^4/4)-(gr5^2/2))+((gr6^4/4)-(gr6^2/2))+((gr7^4/4)-(gr7^2/2))+((gr8^4/4)-(gr8^2/2))+gamma_asymm*(gr0^2*gr1^2+gr0^2*gr2^2+gr0^2*gr3^2+gr0^2*gr4^2+gr0^2*gr5^2+gr0^2*gr6^2+gr0^2*gr7^2+gr0^2*gr8^2+gr1^2*gr2^2+gr1^2*gr3^2+gr1^2*gr4^2+gr1^2*gr5^2+gr1^2*gr6^2+gr1^2*gr7^2+gr1^2*gr8^2+gr2^2*gr3^2+gr2^2*gr4^2+gr2^2*gr5^2+gr2^2*gr6^2+gr2^2*gr7^2+gr2^2*gr8^2+gr3^2*gr4^2+gr3^2*gr5^2+gr3^2*gr6^2+gr3^2*gr7^2+gr3^2*gr8^2+gr4^2*gr5^2+gr4^2*gr6^2+gr4^2*gr7^2+gr4^2*gr8^2+gr5^2*gr6^2+gr5^2*gr7^2+gr5^2*gr8^2+gr6^2*gr7^2+gr6^2*gr8^2+gr7^2*gr8^2)+0.25)'
  #   derivative_order = 2
  # [../]

# ============================================================
# Free energy material for fcc
# ============================================================

  [./free_energy_fcc]
    type = DerivativeParsedMaterial
    block = 0
    property_name = Ffcc
    expression='FCC_C:=Yc_fcc; FCC_VA:=1-Yc_fcc; scale*8.3145*T*(1.0*if(FCC_FE > 1.0e-15,FCC_FE*log(FCC_FE),0) + 1.0*if(FCC_C > 1.0e-15,FCC_C*log(FCC_C),0) + 1.0*if(FCC_VA > 1.0e-15,FCC_VA*log(FCC_VA),0))/(1.0*FCC_FE + 1.0*FCC_C) + 8.3145*T*if(T < 67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA,-0.00174491240370669*T^15/(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^15 - 0.00775516623869638*T^9/(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^9 - 0.174491240370669*T^3/(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^3 + 1 - 0.86033875460538*(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)/T,if(T < -201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA,-0.00174491240370669*T^15/(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^15 - 0.00775516623869638*T^9/(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^9 - 0.174491240370669*T^3/(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^3 + 1 - 0.86033875460538*(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)/T,if(T > -201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA & 201.0*FCC_FE*FCC_C + 201.0*FCC_FE*FCC_VA < 0,-0.0426902268107986*(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^5/T^5 - 0.00135524529558091*(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^15/T^15 - 0.000284601512071991*(-201.0*FCC_FE*FCC_C - 201.0*FCC_FE*FCC_VA + 1.0e-9)^25/T^25,if(T > 67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA & 201.0*FCC_FE*FCC_C + 201.0*FCC_FE*FCC_VA > 0,-0.0426902268107986*(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^5/T^5 - 0.00135524529558091*(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^15/T^15 - 0.000284601512071991*(67.0*FCC_FE*FCC_C + 67.0*FCC_FE*FCC_VA + 1.0e-9)^25/T^25,0))))*log((-2.1*FCC_FE*FCC_C - 2.1*FCC_FE*FCC_VA)*if(2.1*FCC_FE*FCC_C + 2.1*FCC_FE*FCC_VA >= 0,-0.333333333333333,1.0) + 1)/(1.0*FCC_FE + 1.0*FCC_C) - 34671.0*FCC_FE*FCC_C*FCC_VA/(1.0*FCC_FE + 1.0*FCC_C) + (FCC_FE*FCC_C*(81739319.0569838*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 157733333.327276*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) - 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*T^(-1.0))) - 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*T^(-1.0))) + 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*2.71828182845905^(0.5*(0.5*2.71828182845905^(0.5*log(4.000000002e-11*101325 + 1.0)) + 0.5)^(-2.0) - 2.0*(0.5*2.71828182845905^(0.5*log(4.000000002e-11*101325 + 1.0)) + 0.5)^(-1.0) + 1.5)*T^(-1.0))) + 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*2.71828182845905^(0.5875*(0.5*2.71828182845905^(0.5*log(4.84848485090909e-11*101325 + 1.0)) + 0.5)^(-2.0) - 2.35*(0.5*2.71828182845905^(0.5*log(4.84848485090909e-11*101325 + 1.0)) + 0.5)^(-1.0) + 1.7625)*T^(-1.0))) - 15.877*T + 101580260.982824*(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315)^(-2.0) - 765093880.550632*(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315)^(-1.0) + 473161111.115*(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822)^(-2.0) - 3209266666.669*(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822)^(-1.0) - 0.25*(-3.0*2.71828182845905^(0.333333333333333 - 0.333333333333333*2.71828182845905^(0.5*log(4.0e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-0.333333333333333*2.71828182845905^(0.5*log(4.0e-10*101325 + 1.0)) + 0.5*log(4.0e-10*101325 + 1.0) + 0.333333333333333) + 4.0)*if(T < 298.15,-0.00279025825463012*T^2.0,-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*T^(-1.0))) + 4.18320906199373*T + if(T < 1811.0,-1.15*T*log(T) + 8.282*T + 0.00064*T^2.0 + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1462.4,4.9251e+30*T^(-9.0) + 0.940009*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1713.815) - 2975.31547683951) - 0.5*(-1.0*2.71828182845905^(1.0 - 1.0*2.71828182845905^(0.5*log(1.1e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-1.0*2.71828182845905^(0.5*log(1.1e-10*101325 + 1.0)) + 0.5*log(1.1e-10*101325 + 1.0) + 1.0) + 2.0)*(-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*T^(-1.0))) + 0.09009204*T + if(T < 103.0,-0.09009204*T - 2.75e-5*T^3.0 - 1049.14084,if(T < 350.0,1.76583*T*log(T) - 7.39898691*T - 0.01706952*T^2.0 - 988.25091,12000000000.0*T^(-3.0) - 264300000.0*T^(-2.0) + 2562600.0*T^(-1.0) - 24.3*T*log(T) + 170.73*T - 0.0004723*T^2.0 - 17368.441)) + 1049.14084) + if(T < 1811.0,-1.15*T*log(T) + 8.282*T + 0.00064*T^2.0 + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1462.4,4.9251e+30*T^(-9.0) + 0.940009*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1713.815) + if(T < 103.0,-0.09009204*T - 2.75e-5*T^3.0 - 1049.14084,if(T < 350.0,1.76583*T*log(T) - 7.39898691*T - 0.01706952*T^2.0 - 988.25091,12000000000.0*T^(-3.0) - 264300000.0*T^(-2.0) + 2562600.0*T^(-1.0) - 24.3*T*log(T) + 170.73*T - 0.0004723*T^2.0 - 17368.441)) - 960436998.988132*log(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315) - 3627866666.669*log(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822) + 3160223729.73755) + FCC_FE*FCC_VA*(81739319.0569838*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) - 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*T^(-1.0))) + 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*2.71828182845905^(0.5875*(0.5*2.71828182845905^(0.5*log(4.84848485090909e-11*101325 + 1.0)) + 0.5)^(-2.0) - 2.35*(0.5*2.71828182845905^(0.5*log(4.84848485090909e-11*101325 + 1.0)) + 0.5)^(-1.0) + 1.7625)*T^(-1.0))) + 101580260.982824*(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315)^(-2.0) - 765093880.550632*(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315)^(-1.0) - 0.25*(-3.0*2.71828182845905^(0.333333333333333 - 0.333333333333333*2.71828182845905^(0.5*log(4.0e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-0.333333333333333*2.71828182845905^(0.5*log(4.0e-10*101325 + 1.0)) + 0.5*log(4.0e-10*101325 + 1.0) + 0.333333333333333) + 4.0)*if(T < 298.15,-0.00279025825463012*T^2.0,-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-250.0*T^(-1.0))) + 4.18320906199373*T + if(T < 1811.0,-1.15*T*log(T) + 8.282*T + 0.00064*T^2.0 + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1462.4,4.9251e+30*T^(-9.0) + 0.940009*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1713.815) - 2975.31547683951) + if(T < 1811.0,-1.15*T*log(T) + 8.282*T + 0.00064*T^2.0 + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1462.4,4.9251e+30*T^(-9.0) + 0.940009*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1713.815) - 960436998.988132*log(0.203389830494685*2.71828182845905^(0.25*log(3.97306396340067e-11*101325 + 1.0)) + 0.796610169505315) + 581774300.510824))/(1.0*FCC_FE + 1.0*FCC_C)'
    derivative_order=2
    coupled_variables = 'Yc_fcc T'
    constant_names = 'FCC_FE scale'
    constant_expressions = '1 1'
    # constant_names = 'eV_J d'
    # constant_expressions = '6.24150934e+18 1e-21'
    # constant_expressions = '6.24150934e+18 1e-21' # '6.24150934e+18 1e-23'
    # constant_expressions = '1 1e-27' # '6.24150934e+18 1e-23'
    enable_jit = true
  [../]
# ============================================================
# Free energy material for bcc
# ============================================================
  [./free_energy_bcc]
    type = DerivativeParsedMaterial
    block = 0
    property_name = Fbcc
    expression='BCC_C:=Yc_bcc; BCC_VA:=1-Yc_bcc; scale*8.3145*T*(1.0*if(BCC_FE > 1.0e-15,BCC_FE*log(BCC_FE),0) + 3.0*if(BCC_C > 1.0e-15,BCC_C*log(BCC_C),0) + 3.0*if(BCC_VA > 1.0e-15,BCC_VA*log(BCC_VA),0))/(1.0*BCC_FE + 3.0*BCC_C) + 8.3145*T*if(T < -1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA,-0.00153008346379333*T^15/(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^15 - 0.00680037095019259*T^9/(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^9 - 0.153008346379333*T^3/(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^3 + 1 - 0.905299382744389*(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)/T,if(T < 1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA,-0.00153008346379333*T^15/(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^15 - 0.00680037095019259*T^9/(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^9 - 0.153008346379333*T^3/(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^3 + 1 - 0.905299382744389*(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)/T,if(T > 1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA & 1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA > 0,-0.0641731208021339*(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^5/T^5 - 0.00203724193022647*(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^15/T^15 - 0.000427820805347559*(1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA + 1.0e-9)^25/T^25,if(T > -1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA & 1043.0*BCC_FE*BCC_C + 1043.0*BCC_FE*BCC_VA < 0,-0.0641731208021339*(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^5/T^5 - 0.00203724193022647*(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^15/T^15 - 0.000427820805347559*(-1043.0*BCC_FE*BCC_C - 1043.0*BCC_FE*BCC_VA + 1.0e-9)^25/T^25,0))))*log((2.22*BCC_FE*BCC_C + 2.22*BCC_FE*BCC_VA)*if(2.22*BCC_FE*BCC_C + 2.22*BCC_FE*BCC_VA <= 0,-1.0,1.0) + 1)/(1.0*BCC_FE + 3.0*BCC_C) + BCC_FE*BCC_C*BCC_VA*if(T < 3900.0,-190.0*T,0.0)/(1.0*BCC_FE + 3.0*BCC_C) + (BCC_FE*BCC_C*(123749119.992992*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 500239999.980789*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) - 74.8305*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*T^(-1.0))) - 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*T^(-1.0))) + 74.8305*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*2.71828182845905^(0.5*(0.5*2.71828182845905^(0.5*log(4.000000002e-11*101325 + 1.0)) + 0.5)^(-2.0) - 2.0*(0.5*2.71828182845905^(0.5*log(4.000000002e-11*101325 + 1.0)) + 0.5)^(-1.0) + 1.5)*T^(-1.0))) + 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*2.71828182845905^(0.332142857142857*(0.5*2.71828182845905^(0.5*log(5.49019608117647e-11*101325 + 1.0)) + 0.5)^(-2.0) - 1.32857142857143*(0.5*2.71828182845905^(0.5*log(5.49019608117647e-11*101325 + 1.0)) + 0.5)^(-1.0) + 0.996428571428571)*T^(-1.0))) + 75.667*T + 207216388.673325*(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514)^(-2.0) - 1498770592.00399*(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514)^(-1.0) + 1500596666.679*(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822)^(-2.0) - 10177960000.0074*(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822)^(-1.0) - 0.25*(-3.0*2.71828182845905^(0.333333333333333 - 0.333333333333333*2.71828182845905^(0.5*log(5.64705882352941e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-0.333333333333333*2.71828182845905^(0.5*log(5.64705882352941e-10*101325 + 1.0)) + 0.5*log(5.64705882352941e-10*101325 + 1.0) + 0.333333333333333) + 4.0)*(-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*T^(-1.0))) + 9.929623*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 4701.3908) - 1.5*(-1.0*2.71828182845905^(1.0 - 1.0*2.71828182845905^(0.5*log(1.1e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-1.0*2.71828182845905^(0.5*log(1.1e-10*101325 + 1.0)) + 0.5*log(1.1e-10*101325 + 1.0) + 1.0) + 2.0)*(-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-1350.0*T^(-1.0))) + 0.09009204*T + if(T < 103.0,-0.09009204*T - 2.75e-5*T^3.0 - 1049.14084,if(T < 350.0,1.76583*T*log(T) - 7.39898691*T - 0.01706952*T^2.0 - 988.25091,12000000000.0*T^(-3.0) - 264300000.0*T^(-2.0) + 2562600.0*T^(-1.0) - 24.3*T*log(T) + 170.73*T - 0.0004723*T^2.0 - 17368.441)) + 1049.14084) + 3.0*if(T < 103.0,-0.09009204*T - 2.75e-5*T^3.0 - 1049.14084,if(T < 350.0,1.76583*T*log(T) - 7.39898691*T - 0.01706952*T^2.0 - 988.25091,12000000000.0*T^(-3.0) - 264300000.0*T^(-2.0) + 2562600.0*T^(-1.0) - 24.3*T*log(T) + 170.73*T - 0.0004723*T^2.0 - 17368.441)) + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1806737152.00399*log(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514) - 11505520000.0074*log(0.115384615380178*2.71828182845905^(0.25*log(5.77777776355555e-11*101325 + 1.0)) + 0.884615384619822) + 9345250466.68529) + BCC_FE*BCC_VA*(123749119.992992*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) - 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*T^(-1.0))) + 24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*2.71828182845905^(0.332142857142857*(0.5*2.71828182845905^(0.5*log(5.49019608117647e-11*101325 + 1.0)) + 0.5)^(-2.0) - 1.32857142857143*(0.5*2.71828182845905^(0.5*log(5.49019608117647e-11*101325 + 1.0)) + 0.5)^(-1.0) + 0.996428571428571)*T^(-1.0))) + 207216388.673325*(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514)^(-2.0) - 1498770592.00399*(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514)^(-1.0) - 0.25*(-3.0*2.71828182845905^(0.333333333333333 - 0.333333333333333*2.71828182845905^(0.5*log(5.64705882352941e-10*101325 + 1.0))) - 1.0*2.71828182845905^(-0.333333333333333*2.71828182845905^(0.5*log(5.64705882352941e-10*101325 + 1.0)) + 0.5*log(5.64705882352941e-10*101325 + 1.0) + 0.333333333333333) + 4.0)*(-24.9435*T*log(1.0 - 1.0*2.71828182845905^(-300.0*T^(-1.0))) + 9.929623*T + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 4701.3908) + if(T < 50.0,-9.929623*T - 1.975e-6*T^4.0 + 4701.3908,if(T < 170.0,8.75517*T*log(T) - 41.511332*T - 0.13106*T^2.0 + 9.2615e-5*T^3.0 + 4871.6848,if(T < 1811.0,77359.0*T^(-1.0) - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.8927e-8*T^3.0 + 1225.7,2.29603e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581))) - 1806737152.00399*log(0.170454545444861*2.71828182845905^(0.25*log(4.60130717830065e-11*101325 + 1.0)) + 0.82954545455514) + 1167805083.33768))/(1.0*BCC_FE + 3.0*BCC_C)'
    coupled_variables = 'Yc_bcc T'
    derivative_order=2
    constant_names = 'BCC_FE scale'
    constant_expressions = '1 1'
    
    # constant_expressions = '6.24150934e+18 1e-21'
    # constant_expressions = '6.24150934e+18 1e-21' # '6.24150934e+18 1e-23'
    # constant_expressions = '1 1e-27' # '6.24150934e+18 1e-23'
    enable_jit = true
  [../]
        
# ============================================================
# Free energy material for two-phase
# ============================================================
  [./free_energy_chem]
    type = DerivativeTwoPhaseMaterial
    block = 0
    property_name = Ftotal_chem    # Name of the global free energy function (use this in the Parsed Function Kernels)
    fa_name = Ffcc  # f_name of the phase A free energy function
    fb_name = Fbcc  # f_name of the phase B free energy function
    coupled_variables = 'Yc_fcc Yc_bcc T'
    eta = eta     # order parameter that switches between A and B phase
    derivative_order = 2
    outputs = exodus
  [../]

# ============================================================
# Free energy material for multiphase
# ============================================================
  [./free_energy_chem_multiphase]
    type = DerivativeMultiPhaseMaterial
    property_name = Fchem_multiphase
    fi_names = 'Ffcc Ffcc Ffcc Ffcc Ffcc Fbcc Fbcc Fbcc Fbcc'
    hi_names = 'h0 h1 h2 h3 h4 h5 h6 h7 h8'
    etas     = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    coupled_variables = 'Yc_fcc Yc_bcc T'
    derivative_order = 2
  [../]

# ============================================================
# Nucleation material for eta
# ============================================================
# Probability material
# ============================================================
  [./probability]
    type = ParsedMaterial
    property_name = P
    coupled_variables = 'T'
    # Nucleation active between Ms and Mf
    expression = 'if(T <= 623 & T >= 393.15, 4e-3, 0)' # Ms=350 C for martensite formation # mf = 120 C for martensite finish
    outputs = exodus
  [../]

  [./nucleation_eta]
        # step 4
        # The nucleation material is configured to insert nuclei into the free energy
        # tht force the concentration to go to 0.95, and holds this enforcement for 500
        # time units.
        type = DiscreteNucleation
        property_name = Fn
        op_names  = eta
        op_values = 1. # If a point on the mesh is eligible for nucleation, then suddenly change the value of the variable c at that point to op_values (e.g. 1).”
        penalty = 1
        penalty_mode = MATCH
        map = map
        outputs = exodus
  [../]
# ============================================================
# Nucleation material for gr8
# ============================================================
  [./nucleation_gr8]
    # step 4
    # The nucleation material is configured to insert nuclei into the free energy
    # tht force the concentration to go to 0.95, and holds this enforcement for 500
    # time units.
    type = DiscreteNucleation
    property_name = Fn_gr
    op_names  = 'gr8' # gr7 gr6 gr5
    op_values = '1.' # 1. 1. 1. 1. # If a point on the mesh is eligible for nucleation, then suddenly change the value of the variable c at that point to op_values (e.g. 1).”
    penalty = 100
    penalty_mode = MATCH
    map = map
    outputs = exodus
  [../]

# ============================================================
# Multibarrier free energy material
# ============================================================
  [./free_energy_multibarrier]
    type = DerivativeParsedMaterial
    property_name = Fmultibarrier
    
    # coupled_variables = 'gr0 gr1 gr2'
    # material_property_names = 'h0 h1 h2 h3 h4 h5 h6 h7' 
    # expression = '-10*(gr2^2*(1-gr2)^2+gr0^2*(1-gr0)^2+gr1^2*(1-gr1)^2)'
    
    coupled_variables = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    # previously = -1 
    expression = '1*(gr2^2*(1-gr2)^2+gr0^2*(1-gr0)^2+gr1^2*(1-gr1)^2+gr3^2*(1-gr3)^2+gr4^2*(1-gr4)^2+gr5^2*(1-gr5)^2+gr6^2*(1-gr6)^2+gr7^2*(1-gr7)^2+gr8^2*(1-gr8)^2)'
    derivative_order = 2
    enable_jit = true
  [../]
  #[./barrier]
  #  type = MultiBarrierFunctionMaterial
  #  etas = 'eta1 eta2 eta3'
  #[../]

  # =========================
  # Phase–grain partition penalty: enforce
  #   Sfcc ≈ (1 - h), Sbcc ≈ h
  # =========================
  [./phase_grain_penalty]
    type = DerivativeParsedMaterial
    property_name = F_phase_grain
    coupled_variables = 'gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    material_property_names = 'h'
    constant_names = 'Wpg'
    constant_expressions = '1e3'
    expression = '
      Sfcc := gr0^2+gr1+gr2+gr3+gr4;
      Sbcc := gr5+gr6+gr7+gr8;
      Wpg*(Sbcc - h)^2 + Wpg*(Sfcc - (1 - h))^2
    '
    derivative_order = 2
    enable_jit = true
  [../]


# ============================================================
# Total free energy materials
# ============================================================
  [./total_free_energy_eta]
    type = DerivativeSumMaterial
    derivative_order = 2
    coupled_variables = 'eta Yc_fcc Yc_bcc T'
    sum_materials = 'Ftotal_chem Fn'
    property_name = Ftotal_eta
  [../]
        
  [./total_free_energy_gr]
    type = DerivativeSumMaterial
    derivative_order = 2
    coupled_variables = 'eta T gr0 gr1 gr2 gr3 gr4 gr5 gr6 gr7 gr8'
    sum_materials = 'Ftotal_loc F_phase_grain'
    property_name = Ftotal_gr
  [../]
[]

# =========================
# Postprocessors (scalar monitoring)
# =========================
[Postprocessors]
  [dt]
    type = TimestepSize
  []

  [./eta_average]
    type = ElementAverageValue
    variable = eta
  []

  # [./c_average]
  #   type = ElementAverageValue
  #   variable = c
  # []

  [./temperature]
    type = ElementExtremeValue
    variable = T
  []

  [./Ffcc_average]
    type = ElementAverageValue
    variable = Ffcc_total
  []

  [./Fbcc_average]
    type = ElementAverageValue
    variable = Fbcc_total
  []

  [./Ftotal_eta_average]
    type = ElementAverageValue
    variable = Ftotal_eta_aux
  []
  [./Ftotal_gr_average]
    type = ElementAverageValue
    variable = Ftotal_gr_aux
  []
  [./Ftotal_loc_average]
    type = ElementAverageValue
    variable = Ftotal_loc_aux
  []

  [./F_phase_grain_average]
    type = ElementAverageValue
    variable = F_phase_grain_aux
  []

  [./Sfcc_avg]
    type = ElementAverageValue
    variable = Sfcc
  []

  [./Sbcc_avg]
    type = ElementAverageValue
    variable = Sbcc
  []
  [./dG_avg]
    type = ParsedPostprocessor
    expression = 'Fbcc_average - Ffcc_average'
    pp_names = 'Fbcc_average Ffcc_average'
  []
  [./average_gr8]
    type = ElementAverageValue
    variable = gr8
  [../]
    [./average_gr0]
      type = ElementAverageValue
      variable = gr0
    [../]

  [./average_gr1]
    type = ElementAverageValue
    variable = gr1
  []

  [./average_gr2]
    type = ElementAverageValue
    variable = gr2
  []

  [./average_gr3]
    type = ElementAverageValue
    variable = gr3
  []

  [./average_gr4]
    type = ElementAverageValue
    variable = gr4
  []

  [./average_gr5]
    type = ElementAverageValue
    variable = gr5
  []

  [./average_gr6]
    type = ElementAverageValue
    variable = gr6
  []

  [./average_gr7]
    type = ElementAverageValue
    variable = gr7
  []
  [./average_Fn_gr]
    type = ElementAverageValue
    variable = Fn_gr
  []
  [./average_Fn]
    type = ElementAverageValue
    variable = Fn
  []
  [./average_L_fcc]
    type = ElementAverageValue
    variable = L_fcc_aux
  []
  [./average_L_bcc]
    type = ElementAverageValue
    variable = L_bcc_aux
  []
[]

# =========================
# Executioner
# =========================
[Executioner]
  type = Transient
  scheme = bdf2

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'

  l_max_its = 50
  l_tol = 1e-4
  nl_max_its = 10

  start_time = 20
  # dtmin = 1e-3
  dtmax = 0.5
  end_time = 45 # 33.85

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.1
    optimal_iterations = 10
    iteration_window = 2
    growth_factor = 1.5
    cutback_factor = 0.5
  [../]

  # [TimeStepper]
  #   type = IterationAdaptiveDT
  #   dt = 0.01
  #   optimal_iterations = 6
  # []
[]

# =========================
# Outputs
# =========================
[Outputs]
  exodus = true
  csv = true
[]

`

[GiudGiud]

Hello

should you be using the reserve_op parameter of the GrainTracker?
@permcody do you know about the grain remapping? They want to use specific variables for new grains and remapping grains

[permcody]

The current implementation in the GrainTracker doesn't support choosing a subset of order parameter variables for remapping. As you've pointed out, the reserved order parameter is intended for nucleation so it's the right feature to use for that portion of the problem. You can have multiple sets of order parameters in a simulation but that may cause different issues depending on how you are handling your physics kernels. You could attempt to modify the GrainTracker to limit which order parameters are remapped to, but you'd have to be careful to avoid unintended side effects like having grains merge on order parameters that are included or excluded from your list.