When calculating error norms, skip points that are set by radial bc

For Dirichlet bc this probably does not matter (as the boundary points
should have zero error anyway), but when using Neumann boundary
conditions the radial boundaries may give the largest timestep errors
even though they are not actually evolved values - ignore these to avoid
unnecessarily restricting the adaptive timestep.

Author: johnomotani

moment_kinetics/src/boundary_conditions.jl

@@ -10,7 +10,7 @@ export enforce_neutral_boundary_conditions!
 using LinearAlgebra: dot
 using SpecialFunctions: erfc
 
-using ..array_allocation: allocate_shared_float
+using ..array_allocation: allocate_shared_float, allocate_shared_bool
 using ..calculus: integral, reconcile_element_boundaries_MPI!
 using ..communication
 using ..coordinates: coordinate
@@ -219,8 +219,83 @@ function create_r_boundary_info(input_dict, pdf, moments, r, z, vperp, vpa, vzet
         end
     end
 
+    ion_error_norm_skip_inner = allocate_shared_bool(z.n)
+    ion_error_norm_skip_outer = allocate_shared_bool(z.n)
+    electron_error_norm_skip_inner = allocate_shared_bool(z.n)
+    electron_error_norm_skip_outer = allocate_shared_bool(z.n)
+    neutral_error_norm_skip_inner = allocate_shared_bool(z.n)
+    neutral_error_norm_skip_outer = allocate_shared_bool(z.n)
+    @begin_serial_region()
+    @serial_region begin
+        if r.irank == 0
+            for s ∈ inner_sections
+                if (isa(s.ion, ion_r_boundary_section_Neumann)
+                    || isa(s.ion, ion_r_boundary_section_Dirichlet))
+
+                    ion_error_norm_skip_inner[s.z_range] .= true
+                else
+                    ion_error_norm_skip_inner[s.z_range] .= false
+                end
+                if (isa(s.electron, electron_r_boundary_section_Neumann)
+                    || isa(s.electron, electron_r_boundary_section_Dirichlet))
+
+                    electron_error_norm_skip_inner[s.z_range] .= true
+                else
+                    electron_error_norm_skip_inner[s.z_range] .= false
+                end
+                if (isa(s.neutral, neutral_r_boundary_section_Neumann)
+                    || isa(s.neutral, neutral_r_boundary_section_Dirichlet))
+
+                    neutral_error_norm_skip_inner[s.z_range] .= true
+                else
+                    neutral_error_norm_skip_inner[s.z_range] .= false
+                end
+            end
+        else
+            ion_error_norm_skip_inner .= false
+            electron_error_norm_skip_inner .= false
+            neutral_error_norm_skip_inner .= false
+        end
+
+        if r.irank == r.nrank - 1
+            for s ∈ outer_sections
+                if (isa(s.ion, ion_r_boundary_section_Neumann)
+                    || isa(s.ion, ion_r_boundary_section_Dirichlet))
+
+                    ion_error_norm_skip_outer[s.z_range] .= true
+                else
+                    ion_error_norm_skip_outer[s.z_range] .= false
+                end
+                if (isa(s.electron, electron_r_boundary_section_Neumann)
+                    || isa(s.electron, electron_r_boundary_section_Dirichlet))
+
+                    electron_error_norm_skip_outer[s.z_range] .= true
+                else
+                    electron_error_norm_skip_outer[s.z_range] .= false
+                end
+                if (isa(s.neutral, neutral_r_boundary_section_Neumann)
+                    || isa(s.neutral, neutral_r_boundary_section_Dirichlet))
+
+                    neutral_error_norm_skip_outer[s.z_range] .= true
+                else
+                    neutral_error_norm_skip_outer[s.z_range] .= false
+                end
+            end
+        else
+            ion_error_norm_skip_outer .= false
+            electron_error_norm_skip_outer .= false
+            neutral_error_norm_skip_outer .= false
+        end
+    end
+
     r_boundaries = r_boundary_info(Tuple(s for s ∈ inner_sections if s !== nothing),
-                                   Tuple(s for s ∈ outer_sections if s !== nothing))
+                                   Tuple(s for s ∈ outer_sections if s !== nothing),
+                                   ion_error_norm_skip_inner,
+                                   electron_error_norm_skip_inner,
+                                   neutral_error_norm_skip_inner,
+                                   ion_error_norm_skip_outer,
+                                   electron_error_norm_skip_outer,
+                                   neutral_error_norm_skip_outer)
 
     # Check that all boundary points are handled
     if r.irank == 0 && r.n > 1

