Include PIDE integrals from x to the upper boundary of the domain

src/discretization/generate_finite_difference_rules.jl

@@ -77,24 +77,29 @@ function generate_finite_difference_rules(
         # Spherical diffusion scheme
         spherical_diffusion_rules = generate_spherical_diffusion_rules(
             II, s, depvars, derivweights, bmap, indexmap, split_additive_terms(pde))
-        integration_rules = vec(generate_euler_integration_rules(
-            II, s, depvars, indexmap, terms))
+
+        # Integration schemes for boundaries: [lower-bound, x] and [x, upper-bound]
+        integration_rules_from_boundary_to_x = generate_euler_integration_rules(
+            II, s, depvars, indexmap, terms)
     else
         central_deriv_rules_cartesian = []
         advection_rules = []
         nonlinlap_rules = []
         spherical_diffusion_rules = []
         mixed_deriv_rules_cartesian = []
-        integration_rules = []
+        integration_rules_from_boundary_to_x = []
     end
 
     cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
 
-    integration_rules = vcat(integration_rules,
-        vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms)))
+    # Integration schemes for boundaries: [lower-bound, upper-bound]
+    whole_domain_integration_rules = generate_whole_domain_integration_rules(
+        II, s, depvars, indexmap, terms)
+
     return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules),
         vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian),
-        vec(advection_rules), integration_rules)
+        vec(advection_rules), 
+        vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules))
 end
 
 function generate_finite_difference_rules(
@@ -112,10 +117,12 @@ function generate_finite_difference_rules(
     advection_rules = []
     nonlinlap_rules = []
     spherical_diffusion_rules = []
-    integration_rules = []
+    integration_rules_from_boundary_to_x = []
+
+    whole_domain_integration_rules = generate_whole_domain_integration_rules(
+        II, s, depvars, indexmap, terms)
 
-    integration_rules = vcat(integration_rules,
-        vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms)))
     return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules),
-        vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules)
+        vec(central_deriv_rules_cartesian), vec(advection_rules), 
+        vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules))
 end

src/discretization/schemes/integral_expansion/integral_expansion.jl

@@ -1,66 +1,90 @@
 # use the trapezoid rule
-function _euler_integral(II, s, jx, u, ufunc, dx::Number) #where {T,N,Wind,DX<:Number}
+function _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number}
     j, x = jx
-    if II[j] == 1
+    if II[j] == 1 # recursively arrived at lower end of the domain
         return Num(0)
     end
     # unit index in direction of the derivative
     I1 = unitindex(ndims(u, s), j)
     # dx for multiplication
     Itap = [II - I1, II]
-    weights = [dx / 2, dx / 2]
 
+    if dx isa Number # Uniform dx
+        weights = [dx / 2, dx / 2]
+    elseif dx isa AbstractVector # Nonuniform dx
+        weights = fill(dx[II[j] - 1] / 2, 2)
+    else
+        error("Unsupported type of dx: $(type(dx))")
+    end
+
+    # sym_do computes from II to II - I1, 
+    # and recursive call computes from II - I1 to lower end of domain
     return sym_dot(weights, ufunc(u, Itap, x)) +
-           _euler_integral(II - I1, s, jx, u, ufunc, dx)
+           _euler_integral_lower_to_II(II - I1, s, jx, u, ufunc, dx)
 end
-
-# Nonuniform dx
-function _euler_integral(II, s, jx, u, ufunc, dx::AbstractVector) #where {T,N,Wind,DX<:Number}
+function _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number}
     j, x = jx
-    if II[j] == 1
+    if II[j] == length(s, x) # recursively arrived at upper end of the domain
         return Num(0)
     end
     # unit index in direction of the derivative
     I1 = unitindex(ndims(u, s), j)
     # dx for multiplication
-    Itap = [II - I1, II]
-    weights = fill(dx[II[j] - 1] / 2, 2)
+    Itap = [II, II + I1]
 
+    if dx isa Number # Uniform dx
+        weights = [dx / 2, dx / 2]
+    elseif dx isa AbstractVector # Nonuniform dx
+        weights = fill(dx[II[j] + 1] / 2, 2)
+    else
+        error("Unsupported type of dx: $(type(dx))")
+    end
+
+    # sym_do computes from II to II + I1, 
+    # and recursive call computes from II + I1 to upper end of domain
     return sym_dot(weights, ufunc(u, Itap, x)) +
