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"
@@ -34,7 +33,6 @@ 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"

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(x))),eltype(x),1}.(u, u)
+            cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(u, u)
         else
             cache1 = similar(u)
             cache2 = similar(u)

src/MOLFiniteDifference/MOL_discretization.jl

@@ -7,8 +7,8 @@ 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, v)
+cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(x, v)
 @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