moment_kinetics/src/moment_kinetics_structs.jl

@@ -550,6 +550,12 @@ struct r_boundary_info{Tinner <: NTuple{M,r_boundary_section} where M,
                        Touter <: NTuple{N,r_boundary_section} where N}
     inner_sections::Tinner
     outer_sections::Touter
+    ion_error_norm_skip_inner::MPISharedArray{Bool,1}
+    ion_error_norm_skip_outer::MPISharedArray{Bool,1}
+    electron_error_norm_skip_inner::MPISharedArray{Bool,1}
+    electron_error_norm_skip_outer::MPISharedArray{Bool,1}
+    neutral_error_norm_skip_inner::MPISharedArray{Bool,1}
+    neutral_error_norm_skip_outer::MPISharedArray{Bool,1}
 end
 
 struct z_boundary_info

moment_kinetics/src/runge_kutta.jl

@@ -935,50 +935,56 @@ function local_error_norm end
 
 function local_error_norm(f_loworder::MPISharedArray{mk_float,2},
                           f::MPISharedArray{mk_float,2}, rtol, atol; method="Linf",
-                          skip_r_inner=false, skip_z_lower=false, error_sum_zero=0.0)
+                          skip_r_inner, skip_r_outer, skip_z_lower=false,
+                          error_sum_zero=0.0)
+    nr = size(f, 2)
     if method == "Linf"
         f_max = -Inf
         @loop_r_z ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                continue
+            end
             error_norm = abs(f_loworder[iz,ir] - f[iz,ir]) / (rtol*abs(f[iz,ir]) + atol)
             f_max = max(f_max, error_norm)
         end
         return f_max
     elseif method == "L2"
         L2sum = error_sum_zero
         @loop_r_z ir iz begin
-            if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                    (skip_z_lower && iz == 1)
                 continue
             end
             error_norm = ((f_loworder[iz,ir] - f[iz,ir]) / (rtol*abs(f[iz,ir]) + atol))^2
             L2sum += error_norm
         end
         # Will sum results from different processes in shared memory block after returning
         # from this function.
-        nz, nr = size(f_loworder)
-        if skip_r_inner
-            nr -= 1
-        end
-        if skip_z_lower
-            nz -= 1
-        end
         return L2sum
     else
         error("Unrecognized method '$method'")
     end
 end
 function local_error_norm(f_loworder::MPISharedArray{mk_float,3},
                           f::MPISharedArray{mk_float,3}, rtol, atol, neutral=false;
-                          method="Linf", skip_r_inner=false, skip_z_lower=false,
+                          method="Linf", skip_r_inner, skip_r_outer, skip_z_lower=false,
                           error_sum_zero=0.0)
+    nr = size(f, 2)
     if method == "Linf"
         f_max = -Inf
         if neutral
             @loop_sn_r_z isn ir iz begin
+                if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                    continue
+                end
                 error_norm = abs(f_loworder[iz,ir,isn] - f[iz,ir,isn]) / (rtol*abs(f[iz,ir,isn]) + atol)
                 f_max = max(f_max, error_norm)
             end
         else
             @loop_s_r_z is ir iz begin
+                if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                    continue
+                end
                 error_norm = abs(f_loworder[iz,ir,is] - f[iz,ir,is]) / (rtol*abs(f[iz,ir,is]) + atol)
                 f_max = max(f_max, error_norm)
             end
