From 83f198901035ec1cd23dbc0fb1addb94752a74b6 Mon Sep 17 00:00:00 2001 From: Alex Quinlan Date: Tue, 14 Jun 2022 16:07:53 +0200 Subject: [PATCH 1/2] Convergence testing --- src/convergence_test.jl | 43 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 src/convergence_test.jl diff --git a/src/convergence_test.jl b/src/convergence_test.jl new file mode 100644 index 0000000..b4050e5 --- /dev/null +++ b/src/convergence_test.jl @@ -0,0 +1,43 @@ +include("Beam1D.jl") +include("Diagnostics.jl") +include("analytic_testcases.jl") +import Plots +import LinearAlgebra + +L = 1.0 +q0 = 3 +M_0 = 9 +a = 1 +EI_const = 1 +#analytic solution with BCs and q function +ana_sol, BCs, q_func = case_7_constant_load(L , q0 ,EI_const) +""" just copy and paste to run example +ana_sol, BCs, q_func = case_1_supported_beam_(L , a, q0 ,EI_const) +#case_7_constant_load(L , q0 ,EI_const) +#case_8_partly_constant_load(L, q0, a, EI_const) +#case_9_decreasing_load(L, q0, EI_const) +#case_10_free_momentum_at_a(L, a, M_0, EI_const) +""" + +#parameters +pars = (mu=x->1,EI=x->EI_const,q=q_func) + +resolutions = [3, 5, 11, 21, 51, 101, 201] +errors = [] + +for resolution in resolutions + grid = collect(LinRange(0,L,resolution)) + sys = Beam1D.System(Beam1D.Problem(pars,BCs,grid)) + sol = Beam1D.solve_st(sys) + error = LinearAlgebra.norm(ana_sol(grid)-sol(grid)) / resolution + append!(errors, [error]) +end +#plot +plt = Plots.plot(resolutions, errors, xaxis=:log, yaxis=:log, seriestype=:scatter, label = "Convergence") +Plots.ylabel!("error") +Plots.xlabel!("resolution") +println("Plotting convergence. Press enter to continue.") +Plots.gui(plt) +readline() + +# NOTE: This should converge but it currently fails to. From 86adef9a2d2d9870a5995d102bc64af3ac58d53e Mon Sep 17 00:00:00 2001 From: Alex Quinlan Date: Tue, 28 Jun 2022 14:47:21 +0200 Subject: [PATCH 2/2] Working convergence at roughly first order --- Manifest.toml | 24 ++++++++++++++++-- Project.toml | 2 ++ src/convergence_test.jl | 54 ++++++++++++++++++++++------------------- 3 files changed, 53 insertions(+), 27 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 13694c4..0bc4626 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -14,6 +14,18 @@ version = "2.3.0" [[ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +[[Arpack]] +deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"] +git-tree-sha1 = "91ca22c4b8437da89b030f08d71db55a379ce958" +uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +version = "0.5.3" + +[[Arpack_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"] +git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e" +uuid = "68821587-b530-5797-8361-c406ea357684" +version = "3.5.1+1" + [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -410,7 +422,7 @@ uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" version = "2.36.0+0" [[LinearAlgebra]] -deps = ["Libdl"] +deps = ["Libdl", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[LogExpFunctions]] @@ -473,6 +485,10 @@ git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" version = "1.3.5+1" +[[OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" + [[OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" @@ -561,7 +577,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[Random]] -deps = ["Serialization"] +deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[RecipesBase]] @@ -880,6 +896,10 @@ git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" version = "0.15.1+0" +[[libblastrampoline_jll]] +deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" + [[libfdk_aac_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" diff --git a/Project.toml b/Project.toml index 6602882..50b71db 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,10 @@ authors = ["Alex Quinlan", "Henry Jacobson", "Pia Callmer", "Nicola Sabbadini"] version = "0.1.0" [deps] +Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" CubicHermiteSpline = "2338ea46-7d83-48cf-b44e-cf8e0b5b9cc1" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/src/convergence_test.jl b/src/convergence_test.jl index b4050e5..3f80303 100644 --- a/src/convergence_test.jl +++ b/src/convergence_test.jl @@ -4,40 +4,44 @@ include("analytic_testcases.jl") import Plots import LinearAlgebra -L = 1.0 -q0 = 3 -M_0 = 9 -a = 1 -EI_const = 1 -#analytic solution with BCs and q function -ana_sol, BCs, q_func = case_7_constant_load(L , q0 ,EI_const) -""" just copy and paste to run example -ana_sol, BCs, q_func = case_1_supported_beam_(L , a, q0 ,EI_const) -#case_7_constant_load(L , q0 ,EI_const) -#case_8_partly_constant_load(L, q0, a, EI_const) -#case_9_decreasing_load(L, q0, EI_const) -#case_10_free_momentum_at_a(L, a, M_0, EI_const) -""" - -#parameters -pars = (mu=x->1,EI=x->EI_const,q=q_func) - -resolutions = [3, 5, 11, 21, 51, 101, 201] -errors = [] +pars = (mu=x->1, EI=x->1, q=x->-(x<0.5)) + +BCs = Dict((0,'H')=>1, + (0,'G')=>2, + (1,'M')=>3, + (1,'Q')=>4) + + +resolutions = [10, 20, 50, 100, 200, 400, 800] + +best_resolution = last(resolutions)*10 +grid = collect(LinRange(0, 1, best_resolution)) +sys = Beam1D.System(Beam1D.Problem(pars,BCs,grid)) +best = Beam1D.solve_st(sys) + +errors = [] for resolution in resolutions - grid = collect(LinRange(0,L,resolution)) + grid = collect(LinRange(0,1,resolution)) sys = Beam1D.System(Beam1D.Problem(pars,BCs,grid)) sol = Beam1D.solve_st(sys) - error = LinearAlgebra.norm(ana_sol(grid)-sol(grid)) / resolution + error = LinearAlgebra.norm(best(grid) - sol(grid)) / resolution append!(errors, [error]) end #plot -plt = Plots.plot(resolutions, errors, xaxis=:log, yaxis=:log, seriestype=:scatter, label = "Convergence") +plt = Plots.plot(resolutions, errors, xaxis=:log, yaxis=:log, seriestype=:scatter, label = "Convergence", legend=:bottomleft) Plots.ylabel!("error") Plots.xlabel!("resolution") + +ref1 = errors[1]*resolutions[1]*resolutions.^-1 +ref1 = reshape(ref1, length(errors), 1) +ref2 = errors[1]*resolutions[1]^2*resolutions.^-2 +ref2 = reshape(ref2, length(errors), 1) + +Plots.plot!(resolutions, ref1, xaxis=:log, yaxis=:log, label = "Reference 1st Order") +Plots.plot!(resolutions, ref2, xaxis=:log, yaxis=:log, label = "Reference 2nd Order") + println("Plotting convergence. Press enter to continue.") + Plots.gui(plt) readline() - -# NOTE: This should converge but it currently fails to.