GitHub

DifferentialEquations.jl Integration

Overview

PhaseFields.jl integrates with DifferentialEquations.jl using the Method of Lines approach:

  1. Discretize space with finite differences
  2. Convert PDE to system of ODEs
  3. Solve with adaptive time stepping

Grid and Boundary Conditions

Uniform Grid

using PhaseFields

# Create 1D uniform grid
grid = UniformGrid1D(N=101, L=1.0)  # 101 points, length 1.0

# Access properties
grid.N    # Number of points
grid.L    # Domain length
grid.dx   # Grid spacing
grid.x    # Coordinate array

Boundary Conditions

# Neumann (zero-flux)
bc = NeumannBC()

# Dirichlet (fixed values)
bc = DirichletBC(left=0.0, right=1.0)

# Periodic
bc = PeriodicBC()

Laplacian Operator

# In-place computation
∇²φ = zeros(grid.N)
laplacian_1d!(∇²φ, φ, grid.dx, bc)

# Allocating version
∇²φ = laplacian_1d(φ, grid.dx, bc)

# Sparse matrix form (for implicit solvers)
L = laplacian_matrix_1d(grid.N, grid.dx, bc)

Callbacks

Interface Tracking

using PhaseFields
using OrdinaryDiffEq

# Track interface position during simulation
cb, saved = create_interface_saving_callback(grid; saveat=0.1)

sol = solve(prob, Tsit5(); callback=cb)

# Access results
times = saved.t
positions = saved.saveval

Solid Fraction Tracking

cb, saved = create_solid_fraction_callback(; threshold=0.5, saveat=0.1)

sol = solve(prob, Tsit5(); callback=cb)
fractions = saved.saveval

Steady State Termination

# Stop when solution reaches steady state
cb = create_steady_state_callback(abstol=1e-6, reltol=1e-4, min_t=1.0)

sol = solve(prob, Tsit5(); callback=cb)
# sol.t[end] is when steady state was reached

Combined Callbacks

result = create_phase_field_callbacks(grid;
    track_interface=true,
    track_solid_fraction=true,
    terminate_steady_state=true,
    saveat=0.1)

sol = solve(prob, Tsit5(); callback=result.callback)

# Access tracked data
result.interface_data.t
result.interface_data.saveval
result.solid_fraction_data.saveval

Solver Selection

Explicit Solvers (Non-stiff)

# Tsit5: Good general-purpose choice
sol = solve(prob, Tsit5())

# Vern7: Higher accuracy
sol = solve(prob, Vern7())

Implicit Solvers (Stiff)

For thermal coupling or fine grids:

# QNDF: Quasi-constant step BDF (recommended)
sol = solve(prob, QNDF(autodiff=false))

# Rodas5: Rosenbrock method
sol = solve(prob, Rodas5(autodiff=false))

Note: Use autodiff=false with pre-allocated workspaces.

API Reference

Grid and Boundary Conditions

PhaseFields.UniformGrid1D — Type
UniformGrid1D

1D uniform grid for finite difference discretization.

Fields

  • N::Int: Number of grid points
  • L::Float64: Domain length
  • dx::Float64: Grid spacing (computed)
  • x::Vector{Float64}: Grid point coordinates

Example

grid = UniformGrid1D(N=100, L=1.0)
grid isa AbstractDomain{1}  # true
source
PhaseFields.DirichletBC — Type
DirichletBC

Fixed value (Dirichlet) boundary condition for 1D.

Fields

  • left::Float64: Value at left boundary
  • right::Float64: Value at right boundary
source
PhaseFields.laplacian_1d! — Function
laplacian_1d!(∇²u, u, dx, bc::BoundaryCondition)

Compute 1D Laplacian in-place using second-order central differences.

Arguments

  • ∇²u: Output array for Laplacian values (modified in-place)
  • u: Input field values
  • dx: Grid spacing
  • bc: Boundary condition
source
PhaseFields.laplacian_matrix_1d — Function
laplacian_matrix_1d(N, dx, bc::BoundaryCondition)

Create sparse Laplacian matrix for 1D finite differences.

Returns sparse matrix L such that L * u ≈ ∇²u.

source

Callbacks