-           _euler_integral(II - I1, s, jx, u, ufunc, dx)
+           _euler_integral_II_to_upper(II + I1, s, jx, u, ufunc, dx)
 end
 
-function euler_integral(II, s, jx, u, ufunc)
+# An integral from II to end of domain
+function euler_integral(method::Symbol, II, s, jx, u, ufunc)
     j, x = jx
     dx = s.dxs[x]
-    return _euler_integral(II, s, jx, u, ufunc, dx)
+    if method == :lower_boundary_to_x
+        return _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx)
+    elseif method == :x_to_upper_boundary
+        return _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx)
+        # error("Method :x_to_upper_boundary is not implemented for euler_integral.")
+    else
+        error("Unsupported method for euler_integral: $method")
+    end
 end
 
 # An integral across the whole domain (xmin .. xmax)
 function whole_domain_integral(II, s, jx, u, ufunc)
     j, x = jx
     dx = s.dxs[x]
-    if II[j] == length(s, x)
-        return _euler_integral(II, s, jx, u, ufunc, dx)
-    end
-
-    dist2max = length(s, x) - II[j]
+    dist2max = length(s, x) - II[j] # distance to the end of the domain
     I1 = unitindex(ndims(u, s), j)
     Imax = II + dist2max * I1
-    return _euler_integral(Imax, s, jx, u, ufunc, dx)
+
+    return _euler_integral_lower_to_II(Imax, s, jx, u, ufunc, dx)
 end
 
 @inline function generate_euler_integration_rules(
         II::CartesianIndex, s::DiscreteSpace, depvars, indexmap, terms)
     ufunc(u, I, x) = s.discvars[u][I]
 
-    eulerrules = reduce(safe_vcat,
-        [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1],
-              Num(x)))(u) => euler_integral(
-              Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc)
-          for x in ivs(u, s)]
-         for u in depvars],
-        init = [])
+    eulerrules = reduce(safe_vcat, [
+        reduce(safe_vcat, [
+        [
+            # integrals from lower domain end to x:
+            Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(:lower_boundary_to_x, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), 
+            # integrals from x to upper domain end:
+            Integral(x in DomainSets.ClosedInterval(Num(x), s.vars.intervals[x][2]))(u) => euler_integral(:x_to_upper_boundary, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc)
+        ]
+    for x in ivs(u, s)]) for u in depvars] )
+
     return eulerrules
 end
 

src/system_parsing/pde_system_transformation.jl

@@ -156,22 +156,24 @@ function descend_to_incompatible(term, v)
             end
         elseif op isa Integral
             if any(isequal(op.domain.variables), v.x̄)
-                euler = isequal(
-                    op.domain.domain.left, v.intervals[op.domain.variables][1]) &&
-                        isequal(op.domain.domain.right, Num(op.domain.variables))
-                whole = isequal(
-                    op.domain.domain.left, v.intervals[op.domain.variables][1]) &&
+                euler = ( # from lower domain boundary to x:
+                        isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) &&
+                        isequal(op.domain.domain.right, Num(op.domain.variables))) ||
+                        (# from x to upper domain boundary:
+                        isequal(op.domain.domain.left, Num(op.domain.variables)) &&
+                        isequal(op.domain.domain.right, v.intervals[op.domain.variables][2]))
+                whole = isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) &&
                         isequal(op.domain.domain.right, v.intervals[op.domain.variables][2])
                 if any([euler, whole])
                     u = arguments(term)[1]
                     out = check_deriv_arg(u, v)
                     @assert out==(nothing, false) "Integral $term must be purely of a variable, got $u. Try wrapping the integral argument with an auxiliary variable."
                     return (nothing, nothing, false)
                 else
-                    throw(ArgumentError("Integration Domain only supported for integrals from start of iterval to the variable, got $(op.domain.domain) in $(term)"))
+                    throw(ArgumentError("Integration Domain only supported for integrals across whole interval, or from an interval boundary (upper or lower) to the independent variable, got $(op.domain.domain) in $(term)"))
                 end
             else
-                throw(ArgumentError("Integral must be with respect to the independent variable in its upper-bound, got $(op.domain.variables) in $(term)"))
+                throw(ArgumentError("Integral must be with respect to the independent variable in its upper- or lower-bound, got $(op.domain.variables) in $(term)"))
             end
         end
         for arg in SU.arguments(term)