Merge pull request #507 from SciML/upstream_duals

Make dual types concrete and strip out the rest of MOL

Author: ChrisRackauckas

Project.toml

@@ -13,7 +13,6 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -23,7 +22,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
-SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
 BandedMatrices = "0.15.11, 0.16"
@@ -34,19 +32,18 @@ ForwardDiff = "0.10"
 LazyArrays = "0.17, 0.18, 0.19, 0.20, 0.21, 0.22"
 LazyBandedMatrices = "0.5, 0.6, 0.7"
 LoopVectorization = "0.12"
-ModelingToolkit = "7, 8"
 NNlib = "0.6, 0.7"
 NonlinearSolve = "0.3.7"
 Requires = "1"
 RuntimeGeneratedFunctions = "0.4, 0.5"
 SciMLBase = "1.11"
 SparseDiffTools = "1.17"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18"
 julia = "1.6"
 
 [extras]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -56,4 +53,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["FillArrays", "OrdinaryDiffEq", "Parameters", "Random", "SafeTestsets", "SpecialFunctions", "Test", "Zygote"]
+test = ["FillArrays", "LinearSolve", "OrdinaryDiffEq", "Parameters", "Random", "SafeTestsets", "SpecialFunctions", "Test", "Zygote"]

docs/src/index.md

@@ -1,12 +1,14 @@
 # DiffEqOperators.jl
 
 DiffEqOperators.jl is a package for finite difference discretization of partial
-differential equations. It serves two purposes:
+differential equations. It is for building fast lazy operators for high order
+non-uniform finite differences.
 
-1. Building fast lazy operators for high order non-uniform finite differences.
-2. Automated finite difference discretization of symbolically-defined PDEs.
+!!! note
+
+    For automated finite difference discretization of symbolically-defined PDEs,
+    see [MethodOfLines.jl](https://github.com/SciML/MethodOfLines.jl).
 
-#### Note: (2) is still a work in progress!
 
 For the operators, both centered and
 [upwind](https://en.wikipedia.org/wiki/Upwind_scheme) operators are provided,

docs/src/symbolic/molfinitedifference.md

@@ -9,7 +9,7 @@ using SciMLBase: AbstractDiffEqLinearOperator, AbstractDiffEqCompositeOperator,
 import SciMLBase: getops
 using SparseDiffTools
 using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices, LoopVectorization
-using LazyBandedMatrices, ModelingToolkit
+using LazyBandedMatrices
 using RuntimeGeneratedFunctions
 using Requires
 RuntimeGeneratedFunctions.init(@__MODULE__)
@@ -56,9 +56,6 @@ 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)
@@ -82,9 +79,7 @@ export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, Mul
        MultiDimDirectionalBC, ComposedMultiDimBC
 export compose, decompose, perpsize, square_norm, square_norm!, dot_product, dot_product!, cross_product,
        cross_product!
-export discretize, symbolic_discretize
 
 export GhostDerivativeOperator
-export MOLFiniteDifference, center_align, edge_align
 export BoundaryConditionError
 end # module

src/DiffEqOperators.jl

@@ -42,8 +42,8 @@ mutable struct JacVecOperator{T,F,T1,T2,uType,P,tType,O} <:
         opnorm = true,
     ) where {T}
         if autodiff
-            cache1 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(u, u)
-            cache2 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(u, u)
+            cache1 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(u))),eltype(u),1}.(u, ForwardDiff.Partials.(tuple.(u)))
+            cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(u))),eltype(u),1}.(u, ForwardDiff.Partials.(tuple.(u)))
         else
             cache1 = similar(u)
             cache2 = similar(u)

src/MOLFiniteDifference/MOL_discretization.jl

@@ -1,5 +1,3 @@
-using SymbolicUtils
-
 """
 A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension.
 """
@@ -22,61 +20,6 @@ function cartesian_to_linear(I::CartesianIndex, s)  #Not sure if there is a buil
     return out
 end
 
-# Counts the Differential operators for given variable x. This is used to determine
-# the order of a PDE.
-function count_differentials(term, x::Symbolics.Symbolic)
-    S = Symbolics
-    SU = SymbolicUtils
-    if !S.istree(term)
-        return 0
-    else
-        op = SU.operation(term)
-        count_children = sum(map(arg -> count_differentials(arg, x), SU.arguments(term)))
-        if op isa Differential && op.x === x
-            return 1 + count_children
-        end
-        return count_children
-    end
-end
-
-# return list of differential orders in the equation
-function differential_order(eq, x::Symbolics.Symbolic)
-    S = Symbolics
-    SU = SymbolicUtils
-    orders = Set{Int}()
-    if S.istree(eq)
-        op = SU.operation(eq)
-        if op isa Differential
-            push!(orders, count_differentials(eq, x))
-        else
-            for o in map(ch -> differential_order(ch, x), SU.arguments(eq))
-                union!(orders, o)
-            end
-        end
-    end
-    return filter(!iszero, orders)
-end
-
-# find all the dependent variables given by depvar_ops in an expression
-function get_depvars(eq,depvar_ops)
-    S = Symbolics
-    SU = SymbolicUtils
-    depvars = Set()
-    if eq isa Num
-       eq = eq.val
-    end
-    if S.istree(eq)
-        if eq isa Term && any(u->isequal(operation(eq),u),depvar_ops)
-              push!(depvars, eq)
-        else
-            for o in map(x->get_depvars(x,depvar_ops), SU.arguments(eq))
-                union!(depvars, o)
-            end
-        end
-    end
-    depvars
-end
-
 add_dims(A::AbstractArray, n::Int; dims::Int = 1) = cat(ndims(A) + n, A, dims = dims)
 
 ""