PhaseFields.interface_position_1d — Function
interface_position_1d(φ, x; contour=0.5)

Find the position where φ crosses the contour value using linear interpolation.

Arguments

  • φ: Phase field array
  • x: Spatial coordinates array

Keyword Arguments

  • contour: Contour value to track (default: 0.5)

Returns

  • Interface position, or NaN if no crossing found
source
PhaseFields.solid_fraction — Function
solid_fraction(φ; threshold=0.5)

Calculate the fraction of the domain in the solid phase (φ > threshold).

Arguments

  • φ: Phase field array

Keyword Arguments

  • threshold: Value above which a point is considered solid (default: 0.5)

Returns

  • Solid fraction (0 to 1)
source
PhaseFields.create_interface_saving_callback — Function
create_interface_saving_callback(grid; contour=0.5, saveat=nothing, save_everystep=true)

Create a SavingCallback that tracks interface position during integration.

Arguments

  • grid: UniformGrid1D object

Keyword Arguments

  • contour: Contour value to track (default: 0.5)
  • saveat: Specific times to save (default: nothing, uses save_everystep)
  • save_everystep: Save at every solver step (default: true if saveat is empty)

Returns

  • Tuple of (callback, saved_values)

Example

cb, saved = create_interface_saving_callback(grid; saveat=0.1)
sol = solve(prob, Tsit5(); callback=cb)
# Access results:
times = saved.t
positions = saved.saveval
source
PhaseFields.create_solid_fraction_callback — Function
create_solid_fraction_callback(; threshold=0.5, saveat=nothing, save_everystep=true)

Create a SavingCallback that tracks solid fraction during integration.

Keyword Arguments

  • threshold: Value above which a point is considered solid (default: 0.5)
  • saveat: Specific times to save (default: nothing)
  • save_everystep: Save at every solver step (default: true if saveat is empty)

Returns

  • Tuple of (callback, saved_values)

Example

cb, saved = create_solid_fraction_callback(; saveat=0.1)
sol = solve(prob, Tsit5(); callback=cb)
# Access results:
times = saved.t
fractions = saved.saveval
source
PhaseFields.create_steady_state_callback — Function
create_steady_state_callback(; abstol=1e-8, reltol=1e-6, min_t=nothing)

Create a TerminateSteadyState callback that stops integration when the solution reaches steady state.

Keyword Arguments

  • abstol: Absolute tolerance for derivatives (default: 1e-8)
  • reltol: Relative tolerance for derivatives (default: 1e-6)
  • min_t: Minimum time before termination is allowed (default: nothing)

Returns

  • TerminateSteadyState callback

Example

cb = create_steady_state_callback(abstol=1e-6, reltol=1e-4, min_t=1.0)
sol = solve(prob, Tsit5(); callback=cb)
# sol.t[end] will be the time when steady state was reached
source
PhaseFields.create_phase_field_callbacks — Function
create_phase_field_callbacks(grid;
    track_interface=true, track_solid_fraction=false,
    terminate_steady_state=false, kwargs...)

Create a combined set of callbacks for phase field simulations.

Arguments

  • grid: UniformGrid1D object (required if track_interface=true)

Keyword Arguments

  • track_interface: Track interface position (default: true)
  • track_solid_fraction: Track solid fraction (default: false)
  • terminate_steady_state: Enable steady state termination (default: false)
  • contour: Contour value for interface tracking (default: 0.5)
  • threshold: Threshold for solid fraction (default: 0.5)
  • saveat: Times to save tracked values (default: nothing)
  • abstol: Steady state absolute tolerance (default: 1e-8)
  • reltol: Steady state relative tolerance (default: 1e-6)
  • min_t: Minimum time before steady state termination (default: nothing)

Returns

  • Named tuple with:
    • callback: Combined CallbackSet
    • interface_data: SavedValues for interface (or nothing)
    • solid_fraction_data: SavedValues for solid fraction (or nothing)

Example

result = create_phase_field_callbacks(grid;
    track_interface=true,
    terminate_steady_state=true,
    saveat=0.1)

sol = solve(prob, Tsit5(); callback=result.callback)

# Access tracked data
interface_times = result.interface_data.t
interface_positions = result.interface_data.saveval
source

See Also