Merge remote-tracking branch 'upstream/master'

Author: valentinsulzer

.github/workflows/CI.yml

@@ -12,7 +12,10 @@ jobs:
     strategy:
       matrix:
         group:
-          - Interface
+          - OperatorInterface
+          - MOLFiniteDifference
+          - Multithreading
+          - Misc
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1

docs/make.jl

@@ -15,6 +15,9 @@ makedocs(
         "Symbolic Tutorials" => [
             "symbolic_tutorials/mol_heat.md"
         ],
+        "Operator Tutorials" => [
+            "operator_tutorials/kdv.md"
+        ],
         "Operators" => [
             "operators/operator_overview.md",
             "operators/derivative_operators.md",
@@ -26,8 +29,8 @@ makedocs(
         ],
         "Nonlinear Derivatives" => [
             "nonlinear_derivatives/nonlinear_diffusion.md"
-        ]
-    ]
+        ]  
+     ]
 )
 
 deploydocs(

docs/src/operator_tutorials/kdv.md

@@ -0,0 +1,50 @@
+# Solving KdV Solitons with Upwinding Operators
+
+The KdV equation is of the form `uₜ + αuuₓ + βuₓₓₓ = 0`. Here we'll use `α = 6`, `β = 1` for 
+simplicity of the true solution expression.
+
+### 1-Soliton solution using Upwind Difference
+
+The analytical expression for the single soliton case takes the form `u(x,t) = (c/2)/cosh²(√c * ξ/2)`.
+
+`c > 0` (wave speed) ; `ξ  = x - c*t` (moving coordinate)
+
+```julia
+using Test
+using DiffEqOperators, OrdinaryDiffEq, LinearAlgebra
+
+# Space domain and grids
+N = 21
+Δx = 1/(N-1)
+c = 1
+x = -10:Δx:10;
+
+# solution of the single forward moving wave
+ϕ(x,t) = (1/2)*sech.((x .- t)/2).^2 
+
+# Discretizing the PDE at t = 0
+u0 = ϕ(x,0);
+du = zeros(size(x)); 
+
+# Declaring the Upwind operator with winding = -1 since the wave travels from left to right 
+A = UpwindDifference{Float64}(1,3,Δx,length(x),-1);
+
+# Defining the ODE problem
+function KdV(du, u, p, t)
+	bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3) 
+	mul!(du,A,bc*u)
+end
+
+single_solition = ODEProblem(KdV, u0, (0.,5.));
+
+# Solving the ODE problem 
+soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
+
+# Plotting the results, comparing Analytical and Numerical solutions 
+using Plots
+plot(ϕ(x,0), title  = "Single forward moving wave", yaxis="u(x,t)", label = "t = 0.0 (Analytic)")
+plot!(ϕ(x,2), label = "t = 2.0 (Analytic)")
+plot!(soln(0.0), label = "t = 0.0 (Numerical)",ls = :dash)
+plot!(soln(2.0), label = "t = 2.0 (Numerical)",ls = :dash)
+```
+![solution_plot](https://user-images.githubusercontent.com/39168576/111033316-a702a180-8436-11eb-9b8e-749ad9206639.jpeg)

docs/src/symbolic_tutorials/mol_heat.md

@@ -25,15 +25,15 @@ bcs = [u(0,x) ~ cos(x),
 
 # Space and time domains
 domains = [t ∈ IntervalDomain(0.0,1.0),
-        x ∈ IntervalDomain(0.0,1.0)]
+           x ∈ IntervalDomain(0.0,1.0)]
 
 # PDE system
-pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
 
 # Method of lines discretization
 dx = 0.1
 order = 2
-discretization = MOLFiniteDifference(dx,order)
+discretization = MOLFiniteDifference([x=>dx],t)
 
 # Convert the PDE problem into an ODE problem
 prob = discretize(pdesys,discretization)
@@ -43,16 +43,18 @@ using OrdinaryDiffEq
 sol = solve(prob,Tsit5(),saveat=0.2)
 
 # Plot results and compare with exact solution
-x = prob.space[2]
+x = (0:dx:1)[2:end-1]
 t = sol.t
 
 using Plots
 plt = plot()
+
 for i in 1:length(t)
-    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    plot!(x,sol.u[i],label="Numerical, t=$(t[i])")
     scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
 end
 display(plt)
+savefig("plot.png")
 ```
 ### Neumann boundary conditions
 
@@ -79,13 +81,13 @@ domains = [t ∈ IntervalDomain(0.0,1.0),
         x ∈ IntervalDomain(0.0,1.0)]
 
 # PDE system
-pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
 
 # Method of lines discretization
 # Need a small dx here for accuracy
 dx = 0.01
 order = 2
-discretization = MOLFiniteDifference(dx,order)
+discretization = MOLFiniteDifference([x=>dx],t)
 
 # Convert the PDE problem into an ODE problem
 prob = discretize(pdesys,discretization)
@@ -95,16 +97,18 @@ using OrdinaryDiffEq
 sol = solve(prob,Tsit5(),saveat=0.2)
 
 # Plot results and compare with exact solution
-x = prob.space[2]
+x = (0:dx:1)[2:end-1]
 t = sol.t
 
 using Plots
 plt = plot()
+
 for i in 1:length(t)
-    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    plot!(x,sol.u[i],label="Numerical, t=$(t[i])",lw=12)
     scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
 end
 display(plt)
+savefig("plot.png")
 ```
 
 ### Robin boundary conditions
@@ -132,13 +136,13 @@ domains = [t ∈ IntervalDomain(0.0,1.0),
         x ∈ IntervalDomain(-1.0,1.0)]
 
 # PDE system
-pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
+pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
 
 # Method of lines discretization
 # Need a small dx here for accuracy
 dx = 0.05
 order = 2
-discretization = MOLFiniteDifference(dx,order)
+discretization = MOLFiniteDifference([x=>dx],t)
 
 # Convert the PDE problem into an ODE problem
 prob = discretize(pdesys,discretization)
@@ -148,14 +152,16 @@ using OrdinaryDiffEq
 sol = solve(prob,Tsit5(),saveat=0.2)
 
 # Plot results and compare with exact solution
-x = prob.space[2]
+x = (0:dx:1)[2:end-1]
 t = sol.t
 
 using Plots
 plt = plot()
+
 for i in 1:length(t)
-    plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]),label="Numerical, t=$(t[i])")
+    plot!(x,sol.u[i],label="Numerical, t=$(t[i])")
     scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")
 end
 display(plt)
