Adaptivity

Adaptivity

Adaptive mesh refinement for the Pyrolysis.jl package.

This module provides:

• Error indicators (gradient, curvature, reaction rate, composite)
• h-refinement (cell splitting)
• h-coarsening (cell merging)
• r-refinement (node relocation via equidistribution)
• Mesh quality metrics (aspect ratio, smoothness, Jacobian)
• AMR controller for orchestrating adaptation
• Cell depletion detection and handling (both surface and interior)

Usage

using Pyrolysis.Adaptivity

# Create error indicator
indicator = GradientIndicator(:T)

# Compute indicators
η = compute_indicator(indicator, mesh)

# Create AMR controller
controller = AMRController(
	error_indicator = indicator,
	θ_refine = 0.7,
	θ_coarsen = 0.1,
	h_min = 1e-6,
	max_level = 4
)

# Adapt mesh
result = adapt_mesh!(mesh, controller)

# Handle depletion of a thin cell by merging it into a neighbor
neighbor = find_merge_neighbor(mesh, cell_id)
if neighbor > 0
	merge_adjacent_cells!(mesh, cell_id, neighbor; material = material)
end

AspectRatioMetric <: AbstractQualityMetric

Limits the ratio of adjacent cell sizes.

Fields

• max_ratio::Float64: Maximum allowed ratio Δzi / Δz{i±1}

CompositeIndicator <: AbstractErrorIndicator

Combines multiple error indicators.

Component indicators carry incommensurable dimensions (a temperature-gradient indicator is in K, a reaction-rate indicator in 1/s, …), so each component is normalized by its own maximum over active cells before combination. The combined η is dimensionless; the component weights act on these normalized scales.

Fields

• indicators::Vector{AbstractErrorIndicator}: Component indicators
• combination::Symbol: How to combine (:max, :l2, :sum)

CompositeMonitor <: MonitorFunction

Combines multiple monitor functions.

Fields

• monitors::Vector{MonitorFunction}
• combination::Symbol: :max, :sum, or :product

CompositeQualityMetric <: AbstractQualityMetric

Combines multiple quality metrics.

Fields

• metrics::Vector{AbstractQualityMetric}

CurvatureIndicator <: AbstractErrorIndicator

Curvature-based error indicator: ηi = hi² · |∇²φ|_i

More sensitive to sharp fronts than gradient indicator.

Fields

• variable::Union{Symbol,Int}: Variable to track — :T for temperature, or a component index j to track ξ[j, :]
• weight::Float64: Weight for combining with other indicators

GradientIndicator <: AbstractErrorIndicator

Gradient-based error indicator: ηi = hi · |∇φ|_i

Refines where gradients are steep (temperature fronts, reaction zones).

Fields

• variable::Union{Symbol,Int}: Variable to track — :T for temperature, or a component index j to track ξ[j, :]
• weight::Float64: Weight for combining with other indicators

GradientMonitor <: MonitorFunction

Monitor based on solution gradient: ω = √(1 + α|∇φ|²)

Fields

• α::Float64: Scaling parameter for gradient sensitivity
• variable::Symbol: Variable to monitor (:T or component index)

InterfaceMonitor <: MonitorFunction

Monitor that concentrates nodes near an interface.

ω = ωbase + ωpeak · exp(-(z - z_interface)² / (2σ²))

Fields

• z_interface::Float64: Interface position
• σ::Float64: Width of concentration zone
• ω_base::Float64: Base monitor value
• ω_peak::Float64: Peak monitor value at interface

InterfaceProximityIndicator <: AbstractErrorIndicator

Interface proximity indicator: ηi = 1 / (1 + dinterface / h_target)

Forces refinement near tracked interfaces.

Fields

• interface_position::Float64: Current interface z-position
• h_target::Float64: Target cell size near interface
• weight::Float64: Weight for combining with other indicators

JacobianMetric <: AbstractQualityMetric

Ensures positive cell volumes (Jacobian > 0).

Fields

• min_jacobian::Float64: Minimum allowed volume (typically 0 or small positive)

MinCellSizeMetric <: AbstractQualityMetric

Ensures cells don't become too small.

Fields

• h_min::Float64: Minimum allowed cell size

MonitorFunction

Abstract type for monitor functions used in r-refinement.

ReactionRateIndicator <: AbstractErrorIndicator

Reaction rate error indicator: ηi = hi · maxj(|rj|i / max(ξj, ε))

Refines where reactions are active.

