SciML · xtalax · Mar 21, 2023 · Mar 2, 2023 · Mar 2, 2023 · Mar 2, 2023
diff --git a/Project.toml b/Project.toml
@@ -29,7 +29,7 @@ Latexify = "0.15"
 ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
 SafeTestsets = "0.0.1"
-SciMLBase = "1"
+SciMLBase = "1.90, 1.91"
 StaticArrays = "1"
 SymbolicUtils = "1"
 Symbolics = "5"

diff --git a/docs/src/advection_schemes.md b/docs/src/advection_schemes.md
@@ -5,7 +5,7 @@ Used as a keyword argument `advection_scheme` to `MOLFiniteDifference`.
 ```julia
 UpwindScheme(approx_order = 1)
 ```
-Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` keyword argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme).
+Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme).
 
 ## WENO Scheme of Jiang and Shu
 ```julia
@@ -21,4 +21,6 @@ Problems with first order derivatives which multiply one another will need to us
 
 Supports only first order derivatives, other odd order derivatives are unsupported with this scheme. At present, does not support Nonuniform grids, though this is a planned feature.
 
-Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/)
+Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/)
+
+## FunctionalScheme
diff --git a/docs/src/api/utils.md b/docs/src/api/utils.md
@@ -6,6 +6,6 @@ CurrentModule = MethodOfLines
 
 ```@docs
 unitindex
-params
+ivs
 Idx
 ```
diff --git a/docs/src/devnotes.md b/docs/src/devnotes.md
@@ -48,7 +48,7 @@ end
 # Note that indexmap is used along with the function `Idx` to create an equivalent index for the discrete form of `u`,
 # which may have a different number of dimensions to `II`
 function generate_central_difference_rules(II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict)
-    rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in params(u, s)] for u in depvars]
+    rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in ivs(u, s)] for u in depvars]
 
     rules = reduce(vcat, rules)
 

diff --git a/docs/src/tutorials/PIDE.md b/docs/src/tutorials/PIDE.md
@@ -56,7 +56,7 @@ To have an integral over the whole domain, be sure to wrap the integral in an au
 Due to a limitation, the whole domain integral needs to have the same arguments as the integrand, but is constant in x. To use it in an equation one dimension lower, use a boundary value like integral(t, 0.0)
 
 ```@example integrals2
-using MethodOfLines, ModelingToolkit, DomainSets
+using MethodOfLines, ModelingToolkit, DomainSets, OrdinaryDiffEq, Plots
 
 @parameters t, x
 @variables integrand(..) integral(..)
@@ -90,5 +90,6 @@ xdisc = sol[x]
 tdisc = sol[t]
 
 integralsol = sol[integral(t, x)]
-exact = [asf(t_) for t_ in tdisc, x_ in xdisc]
+
+heatmap(integralsol)
 ```
diff --git a/docs/src/tutorials/params.md b/docs/src/tutorials/params.md
@@ -2,7 +2,7 @@
 
 We can also build up more complicated systems with multiple dependent variables and parameters as follows
 
-```@example params1
+```@example ivs1
 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
 
 @parameters t x