+savefig("plot.png")
 ```

src/DiffEqOperators.jl

@@ -42,14 +42,14 @@ include("derivative_operators/coefficient_functions.jl")
 
 ### Composite Operators
 include("composite_operators.jl")
-
-include("MOL_discretization.jl")
-
 include("docstrings.jl")
 
 ### Concretizations
 include("derivative_operators/concretization.jl")
 
+### MOL
+include("MOLFiniteDifference/MOL_discretization.jl")
+
 # The (u,p,t) and (du,u,p,t) interface
 for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition, GhostDerivativeOperator]
   (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
@@ -63,6 +63,7 @@ export AbstractDerivativeOperator, DerivativeOperator,
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC
 export compose, decompose, perpsize
+export discretize, symbolic_discretize
 
 export GhostDerivativeOperator
 export MOLFiniteDifference

src/MOLFiniteDifference/MOL_discretization.jl

@@ -0,0 +1,142 @@
+using ModelingToolkit: operation, istree, arguments
+# Method of lines discretization scheme
+struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization
+    dxs::T
+    time::T2
+    upwind_order::Int
+    centered_order::Int
+end
+
+# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
+MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2) =
+    MOLFiniteDifference(dxs, time, upwind_order, centered_order)
+
+function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
+    pdeeqs = pdesys.eqs isa Vector ? pdesys.eqs : [pdesys.eqs]
+    t = discretization.time
+    nottime = filter(x->x.val != t.val,pdesys.indvars)
+
+    # Discretize space
+
+    space = map(nottime) do x
+        xdomain = pdesys.domain[findfirst(d->x.val == d.variables,pdesys.domain)]
+        @assert xdomain.domain isa IntervalDomain
+        dx = discretization.dxs[findfirst(dxs->x.val == dxs[1].val,discretization.dxs)][2]
+        dx isa Number ? (xdomain.domain.lower:dx:xdomain.domain.upper) : dx
+    end
+    tdomain = pdesys.domain[findfirst(d->t.val == d.variables,pdesys.domain)]
+    @assert tdomain.domain isa IntervalDomain
+    tspan = (tdomain.domain.lower,tdomain.domain.upper)
+
+    # Build symbolic variables
+    indices = CartesianIndices(((axes(s)[1] for s in space)...,))
+    depvars = map(pdesys.depvars) do u
+        [Num(Variable{Symbolics.FnType{Tuple{Any}, Real}}(Base.nameof(ModelingToolkit.operation(u.val)),II.I...))(t) for II in indices]
+    end
+    spacevals = map(y->[Pair(nottime[i],space[i][y.I[i]]) for i in 1:length(nottime)],indices)
+
+    # Build symbolic maps
+    edges = reduce(vcat,[[vcat([Colon() for j in 1:i-1],1,[Colon() for j in i+1:length(nottime)]),
+      vcat([Colon() for j in 1:i-1],size(depvars[1],i),[Colon() for j in i+1:length(nottime)])] for i in 1:length(nottime)])
+
+    #edgeindices = [indices[e...] for e in edges]
+    edgevals = reduce(vcat,[[nottime[i]=>first(space[i]),nottime[i]=>last(space[i])] for i in 1:length(space)])
+    edgevars = [[d[e...] for e in edges] for d in depvars]
+    edgemaps = [spacevals[e...] for e in edges]
+    initmaps = substitute.(pdesys.depvars,[t=>tspan[1]])
+
+    depvarmaps = reduce(vcat,[substitute.((pdesys.depvars[i],),edgevals) .=> edgevars[i] for i in 1:length(pdesys.depvars)])
+    if length(nottime) == 1
+        left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][1],space[j][2]])
+        right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][end-1],space[j][end]])
+        central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:length(nottime)]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:length(nottime)]...)]
+        derivars = [[dot(left_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][2]],depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][1]]]),
+                    dot(right_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][1]],depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][2]]])]
+                    for j in 1:length(pdesys.depvars)]
+        depvarderivmaps = reduce(vcat,[substitute.((Differential(nottime[j])(pdesys.depvars[i]),),edgevals) .=> derivars[i]
+                                       for i in 1:length(pdesys.depvars) for j in 1:length(nottime)])
+   else
+       # TODO: Fix Neumann and Robin on higher dimension
+       depvarderivmaps = []
+   end
+
+    # Generate initial conditions and bc equations
+    u0 = []
+    bceqs = []
+    for bc in pdesys.bcs
+        if ModelingToolkit.operation(bc.lhs) isa Sym && t.val ∉ ModelingToolkit.arguments(bc.lhs)
+            # initial condition
+            # Assume in the form `u(...) ~ ...` for now
+            push!(u0,vec(depvars[findfirst(isequal(bc.lhs),initmaps)] .=> substitute.((bc.rhs,),spacevals)))
+        else
+            # Algebraic equations for BCs
+            i = findfirst(x->occursin(x.val,bc.lhs),first.(depvarmaps))
+
+            # TODO: Fix Neumann and Robin on higher dimension
+            lhs = length(nottime) == 1 ? substitute(bc.lhs,depvarderivmaps[i]) : bc.lhs
+
+            lhs = substitute(lhs,depvarmaps[i])
+            rhs = substitute.((bc.rhs,),edgemaps[i])
+            lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
+            push!(bceqs,lhs .~ rhs)
+        end
+    end
+    u0 = reduce(vcat,u0)
+    bceqs = reduce(vcat,bceqs)
+
+    # Generate PDE Equations
+    interior = indices[[2:length(s)-1 for s in space]...]
+    eqs = vec(map(Base.product(interior,pdeeqs)) do p
+        i,eq = p
+
+        # TODO: Number of points in the central_neighbor_idxs should be dependent
+        # on discretization.centered_order
+        # TODO: Generalize central difference handling to allow higher even order derivatives
+        central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:length(nottime)]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:length(nottime)]...)]
+        central_weights(i,j) = DiffEqOperators.calculate_weights(2, 0.0, [space[j][i[j]-1],space[j][i[j]],space[j][i[j]+1]])
+        central_deriv_rules = [(Differential(nottime[j])^2)(pdesys.depvars[k]) => dot(central_weights(i,j),depvars[k][central_neighbor_idxs(i,j)]) for j in 1:length(nottime), k in 1:length(pdesys.depvars)]
+        valrules = vcat([pdesys.depvars[k] => depvars[k][i] for k in 1:length(pdesys.depvars)],
+                        [nottime[k] => space[k][i[k]] for k in 1:length(nottime)])
+
+        # TODO: Use rule matching for nonlinear Laplacian
+
+        # TODO: upwind rules needs interpolation into `@rule`
+        #forward_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]],space[j][i[j]+1]])
+        #reverse_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]-1],space[j][i[j]]])
+        #upwinding_rules = [@rule(*(~~a,(Differential(nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
+        #                        *(~~a..., ~~b..., dot(reverse_weights(i,j),depvars[k][central_neighbor_idxs(i,j)[1:2]])),
+        #                        *(~~a..., ~~b..., dot(forward_weights(i,j),depvars[k][central_neighbor_idxs(i,j)[2:3]]))))
+        #                        for j in 1:length(nottime), k in 1:length(pdesys.depvars)]
+
+        substitute(eq.lhs,vcat(vec(central_deriv_rules),valrules)) ~ substitute(eq.rhs,vcat(vec(central_deriv_rules),valrules))
+    end)
+
+    # Finalize
+    defaults = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : vcat(u0,pdesys.ps)
+    ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : first.(pdesys.ps)
+    sys = ODESystem(vcat(eqs,unique(bceqs)),t,vec(reduce(vcat,vec(depvars))),ps,defaults=Dict(defaults))
+    sys, tspan
+end
+
+function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
+    sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
+    simpsys = structural_simplify(sys)
+    prob = ODEProblem(simpsys,Pair[],tspan)
+end
+
+# Piracy, to be deleted when https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/251
+# merges
+Base.occursin(needle::ModelingToolkit.SymbolicUtils.Symbolic, haystack::ModelingToolkit.SymbolicUtils.Symbolic) = _occursin(needle, haystack)
+Base.occursin(needle, haystack::ModelingToolkit.SymbolicUtils.Symbolic) = _occursin(needle, haystack)
+Base.occursin(needle::ModelingToolkit.SymbolicUtils.Symbolic, haystack) = _occursin(needle, haystack)
+function _occursin(needle, haystack)
+    isequal(needle, haystack) && return true
+
+    if istree(haystack)
+        args = arguments(haystack)
+        for arg in args
+            occursin(needle, arg) && return true
+        end
+    end
+    return false
+end