Requires per-cell reaction rates: compute_indicator(ind, mesh, T, ξ; rates = r) with r sized like ξ for the current mesh. It therefore cannot run inside the adapt_mesh! cycle (which recomputes indicators after topology changes, when any precomputed rates are stale) — use it standalone, or supply gradient/curvature indicators to the controller instead.

Fields

• weight::Float64: Weight for combining with other indicators

SmoothnessMetric <: AbstractQualityMetric

Limits the size jump between adjacent cells.

Fields

• max_size_jump::Float64: Maximum allowed |Δzi - Δz{i±1}| / Δz_i

are_cells_adjacent(mesh::Mesh1D, cell_a_id::Int, cell_b_id::Int) -> Bool

Check if two cells share a face.

assert_amr_invariants(mesh::Mesh1D)

Cheap consistency check on AMR bookkeeping. Throws an AssertionError if n_active_cells does not match the active mask, or if any of the parallel per-cell arrays disagree in length. Intended to be called after every refinement / coarsening step.

check_quality(mesh::Mesh1D, metric::AbstractQualityMetric) -> Bool

Check if mesh satisfies quality metric.

Returns

• true if mesh passes quality check, false otherwise

compact_mesh_after_merge!(mesh::Mesh1D{NC}, removed_cell_id::Int) where NC

Remove an inactive cell and its orphaned node/face from mesh arrays. Works for any cell position (surface, interior, or bottom-adjacent).

This is the generalized version of compact_mesh_after_surface_merge!.

Algorithm

After mergeadjacentcells!, the mesh has:

• One cell marked inactive (the removed cell)
• One orphaned node (not connected to any active cell)
• One orphaned interior face (was between the merged cells)

This function:

1. Identifies orphaned elements
2. Removes elements using SHIFT (O(N) per element, maintains ordering)
3. Renumbers all indices consistently (simple decrement for elements after removed)
4. Updates connectivity and geometry

Note: We use shift instead of swap-and-pop to maintain the invariant that cell/node/face indices correspond to spatial ordering (lower index = lower z). This is important for 1D mesh operations that rely on left/right semantics.

Arguments

• mesh: Mesh to compact (modified in place)
• removed_cell_id: ID of the cell that was merged away (should be inactive)

Returns

Nothing. Modifies mesh in place.

compact_mesh_after_surface_merge!(mesh::Mesh1D{NC}) where NC

Remove the inactive surface cell and orphaned surface node/face from mesh arrays.

This function should be called AFTER merge_surface_cell! to actually remove the inactive cell from the arrays. This is necessary because:

1. merge_surface_cell! only marks the surface cell inactive
2. The ODE state vector actually removes elements
3. The mesh arrays must match the state vector length

Algorithm

After surface merge, the mesh has:

• Surface cell marked inactive (last active becomes new surface)
• Surface node orphaned (not connected to any active cell)
• Interior face between merged cells orphaned

This function:

1. Identifies orphaned elements (cell, node, face)
2. Removes them from arrays
3. Renumbers all indices consistently
4. Updates connectivity
5. Updates geometry caches

Returns

Nothing. Modifies mesh in place.

Note

This is specifically for surface cell removal where the surface cell is always the highest-indexed active cell. For general compaction with arbitrary inactive cells, a more complex algorithm would be needed.

compute_aspect_ratios(mesh::Mesh1D) -> Vector{Float64}

Compute aspect ratio for each cell (ratio to smaller neighbor).

compute_cell_gradient(φ::Vector{Float64}, mesh::Mesh1D, i::Int) -> Float64

Compute gradient of field φ at cell i.

compute_curvature_1d(φ::Vector{Float64}, mesh::Mesh1D, i::Int) -> Float64

Compute second derivative of field φ at cell i.

compute_gradient_1d(φ::Vector{Float64}, mesh::Mesh1D, i::Int) -> Float64

Compute gradient of field φ at cell i using central differences.

compute_indicator(ind::AbstractErrorIndicator, mesh::Mesh1D, 
				  T::Vector{Float64}, ξ::Matrix{Float64}) -> Vector{Float64}

Compute error indicator values for all active cells.

Arguments

• ind: Error indicator type
• mesh: 1D mesh
• T: Temperature at cell centers [K]
• ξ: Mass concentrations [kg/m³], size (NC, n_cells)

Returns

• Vector of indicator values for each cell

compute_indicator(ind::AbstractErrorIndicator, mesh::Mesh1D; kwargs...) -> Vector{Float64}

