From b32ebce78ebc4e352a9d3eb59c26224585b2a042 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 19 Sep 2024 09:52:00 +0100 Subject: [PATCH 1/2] Tests of the 1D ODE solvers provided by the gauss_legendre.jl module. These 1D ODE solvers are useful to show the basics of how the more complex 2D PDE solvers are constructed, and to show how different types of boundary conditions are imposed. --- moment_kinetics/test/calculus_tests.jl | 292 +++++++++++++++++++++++++ 1 file changed, 292 insertions(+) diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl index 58fc00f1f..5ae798d5f 100644 --- a/moment_kinetics/test/calculus_tests.jl +++ b/moment_kinetics/test/calculus_tests.jl @@ -8,6 +8,7 @@ using moment_kinetics.calculus: laplacian_derivative! using MPI using Random +using LinearAlgebra: mul!, ldiv! function runtests() @testset "calculus" verbose=use_verbose begin @@ -1565,6 +1566,297 @@ function runtests() norm=maxabs_norm) end end + + @testset "GaussLegendre pseudospectral cylindrical laplacian ODE solve, zero" verbose=false begin + @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ + ( + (1, 8, 5.e-2), + (1, 9, 3.e-2), + (1, 10, 4.e-3), + (1, 11, 3.e-3), + (1, 12, 6.e-4), + (1, 13, 4.e-4), + (1, 14, 2.e-4), + (1, 15, 5.e-5), + (1, 16, 3.e-5), + (1, 17, 9.e-6), + + (2, 6, 8.e-3), + (2, 7, 4.e-3), + (2, 8, 4.e-4), + (2, 9, 2.e-4), + (2, 10, 5.e-5), + (2, 11, 1.e-5), + (2, 12, 5.e-6), + (2, 13, 4.e-7), + (2, 14, 4.e-7), + (2, 15, 5.e-8), + (2, 16, 2.e-8), + (2, 17, 5.e-9), + + (3, 6, 2.e-3), + (3, 7, 2.e-4), + (3, 8, 5.e-5), + (3, 9, 3.e-6), + (3, 10, 2.e-6), + (3, 11, 3.e-7), + (3, 12, 6.e-8), + (3, 13, 9.e-9), + (3, 14, 2.e-9), + (3, 15, 4.e-10), + (3, 16, 3.e-11), + (3, 17, 1.e-11), + + (4, 5, 1.e-3), + (4, 6, 8.e-5), + (4, 7, 3.e-5), + (4, 8, 3.e-6), + (4, 9, 6.e-7), + (4, 10, 8.e-8), + (4, 11, 2.e-8), + (4, 12, 3.e-9), + (4, 13, 3.e-10), + (4, 14, 5.e-11), + (4, 15, 3.e-12), + (4, 16, 2.e-12), + (4, 17, 1.e-12), + + (5, 5, 4.e-4), + (5, 6, 3.e-5), + (5, 7, 6.e-6), + (5, 8, 4.e-7), + (5, 9, 9.e-8), + (5, 10, 4.e-9), + (5, 11, 2.e-9), + (5, 12, 8.e-11), + (5, 13, 3.e-11), + (5, 14, 1.e-12), + (5, 15, 4.e-13), + (5, 16, 4.e-13), + (5, 17, 2.e-12), + ) + + # define inputs needed for the test + L = 6.0 + bc = "zero" + element_spacing_option = "uniform" + # create the coordinate struct 'x' + # This test runs effectively in serial, so implicitly uses + # `ignore_MPI=true` to avoid errors due to communicators not being fully + # set up. + x, spectral = define_test_coordinate("vperp"; ngrid=ngrid, + nelement=nelement, L=L, + discretization="gausslegendre_pseudospectral", + bc=bc, + element_spacing_option=element_spacing_option, + collision_operator_dim=false) + expected_f = @. exp(-x.grid^2) + d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2) + # create array for the numerical solution + f = similar(expected_f) + # create array for RHS vector b + b = similar(expected_f) + # solve for f + mul!(b,spectral.mass_matrix,d2f) + # Dirichlet zero BC at upper endpoint + b[end] = 0.0 + # solve ODE + ldiv!(f,spectral.L_matrix_lu,b) + #err = maximum(abs.(f.-expected_f)) + #println("$nelement $ngrid $err") + @test isapprox(f, expected_f, rtol=rtol, atol=1.e-10, + norm=maxabs_norm) + end + end + + @testset "GaussLegendre pseudospectral second derivative ODE solve, zero" verbose=false begin + @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ + ( + (1, 23, 1.e-2), + (1, 24, 5.e-3), + (1, 25, 3.e-3), + (1, 26, 2.e-3), + (1, 27, 1.e-3), + (1, 28, 6.e-4), + (1, 29, 5.e-4), + (1, 30, 3.e-4), + (1, 31, 2.e-4), + (1, 32, 8.e-5), + + (2, 9, 3.e-2), + (2, 10, 2.e-2), + (2, 11, 4.e-3), + (2, 12, 1.e-3), + (2, 13, 4.e-4), + (2, 14, 2.e-4), + (2, 15, 9.e-5), + (2, 16, 3.e-5), + (2, 17, 7.e-6), + + (3, 9, 2.e-2), + (3, 10, 2.e-3), + (3, 11, 6.e-4), + (3, 12, 2.e-4), + (3, 13, 7.e-5), + (3, 14, 1.e-5), + (3, 15, 9.e-6), + (3, 16, 9.e-7), + (3, 17, 1.e-6), + + (4, 7, 3.e-3), + (4, 8, 8.e-4), + (4, 9, 2.e-4), + (4, 10, 4.e-5), + (4, 11, 2.e-5), + (4, 12, 4.e-6), + (4, 13, 9.e-7), + (4, 14, 4.e-7), + (4, 15, 3.e-8), + (4, 16, 3.e-8), + (4, 17, 4.e-9), + + (5, 8, 2.e-4), + (5, 9, 7.e-5), + (5, 10, 7.e-6), + (5, 11, 4.e-6), + (5, 12, 3.e-7), + (5, 13, 3.e-7), + (5, 14, 9.e-9), + (5, 15, 1.e-8), + (5, 16, 3.e-10), + (5, 17, 4.e-10), + ) + + # define inputs needed for the test + L = 12.0 + bc = "zero" + element_spacing_option = "uniform" + # create the coordinate struct 'x' + # This test runs effectively in serial, so implicitly uses + # `ignore_MPI=true` to avoid errors due to communicators not being fully + # set up. + x, spectral = define_test_coordinate("vpa"; ngrid=ngrid, + nelement=nelement, L=L, + discretization="gausslegendre_pseudospectral", + bc=bc, + element_spacing_option=element_spacing_option, + collision_operator_dim=false) + expected_f = @. exp(-x.grid^2) + d2f = @. (4.0*x.grid^2 - 2.0)*exp(-x.grid^2) + # create array for the numerical solution + f = similar(expected_f) + # create array for RHS vector b + b = similar(expected_f) + # solve for f + mul!(b,spectral.mass_matrix,d2f) + # Dirichlet zero BC at lower and upper endpoint + b[1] = 0.0 + b[end] = 0.0 + # solve ODE + ldiv!(f,spectral.L_matrix_lu,b) + #err = maximum(abs.(f.-expected_f)) + #maxfe = maximum(expected_f) + #maxf = maximum(f) + #minf = minimum(f) + #println("$nelement $ngrid $err $maxfe $maxf $minf") + @test isapprox(f, expected_f, rtol=rtol, atol=1.e-10, + norm=maxabs_norm) + end + end + + @testset "GaussLegendre pseudospectral second derivative ODE solve, periodic" verbose=false begin + @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ + ( + (1, 10, 6.e-6), + (1, 11, 2.e-6), + (1, 12, 8.e-8), + (1, 13, 2.e-8), + (1, 14, 1.e-9), + (1, 15, 2.e-10), + (1, 16, 9.e-12), + (1, 17, 2.e-12), + + (2, 9, 5.e-8), + (2, 10, 6.e-9), + (2, 11, 3.e-10), + (2, 12, 3.e-11), + (2, 13, 7.e-13), + (2, 14, 1.e-13), + (2, 15, 1.e-13), + (2, 16, 1.e-13), + (2, 17, 1.e-13), + + (3, 9, 2.e-9), + (3, 10, 7.e-11), + (3, 11, 4.e-12), + (3, 12, 3.e-13), + (3, 13, 3.e-13), + (3, 14, 1.e-13), + (3, 15, 1.e-13), + (3, 16, 3.e-13), + (3, 17, 1.e-13), + + (4, 7, 5.e-8), + (4, 8, 3.e-9), + (4, 9, 8.e-11), + (4, 10, 4.e-12), + (4, 11, 3.e-13), + (4, 12, 1.e-13), + (4, 13, 3.e-13), + (4, 14, 1.e-13), + (4, 15, 1.e-13), + (4, 16, 3.e-13), + (4, 17, 3.e-13), + + (5, 8, 5.e-10), + (5, 9, 9.e-12), + (5, 10, 5.e-13), + (5, 11, 3.e-13), + (5, 12, 1.e-13), + (5, 13, 5.e-13), + (5, 14, 1.e-13), + (5, 15, 2.e-13), + (5, 16, 5.e-13), + (5, 17, 5.e-13), + ) + + # define inputs needed for the test + L = 1.0 + bc = "periodic" + element_spacing_option = "uniform" + phase = pi/3.0 + # create the coordinate struct 'x' + # This test runs effectively in serial, so implicitly uses + # `ignore_MPI=true` to avoid errors due to communicators not being fully + # set up. + x, spectral = define_test_coordinate("vpa"; ngrid=ngrid, + nelement=nelement, L=L, + discretization="gausslegendre_pseudospectral", + bc=bc, + element_spacing_option=element_spacing_option, + collision_operator_dim=false) + expected_f = @. sin((2.0*pi*x.grid/x.L) + phase) + d2f = @. -((2.0*pi/x.L)^2)*sin((2.0*pi*x.grid/x.L)+phase) + # create array for the numerical solution + f = similar(expected_f) + # create array for RHS vector b + b = similar(expected_f) + # solve for f + mul!(b,spectral.mass_matrix,d2f) + # Dirichlet zero BC at lower endpoint, periodic bc at upper endpoint + b[1] = sin((2.0*pi*x.grid[1]/x.L) + phase) # fixes constant piece of solution + b[end] = 0.0 # makes sure periodicity is enforced + # solve ODE + ldiv!(f,spectral.L_matrix_lu,b) + #err = maximum(abs.(f.-expected_f)) + #maxfe = maximum(expected_f) + #maxf = maximum(f) + #minf = minimum(f) + #println("$nelement $ngrid $err $maxfe $maxf $minf") + @test isapprox(f, expected_f, rtol=rtol, atol=1.e-10, + norm=maxabs_norm) + end + end end end From 8f2843a01fbd8aeed892b917e5f4db7832167b33 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 21 Sep 2024 14:18:16 +0100 Subject: [PATCH 2/2] Add a few comments to 1D ODE solver tests --- moment_kinetics/test/calculus_tests.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl index 5ae798d5f..23a77e71a 100644 --- a/moment_kinetics/test/calculus_tests.jl +++ b/moment_kinetics/test/calculus_tests.jl @@ -1651,13 +1651,15 @@ function runtests() element_spacing_option=element_spacing_option, collision_operator_dim=false) expected_f = @. exp(-x.grid^2) - d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2) + # Test solver for + # Laplacian_f = 1/x * d/dx(x*df/dx) + Laplacian_f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2) # create array for the numerical solution f = similar(expected_f) # create array for RHS vector b b = similar(expected_f) # solve for f - mul!(b,spectral.mass_matrix,d2f) + mul!(b,spectral.mass_matrix,Laplacian_f) # Dirichlet zero BC at upper endpoint b[end] = 0.0 # solve ODE @@ -1742,6 +1744,8 @@ function runtests() element_spacing_option=element_spacing_option, collision_operator_dim=false) expected_f = @. exp(-x.grid^2) + # Test solver for + # d2f = d^2f/dx^2 d2f = @. (4.0*x.grid^2 - 2.0)*exp(-x.grid^2) # create array for the numerical solution f = similar(expected_f) @@ -1836,6 +1840,8 @@ function runtests() element_spacing_option=element_spacing_option, collision_operator_dim=false) expected_f = @. sin((2.0*pi*x.grid/x.L) + phase) + # Test solver for + # d2f = d^2f/dx^2 d2f = @. -((2.0*pi/x.L)^2)*sin((2.0*pi*x.grid/x.L)+phase) # create array for the numerical solution f = similar(expected_f)