src/jacvec_operators.jl

@@ -1,14 +1,14 @@
 using DiffEqBase,
-    DiffEqOperators, ForwardDiff, LinearAlgebra, SparseDiffTools, Test
+    DiffEqOperators, ForwardDiff, LinearAlgebra, SparseDiffTools, Test, LinearSolve
 const A = rand(300, 300)
 f(du, u) = mul!(du, A, u)
 f(u) = A * u
 x = rand(300)
 v = rand(300)
 du = similar(x)
 
-cache1 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(x, v)
-cache2 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(x, v)
+cache1 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
+cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
 @test num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f, similar(x), x) * v rtol =
     1e-6
 @test num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈
@@ -81,13 +81,11 @@ ff2 = ODEFunction(
     jac_prototype = JacVecOperator{Float64}(lorenz, u0, autodiff = false),
 )
 
-
 for ff in [ff1, ff2]
     prob = ODEProblem(ff, u0, tspan)
     @test solve(prob, TRBDF2()).retcode == :Success
-    @test solve(prob, TRBDF2(linsolve = LinSolveGMRES())).retcode == :Success
+    @test solve(prob, TRBDF2(linsolve = KrylovJL_GMRES())).retcode == :Success
     @test solve(prob, Exprb32()).retcode == :Success
-    @test_broken sol = solve(prob, Rosenbrock23())
-    @test_broken sol =
-        solve(prob, Rosenbrock23(linsolve = LinSolveGMRES(tol = 1e-10)))
+    @test sol = solve(prob, Rosenbrock23()).retcode == :Success
+    @test sol = solve(prob, Rosenbrock23(linsolve = KrylovJL_GMRES())).retcode == :Success
 end

src/utils.jl

@@ -1,83 +1,9 @@
 using Test, LinearAlgebra
 using DiffEqOperators
-using ModelingToolkit
 
 @testset "utility functions" begin
     @test DiffEqOperators.unit_indices(2) == (CartesianIndex(1,0), CartesianIndex(0,1))
     @test DiffEqOperators.add_dims(zeros(2,2), ndims(zeros(2,2)) + 2) == [6. 6.; 0. 0.; 0. 0.]
     @test DiffEqOperators.perpindex(collect(1:5), 3) == [1, 2, 4, 5]
     @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
 end
-
-@testset "count differentials 1D" begin
-    @parameters t x
-    @variables u(..)
-    Dt = Differential(t)
-
-    Dx = Differential(x)
-    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
-    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 1
-    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
-    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
-    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
-
-    Dxx = Differential(x)^2
-    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
-    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
-    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
-    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
-    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
-
-    Dxxxx = Differential(x)^4
-    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
-    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 4
-    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
-    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
-    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
-end
-
-@testset "count differentials 2D" begin
-    @parameters t x y
-    @variables u(..)
-    Dxx = Differential(x)^2
-    Dyy = Differential(y)^2
-    Dt = Differential(t)
-
-    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
-    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
-    @test first(DiffEqOperators.differential_order(eq.rhs, y.val)) == 2
-    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
-    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
-    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
-    @test isempty(DiffEqOperators.differential_order(eq.lhs, y.val))
-end
-
-@testset "count with mixed terms" begin
-    @parameters t x y
-    @variables u(..)
-    Dxx = Differential(x)^2
-    Dyy = Differential(y)^2
-    Dx = Differential(x)
-    Dy = Differential(y)
-    Dt = Differential(t)
-
-    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
-    @test DiffEqOperators.differential_order(eq.rhs, x.val) == Set([2, 1])
-    @test DiffEqOperators.differential_order(eq.rhs, y.val) == Set([2, 1])
-end
-
-@testset "Kuramoto–Sivashinsky equation" begin
-    @parameters x, t
-    @variables u(..)
-    Dt = Differential(t)
-    Dx = Differential(x)
-    Dx2 = Differential(x)^2
-    Dx3 = Differential(x)^3
-    Dx4 = Differential(x)^4
-
-    α = 1
-    β = 4
-    γ = 1
-    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
-    @test DiffEqOperators.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
-end