From 5b4c38195d7039bd804e9f5aa280de2f377eeb7f Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Tue, 31 Mar 2026 11:25:46 +0100 Subject: [PATCH 1/3] Introduce flags `impose_upper_boundary_condition` and `impose_lower_boundary_condition` to permit cylindrical coordinates (r,z) where the r=0 is included on the grid but no boundary condition is imposed at r=0. --- src/FiniteElementAssembly.jl | 40 +++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl index 0f215ae..f0fba05 100644 --- a/src/FiniteElementAssembly.jl +++ b/src/FiniteElementAssembly.jl @@ -265,6 +265,10 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} include_lower_boundary::Bool # flag to determine if upper endpoint is on the grid include_upper_boundary::Bool + # flag to determine if bc is imposed on lower boundary + impose_lower_boundary_condition::Bool + # flag to determine if bc is imposed on upper boundary + impose_upper_boundary_condition::Bool """ This internal constructor for `FiniteElementCoordinate` takes `scalar_input = ScalarCoordinateInputs(ngrid,nelement,coord_min,coord_max)` to construct @@ -279,7 +283,9 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} # appears in 1D integrals as \int (.) rho(x) d x weight_function::TF=((x)-> 1.0 ), # which boundary condition to store and use for first derivatives - bc::Tbc=NaturalBC() + bc::Tbc=NaturalBC(), + impose_lower_boundary_condition::Bool=true, + impose_upper_boundary_condition::Bool=true ) where {TF <: Function, Tbc <: AbstractBoundaryCondition} ngrid = scalar_input.ngrid nelement = scalar_input.nelement @@ -305,7 +311,9 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} else element_data = nothing end - return FiniteElementCoordinate(name, element_data, weight_function=weight_function, bc=bc) + return FiniteElementCoordinate(name, element_data, weight_function=weight_function, bc=bc, + impose_lower_boundary_condition=impose_lower_boundary_condition, + impose_upper_boundary_condition=impose_upper_boundary_condition) end """ This is the fundamental internal constructor @@ -327,7 +335,9 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} # appears in 1D integrals as \int (.) rho(x) d x weight_function::TF=((x)-> 1.0 ), # which boundary condition to store and use in first derivatives - bc::Tbc=NaturalBC() + bc::Tbc=NaturalBC(), + impose_lower_boundary_condition::Bool=true, + impose_upper_boundary_condition::Bool=true ) where {TF <: Function, Tbc <: AbstractBoundaryCondition} if typeof(element_data) == Nothing # this is a trivial coordinate of length 1 @@ -432,12 +442,18 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} if typeof(bc) == PeriodicBC && (!include_lower_boundary || !include_upper_boundary) error("Must include upper and lower boundaries on the grid for typeof(bc) = $(typeof(bc))") end + if typeof(bc) == PeriodicBC && (!impose_lower_boundary_condition || !impose_upper_boundary_condition) + error("Must impose upper and lower boundary conditions on the grid for typeof(bc) = $(typeof(bc))") + end + impose_lower_boundary_condition = impose_lower_boundary_condition && include_lower_boundary + impose_upper_boundary_condition = impose_upper_boundary_condition && include_upper_boundary return new{typeof(bc)}(name, n_global, ngrid, nelement, domainLength, grid, igrid, ielement, imin, imax, igrid_full, bc, wgts, element_scale, element_shift, element_boundaries, lpoly_data, element_data, derivative_data, - include_lower_boundary, include_upper_boundary) + include_lower_boundary, include_upper_boundary, + impose_lower_boundary_condition, impose_upper_boundary_condition) end end @@ -646,7 +662,7 @@ function impose_boundary_condition_x(boundary_condition::DirichletBC, x::FiniteElementCoordinate,y::FiniteElementCoordinate, weak_form::TF) where {TF <: Function, TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} - if x.include_lower_boundary + if x.impose_lower_boundary_condition # set row to (1, 0, 0, ..., 0) for iy in 1:y.n ixy = ixy_func(1,iy,x.n) @@ -654,7 +670,7 @@ function impose_boundary_condition_x(boundary_condition::DirichletBC, operator_sparse[ixy,ixy] = 1.0 end end - if x.include_upper_boundary + if x.impose_upper_boundary_condition # set row to (1, 0, 0, ..., 0) for iy in 1:y.n ixy = ixy_func(x.n,iy,x.n) @@ -669,7 +685,7 @@ function impose_boundary_condition_x(boundary_condition::PeriodicBC, x::FiniteElementCoordinate,y::FiniteElementCoordinate, weak_form::TF) where {TF <: Function, TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} - if x.include_lower_boundary + if x.impose_lower_boundary_condition # set row to (1, 0, 0, ..., 0, -1) for iy in 1:y.n ixy = ixy_func(1,iy,x.n) @@ -681,7 +697,7 @@ function impose_boundary_condition_x(boundary_condition::PeriodicBC, else error("impose_boundary_condition_x(): x.include_lower_boundary=true is required for boundary_condition_x=PeriodicBC()") end - if x.include_upper_boundary + if x.impose_upper_boundary_condition # assemble lower row contribution to upper row ielement_x = 1 # x index of row assembled to @@ -726,7 +742,7 @@ function impose_boundary_condition_y(boundary_condition::DirichletBC, x::FiniteElementCoordinate,y::FiniteElementCoordinate, weak_form::TF) where {TF <: Function, TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} - if y.include_lower_boundary + if y.impose_lower_boundary_condition # set row to (1, 0, 0, ..., 0) for ix in 1:x.n ixy = ixy_func(ix,1,x.n) @@ -734,7 +750,7 @@ function impose_boundary_condition_y(boundary_condition::DirichletBC, operator_sparse[ixy,ixy] = 1.0 end end - if y.include_upper_boundary + if y.impose_upper_boundary_condition # set row to (1, 0, 0, ..., 0) for ix in 1:x.n ixy = ixy_func(ix,y.n,x.n) @@ -749,7 +765,7 @@ function impose_boundary_condition_y(boundary_condition::PeriodicBC, x::FiniteElementCoordinate,y::FiniteElementCoordinate, weak_form::TF) where {TF <: Function, TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}} - if y.include_lower_boundary + if y.impose_lower_boundary_condition # set row to (1, 0, 0, ..., 0, -1) for ix in 1:x.n ixy = ixy_func(ix,1,x.n) @@ -761,7 +777,7 @@ function impose_boundary_condition_y(boundary_condition::PeriodicBC, else error("impose_boundary_condition_y(): y.include_lower_boundary=true is required for boundary_condition_y=PeriodicBC()") end - if y.include_upper_boundary + if y.impose_upper_boundary_condition # assemble lower row contribution to upper row ielement_y = 1 # y index of row assembled to From 9cac1c3ac2eb40748b6a9882661ea24cb4a1e235 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Tue, 31 Mar 2026 16:29:20 +0100 Subject: [PATCH 2/3] Test the case where a cylindrical coordinate `radial` includes the origin, but no boundary condition is imposed there. --- test/PoissonSolverRadialThetaTests.jl | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/test/PoissonSolverRadialThetaTests.jl b/test/PoissonSolverRadialThetaTests.jl index 6628485..2268d11 100644 --- a/test/PoissonSolverRadialThetaTests.jl +++ b/test/PoissonSolverRadialThetaTests.jl @@ -46,15 +46,24 @@ struct PoissonThetaRadialMatrices end function runtests(; ngrid_theta=15, nelement_theta = 4, Ltheta = 1.0, - ngrid_radial=7, nelement_radial = 4, Lradial = 1.0, print_to_screen=false) + ngrid_radial=7, nelement_radial = 4, Lradial = 1.0, print_to_screen=false, + test_include_origin=false) # create the coordinates theta = FiniteElementCoordinate("theta", ScalarCoordinateInputs(ngrid_theta, nelement_theta, -0.5*Ltheta,0.5*Ltheta,include_boundary_points)) - radial = FiniteElementCoordinate("radial", ScalarCoordinateInputs(ngrid_radial, + if test_include_origin + radial = FiniteElementCoordinate("radial", ScalarCoordinateInputs(ngrid_radial, + nelement_radial, + 0.0,Lradial,include_boundary_points), + weight_function=((r)-> 2.0*pi*r), + impose_lower_boundary_condition=false) + else + radial = FiniteElementCoordinate("radial", ScalarCoordinateInputs(ngrid_radial, nelement_radial, 0.0,Lradial,exclude_lower_boundary_point), weight_function=((r)-> 2.0*pi*r)) + end # allocate fields, noting order (theta,radial) must be the same as in operator definition phi = zeros(radial.n,theta.n) exact_phi = zeros(radial.n,theta.n) @@ -120,8 +129,8 @@ function runtests(; ngrid_theta=15, nelement_theta = 4, Ltheta = 1.0, return nothing end - @testset "Poisson Solver (Radial, Theta) Tests" begin - println("Poisson Solver (Radial, Theta) Tests") + @testset "Poisson Solver (Radial, Theta) Tests test_include_origin=$test_include_origin" begin + println("Poisson Solver (Radial, Theta) Tests test_include_origin=$test_include_origin") test_Poisson(kk=2, phase = 0.0) test_Poisson(kk=2, phase = pi/4.0) test_Poisson(kk=2, phase = pi/3.0) @@ -142,4 +151,5 @@ end # PoissonSolverRadialThetaTests using .PoissonSolverRadialThetaTests -PoissonSolverRadialThetaTests.runtests() \ No newline at end of file +PoissonSolverRadialThetaTests.runtests() +PoissonSolverRadialThetaTests.runtests(test_include_origin=true) \ No newline at end of file From 621b65bc3da88aac5fd8d78d2ae7b36c34db52a7 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 1 Apr 2026 08:47:31 +0100 Subject: [PATCH 3/3] Test on push. --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9201ccb..24646dd 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,7 +1,7 @@ # Based on example from https://github.com/julia-actions/julia-runtest name: Run tests -on: [pull_request, push, workflow_dispatch] +on: [push, workflow_dispatch] permissions: actions: write