Linearized Interface Grain Growth Model

Introduction

Figure 1: Example of order parameter values (, ) and the linearized interface transformed variables (, ),

In a phase field grain growth model, like that developed by Moelans et al. (2008) and implemented in MOOSE, the order parameters that represent the grains tend to follow a profile across the interface. A fine mesh is required to resolve this highly nonlinear profile, such that at least three elements are typically required across the grain boundaries. This fine mesh required at the grain boundaries makes the computational cost of grain growth simulations high.

The cost of grain growth simulations can be reduced by using a change of variables to linearize the interface profile, sometimes referred to as nonlinear preconditioning and first suggested by Glasner (2001). This change in variables is defined by

where is the original order parameter and is the new variable. While follows a profile across the interface, the new variable is linear across the interface, as illustrated in Figure 1. Due to the linear profile across the interface, a single linear element can resolve the profile without error. This significantly decreases the required resolution, and thus the computational cost.

MOOSE Implementation

Linearized interface has been implemented in the multiphase field grain growth model in MOOSE. The basic model is identical to the standard phase field model in MOOSE, but it solves for the transformed linearized interface variables rather than the standard order parameters. The order parameters are created as AuxVariables, allowing them to be visualized. The grain boundary can be visualized using an AuxVariable equal to , calculated using BndsCalcAux in the usual manner.

The MOOSE objects used in the linearized interface grain growth model are listed, below:

Materials

• LinearizedInterfaceFunction - Defines the order parameters from the linear interface variable values at the quadrature points. It uses the ExpressionBuilder so it also calculates the analytical derivatives.

Kernels

• ChangedVariableTimeDerivative - Version of the TimeDerivative kernel after a change of variable.

• ACGrGrPolyLinearizedInterface - Linearized interface version of the ACGrGrPoly kernel.

• ACInterfaceChangedVariable - Version of the ACInterface kernel after a change of variable.

AuxKernels

• LinearizedInterfaceAux - Defines AuxVariable values of the order parameters from the nonlinear variable values at the nodes for the linearized interface grain growth model.

Bounded Solve

When solving the linearized interface version of the grain growth model, Gong et al. (2018) and Chadwick and Voorhees (2021) found that the solve is unstable unless bounds are placed on the possible values for the transformed variables. Therefore, upper and lower bounds must be defined in the [Bounds] block for each variable using ConstantBounds. A typical bounds range is . For example:

