Add nonuniform grid support for interface boundary conditions with UpwindScheme

src/MOL_discretization.jl

@@ -16,13 +16,54 @@ function PDEBase.interface_errors(
     end
 end
 
+function _interface_physical_coords(b, grid1, grid2)
+    if isupper(b)
+        return grid1[end], grid2[1]
+    else
+        return grid1[1], grid2[end]
+    end
+end
+
+function _interface_coords_aligned(coord1, coord2, grid1, grid2)
+    scale = max(abs(grid1[end] - grid1[1]), abs(grid2[end] - grid2[1]))
+    return isapprox(coord1, coord2; atol = sqrt(eps(float(one(scale)))) * scale)
+end
+
+function PDEBase.check_boundarymap(
+        boundarymap, v::PDEBase.VariableMap, discretization::MOLFiniteDifference
+    )
+    return _check_interface_boundarymap(boundarymap, discretization)
+end
+
+# Kept for backwards compatibility with direct 2-arg calls; the PDEBase
+# discretization pipeline invokes the 3-arg hook above.
 function PDEBase.check_boundarymap(boundarymap, discretization::MOLFiniteDifference)
+    return _check_interface_boundarymap(boundarymap, discretization)
+end
+
+function _check_interface_boundarymap(boundarymap, discretization::MOLFiniteDifference)
     bs = filter_interfaces(flatten_vardict(boundarymap))
+    scheme = discretization.advection_scheme
     for b in bs
         dx1 = discretization.dxs[Num(b.x)]
         dx2 = discretization.dxs[Num(b.x2)]
-        if dx1 != dx2
-            throw(ArgumentError("The step size of the connected variables $(b.x) and $(b.x2) must be the same. If you need nonuniform interface boundaries please post an issue on GitHub."))
+        if dx1 isa AbstractVector && dx2 isa AbstractVector
+            if scheme isa UpwindScheme && scheme.order > 1
+                throw(ArgumentError("UpwindScheme(order=$(scheme.order)) is not yet supported with interface or periodic boundary conditions on nonuniform grids, please use the default first order `UpwindScheme()`."))
+            end
+            isequal(b.x, b.x2) && continue
+            coord1, coord2 = _interface_physical_coords(b, dx1, dx2)
+            if !_interface_coords_aligned(coord1, coord2, dx1, dx2)
+                throw(
+                    ArgumentError(
+                        "The physical coordinates at the interface between $(b.x) and $(b.x2) must match for nonuniform grids, got $coord1 and $coord2 at the interface. Please ensure the grids align at the interface boundary. Note that cross-domain periodic (ring) topologies are not supported on nonuniform grids."
+                    )
+                )
+            end
+        elseif dx1 isa AbstractVector || dx2 isa AbstractVector
+            throw(ArgumentError("The interface between $(b.x) and $(b.x2) mixes a scalar step size with a nonuniform grid vector, please supply the same kind of grid on both sides."))
+        elseif dx1 != dx2
+            throw(ArgumentError("The step size of the connected variables $(b.x) and $(b.x2) must be the same."))
         end
     end
     return
@@ -46,7 +87,7 @@ function get_discrete(pdesys, discretization)
     # Create a map of each variable to their boundary conditions including initial conditions
     boundarymap = PDEBase.parse_bcs(get_bcs(pdesys), v, bcorders)
     # Check that the boundary map is valid
-    PDEBase.check_boundarymap(boundarymap, discretization)
+    PDEBase.check_boundarymap(boundarymap, v, discretization)
 
     # Transform system so that it is compatible with the discretization
     if should_transform(pdesys, discretization, boundarymap)

src/discretization/schemes/upwind_difference/upwind_difference.jl

@@ -27,13 +27,116 @@
     return weights, Itap
 end
 
+# Independent variable for a stencil tap on the far side of an interface;
+# returns b.x2 when I.A matches the coupled variable of boundary b, otherwise x.
+@inline function _upwind_stencil_grid_iv(I::RefCartesianIndex, s, bs, j, x)
+    if I.A !== nothing
+        for b in bs
+            if I.A === s.discvars[depvar(b.u2, s)]
+                return b.x2
+            end
+        end
+    end
+    return x
+end
+
+# Domain length L for same-domain periodic interfaces (x ~ x); nothing for
+# cross-domain joints (x1 ~ x2).
+@inline function _upwind_stencil_periodic_length(s, bs, j, x)
+    for b in bs
+        if b isa InterfaceBoundary && isequal(b.x, x) && isequal(b.x2, x)
+            g = s.grid[x]
+            return g[end] - g[firstindex(g)]
+        end
+    end
+    return nothing
+end
+
+# Physical coordinate of stencil tap I on the nonuniform grid, resolving
+# RefCartesianIndex taps to the correct domain. When `raw_j` is supplied
+# (index before bwrap), a same-domain periodic wrap is detected as
+# tap_j != raw_j and the coordinate is shifted by ±L so Fornberg stencils
+# stay monotone; bwrap also wraps raw_j == 1 (u[1] ~ u[end]), which a
+# bounds test on raw_j alone would miss.
+@inline function _upwind_stencil_coord(I, s, bs, j, x; raw_j = nothing)
+    grid_iv = I isa RefCartesianIndex ? _upwind_stencil_grid_iv(I, s, bs, j, x) : x
+    tap_j = I isa RefCartesianIndex ? I.I[j] : I[j]
+    coord = s.grid[grid_iv][tap_j]
+    if raw_j !== nothing && isequal(grid_iv, x) && tap_j != raw_j
+        L = _upwind_stencil_periodic_length(s, bs, j, x)
+        if L !== nothing
+            coord += raw_j > tap_j ? L : -L
+        end
+    end
+    return coord
+end
+
+# Index offsets of the one-sided upwind stencil along the derivative direction.
+@inline function _nonuniform_upwind_stencil_offsets(D, ispositive)
+    return if !ispositive
+        0:(D.stencil_length - 1)
+    else
+        (-D.stencil_length + 1):0
+    end
+end
+
+# True when the stencil leaves the domain or sits at a boundary point where
+# the standard nonuniform upwind stencil does not apply.
+@inline function _nonuniform_upwind_cross_domain_needed(
+        D::DerivativeOperator, II, s, bs, ispositive, u, jx
+    )
+    length(bs) == 0 && return false
+    j, x = jx
+    I1 = unitindex(ndims(u, s), j)
+    l = length(s, x)
+    ij = II[j]
+    offsets = _nonuniform_upwind_stencil_offsets(D, ispositive)
+    is_crossing = any(
+        i -> begin
+            idx = (II + i * I1)[j]
+            return idx < 1 || idx > l
+        end, offsets
+    )
+    is_boundary_blindspot = ispositive ? ij == l : ij == 1
+    return is_crossing || is_boundary_blindspot
+end
+
+# Build (weights, Itap) for nonuniform upwind differencing at interfaces via
+# bwrap and coordinate-based calculate_weights. Weights are calculated
+# dynamically here, as this is the point where the specific variable's
+# boundary map (bs) and cross-domain coordinates are naturally accessible.
+@inline function _nonuniform_upwind_interface_difference(
+        D::DerivativeOperator{T, N, Wind, DX}, II, s, bs,
+        ispositive, u, jx
+    ) where {T, N, Wind, DX <: AbstractVector}
+    j, x = jx
+    I1 = unitindex(ndims(u, s), j)
+    offsets = _nonuniform_upwind_stencil_offsets(D, ispositive)
+    raw_js = [(II + offsets[k] * I1)[j] for k in eachindex(offsets)]
+    Itap = [bwrap(II + offsets[k] * I1, bs, s, jx) for k in eachindex(offsets)]
+    for It in Itap
+        iv = It isa RefCartesianIndex ? _upwind_stencil_grid_iv(It, s, bs, j, x) : x
+        if !(1 <= It[j] <= length(s, iv))
+            throw(ArgumentError("The upwind stencil for $u along $x extends past a non-interface boundary on a nonuniform grid. Stencils wider than first order upwind (higher approximation or derivative orders) are not yet supported with nonuniform interface boundaries."))
+        end
+    end
+    x0 = _upwind_stencil_coord(II, s, bs, j, x)
+    coords = [
+        _upwind_stencil_coord(Itap[k], s, bs, j, x; raw_j = raw_js[k]) for k in eachindex(Itap)
+    ]
+    weights = calculate_weights(D.derivative_order, x0, coords)
+    return weights, Itap
+end
+
 @inline function _upwind_difference(
         D::DerivativeOperator{T, N, Wind, DX}, II, s, bs,
         ispositive, u, jx
     ) where {T, N, Wind, DX <: AbstractVector}
     j, x = jx
-    @assert length(bs) == 0 "Interface boundary conditions are not yet supported for nonuniform dx dimensions, such as $x, please post an issue to https://github.com/SciML/MethodOfLines.jl if you need this functionality."
     I1 = unitindex(ndims(u, s), j)
+    if _nonuniform_upwind_cross_domain_needed(D, II, s, bs, ispositive, u, jx)
+        return _nonuniform_upwind_interface_difference(D, II, s, bs, ispositive, u, jx)
+    end
     if !ispositive
         @assert D.offside == 0