Convenience method that extracts fields from mesh automatically.

compute_smoothness(mesh::Mesh1D) -> Vector{Float64}

Compute smoothness metric for each cell.

compute_target_positions(mesh::Mesh1D, ω::Vector{Float64}) -> Vector{Float64}

Compute target node positions using equidistribution.

Arguments

• mesh: Current mesh
• ω: Monitor function values at cell centers

Returns

• Vector of target z-positions for interior nodes

enforce_quality!(mesh::Mesh1D, metric::AbstractQualityMetric; 
				 max_iter::Int=10) -> Int

Attempt to fix quality violations through node smoothing.

Arguments

• mesh: Mesh to modify
• metric: Quality metric to satisfy
• max_iter: Maximum smoothing iterations

Returns

• Number of violations remaining after enforcement

evaluate_monitor(mon::MonitorFunction, mesh::Mesh1D, 
				 T::Vector{Float64}, ξ::Matrix{Float64}) -> Vector{Float64}

Evaluate monitor function at cell centers.

extract_fields(mesh::Mesh1D{NC}) -> (T::Vector{Float64}, ξ::Matrix{Float64})

Extract temperature and concentration fields from mesh cell states.

extract_state_arrays(mesh::Mesh1D{NC}) -> (T, ξ)

Extract temperature and concentration arrays from mesh.

find_bottom_cell(mesh::Mesh1D) -> Int

Find the ID of the bottom cell (lowest z, active).

This is the cell at z = 0 (substrate side).

find_exterior_neighbor(mesh::Mesh1D, cell_id::Int) -> Int

Find the neighboring cell toward the exterior (higher z).

For an interior cell, this returns the adjacent cell at higher z. Returns 0 if no exterior neighbor exists (cell is at top boundary).

find_interior_neighbor(mesh::Mesh1D, cell_id::Int) -> Int

Find the neighboring cell toward the interior (lower z).

For a surface cell at z = L, this returns the adjacent cell at lower z. Returns 0 if no interior neighbor exists (cell is at bottom boundary).

find_merge_neighbor(mesh::Mesh1D, cell_id::Int) -> Int

Find the best neighbor to merge with for a thin cell.

Selection strategy:

1. If cell has only one neighbor (boundary cell): use that neighbor
2. If cell has two neighbors: prefer interior neighbor (toward substrate)

The interior neighbor is preferred for thermal stability since the substrate side typically has more stable thermal conditions.

Arguments

• mesh: Current mesh
• cell_id: ID of the thin cell

Returns

Neighbor cell ID, or 0 if no valid neighbor exists

find_shared_face(mesh::Mesh1D, cell_a_id::Int, cell_b_id::Int) -> Int

Find the face shared between two adjacent cells.

Returns

• Face ID of shared face, or 0 if cells are not adjacent

find_surface_cell(mesh::Mesh1D) -> Int

Find the ID of the current surface cell (highest z, active).

In the coordinate convention:

• z = 0 is the bottom (substrate)
• z = L is the top (exposed surface)

The surface cell is the active cell with the highest z center.

merge_adjacent_cells!(mesh::Mesh1D{NC}, cell_to_remove_id::Int, 
					  cell_to_keep_id::Int; material=nothing) where NC -> Int

Merge two adjacent cells, removing one and expanding the other.

This is an extended h-coarsening that works on ANY adjacent pair, not just sibling cells from refinement. The key difference from coarsen_cells! is that we don't require a parent relationship.

Arguments