[Bounds]
  [phi0_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi0
    bound_type = upper
    bound_value = 5.0
  []
  [phi0_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi0
    bound_type = lower
    bound_value = -5.0
  []
  [phi1_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi1
    bound_type = upper
    bound_value = 5.0
  []
  [phi1_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi1
    bound_type = lower
    bound_value = -5.0
  []
[]

Note that in order for these bounds to have an effect, the user has to specify the PETSc options -snes_type vinewtonssls or -snes_type vinewtonrsls. A warning will be generated if neither options are specified. The PETSc manual pages for the vinewtonssls algorithm can be found here while the manual page for vinewtonrsls can be found here.

These solve options are demonstrated here:

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON

  petsc_options_iname = '-pc_type -ksp_type -snes_type'
  petsc_options_value = 'bjacobi gmres vinewtonrsls'

  dt = 0.1
  end_time = 0.6
[]

Simplified MOOSE syntax

As with the standard grain growth model in MOOSE, an Action has been created to simplify the input file syntax needed to define a nonlinear preconditioning grain growth model. The GrainGrowthLinearizedInterfaceAction allows the user to specify the number of order parameters to be used in the model and it automatically generates the

• Linearized interface variables

• AuxVariables and corresponding Auxkernels defining the order parameters

• All of the kernels

• Material properties defining the order parameters

• Upper and lower bounds required for the solve

• Auxvariable and kernel defining the grain boundary for visualization

The custom syntax for the linearized interface grain growth model is nearly identical to that for the standard grain growth model, making it simple to switch back and forth between models.

linearized interface Initial Conditions

Initial conditions have been created to transform the IC values for the standard grain growth model to the linearized interface model.

• SmoothCircleICLinearizedInterface - Takes the output from the SmoothCircleIC and performs the change in variable. Used to model a shrinking circular grain.

• PolycrystalColoringICAction - This custom syntax is used for all polycrystal ICs for both the standard grain growth model and with nolinear preconditioning. The parameter nonlinear_preconditioning = true uses the PolycrystalColoringICLinearizedInterface for each variable.

Usage with GrainTracker

GrainTracker, which allows variables to represent multiple grains and remap them to other variables to avoid coalescence, also works with the linearized interface grain growth model. Several changes are needed in its use:

1. The input parameter bound_value should be set to the value used in the bounded solve. This will switch the variable value when a grain is remapped to -bound_value, erasing the grain from the original variable. This is the same parameter used in GrainGrowthLinearizedInterfaceAction and in the linearized interface ICs, so that it can be set for all of the objects from the GlobalParams block.

2. The threshold parameter must be changed to be slightly above the value of the lower bound (-bound_value). It is recommended to use -bound_value + 2*bound_value/10.

The modified syntax is illustrated here:

[GlobalParams]
  # Parameters used by several kernels that are defined globally to simplify input file
  op_num = 8 # Number of order parameters used
  var_name_base = psi # Base name of grains
  bound_value = 5 # +/- bound value
[]

[UserObjects]
  [grain_tracker]
    type = GrainTracker
    threshold = -4
  []
[]

Example Input Files

An example input file for a polycrystal simulation using the linearized interface grain growth model is available:

[GlobalParams]
  bound_value = 5.0
  op_num = 8
  var_name_base = phi
[]

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmax = 1000
  ymax = 1000
  nx = 100
  ny = 100
  uniform_refine = 1
[]

[Modules]
  [PhaseField]
    [GrainGrowthLinearizedInterface]
      op_name_base = gr
      mobility = L
      kappa = kappa_op
    []
  []
[]

[ICs]
  [PolycrystalICs]
    [PolycrystalColoringIC]
      polycrystal_ic_uo = RandomVoronoi
      nonlinear_preconditioning = true
    []
  []
[]

[UserObjects]
  [RandomVoronoi]
    type = PolycrystalVoronoi
    grain_num = 60
    int_width = 10
    rand_seed = 103838
  []
  [grain_tracker]
    type = GrainTracker
    threshold = -4.0
    compute_halo_maps = true # Only necessary for displaying HALOS
  []
[]

[AuxVariables]
  [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
  [var_indices]
    order = CONSTANT
    family = MONOMIAL
  []
  [halos]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [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'
  []
  [halos]
    type = FeatureFloodCountAux
    variable = halos
    flood_counter = grain_tracker
    field_display = HALOS
    execute_on = 'initial timestep_end'
  []
[]

[Materials]
  [properties]
    type = GenericConstantMaterial
    prop_names = 'gbmob gbenergy gbwidth gamma_asymm'
    prop_values = '100 6 10 1.5'
  []
  [kappa_op]
    type = ParsedMaterial
    material_property_names = 'gbenergy gbwidth'
    property_name = kappa_op
    expression = '3/4*gbenergy*gbwidth'
  []
  [L]
    type = ParsedMaterial
    material_property_names = 'gbmob gbwidth'
    property_name = L
    expression = '4/3*gbmob/gbwidth'
  []
  [mu]
    type = ParsedMaterial
    material_property_names = 'gbenergy gbwidth'
    property_name = mu
    expression = '6*gbenergy/gbwidth'
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = 'initial TIMESTEP_END'
  []
[]

[BCs]
  [Periodic]
    [All]
      auto_direction = 'x y'
    []
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = PJFNK

  petsc_options_iname = '-pc_type -pc_hypre_type -snes_type'
  petsc_options_value = 'hypre    boomeramg      vinewtonrsls'

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

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.02
    optimal_iterations = 6
  []
  end_time = 30
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

The linearized interface tests also serve as good examples for how to use the capability.

The evolution of a shrinking circular grain is illustrated here:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmax = 50
  ymax = 50
  nx = 10
  ny = 10
[]

[Variables]
  [phi0]
  []
  [phi1]
  []
[]

[AuxVariables]
  [gr0_aux]
  []
  [gr1_aux]
  []
  [bounds_dummy]
  []
[]

[AuxKernels]
  [gr0]
    type = LinearizedInterfaceAux
    variable = gr0_aux
    nonlinear_variable = phi0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [gr1]
    type = LinearizedInterfaceAux
    variable = gr1_aux
    nonlinear_variable = phi1
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[ICs]
  [phi0_IC]
    type = SmoothCircleICLinearizedInterface
    variable = phi0
    invalue = 1.0
    outvalue = 0.0
    bound_value = 5.0
    radius = 30
    int_width = 10
    x1 = 0.0
    y1 = 0.0
    profile = TANH
  []
  [phi1_IC]
    type = SmoothCircleICLinearizedInterface
    variable = phi1
    invalue = 0.0
    outvalue = 1.0
    bound_value = 5.0
    radius = 30
    int_width = 10
    x1 = 0.0
    y1 = 0.0
    profile = TANH
  []
[]

[Kernels]
  #phi0 Kernels
  [phi0_dot]
    type = ChangedVariableTimeDerivative
    variable = phi0
    order_parameter = gr0
  []
  [phi0_ACInt]
    type = ACInterfaceChangedVariable
    variable = phi0
    kappa_name = kappa_op
    mob_name = L
    order_parameter = gr0
  []
  [gr0_AC]
    type = ACGrGrPolyLinearizedInterface
    variable = phi0
    mob_name = L
    this_op = gr0
    other_ops = gr1
    v = phi1
  []
  #phi1 Kernels
  [phi1_dot]
    type = ChangedVariableTimeDerivative
    variable = phi1
    order_parameter = gr1
  []
  [phi1_ACInt]
    type = ACInterfaceChangedVariable
    variable = phi1
    kappa_name = kappa_op
    mob_name = L
    order_parameter = gr1
  []
  [gr1_AC]
    type = ACGrGrPolyLinearizedInterface
    variable = phi1
    mob_name = L
    this_op = gr1
    other_ops = gr0
    v = phi0
  []
[]

[Materials]
  [gr0]
    type = LinearizedInterfaceFunction
    f_name = gr0
    phi = phi0
  []
  [gr1]
    type = LinearizedInterfaceFunction
    f_name = gr1
    phi = phi1
  []
  [GBEovlution]
    type = GBEvolution
    GBenergy = 0.97
    GBMobility = 0.6e-6
    T = 300
    wGB = 10
  []
[]

[Bounds]
  [phi0_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi0
    bound_type = upper
    bound_value = 5.0
  []
  [phi0_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi0
    bound_type = lower
    bound_value = -5.0
  []
  [phi1_upper_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi1
    bound_type = upper
    bound_value = 5.0
  []
  [phi1_lower_bound]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = phi1
    bound_type = lower
    bound_value = -5.0
  []
[]

[Postprocessors]
  [grain_area_mat]
    type = ElementIntegralMaterialProperty
    mat_prop = gr0
    execute_on = 'initial TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON

  petsc_options_iname = '-pc_type -ksp_type -snes_type'
  petsc_options_value = 'bjacobi gmres vinewtonrsls'

  dt = 0.1
  end_time = 0.6
[]

[Outputs]
  exodus = true
[]

The same problem but using GrainGrowthLinearizedInterfaceAction is illustrated here:

[GlobalParams]
  bound_value = 5.0
[]

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmax = 50
  ymax = 50
  nx = 10
  ny = 10
[]

[Modules]
  [PhaseField]
    [GrainGrowthLinearizedInterface]
      op_num = 2
      var_name_base = phi
      op_name_base = gr
      mobility = L
      kappa = kappa_op
    []
  []
[]

[ICs]
  [phi0_IC]
    type = SmoothCircleICLinearizedInterface
    variable = phi0
    invalue = 1.0
    outvalue = 0.0
    radius = 30
    int_width = 10
    x1 = 0.0
    y1 = 0.0
    profile = TANH
  []
  [phi1_IC]
    type = SmoothCircleICLinearizedInterface
    variable = phi1
    invalue = 0.0
    outvalue = 1.0
    radius = 30
    int_width = 10
    x1 = 0.0
    y1 = 0.0
    profile = TANH
  []
[]

[Materials]
  [GBEovlution]
    type = GBEvolution
    GBenergy = 0.97
    GBMobility = 0.6e-6
    T = 300
    wGB = 10
  []
[]

[Postprocessors]
  [grain_area_mat]
    type = ElementIntegralMaterialProperty
    mat_prop = gr0
    execute_on = 'initial TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON

  petsc_options_iname = '-pc_type -ksp_type -snes_type'
  petsc_options_value = 'bjacobi gmres vinewtonrsls'

  dt = 0.1
  end_time = 0.6
[]

[Outputs]
  exodus = true
[]

The evolution of a five-grain polycrystal is illustrated here:

[GlobalParams]
  bound_value = 5.0
  op_num = 5
  var_name_base = phi
[]

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmax = 100
  ymax = 100
  nx = 20
  ny = 20
[]

[Modules]
  [PhaseField]
    [GrainGrowthLinearizedInterface]
      op_name_base = gr
      mobility = L
      kappa = kappa_op
    []
  []
[]

[ICs]
  [PolycrystalICs]
    [PolycrystalColoringIC]
      polycrystal_ic_uo = RandomVoronoi
      linearized_interface = true
    []
  []
[]

[UserObjects]
  [RandomVoronoi]
    type = PolycrystalVoronoi
    grain_num = 5
    int_width = 10
    rand_seed = 103838
  []
[]

[Materials]
  [GBEovlution]
    type = GBEvolution
    GBenergy = 0.97
    GBMobility = 0.6e-6
    T = 300
    wGB = 10
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON

  petsc_options_iname = '-pc_type -ksp_type -snes_type'
  petsc_options_value = 'bjacobi gmres vinewtonrsls'

  dt = 0.05
  num_steps = 1
[]

[Outputs]
  exodus = true
[]

References