Merge pull request #343 from xtalax/complex2

Complex2.1

.github/workflows/CI.yml

@@ -35,6 +35,7 @@ jobs:
           - Brusselator
           - Mixed_Derivatives
           - Wave_Eq_Staggered
+          - Complex
         version:
           - '1'
           - '1.6'

Project.toml

@@ -31,7 +31,7 @@ Interpolations = "0.14"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
-PDEBase = "0.1.6"
+PDEBase = "0.1.8"
 PrecompileTools = "1"
 RuntimeGeneratedFunctions = "0.5.9"
 SafeTestsets = "0.0.1"

docs/pages.jl

@@ -8,6 +8,7 @@ pages = [
         "tutorials/sispde.md",
         "tutorials/icbc_sampled.md",
         "tutorials/PIDE.md",
+        "tutorials/schroedinger.md",
     ],
     "Performance Tips" => "performance.md",
     "MOLFiniteDifference" => "MOLFiniteDifference.md",

docs/src/tutorials/schroedinger.md

@@ -0,0 +1,48 @@
+# Schrödinger Equation
+
+MethodOfLines can solve linear complex PDEs like the Scrödinger equation:
+
+```@example schro
+using MethodOfLines, OrdinaryDiffEq, Plots, DomainSets, ModelingToolkit
+
+@parameters t, x
+@variables ψ(..)
+
+Dt = Differential(t)
+Dxx = Differential(x)^2
+
+xmin = 0
+xmax = 1
+
+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)
+
+bcs = [ψ(0,x) ~ ψ0(x), 
+    ψ(t,xmin) ~ 0,
+    ψ(t,xmax) ~ 0]
+
+domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)]
+
+@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)])
+
+disc = MOLFiniteDifference([x => 100], t)
+
+prob = discretize(sys, disc)
+
+sol = solve(prob, TRBDF2(), saveat = 0.01)
+
+discx = sol[x]
+disct = sol[t]
+
+discψ = sol[ψ(t, x)]
+anim = @animate for i in 1:length(disct)
+    u = 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)
+```
+
+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!

src/interface/solution/timedep.jl

@@ -61,7 +61,6 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T},
         end |> Dict
         # Build Interpolations
         interp = build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys, discretespace.vars.replaced_vars)
-        @show metadata.complexmap
 
         return SciMLBase.PDETimeSeriesSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata),
             typeof(sol),typeof(sol.errors),typeof(sol.t),typeof(ivgrid),

test/pde_systems/schroedinger.jl

@@ -0,0 +1,55 @@
+using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test
+
+@parameters t, x
+@variables ψ(..)
+
+Dt = Differential(t)
+Dxx = Differential(x)^2
+
+xmin = 0
+xmax = 1
+
+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(2*pi*x)
+
+bcs = [ψ(0,x) ~ ψ0(x), 
+    ψ(t,xmin) ~ 0,
+    ψ(t,xmax) ~ 0]
+
+domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)]
+
+@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)])
+
+disc = MOLFiniteDifference([x => 100], t)
+
+prob = discretize(sys, disc)
+
+sol = solve(prob, TRBDF2(), saveat = 0.01)
+
+discx = sol[x]
+disct = sol[t]
+
+discψ = sol[ψ(t, x)]
+
+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]
+
+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
+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)

test/runtests.jl

@@ -8,6 +8,12 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
 # Start Test Script
 
 @time begin
+    if GROUP == "All" || GROUP == "Complex"
+        @time @safetestset "MOLFiniteDifference Interface" begin
+            include("pde_systems/schroedinger.jl")
+        end
+    end
+
     if GROUP == "All" || GROUP == "Mixed_Derivatives"
         @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin
             include("pde_systems/MOL_Mixed_Deriv.jl")