Adaptivity

Adaptivity

Mesh adaptation for Pyrolysis.jl.

solve uses the depleted-cell parking: a cell that burns out is merged into a neighbour (find_merge_neighbor, park_cell!, with a conservative temperature merge) and its slot moves above the live surface, so the state vector keeps its size.

The module also holds building blocks for mesh refinement that solve does not use: error indicators (gradient, curvature, reaction rate, interface proximity, composite), r-refinement (node relocation by equidistribution of a monitor function) and mesh-quality metrics. The h-refinement, h-coarsening and AMR-controller routines built on them live in Pyrolysis.Experimental. A mesh they refine cannot be passed to solve: refinement leaves inactive parent cells among the active ones, and the solver's kernels assume the active cells come first.

Usage

using Pyrolysis.Adaptivity

# Park a depleted cell: merge it into its neighbour and move the freed slot
# above the live surface
neighbor = find_merge_neighbor(mesh, cell_id)
park_cell!(mesh, cell_id, neighbor; material = material, ale = true)

# Refinement indicator on a mesh
η = compute_indicator(GradientIndicator(:T), mesh)

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

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

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

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

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.

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

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 sensible enthalpy of every component.

Principle

The conserved functional is

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

summed over every component, gas included — the same ρcₚ the temperature equation stores — and T_merged solves H_merged(T_merged) = H_a(T_a) + H_b(T_b). The integral ∫cₚ dT, not cₚ(T)·T, is used deliberately: 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.

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 sensible enthalpy.

Edge case

If neither cell carries any heat capacity (both empty), the energy balance is degenerate and the volume-weighted average is returned.

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

park_cell!(mesh::Mesh1D{NC}, i_park::Int, i_keep::Int;
           material::Material{NC}, ale::Bool) -> Int

Merge live cell i_park into its live neighbour i_keep and move the freed slot into the ghost range. Returns the new live-cell count.

Conservation: the component masses V·ξⱼ of both cells are summed into the keeper (ξ_keep = (V_p ξ_p + V_k ξ_k) / V_keep, with V_keep = V_p + V_k under ALE and the unchanged V_k on a fixed mesh), the temperature solves the sensible-enthalpy balance over every component (merge_temperatures_conservatively), and mesh.m_dry_solid_initial is summed into the keeper. The ghost slot carries the merged temperature, zero concentrations and — under ALE — zero volume.

Under ALE the cells above i_park shift down one index and the nodes above the shared node shift down one index, so the live geometry is unchanged and the ghost node lands on the live surface. On a fixed mesh (ale = false) only the top live cell may be parked; nothing shifts.

The exposed boundary moves to face n_live + 1 (Face1D/Node1D boundary flags and mesh.boundary_faces[1]), mesh.n_active_cells and mesh.active_cell_mask are updated, and the cell volumes and centres are recomputed (update_geometry!).

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