From 9c83e3e56b7900d4b7976a4e5068a66434f8368a Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Thu, 19 Dec 2024 17:07:19 -0500 Subject: [PATCH] Notebooks -> Scripts --- EC_test.ipynb | 100 ----------------------------------- calculations/R2R_Berggren.jl | 48 +++++++++++++++++ p_space_test.ipynb | 62 ---------------------- simple_berggren.ipynb | 70 ------------------------ test/p_space_basic.jl | 21 ++++++++ test/simple_berggren.jl | 30 +++++++++++ 6 files changed, 99 insertions(+), 232 deletions(-) delete mode 100644 EC_test.ipynb create mode 100644 calculations/R2R_Berggren.jl delete mode 100644 p_space_test.ipynb delete mode 100644 simple_berggren.ipynb create mode 100644 test/p_space_basic.jl create mode 100644 test/simple_berggren.jl diff --git a/EC_test.ipynb b/EC_test.ipynb deleted file mode 100644 index 38087bd..0000000 --- a/EC_test.ipynb +++ /dev/null @@ -1,100 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Plots\n", - "include(\"p_space.jl\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "vertices = [0, 0.4 - 0.15im, 0.8, 6]\n", - "subdivisions = [128, 128, 128]\n", - "p, w = get_mesh(vertices, subdivisions)\n", - "mesh_E = p.*p ./ (2*0.5)\n", - "\n", - "# ResonanceEC: Eq. (20)\n", - "V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "training_points = range(0.75, 0.45, 5)\n", - "training_E = Vector{ComplexF64}(undef, length(training_points))\n", - "EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))\n", - "\n", - "for (j, c) in enumerate(training_points)\n", - " evals, evecs = eigen(get_H_matrix(V_system(c), p, w))\n", - " i = identify_pole_i(p, evals)\n", - " training_E[j] = evals[i]\n", - " EC_basis[:, j] = evecs[:, i]\n", - "end\n", - "\n", - "scatter(real.(training_E), imag.(training_E), label=\"training\")\n", - "plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n", - "xlims!(0,1)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "extrapolate_points = range(0.40, 0.25, 5)\n", - "ref_E = 0.2 - 0.1im\n", - "\n", - "exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n", - "extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n", - "\n", - "EC_basis_w = EC_basis .* w\n", - "N_EC = transpose(EC_basis_w) * EC_basis\n", - "\n", - "for (j, c) in enumerate(extrapolate_points)\n", - " exact_E[j] = quick_pole_E(V_system(c))\n", - "\n", - " EC_basis_w = EC_basis .* w\n", - " H = get_H_matrix(V_system(c), p, w)\n", - " H_EC = transpose(EC_basis_w) * H * EC_basis\n", - " evals = eigvals(H_EC, N_EC)\n", - " i = argmin(abs.(evals .- ref_E))\n", - " ref_E = evals[i]\n", - " extrapolate_E[j] = evals[i]\n", - "end\n", - "\n", - "scatter(real.(training_E), imag.(training_E), label=\"training\")\n", - "scatter!(real.(exact_E), imag.(exact_E), label=\"exact\")\n", - "scatter!(real.(extrapolate_E), imag.(extrapolate_E), label=\"extrapolated\")\n", - "plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n", - "xlims!(0,1)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.10.2", - "language": "julia", - "name": "julia-1.10" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.10.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/calculations/R2R_Berggren.jl b/calculations/R2R_Berggren.jl new file mode 100644 index 0000000..8bc17b2 --- /dev/null +++ b/calculations/R2R_Berggren.jl @@ -0,0 +1,48 @@ +using Plots +include("../p_space.jl") + +# contour +vertices = [0, 0.4 - 0.15im, 0.8, 6] +subdivisions = [128, 128, 128] +p, w = get_mesh(vertices, subdivisions) +mesh_E = p.*p ./ (2*0.5) + +# ResonanceEC: Eq. (20) +V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) + +training_points = range(0.75, 0.45, 5) +training_E = Vector{ComplexF64}(undef, length(training_points)) +EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points)) + +for (j, c) in enumerate(training_points) + evals, evecs = eigen(get_H_matrix(V_system(c), p, w)) + i = identify_pole_i(p, evals) + training_E[j] = evals[i] + EC_basis[:, j] = evecs[:, i] +end + +extrapolate_points = range(0.40, 0.25, 5) +ref_E = 0.2 - 0.1im + +exact_E = Vector{ComplexF64}(undef, length(extrapolate_points)) +extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points)) + +EC_basis_w = EC_basis .* w +N_EC = transpose(EC_basis_w) * EC_basis + +for (j, c) in enumerate(extrapolate_points) + exact_E[j] = quick_pole_E(V_system(c)) + + H = get_H_matrix(V_system(c), p, w) + H_EC = transpose(EC_basis_w) * H * EC_basis + evals = eigvals(H_EC, N_EC) + i = argmin(abs.(evals .- ref_E)) + global ref_E = evals[i] + extrapolate_E[j] = evals[i] +end + +scatter(real.(training_E), imag.(training_E), label="training") +scatter!(real.(exact_E), imag.(exact_E), label="exact") +scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") +plot!(real.(mesh_E), imag.(mesh_E), label="contour") +xlims!(0,1) diff --git a/p_space_test.ipynb b/p_space_test.ipynb deleted file mode 100644 index 84ce8c7..0000000 --- a/p_space_test.ipynb +++ /dev/null @@ -1,62 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Plots\n", - "include(\"helper.jl\")\n", - "include(\"p_space.jl\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "vertices = [0, 1 - 0.5im, 2, 6]\n", - "subdivisions = [64, 64, 64]\n", - "p, w = get_mesh(vertices, subdivisions)\n", - "\n", - "scatter(real.(p), imag.(p))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ComplexScaling-FV: Eq. (54)\n", - "V_pq(p, q) = -10 * g1(1, p, q)\n", - "\n", - "H = get_H_matrix(V_pq, p, w)\n", - "evals = eigen(H).values\n", - "\n", - "E_target = 0.258 - 0.164im\n", - "E = nearest(evals, E_target)\n", - "\n", - "print(\"E = $E\")\n", - "scatter(real.(evals), imag.(evals), xlim = (0,1))" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.9.0", - "language": "julia", - "name": "julia-1.9" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.9.0" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/simple_berggren.ipynb b/simple_berggren.ipynb deleted file mode 100644 index cf5bbe8..0000000 --- a/simple_berggren.ipynb +++ /dev/null @@ -1,70 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Plots\n", - "include(\"helper.jl\")\n", - "include(\"p_space.jl\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "vertices = [0, 0.5 - 0.3im, 1, 6]\n", - "subdivisions = [64, 64, 64]\n", - "p, w = get_mesh(vertices, subdivisions)\n", - "\n", - "# resonance example from my thesis\n", - "V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)\n", - "\n", - "basis_eig = eigen(get_H_matrix(V_basis, p, w))\n", - "basis = basis_eig.vectors .* w\n", - "\n", - "evals = basis_eig.values\n", - "E_target = 0.7\n", - "E = include(\"helper.jl\")\n", - "\n", - "print(\"E = $E\")\n", - "scatter(real.(evals), imag.(evals), xlim = (0,2))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ResonanceEC: Eq. (20)\n", - "V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))\n", - "\n", - "H = get_H_matrix(V_system(0.45), p, w)\n", - "H_berggren = transpose(basis) * H * basis\n", - "\n", - "evals = eigvals(H)\n", - "scatter(real.(evals), imag.(evals), xlim = (0, 0.5))" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.10.2", - "language": "julia", - "name": "julia-1.10" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.10.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/test/p_space_basic.jl b/test/p_space_basic.jl new file mode 100644 index 0000000..86d1481 --- /dev/null +++ b/test/p_space_basic.jl @@ -0,0 +1,21 @@ +using Plots +include("../helper.jl") +include("../p_space.jl") + +vertices = [0, 1 - 0.5im, 2, 6] +subdivisions = [64, 64, 64] +p, w = get_mesh(vertices, subdivisions) + +scatter(real.(p), imag.(p)) + +# ComplexScaling-FV: Eq. (54) +V_pq(p, q) = -10 * g1(1, p, q) + +H = get_H_matrix(V_pq, p, w) +evals = eigen(H).values + +E_target = 0.258 - 0.164im +E = nearest(evals, E_target) + +print("E = $E") +scatter(real.(evals), imag.(evals), xlim = (0,1)) diff --git a/test/simple_berggren.jl b/test/simple_berggren.jl new file mode 100644 index 0000000..32bbb62 --- /dev/null +++ b/test/simple_berggren.jl @@ -0,0 +1,30 @@ +using Plots +include("../helper.jl") +include("../p_space.jl") + +vertices = [0, 0.5 - 0.3im, 1, 6] +subdivisions = [64, 64, 64] +p, w = get_mesh(vertices, subdivisions) + +# resonance example from my thesis +V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q) + +basis_eig = eigen(get_H_matrix(V_basis, p, w)) +basis = basis_eig.vectors .* w + +basis_evals = basis_eig.values +E_target = 0.7 +E = nearest(basis_evals, E_target) + +print("pole E = $E") + +# ResonanceEC: Eq. (20) +V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) + +H = get_H_matrix(V_system(0.45), p, w) +H_berggren = transpose(basis) * H * basis + +evals = eigvals(H) +scatter(real.(evals), imag.(evals), label="E") +scatter!(real.(basis_evals), imag.(basis_evals), label="basis", markershape=:x) +xlims!((0, 1)) \ No newline at end of file