Merge pull request #297 from mabarnes/ADI-solver

'Alternating direction implicit' preconditioner

machines/generic-batch-template/compile_dependencies.sh

@@ -77,10 +77,10 @@ if [[ $BUILDHDF5 == "y" && -d hdf5-build ]]; then
 fi
 
 if [[ $BUILDHDF5 == "y" ]]; then
-  HDF5=hdf5-1.14.3
+  HDF5=hdf5-1.14.5
   # Download and extract the source
-  wget -O ${HDF5}.tar.bz2 https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.14/hdf5-1.14.3/src/hdf5-1.14.3.tar.bz2
-  tar xjf ${HDF5}.tar.bz2
+  wget -O ${HDF5}.tar.gz https://support.hdfgroup.org/releases/hdf5/v1_14/v1_14_5/downloads/hdf5-1.14.5.tar.gz
+  tar xjf ${HDF5}.tar.gz
 
   cd $HDF5
 

machines/generic-pc/compile_dependencies.sh

@@ -77,10 +77,10 @@ else
 fi
 
 if [[ $BUILDHDF5 == "y" ]]; then
-  HDF5=hdf5-1.14.3
+  HDF5=hdf5-1.14.5
   # Download and extract the source
-  wget -O ${HDF5}.tar.bz2 https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.14/hdf5-1.14.3/src/hdf5-1.14.3.tar.bz2
-  tar xjf ${HDF5}.tar.bz2
+  wget -O ${HDF5}.tar.gz https://support.hdfgroup.org/releases/hdf5/v1_14/v1_14_5/downloads/hdf5-1.14.5.tar.gz
+  tar xzf ${HDF5}.tar.gz
 
   cd $HDF5
 

machines/marconi/compile_dependencies.sh

@@ -30,10 +30,10 @@ if [ -d hdf5-build ]; then
 fi
 
 if [ $BUILDHDF5 -eq 0 ]; then
-  HDF5=hdf5-1.14.3
+  HDF5=hdf5-1.14.5
   # Download and extract the source
-  wget -O ${HDF5}.tar.bz2 https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.14/hdf5-1.14.3/src/hdf5-1.14.3.tar.bz2
-  tar xjf ${HDF5}.tar.bz2
+  wget -O ${HDF5}.tar.gz https://support.hdfgroup.org/releases/hdf5/v1_14/v1_14_5/downloads/hdf5-1.14.5.tar.gz
+  tar xjf ${HDF5}.tar.gz
 
   cd $HDF5
 

makie_post_processing/makie_post_processing/src/makie_post_processing.jl

