From 4284d6192099ba94258febca6357d2eaca9e4686 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Feb 2023 17:44:43 +0000 Subject: [PATCH 1/5] no deo bruss (was wrong anyway) --- lib/ODEProblemLibrary/src/brusselator_prob.jl | 44 ++++++++++--------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/lib/ODEProblemLibrary/src/brusselator_prob.jl b/lib/ODEProblemLibrary/src/brusselator_prob.jl index 4af847b..135d7c5 100644 --- a/lib/ODEProblemLibrary/src/brusselator_prob.jl +++ b/lib/ODEProblemLibrary/src/brusselator_prob.jl @@ -101,27 +101,19 @@ const D_brusselator_u = CenteredDifference{Float64}(2, 2, 1 / (N_brusselator_1d 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 -end -function init_brusselator_1d(N) - u = zeros(N, 2) - x = range(0, stop = 1, length = N) - for i in 1:N - u[i, 1] = 1 + sin(2pi * x[i]) - u[i, 2] = 3.0 +const N = 32 + +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 - u end """ @@ -157,4 +149,14 @@ From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff an prob_ode_brusselator_1d = ODEProblem(brusselator_1d, 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))) + +function init_brusselator_1d(N) + u = zeros(N, 2) + x = range(0, stop = 1, length = N) + for i in 1:N + u[i, 1] = 1 + sin(2pi * x[i]) + u[i, 2] = 3.0 + end + u +end From d574424dafbe5e269d2ad03bafb473afc920e474 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Feb 2023 17:46:27 +0000 Subject: [PATCH 2/5] a fix --- lib/ODEProblemLibrary/src/brusselator_prob.jl | 31 ++++++++----------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/lib/ODEProblemLibrary/src/brusselator_prob.jl b/lib/ODEProblemLibrary/src/brusselator_prob.jl index 135d7c5..02cbeb5 100644 --- a/lib/ODEProblemLibrary/src/brusselator_prob.jl +++ b/lib/ODEProblemLibrary/src/brusselator_prob.jl @@ -97,25 +97,30 @@ 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) -const N = 32 function brusselator_1d_loop(du, u, p, t) A, B, alpha, dx = p alpha = alpha / dx^2 - @inbounds for i in 2:N-1 + @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] + 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] + 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) + for i in 1:N + u[i, 1] = 1 + sin(2pi * x[i]) + u[i, 2] = 3.0 + end + u +end + """ 1D Brusselator @@ -150,13 +155,3 @@ prob_ode_brusselator_1d = ODEProblem(brusselator_1d, init_brusselator_1d(N_brusselator_1d), (0.0, 10.0), (1.0, 3.0, 1 / 41, zeros(N_brusselator_1d))) - -function init_brusselator_1d(N) - u = zeros(N, 2) - x = range(0, stop = 1, length = N) - for i in 1:N - u[i, 1] = 1 + sin(2pi * x[i]) - u[i, 2] = 3.0 - end - u -end From fd59884dfe4955b2066c8894f61e39bc106ecbc9 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Feb 2023 17:47:09 +0000 Subject: [PATCH 3/5] remove deo --- lib/ODEProblemLibrary/Project.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/ODEProblemLibrary/Project.toml b/lib/ODEProblemLibrary/Project.toml index 48f3ba2..84a7a68 100644 --- a/lib/ODEProblemLibrary/Project.toml +++ b/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" From e309f9a0b2ee14bc3fe1c03d463a59d95c80e5d9 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Feb 2023 18:05:38 +0000 Subject: [PATCH 4/5] remove deo --- lib/ODEProblemLibrary/src/ODEProblemLibrary.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl index a5c814d..c454769 100644 --- a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl +++ b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl @@ -3,7 +3,6 @@ __precompile__(false) module ODEProblemLibrary using DiffEqBase -using DiffEqOperators using Latexify using ModelingToolkit From 628268b93a2d7086d6c3f8ad714d9231835f2ffa Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Feb 2023 18:12:36 +0000 Subject: [PATCH 5/5] fix --- lib/ODEProblemLibrary/src/brusselator_prob.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ODEProblemLibrary/src/brusselator_prob.jl b/lib/ODEProblemLibrary/src/brusselator_prob.jl index 02cbeb5..a0d3ca2 100644 --- a/lib/ODEProblemLibrary/src/brusselator_prob.jl +++ b/lib/ODEProblemLibrary/src/brusselator_prob.jl @@ -151,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 / 41, zeros(N_brusselator_1d)))