From 90aa7f8107d6766c2aa68b6880f7fc3d2970ae3d Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 7 Nov 2024 19:08:51 +0000 Subject: [PATCH 1/3] fix complex bcs and ics --- docs/src/tutorials/schroedinger.md | 6 +- test/pde_systems/schroedinger.jl | 101 +++++++++++++++++++---------- 2 files changed, 71 insertions(+), 36 deletions(-) diff --git a/docs/src/tutorials/schroedinger.md b/docs/src/tutorials/schroedinger.md index c4232c694..0e834a213 100644 --- a/docs/src/tutorials/schroedinger.md +++ b/docs/src/tutorials/schroedinger.md @@ -18,9 +18,9 @@ V(x) = 0.0 eq = [im * Dt(ψ(t, x)) ~ Dxx(ψ(t, x)) + V(x) * ψ(t, x)] # You must enclose complex equations in a vector, even if there is only one equation -ψ0 = x -> sin(2pi * x) +ψ0 = x -> ((1 + im)/sqrt(2))*sinpi(2*x) -bcs = [ψ(0, x) ~ ψ0(x), +bcs = [ψ(0, x) => ψ0(x), # Initial condition must be marked with a => operator ψ(t, xmin) ~ 0, ψ(t, xmax) ~ 0] @@ -46,4 +46,6 @@ end gif(anim, "schroedinger.gif", fps = 10) ``` +Note that complex initial conditions are supported, but must be marked with a `=>` operator. + This represents the second from ground state of a particle in an infinite quantum well, try changing the potential `V(x)`, initial conditions and BCs, it is extremely interesting to see how the wave function evolves even for nonphysical combinations. Be sure to post interesting results on the discourse! diff --git a/test/pde_systems/schroedinger.jl b/test/pde_systems/schroedinger.jl index 81e53442b..d8e6b9c10 100644 --- a/test/pde_systems/schroedinger.jl +++ b/test/pde_systems/schroedinger.jl @@ -1,55 +1,88 @@ using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test -@parameters t, x -@variables ψ(..) +@testset "Schroedinger" begin + @parameters t, x + @variables ψ(..) -Dt = Differential(t) -Dxx = Differential(x)^2 + Dt = Differential(t) + Dxx = Differential(x)^2 -xmin = 0 -xmax = 1 + xmin = 0 + xmax = 1 -V(x) = 0.0 + V(x) = 0.0 -eq = [im * Dt(ψ(t, x)) ~ (Dxx(ψ(t, x)) + V(x) * ψ(t, x))] # You must enclose complex equations in a vector, even if there is only one equation + eq = [im * Dt(ψ(t, x)) ~ (Dxx(ψ(t, x)) + V(x) * ψ(t, x))] # You must enclose complex equations in a vector, even if there is only one equation -ψ0 = x -> sin(2 * pi * x) + ψ0 = x -> sin(2 * pi * x) -bcs = [ψ(0, x) ~ ψ0(x), - ψ(t, xmin) ~ 0, - ψ(t, xmax) ~ 0] + bcs = [ψ(0, x) ~ ψ0(x), + ψ(t, xmin) ~ 0, + ψ(t, xmax) ~ 0] -domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)] + domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)] -@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t, x)]) + @named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t, x)]) -disc = MOLFiniteDifference([x => 100], t) + disc = MOLFiniteDifference([x => 100], t) -prob = discretize(sys, disc) + prob = discretize(sys, disc) -sol = solve(prob, TRBDF2(), saveat = 0.01) + sol = solve(prob, TRBDF2(), saveat = 0.01) -discx = sol[x] -disct = sol[t] + discx = sol[x] + disct = sol[t] -discψ = sol[ψ(t, x)] + discψ = sol[ψ(t, x)] -analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t) + analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t) -analψ = [analytic(t, x) for t in disct, x in discx] + analψ = [analytic(t, x) for t in disct, x in discx] -for i in 1:length(disct) - u = abs.(analψ[i, :]) .^ 2 - u2 = abs.(discψ[i, :]) .^ 2 + for i in 1:length(disct) + u = abs.(analψ[i, :]) .^ 2 + u2 = abs.(discψ[i, :]) .^ 2 - @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-3 + @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-3 + end + #using Plots + + # anim = @animate for i in 1:length(disct) + # u = analψ[i, :] + # u2 = discψ[i, :] + # plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft) + # end + # gif(anim, "schroedinger.gif", fps = 10) end -#using Plots - -# anim = @animate for i in 1:length(disct) -# u = analψ[i, :] -# u2 = discψ[i, :] -# plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft) -# end -# gif(anim, "schroedinger.gif", fps = 10) +@testset "Schroedinger with complex bcs" begin + @parameters t, x + @variables ψ(..) + + Dt = Differential(t) + Dxx = Differential(x)^2 + + xmin = 0 + xmax = 1 + ϵ = 1e-2 + + V(x) = 0.0 + + eq = [(im*ϵ)*Dt(ψ(t,x)) ~ (-0.5*ϵ^2)Dxx(ψ(t,x)) + V(x)*ψ(t,x)] # You must enclose complex equations in a vector, even if there is only one equation + + ψ0 = x -> exp((im/ϵ)*1e-1*sum(x)) + + bcs = [ψ(0,x) ~ ψ0(x), + ψ(t,xmin) ~ exp((im/ϵ)*(1e-1*sum(xmin) - 0.5e-2*t)), + ψ(t,xmax) ~ exp((im/ϵ)*(1e-1*sum(xmax) - 0.5e-2*t))] + + domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)] + + @named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)]) + + disc = MOLFiniteDifference([x => 300], t) + + prob = discretize(sys, disc) + + sol = solve(prob, TRBDF2(), saveat = 0.01) +end \ No newline at end of file From db9d13015d0297ed7b655e413901940c6173a27f Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 7 Nov 2024 19:20:32 +0000 Subject: [PATCH 2/3] fix tests --- test/pde_systems/schroedinger.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/pde_systems/schroedinger.jl b/test/pde_systems/schroedinger.jl index d8e6b9c10..5787e4e4d 100644 --- a/test/pde_systems/schroedinger.jl +++ b/test/pde_systems/schroedinger.jl @@ -14,9 +14,9 @@ using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test eq = [im * Dt(ψ(t, x)) ~ (Dxx(ψ(t, x)) + V(x) * ψ(t, x))] # You must enclose complex equations in a vector, even if there is only one equation - ψ0 = x -> sin(2 * pi * x) + ψ0 = x -> ((1 + im)/sqrt(2))*sinpi(2*x) - bcs = [ψ(0, x) ~ ψ0(x), + bcs = [ψ(0, x) => ψ0(x), # Initial condition must be marked with a => operator ψ(t, xmin) ~ 0, ψ(t, xmax) ~ 0] @@ -35,7 +35,7 @@ using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test discψ = sol[ψ(t, x)] - analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t) + analytic(t, x) = sqrt(2) * sin(2 * pi * x) * exp(-im * 4 * pi^2 * t) * ((1 + im)/sqrt(2)) analψ = [analytic(t, x) for t in disct, x in discx] @@ -43,7 +43,7 @@ using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test u = abs.(analψ[i, :]) .^ 2 u2 = abs.(discψ[i, :]) .^ 2 - @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-3 + @test u ./ maximum(u)≈u2 ./ maximum(u2) atol=1e-2 end #using Plots @@ -72,7 +72,7 @@ end ψ0 = x -> exp((im/ϵ)*1e-1*sum(x)) - bcs = [ψ(0,x) ~ ψ0(x), + bcs = [ψ(0,x) => ψ0(x), ψ(t,xmin) ~ exp((im/ϵ)*(1e-1*sum(xmin) - 0.5e-2*t)), ψ(t,xmax) ~ exp((im/ϵ)*(1e-1*sum(xmax) - 0.5e-2*t))] @@ -85,4 +85,6 @@ end prob = discretize(sys, disc) sol = solve(prob, TRBDF2(), saveat = 0.01) + + @test SciMLBase.successful_retcode(sol) end \ No newline at end of file From c63092c278bf9931e64076ac7e5e42de7e2ff4fb Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 11 Nov 2024 14:29:19 +0000 Subject: [PATCH 3/3] test on new pdebase --- Project.toml | 2 +- test/pde_systems/schroedinger.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f7561c700..9cd88170c 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ Interpolations = "0.14, 0.15" Latexify = "0.15, 0.16" ModelingToolkit = "9.17" OrdinaryDiffEq = "6" -PDEBase = "0.1.12" +PDEBase = "0.1.17" PrecompileTools = "1" RuntimeGeneratedFunctions = "0.5.9" SafeTestsets = "0.0.1" diff --git a/test/pde_systems/schroedinger.jl b/test/pde_systems/schroedinger.jl index 5787e4e4d..d8573bc3a 100644 --- a/test/pde_systems/schroedinger.jl +++ b/test/pde_systems/schroedinger.jl @@ -87,4 +87,4 @@ end sol = solve(prob, TRBDF2(), saveat = 0.01) @test SciMLBase.successful_retcode(sol) -end \ No newline at end of file +end \ No newline at end of file