From 5346026e6d79050d7680bacadac33d456d2b2017 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 2 Apr 2026 11:32:48 +0100 Subject: [PATCH 1/5] Initial version of advection test, using 1st and 2nd Backward Euler. --- test/AdvectionTests.jl | 146 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 test/AdvectionTests.jl diff --git a/test/AdvectionTests.jl b/test/AdvectionTests.jl new file mode 100644 index 0000000..460ced2 --- /dev/null +++ b/test/AdvectionTests.jl @@ -0,0 +1,146 @@ +""" +Module to test advection in 1D with periodic boundary conditions. +""" +module AdvectionTests + +using Test: @testset, @test +using FiniteElementAssembly +using FiniteElementMatrices +using LinearAlgebra + +struct AdvectionData + Mx::Array{Float64,3} + Sx::Array{Float64,3} +end +function AdvectionData(x::FiniteElementCoordinate) + Mx = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + Sx = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + for ielement in 1:x.nelement + x_data = x.element_data[ielement] + Mx[:,:,ielement] = finite_element_matrix(lagrange_x,lagrange_x,0,x_data) + # integrate advection term v d F / dx by parts in weak form + Sx[:,:,ielement] = -finite_element_matrix(lagrange_x,d_lagrange_dx,0,x_data) + # uncomment to not integrate v d F / dx by parts + #Sx[:,:,ielement] = finite_element_matrix(d_lagrange_dx,lagrange_x,0,x_data) + end + return AdvectionData(Mx,Sx) +end + +function advection_1D_test(; ngrid = 20, nelement = 1, atol = 1.0e-13, + delta_t::Float64=1.0, ntime::Int64=1, vx::Float64=1.0, first_order::Bool=false) + # define inputs needed for the test + coord_min = -1.0 + coord_max = 1.0 + # create the coordinate struct 'x' + x = FiniteElementCoordinate("coord", ScalarCoordinateInputs(ngrid, + nelement, coord_min, coord_max, include_boundary_points), + bc = PeriodicBC()) + ad = AdvectionData(x) + # mass matrix + mass_matrix_1D = assemble_operator(ad.Mx, x, PeriodicBC()) + # 1st-order backward Euler time advance + time_advance_1st_order_weak_form = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + @. time_advance_1st_order_weak_form = ad.Mx + delta_t*vx*ad.Sx + # 2nd-order backward Euler time advance + time_advance_2nd_order_weak_form = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + @. time_advance_2nd_order_weak_form = ad.Mx + (2.0/3.0)*delta_t*vx*ad.Sx + #@. time_advance_2nd_order_weak_form = 1.5*ad.Mx + delta_t*vx*ad.Sx + # assembled time advance matrix + time_advance_1st_order_1D = assemble_operator(time_advance_1st_order_weak_form, x, PeriodicBC()) + lu_time_advance_1st_order = lu(time_advance_1st_order_1D) + time_advance_2nd_order_1D = assemble_operator(time_advance_2nd_order_weak_form, x, PeriodicBC()) + lu_time_advance_2nd_order = lu(time_advance_2nd_order_1D) + + # a periodic function + f0 = zeros(x.n) # f at time level n + f1 = zeros(x.n) # f at time level n-1 + f2 = zeros(x.n) # f at time level n-2 + # the exact solution + fexact = zeros(x.n) + # dummy array for RHS + rhs = zeros(x.n) + # dummy array for error + ferr = zeros(x.n) + + # set the initial condition + # and the expected solution + phase = 0.25 + for ix in 1:x.n + xarg = 2.0*x.grid[ix]/x.L + f0[ix] = cospi(xarg + phase) + xparg = 2.0*(x.grid[ix] - vx*delta_t*ntime)/x.L + fexact[ix] = cospi(xparg + phase) + end + + function backward_Euler_1(f::Vector{Float64}) + mul!(rhs, mass_matrix_1D, f) + ldiv!(f, lu_time_advance_1st_order, rhs) + return nothing + end + function backward_Euler_2(f0::Vector{Float64},f1::Vector{Float64},f2::Vector{Float64}) + # rhs = Mx*( (4/3) f1 - (1/3) f2) + # use f2 as dummy as we do not need it again + @. f2 = (4.0/3.0)*f1 - (1.0/3.0)*f2 + mul!(rhs, mass_matrix_1D, f2) + # find new f0 + ldiv!(f0, lu_time_advance_2nd_order, rhs) + # update history + copyto!(f2,f1) + copyto!(f1,f0) + return nothing + end + + if first_order + ntime_1st_order = ntime + ntime_2nd_order = 0 + else + ntime_1st_order = Int64(ntime > 0) + ntime_2nd_order = ntime + end + # println(ntime_1st_order) + # println(ntime_2nd_order) + time = 0.0 + # backward Euler for step 1 + # store initial condition in history + copyto!(f1,f0) + for it in 1:ntime_1st_order + time += delta_t + backward_Euler_1(f0) + end + # update history to prepare for second order advance + copyto!(f2,f1) + copyto!(f1,f0) + # second order backward Euler for remaining steps + for it in 2:ntime_2nd_order + time += delta_t + backward_Euler_2(f0,f1,f2) + end + + # test f + # println(fexact) + # println(f0) + @. ferr = abs(f0 - fexact) + max_ferr = maximum(ferr) + # check total integral vanishes + @test abs(integral((x)->1.0,f0,x)) < 1.0e-10 + # test max(abs.()) norm + @test max_ferr < atol + # test boundary condition + @test abs(f0[1]-f0[end]) < 1.0e-14 + return nothing +end + +function runtests() + @testset "Advection 1D" begin + println("Advection 1D Tests") + advection_1D_test(; ngrid=15, nelement=1, ntime=10000, vx=1.0, delta_t=0.0001, first_order=false, atol=5.0e-7) + advection_1D_test(; ngrid=15, nelement=1, ntime=10000, vx=1.0, delta_t=0.0001, first_order=true, atol=8e-3) + end + return nothing +end + +end # AdvectionTests + +using .AdvectionTests + +AdvectionTests.runtests() \ No newline at end of file From 8dbf3f9cc34ade55fc286416c9a625e5fbf246aa Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 2 Apr 2026 13:10:12 +0100 Subject: [PATCH 2/5] Include advection tests in `test/runtests.jl`. --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 8d5957c..438a17e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ using Test include(joinpath(@__DIR__, "ElectronIonCollisionsTest.jl")) include(joinpath(@__DIR__, "TestParticleCollisionsVpaVperp.jl")) include(joinpath(@__DIR__, "SlowingDownTest.jl")) + include(joinpath(@__DIR__, "AdvectionTests.jl")) end end # FiniteElementAssemblyTests \ No newline at end of file From d2f6d3ad40ae838c7244c8d473be2d066c61e2f2 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 22 Apr 2026 14:41:37 +0100 Subject: [PATCH 3/5] A test of integration in periodic boundary conditions. --- test/PeriodicBcTests.jl | 75 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/test/PeriodicBcTests.jl b/test/PeriodicBcTests.jl index 4a5e02c..bbc29c3 100644 --- a/test/PeriodicBcTests.jl +++ b/test/PeriodicBcTests.jl @@ -4,7 +4,9 @@ Module to test differentiation in 1D with periodic boundary conditions. module PeriodicBcTests using Test: @testset, @test +using LinearAlgebra using FiniteElementAssembly +using FiniteElementMatrices function periodic_differentiation_1D_test(; ngrid = 3, nelement = 1, atol = 1.0e-13) # define inputs needed for the test @@ -36,12 +38,85 @@ function periodic_differentiation_1D_test(; ngrid = 3, nelement = 1, atol = 1.0e return nothing end +function periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 1, atol = 1.0e-13) + # define inputs needed for the test + coord_min = -1.0 + coord_max = 1.0 + # create the coordinate struct 'x' + x = FiniteElementCoordinate("coord", ScalarCoordinateInputs(ngrid, + nelement, coord_min, coord_max, include_boundary_points), + bc = PeriodicBC()) + + # a test of periodic functions + # solve d^2 phi / d x^2 = f = sin(2pi x/L + phase) + # exact solution phi(x) = A + B z - (L/2pi)^2 sin(2pi x/L + phase ) + f = zeros(x.n) + phi = zeros(x.n) + exact_phi = zeros(x.n) + phase = pi/4.0 + for ix in 1:x.n + xarg = 2.0*x.grid[ix]/x.L + phase + f[ix] = sinpi(xarg) + exact_phi[ix] = - (x.L/(2.0*pi))^2*sinpi(xarg) + end + # dummy array + weak_rhs = zeros(x.n) + + # required weak matrices + M_x = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + K_x = Array{Float64,3}(undef,x.ngrid,x.ngrid,x.nelement) + for ielement_x in 1:x.nelement + x_data = x.element_data[ielement_x] + @views M_x[:,:,ielement_x] = finite_element_matrix(lagrange_x,lagrange_x,0,x_data) + @views K_x[:,:,ielement_x] = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,x_data) + end + # assembled sparse matrices + mass_matrix = assemble_operator(M_x, x, PeriodicBC()) + stiffness_matrix = assemble_operator(K_x, x, PeriodicBC()) + lu_mass_matrix = lu(mass_matrix) + # Dirichlet BC. + stiffness_matrix[1,:] .= 0.0 + stiffness_matrix[1,1] = 1.0 + stiffness_matrix[end,:] .= 0.0 + stiffness_matrix[end,end] = 1.0 + lu_Laplacian = lu(stiffness_matrix) + + # form the weak RHS + mul!(weak_rhs,mass_matrix,f) + # impose the Dirichlet BC to fix the constants A and B + weak_rhs[1] = exact_phi[1] + weak_rhs[end] = exact_phi[end] + # solve the system + ldiv!(phi, lu_Laplacian, weak_rhs) + #println(phi) + #println(exact_phi) + max_err_phi = maximum(abs.(phi .- exact_phi)) + # println(max_err_phi) + + # check forward differentiation + stiffness_matrix_differentiation = assemble_operator(K_x, x, PeriodicBC()) + + numerical_f = zeros(x.n) + mul!(weak_rhs,stiffness_matrix_differentiation,exact_phi) + ldiv!(numerical_f,lu_mass_matrix,weak_rhs) + #println(f) + #println(numerical_f) + max_err_f = maximum(abs.(f .- numerical_f)) + # println(max_err_f) + @test max_err_f < 100*atol + @test max_err_phi < atol + return nothing +end + function runtests() @testset "Periodic differentiation 1D" begin println("Periodic Differentiation 1D Tests") periodic_differentiation_1D_test(; ngrid = 25, nelement = 1) periodic_differentiation_1D_test(; ngrid = 25, nelement = 2) periodic_differentiation_1D_test(; ngrid = 25, nelement = 3, atol = 1.0e-12) + periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 1) + periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 2) + periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 3) end return nothing end From 9acf97f4cd3b40ad6be98d8749c94fa255887644 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 22 Apr 2026 14:50:09 +0100 Subject: [PATCH 4/5] Implement Dirichlet BC for 1D operators, update tests. --- src/FiniteElementAssembly.jl | 14 +++++++++++++- test/PeriodicBcTests.jl | 9 ++++----- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl index f0fba05..246f441 100644 --- a/src/FiniteElementAssembly.jl +++ b/src/FiniteElementAssembly.jl @@ -624,13 +624,25 @@ end """ Methods for imposing boundary conditions on 1D operators """ -function impose_boundary_condition_x(boundary_condition::Union{NaturalBC,DirichletBC}, +function impose_boundary_condition_x(boundary_condition::NaturalBC, operator_sparse::TMatrix, n::Int64, ngrid::Int64, QQ::Array{Float64,3} ) where {TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} # do nothing return nothing end +function impose_boundary_condition_x(boundary_condition::DirichletBC, + operator_sparse::TMatrix, + n::Int64, ngrid::Int64, QQ::Array{Float64,3} + ) where {TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} + # set lower row to (1, 0, 0, ..., 0) + operator_sparse[1,:] .= 0.0 + operator_sparse[1,1] = 1.0 + # set upper row to (0, 0, 0, ..., 1) + operator_sparse[end,:] .= 0.0 + operator_sparse[end,end] = 1.0 + return nothing +end function impose_boundary_condition_x(boundary_condition::PeriodicBC, operator_sparse::TMatrix, n::Int64, ngrid::Int64, QQ::Array{Float64,3} diff --git a/test/PeriodicBcTests.jl b/test/PeriodicBcTests.jl index bbc29c3..437cc72 100644 --- a/test/PeriodicBcTests.jl +++ b/test/PeriodicBcTests.jl @@ -72,13 +72,9 @@ function periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 1, atol = 1.0e-13) end # assembled sparse matrices mass_matrix = assemble_operator(M_x, x, PeriodicBC()) - stiffness_matrix = assemble_operator(K_x, x, PeriodicBC()) + stiffness_matrix = assemble_operator(K_x, x, DirichletBC()) lu_mass_matrix = lu(mass_matrix) # Dirichlet BC. - stiffness_matrix[1,:] .= 0.0 - stiffness_matrix[1,1] = 1.0 - stiffness_matrix[end,:] .= 0.0 - stiffness_matrix[end,end] = 1.0 lu_Laplacian = lu(stiffness_matrix) # form the weak RHS @@ -114,6 +110,9 @@ function runtests() periodic_differentiation_1D_test(; ngrid = 25, nelement = 1) periodic_differentiation_1D_test(; ngrid = 25, nelement = 2) periodic_differentiation_1D_test(; ngrid = 25, nelement = 3, atol = 1.0e-12) + end + @testset "Periodic ODE Integration 1D" begin + println("Periodic ODE Integration 1D Tests") periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 1) periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 2) periodic_ODE_solve_1D_test(; ngrid = 20, nelement = 3) From 2c08db02561cc7002d021a5c44532069ff291db2 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 6 May 2026 16:17:42 +0100 Subject: [PATCH 5/5] Reorganise advection test with periodic boundary conditions. --- ...AdvectionTests.jl => PeriodicBcAdvectionTests.jl} | 12 ++++++------ test/runtests.jl | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) rename test/{AdvectionTests.jl => PeriodicBcAdvectionTests.jl} (95%) diff --git a/test/AdvectionTests.jl b/test/PeriodicBcAdvectionTests.jl similarity index 95% rename from test/AdvectionTests.jl rename to test/PeriodicBcAdvectionTests.jl index 460ced2..458a9db 100644 --- a/test/AdvectionTests.jl +++ b/test/PeriodicBcAdvectionTests.jl @@ -1,7 +1,7 @@ """ Module to test advection in 1D with periodic boundary conditions. """ -module AdvectionTests +module PeriodicBcAdvectionTests using Test: @testset, @test using FiniteElementAssembly @@ -131,16 +131,16 @@ function advection_1D_test(; ngrid = 20, nelement = 1, atol = 1.0e-13, end function runtests() - @testset "Advection 1D" begin - println("Advection 1D Tests") + @testset "Periodic Advection 1D" begin + println("Periodic Advection 1D Tests") advection_1D_test(; ngrid=15, nelement=1, ntime=10000, vx=1.0, delta_t=0.0001, first_order=false, atol=5.0e-7) advection_1D_test(; ngrid=15, nelement=1, ntime=10000, vx=1.0, delta_t=0.0001, first_order=true, atol=8e-3) end return nothing end -end # AdvectionTests +end # PeriodicBcAdvectionTests -using .AdvectionTests +using .PeriodicBcAdvectionTests -AdvectionTests.runtests() \ No newline at end of file +PeriodicBcAdvectionTests.runtests() \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 438a17e..cb1e156 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,10 +10,10 @@ using Test include(joinpath(@__DIR__, "PoissonSolverThetaZedTests.jl")) include(joinpath(@__DIR__, "PoissonSolverRadialThetaTests.jl")) include(joinpath(@__DIR__, "PeriodicBcTests.jl")) + include(joinpath(@__DIR__, "PeriodicBcAdvectionTests.jl")) include(joinpath(@__DIR__, "ElectronIonCollisionsTest.jl")) include(joinpath(@__DIR__, "TestParticleCollisionsVpaVperp.jl")) include(joinpath(@__DIR__, "SlowingDownTest.jl")) - include(joinpath(@__DIR__, "AdvectionTests.jl")) end end # FiniteElementAssemblyTests \ No newline at end of file