From 2f56811960161deae37c5795290d9df21a9afa65 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Sat, 29 Jul 2023 13:13:35 -0400 Subject: [PATCH 1/3] Add Schroedinger Decapode --- examples/quantum/schroedinger.jl | 54 ++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 examples/quantum/schroedinger.jl diff --git a/examples/quantum/schroedinger.jl b/examples/quantum/schroedinger.jl new file mode 100644 index 00000000..9810451b --- /dev/null +++ b/examples/quantum/schroedinger.jl @@ -0,0 +1,54 @@ +using Catlab, CombinatorialSpaces, Decapodes +using LinearAlgebra, MultiScaleArrays, MLStyle, OrdinaryDiffEq +using GeometryBasics: Point2 +Point2D = Point2{Float64} + +Schroedinger = @decapode begin + (i,h,m)::Constant + V::Parameter + Psi::Form0 + + ∂ₜ(Psi) == ((-1 * (h^2)/(2*m))*Δ(Psi) + V * Psi) / (i*h) +end + +infer_types!(Schroedinger, op1_inf_rules_1D, op2_inf_rules_1D) +resolve_overloads!(Schroedinger, op1_res_rules_1D, op2_res_rules_1D) + +function generate(sd, my_symbol; hodge=GeometricHodge()) + op = @match my_symbol begin + :Δ => x -> begin + lap = Δ(0,sd) + lap * x + end + :^ => (x,y) -> x .^ y + x => error("Unmatched operator $my_symbol") + end + return (args...) -> op(args...) +end + +s′ = EmbeddedDeltaSet1D{Bool, Point2D}() +add_vertices!(s′, 30, point=Point2D.(range(-π/2 + π/32, π/2 - π/32, length=30), 0)) +add_edges!(s′, 1:nv(s′)-1, 2:nv(s′)) +orient!(s′) +s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s′) +subdivide_duals!(s, Circumcenter()) + +sim = eval(gensim(Schroedinger, dimension=1)) +fₘ = sim(s, generate) + +Psi = zeros(nv(s)) + +u₀ = construct(PhysicsState, [VectorForm(Psi)], Float64[], [:Psi]) +constants_and_parameters = ( + i = 1, # TODO: Relax assumption of Float64 + V = t -> begin + z = zeros(nv(s)) + z + end, + h = 6.5e-16, # Planck constant in [eV s] + m = 5.49e-4, # mass of electron in [eV] +) + +prob = ODEProblem(fₘ, u₀, (0, 1e-16), constants_and_parameters) +soln = solve(prob, Tsit5()) + From 84f398b57d25b86710b12d837dd2ec13cbfb9395 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Sat, 29 Jul 2023 13:55:28 -0400 Subject: [PATCH 2/3] Assume ComplexF64 --- examples/quantum/schroedinger.jl | 8 ++++---- src/simulation.jl | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/quantum/schroedinger.jl b/examples/quantum/schroedinger.jl index 9810451b..a2db308d 100644 --- a/examples/quantum/schroedinger.jl +++ b/examples/quantum/schroedinger.jl @@ -36,13 +36,13 @@ subdivide_duals!(s, Circumcenter()) sim = eval(gensim(Schroedinger, dimension=1)) fₘ = sim(s, generate) -Psi = zeros(nv(s)) +Psi = zeros(ComplexF64, nv(s)) -u₀ = construct(PhysicsState, [VectorForm(Psi)], Float64[], [:Psi]) +u₀ = construct(PhysicsState, [VectorForm(Psi)], ComplexF64[], [:Psi]) constants_and_parameters = ( - i = 1, # TODO: Relax assumption of Float64 + i = im, # TODO: Relax assumption of Float64 V = t -> begin - z = zeros(nv(s)) + z = zeros(ComplexF64, nv(s)) z end, h = 6.5e-16, # Planck constant in [eV s] diff --git a/src/simulation.jl b/src/simulation.jl index ac151b13..44bbc357 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -98,7 +98,7 @@ Base.Expr(c::AllocVecCall) = begin _ => return :AllocVecCall_Error end - :($(c.name) = Vector{Float64}(undef, nparts(mesh, $(QuoteNode(resolved_form))))) + :($(c.name) = Vector{ComplexF64}(undef, nparts(mesh, $(QuoteNode(resolved_form))))) end #= function get_form_number(d::SummationDecapode, var_id::Int) @@ -225,8 +225,8 @@ function get_vars_code(d::AbstractNamedDecapode, vars::Vector{Symbol}) elseif all(d[incident(d, s, :name) , :type] .== :Parameter) :($s = (p.$s)(t)) elseif all(d[incident(d, s, :name) , :type] .== :Literal) - # TODO: Fix this. We assume that all literals are Float64s. - :($s = $(parse(Float64, String(s)))) + # TODO: Fix this. We assume that all literals are ComplexF64s. + :($s = $(parse(ComplexF64, String(s)))) else # TODO: If names are not unique, then the type is assumed to be a # form for all of the vars sharing a same name. @@ -681,7 +681,7 @@ function closest_point(p1, p2, dims) end end end - Point3{Float64}(p_res...) + Point3{ComplexF64}(p_res...) end function flat_op(s::AbstractDeltaDualComplex2D, X::AbstractVector; dims=[Inf, Inf, Inf]) From 44fd36f4046003038e1621d41d53e6acef1a9848 Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Sat, 29 Jul 2023 18:20:14 -0700 Subject: [PATCH 3/3] Use unicode, add gif, add ICs for Schroedinger --- examples/quantum/schroedinger.jl | 37 ++++++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/examples/quantum/schroedinger.jl b/examples/quantum/schroedinger.jl index a2db308d..a2aaa96d 100644 --- a/examples/quantum/schroedinger.jl +++ b/examples/quantum/schroedinger.jl @@ -1,18 +1,19 @@ using Catlab, CombinatorialSpaces, Decapodes -using LinearAlgebra, MultiScaleArrays, MLStyle, OrdinaryDiffEq +using GLMakie, JLD2, LinearAlgebra, MultiScaleArrays, MLStyle, OrdinaryDiffEq using GeometryBasics: Point2 Point2D = Point2{Float64} Schroedinger = @decapode begin (i,h,m)::Constant V::Parameter - Psi::Form0 + Ψ::Form0 - ∂ₜ(Psi) == ((-1 * (h^2)/(2*m))*Δ(Psi) + V * Psi) / (i*h) + ∂ₜ(Ψ) == ((-1 * (h^2)/(2*m))*Δ(Ψ) + V * Ψ) / (i*h) end infer_types!(Schroedinger, op1_inf_rules_1D, op2_inf_rules_1D) resolve_overloads!(Schroedinger, op1_res_rules_1D, op2_res_rules_1D) +to_graphviz(Schroedinger) function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin @@ -27,7 +28,7 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end s′ = EmbeddedDeltaSet1D{Bool, Point2D}() -add_vertices!(s′, 30, point=Point2D.(range(-π/2 + π/32, π/2 - π/32, length=30), 0)) +add_vertices!(s′, 100, point=Point2D.(range(-1, 1, length=100), 0)) add_edges!(s′, 1:nv(s′)-1, 2:nv(s′)) orient!(s′) s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s′) @@ -36,9 +37,10 @@ subdivide_duals!(s, Circumcenter()) sim = eval(gensim(Schroedinger, dimension=1)) fₘ = sim(s, generate) -Psi = zeros(ComplexF64, nv(s)) +Ψ = zeros(ComplexF64, nv(s)) +Ψ[49] = 1e-16 -u₀ = construct(PhysicsState, [VectorForm(Psi)], ComplexF64[], [:Psi]) +u₀ = construct(PhysicsState, [VectorForm(Ψ)], ComplexF64[], [:Ψ]) constants_and_parameters = ( i = im, # TODO: Relax assumption of Float64 V = t -> begin @@ -49,6 +51,27 @@ constants_and_parameters = ( m = 5.49e-4, # mass of electron in [eV] ) -prob = ODEProblem(fₘ, u₀, (0, 1e-16), constants_and_parameters) +tₑ = 1e12 +prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) soln = solve(prob, Tsit5()) +@save "schroedinger.jld2" soln + +lines(map(x -> x[1], point(s′)), map(x -> x.re, findnode(soln(0.0), :Ψ))) +lines(map(x -> x[1], point(s′)), map(x -> x.re, findnode(soln(tₑ), :Ψ))) + +begin +# Initial frame +frames = 100 +fig = Figure(resolution = (800, 800)) +ax1 = Axis(fig[1,1]) +xlims!(ax1, extrema(map(x -> x[1], point(s′)))) +ylims!(ax1, extrema(map(x -> x.re, findnode(soln(tₑ), :Ψ)))) +Label(fig[1,1,Top()], "Ψ from Schroedinger Wave Equation") +Label(fig[2,1,Top()], "Line plot of real portion of Ψ, every $(tₑ/frames) time units") + +# Animation +record(fig, "schroedinger.gif", range(0.0, tₑ; length=frames); framerate = 15) do t + lines!(fig[1,1], map(x -> x[1], point(s′)), map(x -> x.re, findnode(soln(t), :Ψ))) +end +end