• mesh: The mesh to modify
• cell_to_remove_id: Cell to remove (typically depleted surface cell)
• cell_to_keep_id: Cell to expand (will absorb the removed cell's volume/mass)
• material=nothing: Optional material for energy-conserving temperature merge. If provided, uses ρcₚ-weighted energy conservation. Otherwise, uses simple volume-weighted averaging.

Returns

• ID of the merged cell (same as celltokeep_id)

Algorithm

1. Verify cells are adjacent
2. Conservative merge: combine states with mass/energy conservation
3. Expand the kept cell's geometry to encompass both
4. Update face connectivity
5. Mark removed cell as inactive

Physical Use Case

During pyrolysis, surface cells become depleted as solid decomposes. When the surface cell is mostly gas, we merge it with the interior neighbor to remove the singularity in dT/dt = Q/(ρcₚ).

Conservation

• Mass: ξmerged = (Va * ξa + Vb * ξb) / Vmerged
• Energy: When material is provided, uses ρcₚ-weighted energy conservation. Otherwise, uses simple volume-weighted temperature averaging.

Note

Unlike coarsen_cells!, this function:

• Does NOT require a parent-child relationship
• Can work on any adjacent cell pair
• Is designed for removing depleted cells, not low-error coarsening

merge_surface_cell!(mesh::Mesh1D{NC}; material=nothing) where NC -> Int

Convenience function to merge the depleted surface cell with its interior neighbor.

Arguments

• mesh: The mesh to modify
• material=nothing: Optional material for energy-conserving temperature merge. If provided, uses ρcₚ-weighted energy conservation. Otherwise, uses simple volume-weighted averaging.

Returns

• ID of the merged cell, or 0 if no merge was possible

Usage

using Pyrolysis.Adaptivity

# When surface depletion is detected
if is_surface_depleted(mesh, material, 0.01)
	merged_id = merge_surface_cell!(mesh; material=material)
	if merged_id > 0
		# Update state vector etc.
	end
end

merge_temperatures_conservatively(V_a::Float64, state_a::CellState{NC},
								   V_b::Float64, state_b::CellState{NC},
								   material::Material{NC}) where NC -> Float64

Compute the merged temperature by conserving the condensed-phase sensible enthalpy.

Principle

The conserved functional is

H(T) = Σ_{j ∉ gas} m_j ∫_{T_REF}^{T} cₚ,j(T′) dT′ ,   m_j = V·ξ_j

and T_merged solves H_merged(T_merged) = H_a(T_a) + H_b(T_b). Two choices here are deliberate:

• ∫cₚ dT, not cₚ(T)·T. For temperature-dependent cₚ the product form is not an enthalpy: with cₚ = a + bT it carries the b-term with coefficient 1 instead of ½, a few-K systematic bias per merge when a hot surface cell merges into a cooler interior one. The integral form is exact for any cₚ(T) (4-point Gauss–Legendre per component). The datum T_REF cancels because each component's total mass is identical before and after the merge.
• Gas components are excluded. Formulation A stores thermal energy in the condensed matrix only (the gas is thermally quasi-steady), so gas heat capacity does not belong in the merge weighting.

The solve is a bracketed Newton iteration on [min(T_a,T_b), max(T_a,T_b)]; H is strictly monotone (cₚ > 0) so the root is unique and always inside the bracket.

Arguments

• V_a: Volume of cell a [m³]
• state_a: Cell state for cell a (T, ξ)
• V_b: Volume of cell b [m³]
• state_b: Cell state for cell b (T, ξ)
• material: Material definition with component heat capacity functions

Returns

Merged temperature [K] conserving condensed-phase sensible enthalpy.

Edge case

If neither cell carries condensed heat capacity (both essentially pure gas), the condensed energy balance is degenerate and the volume-weighted average is returned.

merge_temperatures_conservatively_from_values(V_a::Float64, T_a::Float64, ξ_a::NTuple{NC,Float64},
											   V_b::Float64, T_b::Float64, ξ_b::NTuple{NC,Float64},
											   material::Material{NC}) where NC -> Float64

Alternate interface for energy-conserving temperature merge using raw values instead of CellState. Useful when working directly with state vectors.

mesh_quality_summary(mesh::Mesh1D) -> NamedTuple

Compute summary statistics of mesh quality.

Returns

NamedTuple with:

• n_active_cells: Number of active cells
• min_cell_size: Minimum cell size
• max_cell_size: Maximum cell size
• max_aspect_ratio: Maximum aspect ratio
• max_smoothness: Maximum smoothness deviation
• all_positive_volume: Whether all volumes are positive

quality_violations(mesh::Mesh1D, metric::AbstractQualityMetric) -> Vector{Int}

Find cells that violate quality metric.

Returns

• Vector of cell indices that fail the quality check

relocate_nodes!(mesh::Mesh1D{NC}; monitor::MonitorFunction=GradientMonitor(),
				max_iter::Int=10, relaxation::Float64=0.5, 
				tol::Float64=1e-3, h_min::Float64=1e-8) where NC

Relocate interior nodes using equidistribution.

Arguments

• mesh: Mesh to modify (in-place)
• monitor: Monitor function for equidistribution
• max_iter: Maximum iterations
• relaxation: Relaxation factor (0 < β ≤ 1)
• tol: Convergence tolerance (max relative node movement)
• h_min: Minimum allowed cell size

Returns

• Number of iterations performed