@@ -988,15 +994,17 @@ function local_error_norm(f_loworder::MPISharedArray{mk_float,3},
         L2sum = error_sum_zero
         if neutral
             @loop_sn_r_z isn ir iz begin
-                if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+                if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                        (skip_z_lower && iz == 1)
                     continue
                 end
                 error_norm = ((f_loworder[iz,ir,isn] - f[iz,ir,isn]) / (rtol*abs(f[iz,ir,isn]) + atol))^2
                 L2sum += error_norm
             end
         else
             @loop_s_r_z is ir iz begin
-                if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+                if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                        (skip_z_lower && iz == 1)
                     continue
                 end
                 error_norm = ((f_loworder[iz,ir,is] - f[iz,ir,is]) / (rtol*abs(f[iz,ir,is]) + atol))^2
@@ -1005,108 +1013,117 @@ function local_error_norm(f_loworder::MPISharedArray{mk_float,3},
         end
         # Will sum results from different processes in shared memory block after returning
         # from this function.
-        nz, nr, nspecies = size(f_loworder)
-        if skip_r_inner
-            nr -= 1
-        end
-        if skip_z_lower
-            nz -= 1
-        end
         return L2sum
     else
         error("Unrecognized method '$method'")
     end
 end
 function local_error_norm(f_loworder::MPISharedArray{mk_float,4},
                           f::MPISharedArray{mk_float,4}, rtol, atol; method="Linf",
-                          skip_r_inner=false, skip_z_lower=false, error_sum_zero=0.0)
+                          skip_r_inner, skip_r_outer, skip_z_lower=false,
+                          error_sum_zero=0.0)
+    nr = size(f, 4)
     if method == "Linf"
         f_max = -Inf
-        @loop_r_z_vperp_vpa ir iz ivperp ivpa begin
-            error_norm = abs(f_loworder[ivpa,ivperp,iz,ir] - f[ivpa,ivperp,iz,ir]) /
-                         (rtol*abs(f[ivpa,ivperp,iz,ir]) + atol)
-            f_max = max(f_max, error_norm)
+        @loop_r_z ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                continue
+            end
+            @loop_vperp_vpa ivperp ivpa begin
+                error_norm = abs(f_loworder[ivpa,ivperp,iz,ir] - f[ivpa,ivperp,iz,ir]) /
+                             (rtol*abs(f[ivpa,ivperp,iz,ir]) + atol)
+                f_max = max(f_max, error_norm)
+            end
         end
         return f_max
     elseif method == "L2"
         L2sum = error_sum_zero
-        @loop_r_z_vperp_vpa ir iz ivperp ivpa begin
-            if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+        @loop_r_z ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                    (skip_z_lower && iz == 1)
                 continue
             end
-            error_norm = ((f_loworder[ivpa,ivperp,iz,ir] - f[ivpa,ivperp,iz,ir]) /
-                          (rtol*abs(f[ivpa,ivperp,iz,ir]) + atol))^2
-            L2sum += error_norm
+            @loop_vperp_vpa ivperp ivpa begin
+                error_norm = ((f_loworder[ivpa,ivperp,iz,ir] - f[ivpa,ivperp,iz,ir]) /
+                              (rtol*abs(f[ivpa,ivperp,iz,ir]) + atol))^2
+                L2sum += error_norm
+            end
         end
         # Will sum results from different processes in shared memory block after returning
         # from this function.
-        nvpa, nvperp, nz, nr = size(f_loworder)
-        if skip_r_inner
-            nr -= 1
-        end
-        if skip_z_lower
-            nz -= 1
-        end
         return L2sum
     else
         error("Unrecognized method '$method'")
     end
 end
 function local_error_norm(f_loworder::MPISharedArray{mk_float,5},
                           f::MPISharedArray{mk_float,5}, rtol, atol; method="Linf",
-                          skip_r_inner=false, skip_z_lower=false, error_sum_zero=0.0)
+                          skip_r_inner, skip_r_outer, skip_z_lower=false,
+                          error_sum_zero=0.0)
+    nr = size(f, 4)
     if method == "Linf"
         f_max = -Inf
-        @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
-            error_norm = abs(f_loworder[ivpa,ivperp,iz,ir,is] - f[ivpa,ivperp,iz,ir,is]) /
-                         (rtol*abs(f[ivpa,ivperp,iz,ir,is]) + atol)
-            f_max = max(f_max, error_norm)
+        @loop_s_r_z is ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                continue
+            end
+            @loop_vperp_vpa ivperp ivpa begin
+                error_norm = abs(f_loworder[ivpa,ivperp,iz,ir,is] - f[ivpa,ivperp,iz,ir,is]) /
+                             (rtol*abs(f[ivpa,ivperp,iz,ir,is]) + atol)
+                f_max = max(f_max, error_norm)
+            end
         end
         return f_max
     elseif method == "L2"
         L2sum = error_sum_zero
