Merge pull request #103 from xtalax/bruss

no deo bruss (was wrong anyway)

Author: ChrisRackauckas

lib/ODEProblemLibrary/Project.toml

@@ -4,7 +4,6 @@ version = "0.1.4"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
@@ -15,7 +14,6 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 [compat]
 Aqua = "0.5"
 DiffEqBase = "6"
-DiffEqOperators = "4"
 Latexify = "0.15"
 ModelingToolkit = "7,8"
 RuntimeGeneratedFunctions = "0.5"

lib/ODEProblemLibrary/src/ODEProblemLibrary.jl

@@ -3,7 +3,6 @@ __precompile__(false)
 module ODEProblemLibrary
 
 using DiffEqBase
-using DiffEqOperators
 using Latexify
 using ModelingToolkit
 

lib/ODEProblemLibrary/src/brusselator_prob.jl

@@ -97,23 +97,20 @@ prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                       length(xyd_brusselator)))
 
 const N_brusselator_1d = 40
-const D_brusselator_u = CenteredDifference{Float64}(2, 2, 1 / (N_brusselator_1d - 1),
-                                                    N_brusselator_1d)
-const D_brusselator_v = CenteredDifference{Float64}(2, 2, 1 / (N_brusselator_1d - 1),
-                                                    N_brusselator_1d)
-function brusselator_1d(du, u_, p, t)
-    A, B, α, buffer = p
-    u = @view(u_[:, 1])
-    v = @view(u_[:, 2])
-    mul!(buffer, D_brusselator_u, u)
-    Du = buffer
-    @. du[:, 1] = A + u^2 * v - (B + 1) * u + α * Du
-
-    mul!(buffer, D_brusselator_v, v)
-    Dv = buffer
-    @. du[:, 2] = B * u - u^2 * v + α * Dv
-    nothing
+
+function brusselator_1d_loop(du, u, p, t)
+    A, B, alpha, dx = p
+    alpha = alpha / dx^2
+    @inbounds for i in 2:(N - 1)
+        x = xyd_brusselator[i]
+        ip1, im1 = i + 1, i - 1
+        du[i, 1] = alpha * (u[im1, 1] + u[ip1, 1] - 2u[i, 1]) +
+                   A + u[i, 1]^2 * u[i, 2] - (B + 1) * u[i, 1]
+        du[i, 2] = alpha * (u[im1, 2] + u[ip1, 2] - 2u[i, 2]) +
+                   B * u[i, 1] - u[i, 1]^2 * u[i, 2]
+    end
 end
+
 function init_brusselator_1d(N)
     u = zeros(N, 2)
     x = range(0, stop = 1, length = N)
@@ -154,7 +151,7 @@ v(0,t) = v(1,t) = 3
 
 From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 6
 """
-prob_ode_brusselator_1d = ODEProblem(brusselator_1d,
+prob_ode_brusselator_1d = ODEProblem(brusselator_1d_loop,
                                      init_brusselator_1d(N_brusselator_1d),
                                      (0.0, 10.0),
-                                     (1.0, 3.0, 1 / 50, zeros(N_brusselator_1d)))
+                                     (1.0, 3.0, 1 / 41, zeros(N_brusselator_1d)))