@@ -50,7 +50,7 @@ gif(anim, "plot.gif",fps=30)
 
 The system does not need to be re-discretized every time we want to plot with different parameters, the system can be remade with new parameters with `remake`. See the `ModelingToolkit.jl` [docs](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/#Algebraic-relations-and-structural-simplification) for more ways to manipulate a `prob` post discretization.
 
-```@example params2
+```@example ivs2
 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
 
 @parameters t x

diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
@@ -4,8 +4,8 @@ function interface_errors(depvars, indvars, discretization)
     for x in indvars
         @assert haskey(discretization.dxs, Num(x)) || haskey(discretization.dxs, x) "Variable $x has no step size"
     end
-    if !(typeof(discretization.advection_scheme) ∈ [UpwindScheme, WENOScheme])
-        throw(ArgumentError("Only `UpwindScheme()` and `WENOScheme()` are supported advection schemes."))
+    if !any(s -> discretization.advection_scheme isa s,  [UpwindScheme, FunctionalScheme])
+        throw(ArgumentError("Only `UpwindScheme()` and `FunctionalScheme()` are supported advection schemes. Got $(typeof(discretization.advection_scheme))."))
     end
     if !(typeof(discretization.disc_strategy) ∈ [ScalarizedDiscretization])
         throw(ArgumentError("Only `ScalarizedDiscretization()` are supported discretization strategies."))
@@ -107,7 +107,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         eqvar = interiormap.var[pde]
 
         # * Assumes that all variables in the equation have same dimensionality except edgevals
-        args = params(eqvar, s)
+        args = ivs(eqvar, s)
         indexmap = Dict([args[i] => i for i in 1:length(args)])
         if disc_strategy isa ScalarizedDiscretization
             # Generate the equations for the interior points
@@ -147,7 +147,8 @@ function SciMLBase.discretize(pdesys::PDESystem,discretization::MethodOfLines.MO
                     return prob
                 else
                     f = ODEFunction(pdesys, discretization, analytic=analytic, discretization.kwargs..., kwargs...)
-                    return ODEProblem(f, prob.u0, prob.tspan; discretization.kwargs..., kwargs...)
+
+                    return ODEProblem(f, prob.u0, prob.tspan, prob.p; discretization.kwargs..., kwargs...)
                 end
             end
         end

diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl
@@ -64,13 +64,13 @@ end
 
 function get_gridloc(u, s)
     if isequal(operation(u), getindex)
+        # Remember arguments of getindex have u(t) first
         return _get_gridloc(s, arguments(u)...)
     else
         return (operation(u), [])
     end
 end
 
-
 function generate_function_from_gridlocs(analyticmap, gridlocs, s)
     is_t_first_map = Dict(map(s.ū) do u
         operation(u) => (findfirst(x -> isequal(s.time, x), arguments(u)) == 1)
@@ -84,14 +84,14 @@ function generate_function_from_gridlocs(analyticmap, gridlocs, s)
         is_t_first = is_t_first_map[uop]
         _f = analyticmap[opsmap[uop]]
         if is_t_first
-            return t -> _f(t, x̄...)
+            return (p, t) -> _f(p, t, x̄...)
         else
-            return t -> _f(x̄..., t)
+            return (p, t) -> _f(p, x̄..., t)
         end
     end
 
     f = (u0, p, t) -> map(fs_) do f_
-        f_(t)
+        f_(p, t)
     end
 
     return f
@@ -100,7 +100,7 @@ end
 function newindex(u_, II, s, indexmap)
     u = depvar(u_, s)
     args_ = remove(arguments(u_), s.time)
-    args = params(u, s)
+    args = ivs(u, s)
     is = map(enumerate(args_)) do (j, x)
         if haskey(indexmap, x)
             II[indexmap[x]]
@@ -135,8 +135,14 @@ end
 function chebyspace(N, dom)
     interval = dom.domain
     a, b = DomainSets.infimum(interval), DomainSets.supremum(interval)
-    x = reverse([(a + b) / 2 + (b - a) / 2 * cos(π * (2k - 1) / (2N)) for k in 1:N])
+    x = reverse([(a + b) / 2 + (b - a) / 2 * cospi((2k - 1) / (2N)) for k in 1:N])
     x[1] = a
     x[end] = b
     return dom.variables => x
 end
+
+@inline function sym_dot(a, b)
+    mapreduce((+), zip(a, b)) do (a_, b_)
+        a_ * b_
+    end
+end
diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl
@@ -3,7 +3,8 @@ using LinearAlgebra
 using SciMLBase
 using DiffEqBase
 using ModelingToolkit
-using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states
+using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states,
+parameters, defaults, varmap_to_vars
 using SymbolicUtils, Symbolics
 using Symbolics: unwrap, solve_for, expand_derivatives, diff2term, setname, rename, similarterm
 using SymbolicUtils: operation, arguments
@@ -60,6 +61,7 @@ include("system_parsing/pde_system_transformation.jl")
 include("discretization/interface_boundary.jl")
 
 # Schemes
+include("discretization/schemes/function_scheme/function_scheme.jl")
 include("discretization/schemes/centered_difference/centered_difference.jl")
 include("discretization/schemes/upwind_difference/upwind_difference.jl")
 include("discretization/schemes/half_offset_centred_difference.jl")
@@ -79,6 +81,6 @@ include("scalar_discretization.jl")
 include("MOL_discretization.jl")
 
 export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace
-export UpwindScheme, WENOScheme
+export UpwindScheme, WENOScheme, FunctionalScheme
 
 end
diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl
@@ -136,41 +136,43 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
         end
         if t === nothing
             uaxes = collect(axes(grid[x])[1] for x in arguments(u))
-            u => collect(first(@variables $sym[uaxes...]))
+            u => unwrap.(collect(first(@variables $sym[uaxes...])))
         elseif isequal(SymbolicUtils.arguments(u), [t])
-            u => fill(u, ()) #Create a 0-dimensional array
+            u => fill(safe_unwrap(u), ()) #Create a 0-dimensional array
         else
             uaxes = collect(axes(grid[x])[1] for x in remove(arguments(u), t))
-            u => collect(first(@variables $sym(t)[uaxes...]))
+            u => unwrap.(collect(first(@variables $sym(t)[uaxes...])))
         end
     end
 
     return DiscreteSpace{nspace,length(depvars),G}(vars, Dict(depvarsdisc), axies, grid, Dict(dxs), Dict(Iaxies), Dict(Igrid))
 end
 
 function Base.getproperty(s::DiscreteSpace, p::Symbol)
-    if p in [:ū, :x̄, :time, :args, :x2i, :i2x]
+    if p in [:ū, :x̄, :ps, :time, :args, :x2i, :i2x]
         getfield(s.vars, p)
     else
         getfield(s, p)
     end
 end
 
+params(s::DiscreteSpace) = s.ps
+
 get_grid_type(::DiscreteSpace{N,M,G}) where {N,M,G} = G
 
 prepare_dx(dx::Integer, xdomain, ::CenterAlignedGrid) = (xdomain[2] - xdomain[1])/(dx - 1)
 prepare_dx(dx::Integer, xdomain, ::EdgeAlignedGrid) = (xdomain[2] - xdomain[1])/dx
 prepare_dx(dx, xdomain, ::AbstractGrid) = dx
 
-nparams(::DiscreteSpace{N,M}) where {N,M} = N
+nivs(::DiscreteSpace{N,M}) where {N,M} = N
 nvars(::DiscreteSpace{N,M}) where {N,M} = M
 
 """
-    params(u, s::DiscreteSpace)
+    ivs(u, s::DiscreteSpace)
 
 Filter out the time variable and get the spatial variables of `u` in `s`.
 """
-params(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time)
+ivs(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time)
 Base.ndims(u, s::DiscreteSpace) = ndims(s.discvars[depvar(u, s)])
 
 Base.length(s::DiscreteSpace, x) = length(s.grid[x])
@@ -185,9 +187,9 @@ of `II` that corresponds to only the spatial arguments of `u`.
 """
 function Idx(II::CartesianIndex, s::DiscreteSpace, u, indexmap)
     # We need to construct a new index as indices may be of different size
-    length(params(u, s)) == 0 && return CartesianIndex()
-    !all(x -> haskey(indexmap, x), params(u, s)) && return II
-    is = [II[indexmap[x]] for x in params(u, s)]
+    length(ivs(u, s)) == 0 && return CartesianIndex()
+    !all(x -> haskey(indexmap, x), ivs(u, s)) && return II
+    is = [II[indexmap[x]] for x in ivs(u, s)]
 
     II = CartesianIndex(is...)
     return II
@@ -198,7 +200,7 @@ A function that returns what to replace independent variables with in boundary e
 """
 @inline function axiesvals(s::DiscreteSpace{N,M,G}, u_, x_, I) where {N,M,G}
     u = depvar(u_, s)
-    map(params(u, s)) do x
+    map(ivs(u, s)) do x
         if isequal(x, x_)
             x => (I[x2i(s, u, x)] == 1 ? first(s.axies[x]) : last(s.axies[x]))
         else
@@ -207,8 +209,8 @@ A function that returns what to replace independent variables with in boundary e
     end
 end
 
-gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in params(u, s)], s.Igrid[u])
-gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in params(u, s)]
+gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in ivs(u, s)], s.Igrid[u])
+gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in ivs(u, s)]
 
 
 varmaps(s::DiscreteSpace, depvars, II, indexmap) = [u => s.discvars[u][Idx(II, s, u, indexmap)] for u in depvars]

diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl
@@ -1,5 +1,5 @@
 @inline function generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary::AbstractTruncatingBoundary)
-    args = params(depvar(boundary.u, s), s)
+    args = ivs(depvar(boundary.u, s), s)
     indexmap = Dict([args[i]=>i for i in 1:length(args)])
     push!(bceqs, generate_bc_eqs(s, boundaryvalfuncs, boundary, interiormap, indexmap))
 end
@@ -48,7 +48,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights
     # replace u(t,0) with u₁, etc
 
     u = depvar(u_, s)
-    args = params(u, s)
+    args = ivs(u, s)
     j = findfirst(isequal(x_), args)
     IIold = II
     # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
@@ -99,7 +99,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights
 
     # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary
     u = depvar(u_, s)
-    args = params(u, s)
+    args = ivs(u, s)
     j = findfirst(isequal(x_), args)
     IIold = II
     # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
@@ -155,7 +155,7 @@ Pads the boundaries with extrapolation equations, extrapolated with 6th order la
 Reuses `central_difference` as this already dispatches the correct stencil, given a `DerivativeOperator` which contains the correct weights.
 """
 function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
-    args = params(u, s)
+    args = ivs(u, s)
     length(args) == 0 && return
 
     lowerextents, upperextents = interiormap.stencil_extents[pde]
@@ -167,7 +167,12 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
     for (j, x) in enumerate(args)
         ninterp = lowerextents[j] - vlower[j]
         I1 = unitindex(length(args), j)
+        bs = bcmap[operation(u)][x]
+        haslower, hasupper = haslowerupper(bs, x)
         while ninterp >= vlower[j]
+            if haslower
+                break
+            end
             for Il in (edge(interiormap, s, u, j, true) .+ (ninterp * I1,))
                 expr = central_difference(derivweights.boundary[x], Il, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc)
                 push!(eqmap[Il], expr)
@@ -176,6 +181,9 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
         end
         ninterp = upperextents[j] - vupper[j]
         while ninterp >= vupper[j]
+            if hasupper
+                break
+            end
             for Iu in (edge(interiormap, s, u, j, false) .- (ninterp * I1,))
                 expr = central_difference(derivweights.boundary[x], Iu, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc)
                 push!(eqmap[Iu], expr)
@@ -201,7 +209,7 @@ end
 
 @inline function generate_corner_eqs!(bceqs, s, interiormap, N, u)
     interior = interiormap.I[interiormap.pde[u]]
-    ndims(u, s) == 0 && return
+    N == 0 || ndims(u, s) == 0 && return
     sd(i, j) = selectdim(interior, j, i)
     domain = setdiff(s.Igrid[u], interior)
     II1 = unitindices(N)

diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl
@@ -54,10 +54,12 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace,
         # Advection rules
         if derivweights.advection_scheme isa UpwindScheme
             advection_rules = generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
-        elseif derivweights.advection_scheme isa WENOScheme
-            advection_rules = generate_WENO_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+        elseif derivweights.advection_scheme isa FunctionalScheme
+            advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s,
+                                                       depvars, derivweights, bmap, indexmap, terms)
             advection_rules = vcat(advection_rules,
-                                generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms; skip = [1]))
+                                generate_winding_rules(II, s, depvars, derivweights, bmap,
+                                                       indexmap, terms; skip = [1]))
         else
             error("Unsupported advection scheme $(derivweights.advection_scheme) encountered.")
         end

diff --git a/src/discretization/generate_ic_defaults.jl b/src/discretization/generate_ic_defaults.jl
@@ -6,7 +6,7 @@ function generate_ic_defaults(tconds, s, ::ScalarizedDiscretization)
                 throw(ArgumentError("Upper boundary condition $(ic.eq) on time variable is not supported, please use a change of variables `t => -τ` to make this an initial condition."))
             end
 
-            args = params(depvar(ic.u, s), s)
+            args = ivs(depvar(ic.u, s), s)
             indexmap = Dict([args[i]=>i for i in 1:length(args)])
             D = ic.order == 0 ? identity : (Differential(t)^ic.order)
             defaultvars = D.(s.discvars[depvar(ic.u, s)])
-Original file line number
+Diff line change
@@ Expand Up / @@ -6,6 +6,6 @@ CurrentModule = MethodOfLines @@
     ```@docs
     unitindex
-    params
+    ivs
     Idx
     ```