@@ -8079,8 +8079,10 @@ function timestep_diagnostics(run_info, run_info_dfns; plot_prefix=nothing, it=n
                     for p ∈ nl_prefixes
                         nonlinear_iterations = get_variable(ri, "$(p)_nonlinear_iterations_per_solve")
                         linear_iterations = get_variable(ri, "$(p)_linear_iterations_per_nonlinear_iteration")
+                        precon_iterations = get_variable(ri, "$(p)_precon_iterations_per_linear_iteration")
                         plot_1d(time, nonlinear_iterations, label=prefix * " " * p * " NL per solve", ax=ax)
                         plot_1d(time, linear_iterations, label=prefix * " " * p * " L per NL", ax=ax)
+                        plot_1d(time, precon_iterations, label=prefix * " " * p * " P per L", ax=ax)
                     end
                 end
             end

moment_kinetics/src/analysis.jl

@@ -15,7 +15,7 @@ using ..interpolation: interpolate_to_grid_1d
 using ..load_data: open_readonly_output_file, get_nranks, load_pdf_data, load_rank_data
 using ..load_data: load_distributed_ion_pdf_slice
 using ..looping
-using ..type_definitions: mk_int
+using ..type_definitions: mk_int, mk_float
 using ..velocity_moments: integrate_over_vspace
 
 using FFTW
@@ -595,7 +595,8 @@ const default_epsilon = 1.0e-4
 
 """
     steady_state_residuals(variable, variable_at_previous_time, dt;
-                           epsilon=$default_epsilon, use_mpi=false)
+                           epsilon=$default_epsilon, use_mpi=false,
+                           only_max_abs=false)
 
 Calculate how close a variable is to steady state.
 
@@ -630,33 +631,36 @@ initialised, and that `variable` has r and z dimensions but no species dimension
 distributed-memory MPI, this routine will double-count the points on block boundaries.
 
 If `only_max_abs=true` is passed, then only calculate the 'maxium absolute residual'. In
-this case the OrderedDict returned will have only one entry, for `"max absolute
-residual"`.
+this case just returns the "max absolute residual", not an OrderedDict.
 """
 function steady_state_residuals(variable, variable_at_previous_time, dt;
                                 epsilon=default_epsilon, use_mpi=false,
                                 only_max_abs=false)
+    return steady_state_residuals(variable, variable_at_previous_time, dt, use_mpi,
+                                  only_max_abs, epsilon)
+end
+function steady_state_residuals(variable, variable_at_previous_time, dt, use_mpi,
+                                only_max_abs=false, epsilon=default_epsilon)
     square_residual_norms =
-        steady_state_square_residuals(variable, variable_at_previous_time, dt;
-                                      epsilon=epsilon, use_mpi=use_mpi,
-                                      only_max_abs=only_max_abs)
+        steady_state_square_residuals(variable, variable_at_previous_time, dt, nothing,
+                                      use_mpi, only_max_abs, epsilon)
     if global_rank[] == 0
         if only_max_abs
             # In this case as an optimisation the residual was not squared, so do not need
             # to square-root here
             return square_residual_norms
         else
-            return OrderedDict(k=>sqrt.(v) for (k,v) ∈ square_residual_norms)
+            return OrderedDict{String,Vector{mk_float}}(k=>sqrt.(v) for (k,v) ∈ square_residual_norms)
         end
     else
         return nothing
     end
 end
 
 """
-    steady_state_square_residuals(variable, variable_at_previous_time, dt;
-                                  variable_max=nothing, epsilon=1.0e-4,
-                                  use_mpi=false, only_max_abs=false)
+    steady_state_square_residuals(variable, variable_at_previous_time, dt,
+                                  variable_max=nothing, use_mpi=false,
+                                  only_max_abs=false, epsilon=$default_epsilon)
 
 Used to calculate the mean square residual for [`steady_state_residuals`](@ref).
 
@@ -668,9 +672,9 @@ See [`steady_state_residuals`](@ref) for documenation of the other arguments. Th
 values of [`steady_state_residuals`](@ref) are the square-root of the return values of
 this function.
 """
-function steady_state_square_residuals(variable, variable_at_previous_time, dt;
-                                       variable_max=nothing, epsilon=default_epsilon,
-                                       use_mpi=false, only_max_abs=false)
+function steady_state_square_residuals(variable, variable_at_previous_time, dt,
+                                       variable_max=nothing, use_mpi=false,
+                                       only_max_abs=false, epsilon=default_epsilon)
     if ndims(dt) == 0
         t_dim = ndims(variable) + 1
     else
@@ -797,10 +801,9 @@ function steady_state_square_residuals(variable, variable_at_previous_time, dt;
                                              (size(packed_results)..., n_blocks[]))
 
             if only_max_abs
-                return OrderedDict(
-                           "max absolute residual"=>maximum(gathered_block_results, dims=2))
+                return maximum(gathered_block_results, dims=2)
             else
-                return OrderedDict(
+                return OrderedDict{String,mk_float}(
                            "RMS absolute residual"=>mean(@view(gathered_block_results[:,1,:]), dims=2),
                            "max absolute residual"=>maximum(@view(gathered_block_results[:,2,:]), dims=2),
                            "RMS relative residual"=>mean(@view(gathered_block_results[:,3,:]), dims=2),
@@ -817,21 +820,19 @@ function steady_state_square_residuals(variable, variable_at_previous_time, dt;
 
         if only_max_abs
             absolute_residual =
-                _steady_state_residual(variable, variable_at_previous_time, reshaped_dt)
+                _steady_state_absolute_residual(variable, variable_at_previous_time, reshaped_dt)
             # Need to wrap the maximum(...) in a call to vec(...) so that we return a
             # Vector, not an N-dimensional array where the first (N-1) dimensions all have
             # size 1.
-            return OrderedDict(
-                       "max absolute residual"=>vec(maximum(absolute_residual;
-                                                            dims=tuple((1:t_dim-1)...))))
+            return vec(maximum(absolute_residual; dims=tuple((1:t_dim-1)...)))
         else
             absolute_square_residual, relative_square_residual =
                 _steady_state_square_residual(variable, variable_at_previous_time,
                                               reshaped_dt, epsilon, variable_max)
             # Need to wrap the mean(...) or maximum(...) in a call to vec(...) so that we
             # return a Vector, not an N-dimensional array where the first (N-1) dimensions all
             # have size 1.
-            return OrderedDict(
+            return OrderedDict{String,Vector{mk_float}}(
                        "RMS absolute residual"=>vec(mean(absolute_square_residual;
                                                          dims=tuple((1:t_dim-1)...))),
                        "max absolute residual"=>vec(maximum(absolute_square_residual;

moment_kinetics/src/boundary_conditions.jl

@@ -1029,7 +1029,7 @@ function enforce_v_boundary_condition_local!(f, bc, speed, v_diffusion, v, v_spe
 
         D0 = v_spectral.lobatto.Dmat[end,:]
         # adjust F(vpa = L/2) so that d F / d vpa = 0 at vpa = L/2
-        f[end] = -sum(D0[1:ngrid-1].*f[end-v.ngrid+1:end-1])/D0[v.ngrid]
+        f[end] = -sum(D0[1:v.ngrid-1].*f[end-v.ngrid+1:end-1])/D0[v.ngrid]
     elseif bc == "periodic"
         f[1] = 0.5*(f[1]+f[end])
         f[end] = f[1]
@@ -1038,6 +1038,7 @@ function enforce_v_boundary_condition_local!(f, bc, speed, v_diffusion, v, v_spe
     else
         error("Unsupported boundary condition option '$bc' for $(v.name)")
     end
+    return nothing
 end
 
 """

moment_kinetics/src/coordinates.jl

@@ -381,9 +381,9 @@ function define_coordinate(coord_input::NamedTuple; parallel_io::Bool=false,
     end
 
     coord = coordinate(coord_input.name, n_global, n_local, coord_input.ngrid,
-        coord_input.nelement, coord_input.nelement_local, nrank, irank, coord_input.L,
-        grid, cell_width, igrid, ielement, imin, imax, igrid_full,
-        coord_input.discretization, coord_input.finite_difference_option,
+        coord_input.nelement, coord_input.nelement_local, nrank, irank,
+        mk_float(coord_input.L), grid, cell_width, igrid, ielement, imin, imax,
+        igrid_full, coord_input.discretization, coord_input.finite_difference_option,
         coord_input.cheb_option, coord_input.bc, coord_input.boundary_parameters, wgts,
         uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         copy(scratch), copy(scratch), copy(scratch), copy(scratch), copy(scratch),