-        @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
-            if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+        @loop_s_r_z is ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                    (skip_z_lower && iz == 1)
                 continue
             end
-            error_norm = ((f_loworder[ivpa,ivperp,iz,ir,is] - f[ivpa,ivperp,iz,ir,is]) /
-                          (rtol*abs(f[ivpa,ivperp,iz,ir,is]) + atol))^2
-            L2sum += error_norm
+            @loop_vperp_vpa ivperp ivpa begin
+                error_norm = ((f_loworder[ivpa,ivperp,iz,ir,is] - f[ivpa,ivperp,iz,ir,is]) /
+                              (rtol*abs(f[ivpa,ivperp,iz,ir,is]) + atol))^2
+                L2sum += error_norm
+            end
         end
         # Will sum results from different processes in shared memory block after returning
         # from this function.
-        nvpa, nvperp, nz, nr, nspecies = size(f_loworder)
-        if skip_r_inner
-            nr -= 1
-        end
-        if skip_z_lower
-            nz -= 1
-        end
         return L2sum
     else
         error("Unrecognized method '$method'")
     end
 end
 function local_error_norm(f_loworder::MPISharedArray{mk_float,6},
                           f::MPISharedArray{mk_float,6}, rtol, atol; method="Linf",
-                          skip_r_inner=false, skip_z_lower=false, error_sum_zero=0.0)
+                          skip_r_inner, skip_r_outer, skip_z_lower=false,
+                          error_sum_zero=0.0)
+    nr = size(f, 5)
     if method == "Linf"
         f_max = -Inf
-        @loop_sn_r_z_vzeta_vr_vz isn ir iz ivzeta ivr ivz begin
-            error_norm = abs(f_loworder[ivz,ivr,ivzeta,iz,ir,isn] - f[ivz,ivr,ivzeta,iz,ir,isn]) /
-                         (rtol*abs(f[ivz,ivr,ivzeta,iz,ir,isn]) + atol)
-            f_max = max(f_max, error_norm)
+        @loop_sn_r_z isn ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr)
+                continue
+            end
+            @loop_vzeta_vr_vz ivzeta ivr ivz begin
+                error_norm = abs(f_loworder[ivz,ivr,ivzeta,iz,ir,isn] - f[ivz,ivr,ivzeta,iz,ir,isn]) /
+                             (rtol*abs(f[ivz,ivr,ivzeta,iz,ir,isn]) + atol)
+                f_max = max(f_max, error_norm)
+            end
         end
         return f_max
     elseif method == "L2"
         L2sum = error_sum_zero
-        @loop_sn_r_z_vzeta_vr_vz isn ir iz ivzeta ivr ivz begin
-            if (skip_r_inner && ir == 1) || (skip_z_lower && iz == 1)
+        @loop_sn_r_z isn ir iz begin
+            if (skip_r_inner[iz] && ir == 1) || (skip_r_outer[iz] && ir == nr) ||
+                    (skip_z_lower && iz == 1)
                 continue
             end
-            error_norm = ((f_loworder[ivz,ivr,ivzeta,iz,ir,isn] - f_loworder[ivz,ivr,ivzeta,iz,ir,isn]) /
-                          (rtol*abs(f[ivz,ivr,ivzeta,iz,ir,isn]) + atol))^2
-            L2sum += error_norm
+            @loop_vzeta_vr_vz ivzeta ivr ivz begin
+                error_norm = ((f_loworder[ivz,ivr,ivzeta,iz,ir,isn] - f_loworder[ivz,ivr,ivzeta,iz,ir,isn]) /
+                              (rtol*abs(f[ivz,ivr,ivzeta,iz,ir,isn]) + atol))^2
+                L2sum += error_norm
+            end
         end
         # Will sum results from different processes in shared memory block after returning
         # from this function.