# 18.1. Reaction Diffusion Solver

The phase-field model solves the phase-field evolution using a reaction-diffusion solver implemented in Sierra/SM. The reaction-diffusion command block is as follows:

    BEGIN REACTION DIFFUSION rxndiffname
      BLOCK = <string list>block_names
      REMOVE BLOCK = <string list>remove_block_names
      INCLUDE ALL BLOCKS

      ASSEMBLY = <string list>assembly_names
      REMOVE ASSEMBLY = <string list>remove_assembly_names

      ACTIVE PERIODS = <string list> active_period_names
      INACTIVE PERIODS = <string list> inactive_period_names

      PHASE FIELD BOUND CONSTRAINTS

      BEGIN PRESCRIBED FIELD
        SURFACE = <string>surf_name
        NODE SET|NODESET = <string>node_set_name
        FUNCTION = <string>func_name
      END
      BEGIN PRESCRIBED FLUX
        SURFACE = <string>surf_name
        NODE SET|NODESET = <string>node_set_name
        FUNCTION = <string>func_name
      END

      GRADIENT CONFIGURATION = MODEL|CURRENT(MODEL)
      USE FINITE ELEMENT MODEL = <string>model_name
          [Model Coordinates Are <string> nodal_variable_name]

      SOLVE AT INITIALIZATION = OFF|ON(OFF)
      SOLVE STEP INCREMENT = <integer>solv_step_incr(1)
      SOLVE TIME INCREMENT = <real>solv_time_incr

      EQUATION SYSTEM = LINEAR|NONLINEAR(LINEAR)
      USE LINEAR SOLVER = FALSE|TRUE(TRUE)

      BEGIN SOLVER OPTIONS
        ACCEPTABLE RESIDUAL = <real>accept_resid(1.0e-14)
        ACCEPTABLE RELATIVE RESIDUAL = <real>accept_rel_resid(1.0e-12)
        TARGET RESIDUAL = <real>target_resid(1.0e-15)
        TARGET RELATIVE RESIDUAL = <real>target_rel_resid(1.0e-13)
        MAXIMUM ITERATIONS = <integer>max_iter(1000)
        LINE SEARCH BACKTRACK [<real>min_step(1.0e-8) <real>reduction(0.1)
            <real>init_step(1.0)
            BISECTION|CUTBACK|SECANT|TRUST_REGION|CRITICAL_POINT(BISECTION)
            <real>dispinc(0.5) <real>rotinc(0.0) ]
        LINEAR SOLVER = <string>linear_solver_name(DEFAULT_FETI_SOLVER)
        PRECONDITIONER = PROBE|IDENTITY(PROBE)
      END

    END

The commands regarding `BLOCK`, `ASSEMBLY`, and `ACTIVE PERIODS` are common to Sierra/SM, and allow the user to specify reaction-diffusion solution only on the blocks or time periods of interest.

The `SOLVE AT INITIALIZATION` command provides a way to solve the phase field before any mechanical loading is applied. This is particularly useful when applying a phase boundary condition (e.g. prescribed field, prescribed flux), as the phase field will not be in equilibrium prior to the mechanics solve in the first time-step.

Explicit only

The commands `SOLVE STEP INCREMENT` and `SOLVE TIME INCREMENT` offer a way of decreasing the frequency of reaction diffusion solves, particularly for the use-case involving implicit reaction diffusion solves with explicit time incrementation of the solid mechanics solution. Individually, the two commands define the period (in steps, or in time) between successive reaction-diffusion solves. It is an error to specify both commands simultaneously.

The command `PHASE FIELD BOUND CONSTRAINTS` offers a way of ensuring satisfaction of the phase field bound constraints. Specifically, this command restricts the reaction-diffusion solution to the interval \(\phi \in [0,1]\) and also ensures that the field evolves monotonically \(\Delta \phi \leq 0\). Implementation is done using a bound-constrained conjugate gradient method [[1]] with modifications around the active set selection.

The commands `BEGIN PRESCRIBED FIELD` and `BEGIN PRESCRIBED FLUX` offer a means to specify Dirichlet or Neumann boundary conditions on the phase field. The prescribed field is commonly used to instantiate the model with an initial damage field, such as a phase field pre-crack. Often this provides for better numerical performance than a meshed sharp crack with no initial phase field. The prescribed flux has not been tested, as it lacks a clear physical motivation. It is recommended to use the `SOLVE AT INITIALIZATION` command with these options to ensure that the phase field is in equilibrium before performing any mechanical solution.

The command `GRADIENT CONFIGURATION` provides the option to compute the phase-field gradient in model or current coordinates. Model coordinates are recommended.

The command `USE FINITE ELEMENT MODEL` provides the user the option to specify the finite element model for the phase field. By default, the same finite element model as the Adagio region is used. Using a different finite element mesh for the reaction-diffusion solve has not been tested and is not recommended.

The commands `EQUATION SYSTEM` and `USE LINEAR SYSTEM` are used to determine whether to solve the reaction-diffusion equation as a linear or nonlinear partial differential equation. For the current Phase Field FeFp implementation, the phase-field evolution equation is a linear partial differential equation when the degradation function parameter \(\gamma = 0\) (corresponding to \(\psi_c = \frac{3 G_c}{16 \ell}\)), so the linear reaction-diffusion solve may be used: `LINEAR`, `TRUE`; this uses the default FETI solver to solve the linear system. For other cases (\(\gamma \ne 0\)), it is the phase-field evolution equation is a nonlinear partial differential equation, so it is recommended to use the nonlinear reaction-diffusion solver: `NONLINEAR`, `FALSE`; this uses a preconditioned conjugate gradient method to iteratively solve the system.

The `BEGIN SOLVER OPTIONS` block offers options to customize the preconditioned conjugate gradient solve (`EQUATION SYSTEM = NONLINEAR`, `USE LINEAR SYSTEM = FALSE`), with commands that reflect a limited set of the commands in the Sierra/SM `BEGIN CG` block. Two nodal preconditioners are available: identity and probe. The probe preconditioner defines its entries as the quotient of the phase system residual divided by the nodal phase stiffness; the phase stiffness is computed as the derivative of the phase force-internal at that node with respect to phase, estimated using a forward finite difference (\(\Delta \phi = -1.0\cdot10^{-10}\)).

[1]

Edwin A.H. Vollebregt. The bound-constrained conjugate gradient method for non-negative matrices. Journal of Optimization Theory and Applications, 162(3):931–953, 2014.