Merge pull request #15 from YingboMa/master

WIP: Add more stiff ODE problems

Author: ChrisRackauckas

REQUIRE

@@ -2,3 +2,4 @@ julia 0.6
 ParameterizedFunctions 2.0.0
 DiffEqBase 3.0.3
 DiffEqPDEBase
+DiffEqOperators

src/DiffEqProblemLibrary.jl

@@ -2,9 +2,9 @@ __precompile__()
 
 module DiffEqProblemLibrary
 
-using DiffEqBase, ParameterizedFunctions, DiffEqPDEBase
+using DiffEqBase, ParameterizedFunctions, DiffEqPDEBase, DiffEqOperators
 
-include("ode_premade_problems.jl")
+include("ode/ode_premade_problems.jl")
 include("dae_premade_problems.jl")
 include("dde_premade_problems.jl")
 include("sde_premade_problems.jl")
@@ -13,7 +13,10 @@ include("sde_premade_problems.jl")
 export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
        prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
        prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
-       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades
+       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
+       prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
+       prob_ode_hires, prob_ode_pollution, prob_ode_filament,
+       SolverDiffEq, filament_prob
 
  #SDE Example Problems
  export prob_sde_wave, prob_sde_linear, prob_sde_cubic, prob_sde_2Dlinear, prob_sde_lorenz,

src/ode/brusselator_prob.jl

@@ -0,0 +1,162 @@
+"""
+2D Brusselator
+
+```math
+\\begin{align}
+\\frac{\\partial u}{\\partial t} &= 1 + u^2v - 4.4u + \\alpha(\frac{\\partial^2 u}{\\partial x^2} + \frac{\\partial^2 u}{\\partial y^2}) + f(x, y, t)
+\\frac{\\partial v}{\\partial t} &= 3.4u - u^2v + \\alpha(\frac{\\partial^2 u}{\\partial x^2} + \frac{\\partial^2 u}{\\partial y^2})
+\\end{align}
+```
+
+where
+
+```math
+f(x, y, t) = \\begin{cases}
+5 & \\quad \\text{if } (x-0.3)^2+(y-0.6)^2 ≤ 0.1^2 \\text{ and } t ≥ 1.1 \\\\
+0 & \\quad \\text{else}
+\\end{cases}
+```
+
+and the initial conditions are
+
+```math
+\\begin{align}
+u(x, y, 0) &= 22\\cdot y(1-y)^{3/2} \\\\
+v(x, y, 0) &= 27\\cdot x(1-x)^{3/2}
+\\end{align}
+```
+
+with the periodic boundary condition
+
+```math
+\\begin{align}
+u(x+1,y,t) &= u(x,y,t) \\\\
+u(x,y+1,t) &= u(x,y,t)
+\\end{align}
+```
+
+From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 152
+"""
+brusselator_f(x, y, t) = ifelse((((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
+                                (t >= 1.1), 5., 0.)
+function limit(a, N)
+  if a == N+1
+    return 1
+  elseif a == 0
+    return N
+  else
+    return a
+  end
+end
+function brusselator_2d_loop(du, u, p, t)
+  @inbounds begin
+    A, B, α, xyd, dx, N = p
+    α = α/dx^2
+    for I in CartesianRange((N, N))
+      x = xyd[I[1]]
+      y = xyd[I[2]]
+      i = I[1]
+      j = I[2]
+      ip1 = limit(i+1, N)
+      im1 = limit(i-1, N)
+      jp1 = limit(j+1, N)
+      jm1 = limit(j-1, N)
+      du[i,j,1] = α*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
+      B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
+    end
+    for I in CartesianRange((N, N))
+      i = I[1]
+      j = I[2]
+      ip1 = limit(i+1, N)
+      im1 = limit(i-1, N)
+      jp1 = limit(j+1, N)
+      jm1 = limit(j-1, N)
+      du[i,j,2] = α*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
+      A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
+    end
+  end
+end
+function init_brusselator_2d(xyd)
+  N = length(xyd)
+  u = zeros(N, N, 2)
+  for I in CartesianRange((N, N))
+    x = xyd[I[1]]
+    y = xyd[I[2]]
+    u[I,1] = 22*(y*(1-y))^(3/2)
+    u[I,2] = 27*(x*(1-x))^(3/2)
+  end
+  u
+end
+xyd_brusselator = linspace(0,1,32)
+prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
+                                     init_brusselator_2d(xyd_brusselator),
+                                     (0.,11.5),
+                                     (3.4, 1., 10.,
+                                      xyd_brusselator, step(xyd_brusselator),
+                                      length(xyd_brusselator)))
+
+"""
+1D Brusselator
+
+```math
+\\begin{align}
+\\frac{\\partial u}{\\partial t} &= A + u^2v - (B+1)u + \\alpha\frac{\\partial^2 u}{\\partial x^2}
+\\frac{\\partial v}{\\partial t} &= Bu - u^2v + \\alpha\frac{\\partial^2 u}{\\partial x^2}
+\\end{align}
+```
+
+and the initial conditions are
+
+```math
+\\begin{align}
+u(x,0) &= 1+\\sin(2π x) \\\\
+v(x,0) &= 3
+\\end{align}
+```
+
+with the boundary condition
+
+```math
+\\begin{align}
+u(0,t) &= u(1,t) = 1 \\\\
+v(0,t) &= v(1,t) = 3
+\\end{align}
+```
+
+From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 6
+"""
+const N_brusselator_1d = 40
+const D_brusselator_u = DerivativeOperator{Float64}(2,2,1/(N_brusselator_1d-1),
+                                                    N_brusselator_1d,
+                                                    :Dirichlet,:Dirichlet;
+                                                    BC=(1.,1.))
+const D_brusselator_v = DerivativeOperator{Float64}(2,2,1/(N_brusselator_1d-1),
+                                                    N_brusselator_1d,
+                                                    :Dirichlet,:Dirichlet;
+                                                    BC=(3.,3.))
+function brusselator_1d(du, u_, p, t)
+    A, B, α, buffer = p
+    u = @view(u_[:, 1])
+    v = @view(u_[:, 2])
+    A_mul_B!(buffer, D_brusselator_u, u)
+    Du = buffer
+    @. du[:, 1] = A + u^2*v - (B+1)*u + α*Du
+
+    A_mul_B!(buffer, D_brusselator_v, v)
+    Dv = buffer
+    @. du[:, 2] = B*u - u^2*v + α*Dv
+    nothing
+end
+function init_brusselator_1d(N)
+  u = zeros(N, 2)
+  x = linspace(0, 1, N)
+  for i in 1:N
+    u[i, 1] = 1 + sin(2pi*x[i])
+    u[i, 2] = 3.
+  end
+  u
+end
+prob_ode_brusselator_1d = ODEProblem(brusselator_1d,
+                                    init_brusselator_1d(N_brusselator_1d),
+                                    (0.,10.),
+                                    (1., 3., 1/50, zeros(N_brusselator_1d)))