From 8f120810bc63c8e530afb1b6cc07792f6d877fff Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 8 Sep 2023 17:03:44 +0100 Subject: [PATCH 01/23] Attempt to start making a nonuniformly spaced elemental grid --- src/coordinates.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/coordinates.jl b/src/coordinates.jl index 704f67023..7e4cb04f5 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -4,6 +4,7 @@ module coordinates export define_coordinate, write_coordinate export equally_spaced_grid +export set_element_boundaries using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float, allocate_int @@ -84,6 +85,12 @@ struct coordinate local_io_range::UnitRange{Int64} # global range to write into in output file global_io_range::UnitRange{Int64} + # scale for each element + element_scale::Array{mk_float,1} + # shift for each element + element_shift::Array{mk_float,1} + # boundaries for each element + element_boundaries::Array{mk_float,1} end """ @@ -105,6 +112,14 @@ function define_coordinate(input, parallel_io::Bool=false) # obtain (local) index mapping from the grid within each element # to the full grid imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local) + # initialise the data used to construct the grid + # boundaries for each element + element_boundaries = allocate_float(nelement_global) + element_boundaries_local = allocate_float(nelement_local) + + element_scale = allocate_float(nelement_local) + # shift for each element + element_shift = allocate_float(nelement_local) # initialize the grid and the integration weights associated with the grid # also obtain the Chebyshev theta grid and spacing if chosen as discretization option grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global, @@ -126,7 +141,6 @@ function define_coordinate(input, parallel_io::Bool=false) # endpoints, so only two pieces of information must be shared send_buffer = allocate_float(1) receive_buffer = allocate_float(1) - # Add some ranges to support parallel file io if !parallel_io # No parallel io, just write everything @@ -149,7 +163,7 @@ function define_coordinate(input, parallel_io::Bool=false) cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option, input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm, - local_io_range, global_io_range) + local_io_range, global_io_range, element_scale, element_shift, element_boundaries) if input.discretization == "chebyshev_pseudospectral" && coord.n > 1 # create arrays needed for explicit Chebyshev pseudospectral treatment in this @@ -167,6 +181,40 @@ function define_coordinate(input, parallel_io::Bool=false) return coord, spectral end +function set_element_boundaries(nelement_global, L, element_spacing_option) + # set global element boundaries + element_boundaries = allocate_float(nelement_global+1) + if element_spacing_option == "sqrt" + # number of boundaries of sqrt grid + nsqrt = floor(mk_int,(nelement_global+1)/2) + println("nsqrt",nsqrt) + # number of boundaries of uniform grid + nuniform = nelement_global + 3 - 2*nsqrt + println("nuniform",nuniform) + DL = L/6.0 # 1/3 of the domain is uniformly spaced + delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section + for j in 1:nsqrt + element_boundaries[j] = -(L/2.0) + ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2 + end + println(element_boundaries) + for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt + element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) + end + println(element_boundaries) + for j in 1:nsqrt + element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2 + end + println(element_boundaries) + elseif element_spacing_option == "uniform" # uniform spacing + for j in 1:nelement_global+1 + element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5) + end + else + println("ERROR: element_spacing_option: ",element_spacing_option, " not supported") + end + return element_boundaries +end + """ setup a grid with n_global grid points on the interval [-L/2,L/2] """ From 5475949b8c4602e2d0ed9c70d0267e2699fac32d Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 8 Sep 2023 20:17:35 +0100 Subject: [PATCH 02/23] Version of element boundary routine working for nelement >= 2 --- src/coordinates.jl | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/coordinates.jl b/src/coordinates.jl index 7e4cb04f5..d38a8e8af 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -114,12 +114,12 @@ function define_coordinate(input, parallel_io::Bool=false) imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local) # initialise the data used to construct the grid # boundaries for each element - element_boundaries = allocate_float(nelement_global) - element_boundaries_local = allocate_float(nelement_local) + element_boundaries = allocate_float(input.nelement_global) + element_boundaries_local = allocate_float(input.nelement_local) - element_scale = allocate_float(nelement_local) + element_scale = allocate_float(input.nelement_local) # shift for each element - element_shift = allocate_float(nelement_local) + element_shift = allocate_float(input.nelement_local) # initialize the grid and the integration weights associated with the grid # also obtain the Chebyshev theta grid and spacing if chosen as discretization option grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global, @@ -186,23 +186,39 @@ function set_element_boundaries(nelement_global, L, element_spacing_option) element_boundaries = allocate_float(nelement_global+1) if element_spacing_option == "sqrt" # number of boundaries of sqrt grid - nsqrt = floor(mk_int,(nelement_global+1)/2) + nsqrt = floor(mk_int,(nelement_global)/2) + 1 + if nelement_global%2 > 0 # odd + #fac = 0.0 + #for j in 2:nsqrt + # fac += 2.0*(((j-2)/(nsqrt-1))^2 - ((j-1)/(nsqrt-1))^2 ) + #end + x = ((nsqrt-2)/(nsqrt-1))^2 + if nsqrt < 3 + fac = 2.0/3.0 + else + #fac = 0.5*(sqrt( 1.0 + 4.0*x) - 1.0)/x + fac = 1.0/( 3.0/2.0 - 0.5*((nsqrt-2)/(nsqrt-1))^2) + end + else + fac = 1.0 + end + println("nsqrt",nsqrt) # number of boundaries of uniform grid - nuniform = nelement_global + 3 - 2*nsqrt - println("nuniform",nuniform) - DL = L/6.0 # 1/3 of the domain is uniformly spaced - delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section + #nuniform = nelement_global + 3 - 2*nsqrt + #println("nuniform",nuniform) + #DL = L/6.0 # 1/3 of the domain is uniformly spaced + #delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section for j in 1:nsqrt - element_boundaries[j] = -(L/2.0) + ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2 + element_boundaries[j] = -(L/2.0) + fac*(L/2.0)*((j-1)/(nsqrt-1))^2 end println(element_boundaries) - for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt - element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) - end + #for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt + # element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) + #end println(element_boundaries) for j in 1:nsqrt - element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2 + element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2 end println(element_boundaries) elseif element_spacing_option == "uniform" # uniform spacing From e7788087bc52a346916e2386577c47c0c80ed53a Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 11 Sep 2023 09:22:00 +0100 Subject: [PATCH 03/23] Merge the new method of element boundary specification into the latest master branch. Use the feature by specifying the option z_element_spacing_option = "sqrt". The default value is "uniform", which reproduces previous behaviour. The new method of specifying element boundaries gives the potential to optimize the element spacing for new applications in the future. --- src/chebyshev.jl | 30 ++++++++------- src/coordinates.jl | 56 ++++++++++++++-------------- src/file_io.jl | 4 ++ src/input_structs.jl | 4 ++ src/load_data.jl | 4 +- src/moment_kinetics_input.jl | 59 +++++++++++++++++++----------- src/post_processing.jl | 5 ++- test/calculus_tests.jl | 56 +++++++++++++++++++--------- test/interpolation_tests.jl | 6 ++- test/nonlinear_sound_wave_tests.jl | 7 ++-- test/wall_bc_tests.jl | 5 ++- 11 files changed, 145 insertions(+), 91 deletions(-) diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 0908ac98e..854cdf7a5 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -42,8 +42,8 @@ end """ initialize chebyshev grid scaled to interval [-box_length/2, box_length/2] """ -function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n, - irank, box_length, imin, imax) +function scaled_chebyshev_grid(ngrid, nelement_local, n, + element_scale, element_shift, imin, imax) # initialize chebyshev grid defined on [1,-1] # with n grid points chosen to facilitate # the fast Chebyshev transform (aka the discrete cosine transform) @@ -55,7 +55,7 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n, # setup the scale factor by which the Chebyshev grid on [-1,1] # is to be multiplied to account for the full domain [-L/2,L/2] # and the splitting into nelement elements with ngrid grid points - scale_factor = 0.5*box_length/float(nelement_global) + #scale_factor = 0.5*box_length/float(nelement_global) # account for the fact that the minimum index needed for the chebyshev_grid # within each element changes from 1 to 2 in going from the first element # to the remaining elements @@ -63,8 +63,10 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n, @inbounds for j ∈ 1:nelement_local #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor # amount by which to shift the centre of this element from zero - iel_global = j + irank*nelement_local - shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5) + #iel_global = j + irank*nelement_local + #shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5) + scale_factor = element_scale[j] + shift = element_shift[j] # reverse the order of the original chebyshev_grid (ran from [1,-1]) # and apply the scale factor and shift grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift @@ -72,7 +74,7 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n, # to avoid double-counting boundary element k = 2 end - wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor) + wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale) return grid, wgts end @@ -90,7 +92,7 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info) @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df)) # note that one must multiply by 2*nelement/L to get derivative # in scaled coordinate - scale_factor = 2.0*float(coord.nelement_global)/coord.L + #scale_factor = 2.0*float(coord.nelement_global)/coord.L # scale factor is (length of a single element/2)^{-1} # variable k will be used to avoid double counting of overlapping point @@ -112,7 +114,7 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info) # and multiply by scaling factor needed to go # from Chebyshev z coordinate to actual z for i ∈ 1:coord.ngrid - df[i,j] *= scale_factor + df[i,j] /= coord.element_scale[j] end k = 1 end @@ -323,23 +325,25 @@ end returns wgts array containing the integration weights associated with all grid points for Clenshaw-Curtis quadrature """ -function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor) +function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale) # create array containing the integration weights wgts = zeros(mk_float, n) # calculate the modified Chebshev moments of the first kind μ = chebyshevmoments(ngrid) + # calculate the raw weights for a normalised grid on [-1,1] + w = clenshawcurtisweights(μ) @inbounds begin # calculate the weights within a single element and # scale to account for modified domain (not [-1,1]) - wgts[1:ngrid] = clenshawcurtisweights(μ)*scale_factor + wgts[1:ngrid] = w*element_scale[1] if nelement_local > 1 # account for double-counting of points at inner element boundaries - wgts[ngrid] *= 2 + #wgts[ngrid] *= 2 for j ∈ 2:nelement_local - wgts[imin[j]:imax[j]] .= wgts[2:ngrid] + wgts[imin[j]-1:imax[j]] .+= w*element_scale[j] end # remove double-counting of outer element boundary for last element - wgts[n] *= 0.5 + #wgts[n] *= 0.5 end end return wgts diff --git a/src/coordinates.jl b/src/coordinates.jl index d38a8e8af..b44a7b2ec 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -4,7 +4,9 @@ module coordinates export define_coordinate, write_coordinate export equally_spaced_grid +# testing export set_element_boundaries +export init_grid using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float, allocate_int @@ -90,7 +92,9 @@ struct coordinate # shift for each element element_shift::Array{mk_float,1} # boundaries for each element - element_boundaries::Array{mk_float,1} + #element_boundaries::Array{mk_float,1} + # option used to set up element spacing + element_spacing_option::String end """ @@ -114,16 +118,12 @@ function define_coordinate(input, parallel_io::Bool=false) imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local) # initialise the data used to construct the grid # boundaries for each element - element_boundaries = allocate_float(input.nelement_global) - element_boundaries_local = allocate_float(input.nelement_local) - - element_scale = allocate_float(input.nelement_local) - # shift for each element - element_shift = allocate_float(input.nelement_local) + element_boundaries = set_element_boundaries(input.nelement_global, input.L, input.element_spacing_option) + # shift and scale factors for each local element + element_scale, element_shift = set_element_scale_and_shift(input.nelement_global, input.nelement_local, input.irank, element_boundaries) # initialize the grid and the integration weights associated with the grid # also obtain the Chebyshev theta grid and spacing if chosen as discretization option - grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global, - input.nelement_local, n_global, n_local, input.irank, input.L, + grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_local, n_global, n_local, input.irank, input.L, element_scale, element_shift, imin, imax, igrid, input.discretization, input.name) # calculate the widths of the cells between neighboring grid points cell_width = grid_spacing(grid, n_local) @@ -163,7 +163,7 @@ function define_coordinate(input, parallel_io::Bool=false) cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option, input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm, - local_io_range, global_io_range, element_scale, element_shift, element_boundaries) + local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)#, element_boundaries) if input.discretization == "chebyshev_pseudospectral" && coord.n > 1 # create arrays needed for explicit Chebyshev pseudospectral treatment in this @@ -184,7 +184,7 @@ end function set_element_boundaries(nelement_global, L, element_spacing_option) # set global element boundaries element_boundaries = allocate_float(nelement_global+1) - if element_spacing_option == "sqrt" + if element_spacing_option == "sqrt" && nelement_global > 3 # number of boundaries of sqrt grid nsqrt = floor(mk_int,(nelement_global)/2) + 1 if nelement_global%2 > 0 # odd @@ -203,25 +203,14 @@ function set_element_boundaries(nelement_global, L, element_spacing_option) fac = 1.0 end - println("nsqrt",nsqrt) - # number of boundaries of uniform grid - #nuniform = nelement_global + 3 - 2*nsqrt - #println("nuniform",nuniform) - #DL = L/6.0 # 1/3 of the domain is uniformly spaced - #delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section for j in 1:nsqrt element_boundaries[j] = -(L/2.0) + fac*(L/2.0)*((j-1)/(nsqrt-1))^2 end - println(element_boundaries) - #for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt - # element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) - #end - println(element_boundaries) for j in 1:nsqrt element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2 end - println(element_boundaries) - elseif element_spacing_option == "uniform" # uniform spacing + + elseif element_spacing_option == "uniform" || nelement_global < 4 # uniform spacing for j in 1:nelement_global+1 element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5) end @@ -231,10 +220,23 @@ function set_element_boundaries(nelement_global, L, element_spacing_option) return element_boundaries end +function set_element_scale_and_shift(nelement_global, nelement_local, irank, element_boundaries) + element_scale = allocate_float(nelement_local) + element_shift = allocate_float(nelement_local) + + for j in 1:nelement_local + iel_global = j + irank*nelement_local + upper_boundary = element_boundaries[iel_global+1] + lower_boundary = element_boundaries[iel_global] + element_scale[j] = 0.5*(upper_boundary-lower_boundary) + element_shift[j] = 0.5*(upper_boundary+lower_boundary) + end + return element_scale, element_shift +end """ setup a grid with n_global grid points on the interval [-L/2,L/2] """ -function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, irank, L, +function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_scale, element_shift, imin, imax, igrid, discretization, name) uniform_grid = equally_spaced_grid(n_global, n_local, irank, L) uniform_grid_shifted = equally_spaced_grid_shifted(n_global, n_local, irank, L) @@ -251,7 +253,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir elseif discretization == "chebyshev_pseudospectral" if name == "vperp" # initialize chebyshev grid defined on [-L/2,L/2] - grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax) + grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax) grid .= grid .+ L/2.0 # shift to [0,L] appropriate to vperp variable wgts = 2.0 .* wgts .* grid # to include 2 vperp in jacobian of integral # see note above on normalisation @@ -262,7 +264,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir # needed to obtain Chebyshev spectral coefficients # 'wgts' are the integration weights attached to each grid points # that are those associated with Clenshaw-Curtis quadrature - grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax) + grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax) end elseif discretization == "finite_difference" if name == "vperp" diff --git a/src/file_io.jl b/src/file_io.jl index bb53a90f2..5dd229843 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -521,6 +521,10 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_ write_single_value!(group, "bc", coord.bc; parallel_io=parallel_io, description="boundary condition for $coord_name") + # write the boundary condition for the coordinate + write_single_value!(group, "element_spacing_option", coord.element_spacing_option; parallel_io=parallel_io, + description="element_spacing_option for $coord_name") + return group end diff --git a/src/input_structs.jl b/src/input_structs.jl index b013f703b..a08168e58 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -127,6 +127,8 @@ mutable struct grid_input_mutable bc::String # mutable struct containing advection speed options advection::advection_input_mutable + # string option determining boundary spacing + element_spacing_option::String end """ @@ -156,6 +158,8 @@ struct grid_input advection::advection_input # MPI communicator comm::MPI.Comm + # string option determining boundary spacing + element_spacing_option::String end """ diff --git a/src/load_data.jl b/src/load_data.jl index 206873bbb..5507965b1 100644 --- a/src/load_data.jl +++ b/src/load_data.jl @@ -213,6 +213,7 @@ function load_coordinate_data(fid, name; printout=false) discretization = load_variable(coord_group, "discretization") fd_option = load_variable(coord_group, "fd_option") bc = load_variable(coord_group, "bc") + element_spacing_option = load_variable(coord_group, "element_spacing_option") nelement_local = nothing if n_local == 1 && ngrid == 1 @@ -225,11 +226,10 @@ function load_coordinate_data(fid, name; printout=false) else nelement_global = (n_global-1) ÷ (ngrid-1) end - # Define input to create coordinate struct input = grid_input(name, ngrid, nelement_global, nelement_local, nrank, irank, L, discretization, fd_option, bc, advection_input("", 0.0, 0.0, 0.0), - MPI.COMM_NULL) + MPI.COMM_NULL, element_spacing_option) coord, spectral = define_coordinate(input) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 6663d6e11..cf1cf4546 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -340,7 +340,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the boundary condition to impose in r # supported options are "periodic" and "Dirichlet" r.bc = get(scan_input, "r_bc", "periodic") - + r.element_spacing_option = get(scan_input, "r_element_spacing_option", "uniform") + # overwrite some default parameters related to the z grid # ngrid is number of grid points per element z.ngrid = get(scan_input, "z_ngrid", 9) @@ -357,7 +358,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the boundary condition to impose in z # supported options are "constant", "periodic" and "wall" z.bc = get(scan_input, "z_bc", "wall") - + z.element_spacing_option = get(scan_input, "z_element_spacing_option", "uniform") + # overwrite some default parameters related to the vpa grid # ngrid is the number of grid points per element vpa.ngrid = get(scan_input, "vpa_ngrid", 17) @@ -374,7 +376,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # supported options are "chebyshev_pseudospectral" and "finite_difference" vpa.discretization = get(scan_input, "vpa_discretization", "chebyshev_pseudospectral") vpa.fd_option = get(scan_input, "vpa_finite_difference_option", "third_order_upwind") - + vpa.element_spacing_option = get(scan_input, "vpa_element_spacing_option", "uniform") + num_diss_params = setup_numerical_dissipation( get(scan_input, "numerical_dissipation", Dict{String,Any}()), true) @@ -395,7 +398,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the discretization option for the vperp grid # supported options are "finite_difference_vperp" "chebyshev_pseudospectral_vperp" vperp.discretization = get(scan_input, "vperp_discretization", "chebyshev_pseudospectral_vperp") - + vperp.element_spacing_option = get(scan_input, "vperp_element_spacing_option", "uniform") + # overwrite some default parameters related to the gyrophase grid # ngrid is the number of grid points per element gyrophase.ngrid = get(scan_input, "gyrophase_ngrid", 17) @@ -421,7 +425,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the discretization option for the vz grid # supported options are "chebyshev_pseudospectral" and "finite_difference" vz.discretization = get(scan_input, "vz_discretization", vpa.discretization) - + vz.element_spacing_option = get(scan_input, "vz_element_spacing_option", "uniform") + # overwrite some default parameters related to the vr grid # ngrid is the number of grid points per element vr.ngrid = get(scan_input, "vr_ngrid", 1) @@ -437,7 +442,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the discretization option for the vr grid # supported options are "chebyshev_pseudospectral" and "finite_difference" vr.discretization = get(scan_input, "vr_discretization", "chebyshev_pseudospectral") - + vr.element_spacing_option = get(scan_input, "vr_element_spacing_option", "uniform") + # overwrite some default parameters related to the vzeta grid # ngrid is the number of grid points per element vzeta.ngrid = get(scan_input, "vzeta_ngrid", 1) @@ -453,6 +459,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # determine the discretization option for the vzeta grid # supported options are "chebyshev_pseudospectral" and "finite_difference" vzeta.discretization = get(scan_input, "vzeta_discretization", "chebyshev_pseudospectral") + vzeta.element_spacing_option = get(scan_input, "vzeta_element_spacing_option", "uniform") end is_1V = (vperp.ngrid == vperp.nelement_global == 1 && vzeta.ngrid == @@ -489,37 +496,37 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) z_advection_immutable = advection_input(z.advection.option, z.advection.constant_speed, z.advection.frequency, z.advection.oscillation_amplitude) z_immutable = grid_input("z", z.ngrid, z.nelement_global, z.nelement_local, nrank_z, irank_z, z.L, - z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z) + z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z, z.element_spacing_option) r_advection_immutable = advection_input(r.advection.option, r.advection.constant_speed, r.advection.frequency, r.advection.oscillation_amplitude) r_immutable = grid_input("r", r.ngrid, r.nelement_global, r.nelement_local, nrank_r, irank_r, r.L, - r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r) + r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r, r.element_spacing_option) # for dimensions below which do not currently use distributed-memory MPI # assign dummy values to nrank, irank and comm of coord struct vpa_advection_immutable = advection_input(vpa.advection.option, vpa.advection.constant_speed, vpa.advection.frequency, vpa.advection.oscillation_amplitude) vpa_immutable = grid_input("vpa", vpa.ngrid, vpa.nelement_global, vpa.nelement_local, 1, 0, vpa.L, - vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL) + vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL, vpa.element_spacing_option) vperp_advection_immutable = advection_input(vperp.advection.option, vperp.advection.constant_speed, vperp.advection.frequency, vperp.advection.oscillation_amplitude) vperp_immutable = grid_input("vperp", vperp.ngrid, vperp.nelement_global, vperp.nelement_local, 1, 0, vperp.L, - vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL) + vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL, vperp.element_spacing_option) gyrophase_advection_immutable = advection_input(gyrophase.advection.option, gyrophase.advection.constant_speed, gyrophase.advection.frequency, gyrophase.advection.oscillation_amplitude) gyrophase_immutable = grid_input("gyrophase", gyrophase.ngrid, gyrophase.nelement_global, gyrophase.nelement_local, 1, 0, gyrophase.L, - gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL) + gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL, gyrophase.element_spacing_option) vz_advection_immutable = advection_input(vz.advection.option, vz.advection.constant_speed, vz.advection.frequency, vz.advection.oscillation_amplitude) vz_immutable = grid_input("vz", vz.ngrid, vz.nelement_global, vz.nelement_local, 1, 0, vz.L, - vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL) + vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL, vz.element_spacing_option) vr_advection_immutable = advection_input(vr.advection.option, vr.advection.constant_speed, vr.advection.frequency, vr.advection.oscillation_amplitude) vr_immutable = grid_input("vr", vr.ngrid, vr.nelement_global, vr.nelement_local, 1, 0, vr.L, - vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL) + vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL, vr.element_spacing_option) vzeta_advection_immutable = advection_input(vzeta.advection.option, vzeta.advection.constant_speed, vzeta.advection.frequency, vzeta.advection.oscillation_amplitude) vzeta_immutable = grid_input("vzeta", vzeta.ngrid, vzeta.nelement_global, vzeta.nelement_local, 1, 0, vzeta.L, - vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL) + vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL, vzeta.element_spacing_option) species_charged_immutable = Array{species_parameters,1}(undef,n_ion_species) species_neutral_immutable = Array{species_parameters,1}(undef,n_neutral_species) @@ -660,10 +667,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_z = advection_input_mutable(advection_option_z, advection_speed_z, frequency_z, oscillation_amplitude_z) + element_spacing_option_z = "uniform" # create a mutable structure containing the input info related to the z grid z = grid_input_mutable("z", ngrid_z, nelement_global_z, nelement_local_z, L_z, discretization_option_z, finite_difference_option_z, boundary_option_z, - advection_z) + advection_z, element_spacing_option_z) #################### parameters related to the r grid ###################### # ngrid_r is number of grid points per element ngrid_r = 1 @@ -700,10 +708,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for r advection_r = advection_input_mutable(advection_option_r, advection_speed_r, frequency_r, oscillation_amplitude_r) + element_spacing_option_r = "uniform" # create a mutable structure containing the input info related to the r grid r = grid_input_mutable("r", ngrid_r, nelement_global_r, nelement_local_r, L_r, discretization_option_r, finite_difference_option_r, boundary_option_r, - advection_r) + advection_r, element_spacing_option_r) ############################################################################ ################### parameters related to the vpa grid ##################### # ngrid_vpa is the number of grid points per element @@ -738,10 +747,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_vpa = advection_input_mutable(advection_option_vpa, advection_speed_vpa, frequency_vpa, oscillation_amplitude_vpa) + element_spacing_option_vpa = "uniform" # create a mutable structure containing the input info related to the vpa grid vpa = grid_input_mutable("vpa", ngrid_vpa, nelement_vpa, nelement_vpa, L_vpa, discretization_option_vpa, finite_difference_option_vpa, boundary_option_vpa, - advection_vpa) + advection_vpa, element_spacing_option_vpa) ############################################################################ ################### parameters related to the vperp grid ##################### # ngrid_vperp is the number of grid points per element @@ -775,10 +785,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_vperp = advection_input_mutable(advection_option_vperp, advection_speed_vperp, frequency_vperp, oscillation_amplitude_vperp) + element_spacing_option_vperp = "uniform" # create a mutable structure containing the input info related to the vperp grid vperp = grid_input_mutable("vperp", ngrid_vperp, nelement_vperp, nelement_vperp, L_vperp, discretization_option_vperp, finite_difference_option_vperp, boundary_option_vperp, - advection_vperp) + advection_vperp, element_spacing_option_vperp) ############################################################################ ################### parameters related to the gyrophase grid ##################### # ngrid_gyrophase is the number of grid points per element @@ -798,10 +809,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) oscillation_amplitude_gyrophase = 1.0 advection_gyrophase = advection_input_mutable(advection_option_gyrophase, advection_speed_gyrophase, frequency_gyrophase, oscillation_amplitude_gyrophase) + element_spacing_option_gyrophase = "uniform" # create a mutable structure containing the input info related to the gyrophase grid gyrophase = grid_input_mutable("gyrophase", ngrid_gyrophase, nelement_gyrophase, nelement_gyrophase, L_gyrophase, discretization_option_gyrophase, finite_difference_option_gyrophase, boundary_option_gyrophase, - advection_gyrophase) + advection_gyrophase, element_spacing_option_gyrophase) ############################################################################ ################### parameters related to the vr grid ##################### # ngrid_vr is the number of grid points per element @@ -833,10 +845,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_vr = advection_input_mutable(advection_option_vr, advection_speed_vr, frequency_vr, oscillation_amplitude_vr) + element_spacing_option_vr = "uniform" # create a mutable structure containing the input info related to the vr grid vr = grid_input_mutable("vr", ngrid_vr, nelement_vr, nelement_vr, L_vr, discretization_option_vr, finite_difference_option_vr, boundary_option_vr, - advection_vr) + advection_vr, element_spacing_option_vr) ############################################################################ ################### parameters related to the vz grid ##################### # ngrid_vz is the number of grid points per element @@ -868,10 +881,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_vz = advection_input_mutable(advection_option_vz, advection_speed_vz, frequency_vz, oscillation_amplitude_vz) + element_spacing_option_vz = "uniform" # create a mutable structure containing the input info related to the vz grid vz = grid_input_mutable("vz", ngrid_vz, nelement_vz, nelement_vz, L_vz, discretization_option_vz, finite_difference_option_vz, boundary_option_vz, - advection_vz) + advection_vz, element_spacing_option_vz) ############################################################################ ################### parameters related to the vzeta grid ##################### # ngrid_vzeta is the number of grid points per element @@ -903,10 +917,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # mutable struct containing advection speed options/inputs for z advection_vzeta = advection_input_mutable(advection_option_vzeta, advection_speed_vzeta, frequency_vzeta, oscillation_amplitude_vzeta) + element_spacing_option_vzeta = "uniform" # create a mutable structure containing the input info related to the vzeta grid vzeta = grid_input_mutable("vzeta", ngrid_vzeta, nelement_vzeta, nelement_vzeta, L_vzeta, discretization_option_vzeta, finite_difference_option_vzeta, boundary_option_vzeta, - advection_vzeta) + advection_vzeta, element_spacing_option_vzeta) ############################################################################# # define default values and create corresponding mutable structs holding # information about the composition of the species and their initial conditions diff --git a/src/post_processing.jl b/src/post_processing.jl index 69f1021f3..d2d1f7ef3 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -196,10 +196,11 @@ end function construct_global_zr_coords(r_local, z_local) function make_global_input(coord_local) + #element_spacing_option = "" # dummy value return grid_input(coord_local.name, coord_local.ngrid, coord_local.nelement_global, coord_local.nelement_global, 1, 0, coord_local.L, coord_local.discretization, coord_local.fd_option, coord_local.bc, - coord_local.advection, MPI.COMM_NULL) + coord_local.advection, MPI.COMM_NULL, coord_local.element_spacing_option) end r_global, r_global_spectral = define_coordinate(make_global_input(r_local)) @@ -1054,7 +1055,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) composition, species, collisions, geometry, drive_input, num_diss_params, manufactured_solns_input = input - if !is_1D1V + if !is_1D1V || true # make plots and animations of the phi, Ez and Er plot_charged_moments_2D(density, parallel_flow, parallel_pressure, time, z_global.grid, r_global.grid, iz0, ir0, n_ion_species, diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl index c1bbdd483..eab6e882e 100644 --- a/test/calculus_tests.jl +++ b/test/calculus_tests.jl @@ -39,9 +39,11 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - discretization, fd_option, bc, adv_input, comm) + discretization, fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' x, spectral = define_coordinate(input) # create array for the function f(x) to be differentiated/integrated @@ -86,9 +88,11 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - input = grid_input("coord", ngrid, nelement, + element_spacing_option = "uniform" # dummy value + input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "finite_difference", fd_option, bc, adv_input, comm) + "finite_difference", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' x, spectral = define_coordinate(input) @@ -132,10 +136,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "finite_difference", fd_option, bc, adv_input, comm) + "finite_difference", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' x, spectral = define_coordinate(input) @@ -175,10 +181,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "finite_difference", fd_option, bc, adv_input, comm) + "finite_difference", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' x, spectral = define_coordinate(input) @@ -226,10 +234,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "finite_difference", fd_option, bc, adv_input, comm) + "finite_difference", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' x, spectral = define_coordinate(input) @@ -440,10 +450,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "chebyshev_pseudospectral", fd_option, bc, adv_input, comm) + "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. x, spectral = define_coordinate(input) @@ -634,10 +646,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "chebyshev_pseudospectral", fd_option, bc, adv_input, comm) + "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. x, spectral = define_coordinate(input) @@ -676,10 +690,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "chebyshev_pseudospectral", fd_option, bc, adv_input, comm) + "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. x, spectral = define_coordinate(input) @@ -726,10 +742,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "chebyshev_pseudospectral", fd_option, bc, adv_input, comm) + "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. x, spectral = define_coordinate(input) @@ -940,9 +958,11 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - "chebyshev_pseudospectral", fd_option, bc, adv_input, comm) + "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. x, spectral = define_coordinate(input) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index b40a55ff7..33cd0fbaa 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -35,10 +35,12 @@ function runtests() nelement_local = nelement nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, - discretization, fd_option, bc, adv_input, comm) + discretization, fd_option, bc, adv_input, comm, + element_spacing_option) # create the coordinate struct 'z' z, spectral = define_coordinate(input) diff --git a/test/nonlinear_sound_wave_tests.jl b/test/nonlinear_sound_wave_tests.jl index 2742760cd..703c71b21 100644 --- a/test/nonlinear_sound_wave_tests.jl +++ b/test/nonlinear_sound_wave_tests.jl @@ -348,17 +348,18 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) adv_input = advection_input("default", 1.0, 0.0, 0.0) nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], test_input["z_nelement"], nrank_per_block, irank, z_L, test_input["z_discretization"], "", "periodic", #test_input["z_bc"], - adv_input,comm) + adv_input,comm, element_spacing_option) z, z_spectral = define_coordinate(input) input = grid_input("coord", test_input["vpa_ngrid"], test_input["vpa_nelement"], test_input["vpa_nelement"], nrank_per_block, irank, vpa_L, test_input["vpa_discretization"], "", - test_input["vpa_bc"], adv_input, comm) + test_input["vpa_bc"], adv_input, comm, element_spacing_option) vpa, vpa_spectral = define_coordinate(input) # Test against values interpolated onto 'expected' grid which is fairly coarse no we diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl index 9a5708921..9b51063e2 100644 --- a/test/wall_bc_tests.jl +++ b/test/wall_bc_tests.jl @@ -191,11 +191,12 @@ function run_test(test_input, expected_phi, tolerance; args...) adv_input = advection_input("default", 1.0, 0.0, 0.0) nrank_per_block = 0 # dummy value irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], test_input["z_nelement"], nrank_per_block, irank, 1.0, test_input["z_discretization"], "", test_input["z_bc"], - adv_input, comm) + adv_input, comm, element_spacing_option) z, z_spectral = define_coordinate(input) # Cross comparison of all discretizations to same benchmark From 4ff15a702fa134687077d3c8b986ad4d9c1d6527 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 11 Sep 2023 19:03:01 +0100 Subject: [PATCH 04/23] Flush stdout after printing progress messages When running as a batch job, status messages often do not get flushed immediately by default, so force this. Also add a message for when setup_moment_kinetics() starts, so if the code hangs we can see if it was during compilation or during the actual initialization. --- src/moment_kinetics.jl | 5 +++++ src/time_advance.jl | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index c6eb39282..d1ddaeaf5 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -347,6 +347,11 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, # Set up MPI initialize_comms!() + if global_rank[] == 0 + println("Starting setup ", Dates.format(now(), dateformat"H:MM:SS")) + flush(stdout) + end + input = mk_input(input_dict; save_inputs_to_txt=true, ignore_MPI=false) # obtain input options from moment_kinetics_input.jl # and check input to catch errors diff --git a/src/time_advance.jl b/src/time_advance.jl index 5b6d473b1..d8f499333 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -736,6 +736,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro @serial_region begin if global_rank[] == 0 println("beginning time advance ", Dates.format(now(), dateformat"H:MM:SS")) + flush(stdout) end end @@ -771,6 +772,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro if isfile(t_input.stopfile) # Stop cleanly if a file called 'stop' was created println("Found 'stop' file $(t_input.stopfile), aborting run") + flush(stdout) finish_now = true end end @@ -829,12 +831,14 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro end if global_rank[] == 0 println(" residuals:", result_string) + flush(stdout) end if t_input.converged_residual_value > 0.0 if global_rank[] == 0 if all(r < t_input.converged_residual_value for r ∈ all_residuals) println("Run converged! All tested residuals less than ", t_input.converged_residual_value) + flush(stdout) finish_now = true end end @@ -843,6 +847,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro else if global_rank[] == 0 println() + flush(stdout) end end @@ -972,6 +977,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro if global_rank[] == 0 println("writing distribution functions at step ", i," ", Dates.format(now(), dateformat"H:MM:SS")) + flush(stdout) end end write_dfns_data_to_binary(pdf.charged.norm, pdf.neutral.norm, moments, fields, From 50449bc99e9b28d87ecf581c8b9a944dc3db5594 Mon Sep 17 00:00:00 2001 From: mrhardman Date: Wed, 13 Sep 2023 09:03:53 +0100 Subject: [PATCH 05/23] Update src/chebyshev.jl Delete commented out code Co-authored-by: John Omotani --- src/chebyshev.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 854cdf7a5..b62eb5601 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -63,8 +63,6 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n, @inbounds for j ∈ 1:nelement_local #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor # amount by which to shift the centre of this element from zero - #iel_global = j + irank*nelement_local - #shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5) scale_factor = element_scale[j] shift = element_shift[j] # reverse the order of the original chebyshev_grid (ran from [1,-1]) From d76373c0f88f61a2cff20466caedb6f2bcec8ac5 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 13 Sep 2023 13:07:24 +0100 Subject: [PATCH 06/23] Address @johnomotani comments regarding commented out lines and documentation. --- src/chebyshev.jl | 32 ++++++++++++++++++++------------ src/coordinates.jl | 14 ++------------ src/file_io.jl | 2 +- src/post_processing.jl | 3 +-- 4 files changed, 24 insertions(+), 27 deletions(-) diff --git a/src/chebyshev.jl b/src/chebyshev.jl index b62eb5601..2a235d39d 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -41,6 +41,22 @@ end """ initialize chebyshev grid scaled to interval [-box_length/2, box_length/2] +we no longer pass the box_length to this function, but instead pass precomputed +arrays element_scale and element_shift that are needed to compute the grid. + +ngrid -- number of points per element (including boundary points) +nelement_local -- number of elements in the local (distributed memory MPI) grid +n -- total number of points in the local grid (excluding duplicate points) +element_scale -- the scale factor in the transform from the coordinates + where the element limits are -1, 1 to the coordinate where + the limits are Aj = coord.grid[imin[j]-1] and Bj = coord.grid[imax[j]] + element_scale = 0.5*(Bj - Aj) +element_shift -- the centre of the element in the extended grid coordinate + element_shift = 0.5*(Aj + Bj) +imin -- the array of minimum indices of each element on the extended grid. + By convention, the duplicated points are not included, so for element index j > 1 + the lower boundary point is actually imin[j] - 1 +imax -- the array of maximum indices of each element on the extended grid. """ function scaled_chebyshev_grid(ngrid, nelement_local, n, element_scale, element_shift, imin, imax) @@ -55,14 +71,12 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n, # setup the scale factor by which the Chebyshev grid on [-1,1] # is to be multiplied to account for the full domain [-L/2,L/2] # and the splitting into nelement elements with ngrid grid points - #scale_factor = 0.5*box_length/float(nelement_global) + # account for the fact that the minimum index needed for the chebyshev_grid # within each element changes from 1 to 2 in going from the first element # to the remaining elements k = 1 @inbounds for j ∈ 1:nelement_local - #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor - # amount by which to shift the centre of this element from zero scale_factor = element_scale[j] shift = element_shift[j] # reverse the order of the original chebyshev_grid (ran from [1,-1]) @@ -88,11 +102,9 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info) # check array bounds @boundscheck nelement == size(chebyshev.f,2) || throw(BoundsError(chebyshev.f)) @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df)) - # note that one must multiply by 2*nelement/L to get derivative - # in scaled coordinate - #scale_factor = 2.0*float(coord.nelement_global)/coord.L - # scale factor is (length of a single element/2)^{-1} - + # note that one must multiply by a coordinate transform factor 1/element_scale[j] + # for each element j to get derivative on the extended grid + # variable k will be used to avoid double counting of overlapping point # at element boundaries (see below for further explanation) k = 0 @@ -335,13 +347,9 @@ function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_s # scale to account for modified domain (not [-1,1]) wgts[1:ngrid] = w*element_scale[1] if nelement_local > 1 - # account for double-counting of points at inner element boundaries - #wgts[ngrid] *= 2 for j ∈ 2:nelement_local wgts[imin[j]-1:imax[j]] .+= w*element_scale[j] end - # remove double-counting of outer element boundary for last element - #wgts[n] *= 0.5 end end return wgts diff --git a/src/coordinates.jl b/src/coordinates.jl index b44a7b2ec..6ad0164b7 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -4,9 +4,7 @@ module coordinates export define_coordinate, write_coordinate export equally_spaced_grid -# testing export set_element_boundaries -export init_grid using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float, allocate_int @@ -91,8 +89,6 @@ struct coordinate element_scale::Array{mk_float,1} # shift for each element element_shift::Array{mk_float,1} - # boundaries for each element - #element_boundaries::Array{mk_float,1} # option used to set up element spacing element_spacing_option::String end @@ -163,7 +159,7 @@ function define_coordinate(input, parallel_io::Bool=false) cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option, input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm, - local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)#, element_boundaries) + local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option) if input.discretization == "chebyshev_pseudospectral" && coord.n > 1 # create arrays needed for explicit Chebyshev pseudospectral treatment in this @@ -188,15 +184,9 @@ function set_element_boundaries(nelement_global, L, element_spacing_option) # number of boundaries of sqrt grid nsqrt = floor(mk_int,(nelement_global)/2) + 1 if nelement_global%2 > 0 # odd - #fac = 0.0 - #for j in 2:nsqrt - # fac += 2.0*(((j-2)/(nsqrt-1))^2 - ((j-1)/(nsqrt-1))^2 ) - #end - x = ((nsqrt-2)/(nsqrt-1))^2 if nsqrt < 3 fac = 2.0/3.0 else - #fac = 0.5*(sqrt( 1.0 + 4.0*x) - 1.0)/x fac = 1.0/( 3.0/2.0 - 0.5*((nsqrt-2)/(nsqrt-1))^2) end else @@ -210,7 +200,7 @@ function set_element_boundaries(nelement_global, L, element_spacing_option) element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2 end - elseif element_spacing_option == "uniform" || nelement_global < 4 # uniform spacing + elseif element_spacing_option == "uniform" || (element_spacing_option == "sqrt" && nelement_global < 4) # uniform spacing for j in 1:nelement_global+1 element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5) end diff --git a/src/file_io.jl b/src/file_io.jl index 5dd229843..d77ae35dd 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -521,7 +521,7 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_ write_single_value!(group, "bc", coord.bc; parallel_io=parallel_io, description="boundary condition for $coord_name") - # write the boundary condition for the coordinate + # write the element spacing option for the coordinate write_single_value!(group, "element_spacing_option", coord.element_spacing_option; parallel_io=parallel_io, description="element_spacing_option for $coord_name") diff --git a/src/post_processing.jl b/src/post_processing.jl index d2d1f7ef3..5b3857526 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -196,7 +196,6 @@ end function construct_global_zr_coords(r_local, z_local) function make_global_input(coord_local) - #element_spacing_option = "" # dummy value return grid_input(coord_local.name, coord_local.ngrid, coord_local.nelement_global, coord_local.nelement_global, 1, 0, coord_local.L, coord_local.discretization, coord_local.fd_option, coord_local.bc, @@ -1055,7 +1054,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) composition, species, collisions, geometry, drive_input, num_diss_params, manufactured_solns_input = input - if !is_1D1V || true + if !is_1D1V # make plots and animations of the phi, Ez and Er plot_charged_moments_2D(density, parallel_flow, parallel_pressure, time, z_global.grid, r_global.grid, iz0, ir0, n_ion_species, From d81cf24dfd872a99ea3483e2917f84d333667f00 Mon Sep 17 00:00:00 2001 From: mrhardman Date: Wed, 13 Sep 2023 16:57:09 +0100 Subject: [PATCH 07/23] Update src/chebyshev.jl Remove out of date comment Co-authored-by: John Omotani --- src/chebyshev.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 2a235d39d..30209ec7b 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -68,9 +68,6 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n, chebyshev_grid = chebyshevpoints(ngrid) # create array for the full grid grid = allocate_float(n) - # setup the scale factor by which the Chebyshev grid on [-1,1] - # is to be multiplied to account for the full domain [-L/2,L/2] - # and the splitting into nelement elements with ngrid grid points # account for the fact that the minimum index needed for the chebyshev_grid # within each element changes from 1 to 2 in going from the first element From 6b35ca3315a9dbf12c509893f9445d04cbae5c66 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 10 Sep 2023 10:51:50 +0100 Subject: [PATCH 08/23] If a NaN is found, stop the simulation When writing output, check for NaNs and Allreduce/Bcast the result so that the simulation can stop cleanly. --- src/time_advance.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/time_advance.jl b/src/time_advance.jl index 3a7ea1dc6..350d80f25 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -8,8 +8,8 @@ export time_advance! using MPI using ..type_definitions: mk_float using ..array_allocation: allocate_float, allocate_shared_float -using ..communication: _block_synchronize, block_rank, global_size, iblock_index, - comm_world, MPISharedArray, global_rank +using ..communication +using ..communication: _block_synchronize using ..debugging using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_data_to_binary, debug_dump using ..looping @@ -785,6 +785,25 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro flush(stdout) finish_now = true end + + # If NaNs are present, they should propagate into every field, so only need to + # check one. Choose phi because it is small (it has no species or velocity + # dimensions). If a NaN is found, stop the simulation. + if block_rank[] == 0 + if any(isnan.(fields.phi)) + println("Found NaN, stopping simulation") + found_nan = 1 + else + found_nan = 0 + end + found_nan = MPI.Allreduce(found_nan, +, comm_inter_block[]) + else + found_nan = 0 + end + found_nan = MPI.Bcast(found_nan, 0, comm_block[]) + if found_nan != 0 + finish_now = true + end end # write moments data to file every nwrite_moments time steps if mod(i,t_input.nwrite_moments) == 0 || finish_now From 081174913bcd25ee07c0338d29cf9d1ba2afae18 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 19 Sep 2023 18:01:17 +0100 Subject: [PATCH 09/23] Simplify restart functionality Move actual implementation of restart functionality inside setup_moment_kinetics!() - this also fixes a bug where `iblock_index[]` was not initialised when `get_backup_filename()` needed to use it. This allows the `restart_moment_kinetics()` function to be removed. It is replaced by a couple of arguments to run_moment_kinetics. The `restart_moment_kinetics.jl` script is removed. Now either pass the `--restart` flag to `run_moment_kinetics()` to restart from the latest file in the run directory, or pass the output file to restart from as the second positional argument. --- machines/README.md | 2 +- machines/archer/jobscript-restart.template | 2 +- machines/marconi/jobscript-restart.template | 2 +- restart_moment_kinetics.jl | 7 - src/command_line_options.jl | 4 + src/moment_kinetics.jl | 177 +++++++------------- src/post_processing.jl | 12 +- submit-restart.sh | 13 +- 8 files changed, 74 insertions(+), 145 deletions(-) delete mode 100644 restart_moment_kinetics.jl diff --git a/machines/README.md b/machines/README.md index ac977fdf3..01031eda2 100644 --- a/machines/README.md +++ b/machines/README.md @@ -67,7 +67,7 @@ $ ./submit-restart.sh .toml ``` will submit a job to run and post-process a restart using input file. The simulation will restart from the last time point of the previous run -(`restart_moment_kinetics.jl` supports more flexibility, but for now you would +(`run_moment_kinetics.jl` supports more flexibility, but for now you would need to write your own submission script to pass the options needed for that). Default parameters for the runs (number of nodes, time limit, etc.) were set up diff --git a/machines/archer/jobscript-restart.template b/machines/archer/jobscript-restart.template index c4018292c..52501d169 100644 --- a/machines/archer/jobscript-restart.template +++ b/machines/archer/jobscript-restart.template @@ -22,6 +22,6 @@ export SRUN_CPUS_PER_TASK=$SLURM_CPUS_PER_TASK echo "running INPUTFILE $(date)" -srun --distribution=block:block --hint=nomultithread --ntasks=$SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no restart_moment_kinetics.jl INPUTFILE RESTARTFROM +srun --distribution=block:block --hint=nomultithread --ntasks=$SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no run_moment_kinetics.jl --restart INPUTFILE RESTARTFROM echo "finished INPUTFILE $(date)" diff --git a/machines/marconi/jobscript-restart.template b/machines/marconi/jobscript-restart.template index 53ba20059..420687c93 100644 --- a/machines/marconi/jobscript-restart.template +++ b/machines/marconi/jobscript-restart.template @@ -18,6 +18,6 @@ source julia.env echo "running INPUTFILE $(date)" -mpirun -np $SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no restart_moment_kinetics.jl INPUTFILE RESTARTFROM +mpirun -np $SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no run_moment_kinetics.jl --restart INPUTFILE RESTARTFROM echo "finished INPUTFILE $(date)" diff --git a/restart_moment_kinetics.jl b/restart_moment_kinetics.jl deleted file mode 100644 index c82916954..000000000 --- a/restart_moment_kinetics.jl +++ /dev/null @@ -1,7 +0,0 @@ -# provide option of running from command line via 'julia run_moment_kinetics.jl' -using Pkg -Pkg.activate(".") - -using moment_kinetics - -restart_moment_kinetics() diff --git a/src/command_line_options.jl b/src/command_line_options.jl index ddee48f9d..8d3404f00 100644 --- a/src/command_line_options.jl +++ b/src/command_line_options.jl @@ -28,6 +28,10 @@ const s = ArgParseSettings() "--long" help = "Include more tests, increasing test run time." action = :store_true + "--restart" + help = "Restart from latest output file in run directory (ignored if " * + "`restartfile` is passed)" + action = :store_true "--restart-time-index" help = "Time index in output file to restart from, defaults to final time point" arg_type = Int diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index d1ddaeaf5..80535db2f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -2,7 +2,7 @@ """ module moment_kinetics -export run_moment_kinetics, restart_moment_kinetics +export run_moment_kinetics using MPI @@ -89,11 +89,13 @@ using .type_definitions: mk_int """ main function that contains all of the content of the program """ -function run_moment_kinetics(to::TimerOutput, input_dict=Dict()) +function run_moment_kinetics(to::TimerOutput, input_dict=Dict(); restart=false, + restart_time_index=-1) mk_state = nothing try # set up all the structs, etc. needed for a run - mk_state = setup_moment_kinetics(input_dict) + mk_state = setup_moment_kinetics(input_dict; restart=restart, + restart_time_index=restart_time_index) # solve the 1+1D kinetic equation to advance f in time by nstep time steps if run_type == performance_test @@ -136,27 +138,38 @@ end """ overload which takes a filename and loads input """ -function run_moment_kinetics(to::TimerOutput, input_filename::String) - return run_moment_kinetics(to, read_input_file(input_filename)) +function run_moment_kinetics(to::TimerOutput, input_filename::String; restart=false, + restart_time_index=-1) + return run_moment_kinetics(to, read_input_file(input_filename); restart=restart, + restart_time_index=restart_time_index) end """ overload with no TimerOutput arguments """ -function run_moment_kinetics(input) - return run_moment_kinetics(TimerOutput(), input) +function run_moment_kinetics(input; restart=false, restart_time_index=-1) + return run_moment_kinetics(TimerOutput(), input; restart=restart, + restart_time_index=restart_time_index) end """ overload which gets the input file name from command line arguments """ function run_moment_kinetics() - inputfile = get_options()["inputfile"] - if inputfile == nothing - run_moment_kinetics(Dict()) + options = get_options() + inputfile = options["inputfile"] + restart = options["restart"] + if options["restart_file"] !== nothing + restart = options["restartfile"] + end + restart_time_index = options["restart_time_index"] + if inputfile === nothing + this_input = Dict() else - run_moment_kinetics(inputfile) + this_input = inputfile end + run_moment_kinetics(this_input; restart=restart, + restart_time_index=restart_time_index) end """ @@ -227,109 +240,6 @@ function get_backup_filename(filename) backup_moments_filename, backup_prefix_iblock end -""" - restart_moment_kinetics(input_filename::String, - restart_filename::Union{String,Nothing}=nothing, - time_index::Int=-1) - -Restart moment kinetics from an existing run. Space/velocity-space resolution in the -input must be the same as for the original run. - -`input_filename` is the input file to use. - -`restart_filename` can be used to pick a particular distribution-functions-output file to -restart from. By default will use the most recent one (the one without the numerical -suffix) in the run directory. - -`time_index` can be passed to select the time index from `restart_filename` to restart -from. By default the latest time point is used. -""" -function restart_moment_kinetics(input_filename::String, - restart_filename::Union{String,Nothing}=nothing, - time_index::Int=-1) - restart_moment_kinetics(read_input_file(input_filename), restart_filename, - time_index) - return nothing -end -function restart_moment_kinetics() - options = get_options() - inputfile = options["inputfile"] - if inputfile === nothing - error("Must pass input file as first argument to restart a run.") - end - restartfile = options["restartfile"] - if restartfile === nothing - error("Must pass output file to restart from as second argument.") - end - time_index = options["restart-time-index"] - - restart_moment_kinetics(inputfile, restartfile, time_index) - - return nothing -end -function restart_moment_kinetics(input_dict::Dict, - restart_filename::Union{String,Nothing}=nothing, - time_index::Int=-1) - - if restart_filename === nothing - run_name = input_dict["run_name"] - base_directory = get(input_dict, "base_directory", "runs") - output_dir = joinpath(base_directory, run_name) - io_settings = get(input_dict, "output", Dict{String,Any}()) - binary_format = get(io_settings, "binary_format", hdf5) - if binary_format === hdf5 - ext = "h5" - elseif binary_format === netcdf - ext = "cdf" - else - error("Unrecognized binary_format '$binary_format'") - end - restart_filename = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1] - end - - try - # Move the output file being restarted from to make sure it doesn't get - # overwritten. - dfns_filename, backup_dfns_filename, parallel_io, moments_filename, - backup_moments_filename, backup_prefix_iblock = - get_backup_filename(restart_filename) - # Ensure every process got the filenames and checked files exist before moving - # files - MPI.Barrier(comm_world) - if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0) - mv(dfns_filename, backup_dfns_filename) - mv(moments_filename, backup_moments_filename) - end - # Ensure files have been moved before any process tries to read from them - MPI.Barrier(comm_world) - - # Set up all the structs, etc. needed for a run. - mk_state = setup_moment_kinetics(input_dict, - restart_prefix_iblock=backup_prefix_iblock, - restart_time_index=time_index) - - try - time_advance!(mk_state...) - finally - # clean up i/o and communications - # last 2 elements of mk_state are `io` and `cdf` - cleanup_moment_kinetics!(mk_state[end-2:end]...) - end - catch e - # Stop code from hanging when running on multiple processes if only one of them - # throws an error - if global_size[] > 1 - println("Abort called on rank $(block_rank[]) due to error. Error message " - * "was:\n", e) - MPI.Abort(comm_world, 1) - end - - rethrow(e) - end - - return nothing -end - """ Perform all the initialization steps for a run. @@ -339,8 +249,8 @@ reload data from time index given by `restart_time_index` for a restart. `debug_loop_type` and `debug_loop_parallel_dims` are used to force specific set ups for parallel loop ranges, and are only used by the tests in `debug_test/`. """ -function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, - restart_time_index=-1, +function setup_moment_kinetics(input_dict::Dict; + restart=false, restart_time_index=-1, debug_loop_type::Union{Nothing,NTuple{N,Symbol} where N}=nothing, debug_loop_parallel_dims::Union{Nothing,NTuple{N,Symbol} where N}=nothing) @@ -402,7 +312,7 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments, collisions, num_diss_params) - if restart_prefix_iblock === nothing + if restart === false restarting = false # initialize f(z,vpa) and the lowest three v-space moments (density(z), upar(z) and ppar(z)), # each of which may be evolved separately depending on input choices. @@ -416,10 +326,41 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, else restarting = true + if restart === true + run_name = input_dict["run_name"] + base_directory = get(input_dict, "base_directory", "runs") + output_dir = joinpath(base_directory, run_name) + io_settings = get(input_dict, "output", Dict{String,Any}()) + binary_format = get(io_settings, "binary_format", hdf5) + if binary_format === hdf5 + ext = "h5" + elseif binary_format === netcdf + ext = "cdf" + else + error("Unrecognized binary_format '$binary_format'") + end + restart = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1] + end + + # Move the output file being restarted from to make sure it doesn't get + # overwritten. + dfns_filename, backup_dfns_filename, parallel_io, moments_filename, + backup_moments_filename, backup_prefix_iblock = + get_backup_filename(restart) + # Ensure every process got the filenames and checked files exist before moving + # files + MPI.Barrier(comm_world) + if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0) + mv(dfns_filename, backup_dfns_filename) + mv(moments_filename, backup_moments_filename) + end + # Ensure files have been moved before any process tries to read from them + MPI.Barrier(comm_world) + # Reload pdf and moments from an existing output file code_time, previous_runs_info, restart_time_index = reload_evolving_fields!(pdf, moments, boundary_distributions, - restart_prefix_iblock, restart_time_index, + backup_prefix_iblock, restart_time_index, composition, r, z, vpa, vperp, vzeta, vr, vz) _block_synchronize() end diff --git a/src/post_processing.jl b/src/post_processing.jl index 9cee0bf5a..1c38015bc 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -373,12 +373,12 @@ are passed, the plots/movies are given names beginning with `compare_` and are c in the `comparison_plots/` subdirectory. By default plots output from all restarts in a directory. To select a single run, pass the -`run_index` argument - the value corresponds to the `_` suffix given to output files by -`restart_moment_kinetics()`. `run_index` can be an integer (which is applied to all -directories in `prefix...`), or a tuple of integers (which should have the same length as -the number of directories passed to `prefix...`). Use `run_index=-1` to get the most -recent run (which does not have a `_` suffix). Note that `run_index` is only used when -a directory (rather than the prefix of a specific output file) is passed to `prefix...` +`run_index` argument - the value corresponds to the `_` suffix given to output files +when restarting. `run_index` can be an integer (which is applied to all directories in +`prefix...`), or a tuple of integers (which should have the same length as the number of +directories passed to `prefix...`). Use `run_index=-1` to get the most recent run (which +does not have a `_` suffix). Note that `run_index` is only used when a directory +(rather than the prefix of a specific output file) is passed to `prefix...` """ function analyze_and_plot_data(prefix...; run_index=nothing) if length(prefix) == 0 diff --git a/submit-restart.sh b/submit-restart.sh index 759c247c4..3186173a4 100755 --- a/submit-restart.sh +++ b/submit-restart.sh @@ -94,19 +94,10 @@ RUNNAME=$(util/get-run-name.jl $INPUTFILE) RUNDIR=runs/$RUNNAME/ mkdir -p $RUNDIR -# Get default file to restart from, which is the latest run in $RUNDIR -if [[ -z $RESTARTFROM ]]; then - # "shopt -s extglob" is needed to let us use the ?() syntax within a script - # (it doesn't seem to be needed in an interactive shell!). See - # https://www.linuxjournal.com/content/pattern-matching-bash - shopt -s extglob - RESTARTFROM=$(ls $RUNDIR/$RUNNAME.dfns*.?(h5|cdf) | head -n 1) -fi - if [[ $POSTPROC -eq 0 ]]; then - echo "Submitting $INPUTFILE for restart from $RESTARTFROM and post-processing..." + echo "Submitting $INPUTFILE for restart from '$RESTARTFROM' and post-processing..." else - echo "Submitting $INPUTFILE for restart from $RESTARTFROM..." + echo "Submitting $INPUTFILE for restart from '$RESTARTFROM'..." fi # Create a submission script for the run From b42c7c7ebe03fdc2f25d8f4ab47635e8e8d8a784 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 20 Sep 2023 10:58:16 +0100 Subject: [PATCH 10/23] Fix typos in option names --- src/moment_kinetics.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 80535db2f..56594757f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -159,10 +159,10 @@ function run_moment_kinetics() options = get_options() inputfile = options["inputfile"] restart = options["restart"] - if options["restart_file"] !== nothing + if options["restartfile"] !== nothing restart = options["restartfile"] end - restart_time_index = options["restart_time_index"] + restart_time_index = options["restart-time-index"] if inputfile === nothing this_input = Dict() else From 828d9bef3ec88c09f3d00f3e72aeb61f704e303b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 20 Sep 2023 11:49:33 +0100 Subject: [PATCH 11/23] Create `restart_filename` variable, don't modify value of `restart` arg --- src/moment_kinetics.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 56594757f..5067b2987 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -250,7 +250,7 @@ reload data from time index given by `restart_time_index` for a restart. parallel loop ranges, and are only used by the tests in `debug_test/`. """ function setup_moment_kinetics(input_dict::Dict; - restart=false, restart_time_index=-1, + restart::Union{Bool,AbstractString}=false, restart_time_index::mk_int=-1, debug_loop_type::Union{Nothing,NTuple{N,Symbol} where N}=nothing, debug_loop_parallel_dims::Union{Nothing,NTuple{N,Symbol} where N}=nothing) @@ -339,14 +339,16 @@ function setup_moment_kinetics(input_dict::Dict; else error("Unrecognized binary_format '$binary_format'") end - restart = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1] + restart_filename = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1] + else + restart_filename = restart end # Move the output file being restarted from to make sure it doesn't get # overwritten. dfns_filename, backup_dfns_filename, parallel_io, moments_filename, backup_moments_filename, backup_prefix_iblock = - get_backup_filename(restart) + get_backup_filename(restart_filename) # Ensure every process got the filenames and checked files exist before moving # files MPI.Barrier(comm_world) From e3689786a4e1531fbfc22f83187186fc2f5cd450 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 4 Sep 2023 09:12:50 +0100 Subject: [PATCH 12/23] If restarting in a different run directory, do not move existing files Sometimes it is useful to restart one run from a separate run (whose output is in a different directory). In this case there is no need to move the output files that are being read from, as they would not be overwritten by the new run. So do not move them, to keep the directory of the run being restarted from in its original state. --- src/moment_kinetics.jl | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 5067b2987..74519f2d2 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -236,8 +236,9 @@ function get_backup_filename(filename) end backup_dfns_filename == "" && error("Failed to find a name for backup file.") backup_prefix_iblock = ("$(basename)_$(counter)", iblock) + original_prefix_iblock = (basename, iblock) return dfns_filename, backup_dfns_filename, parallel_io, moments_filename, - backup_moments_filename, backup_prefix_iblock + backup_moments_filename, backup_prefix_iblock, original_prefix_iblock end """ @@ -326,10 +327,11 @@ function setup_moment_kinetics(input_dict::Dict; else restarting = true + run_name = input_dict["run_name"] + base_directory = get(input_dict, "base_directory", "runs") + output_dir = joinpath(base_directory, run_name) if restart === true run_name = input_dict["run_name"] - base_directory = get(input_dict, "base_directory", "runs") - output_dir = joinpath(base_directory, run_name) io_settings = get(input_dict, "output", Dict{String,Any}()) binary_format = get(io_settings, "binary_format", hdf5) if binary_format === hdf5 @@ -347,15 +349,25 @@ function setup_moment_kinetics(input_dict::Dict; # Move the output file being restarted from to make sure it doesn't get # overwritten. dfns_filename, backup_dfns_filename, parallel_io, moments_filename, - backup_moments_filename, backup_prefix_iblock = + backup_moments_filename, backup_prefix_iblock, original_prefix_iblock = get_backup_filename(restart_filename) + # Ensure every process got the filenames and checked files exist before moving # files MPI.Barrier(comm_world) - if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0) - mv(dfns_filename, backup_dfns_filename) - mv(moments_filename, backup_moments_filename) + + if abspath(output_dir) == abspath(dirname(dfns_filename)) + # Only move the file if it is in our current run directory. Otherwise we are + # restarting from another run, and will not be overwriting the file. + if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0) + mv(dfns_filename, backup_dfns_filename) + mv(moments_filename, backup_moments_filename) + end + else + # Reload from dfns_filename without moving the file + backup_prefix_iblock = original_prefix_iblock end + # Ensure files have been moved before any process tries to read from them MPI.Barrier(comm_world) From e0b70a0aea1ac1eaa6ea08a72815a0f300c1d863 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 4 Sep 2023 09:20:16 +0100 Subject: [PATCH 13/23] More informative error message when distibuted MPI setup fails --- src/communication.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/communication.jl b/src/communication.jl index 8e9033e19..3259a4a06 100644 --- a/src/communication.jl +++ b/src/communication.jl @@ -125,8 +125,12 @@ function setup_distributed_memory_MPI(z_nelement_global,z_nelement_local,r_nelem end # throw an error if user specified information is inconsistent if (nrank_per_zr_block*nblocks < nrank_global) - error("ERROR: You must choose global number of processes to be an integer multiple of the number of \n - nblocks = (r_nelement_global/r_nelement_local)*(z_nelement_global/z_nelement_local)") + error("ERROR: You must choose global number of processes to be an integer " + * "multiple of the number of\n" + * "nblocks($nblocks) = (r_nelement_global($r_nelement_global)/" + * "r_nelement_local($r_nelement_local))*" + * "(z_nelement_global($z_nelement_global)/" + * "z_nelement_local($z_nelement_local))") end # assign information regarding shared-memory blocks From 4c59b94f7e5986b26a643e7bef34696582e0432f Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 4 Sep 2023 21:06:14 +0100 Subject: [PATCH 14/23] More helpful error message when file to restart from not found --- src/moment_kinetics.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 74519f2d2..880284a0f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -341,7 +341,13 @@ function setup_moment_kinetics(input_dict::Dict; else error("Unrecognized binary_format '$binary_format'") end - restart_filename = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1] + restart_filename_pattern = joinpath(output_dir, run_name * ".dfns*." * ext) + restart_filename_glob = glob(restart_filename_pattern) + if length(restart_filename_glob) == 0 + error("No output file to restart from found matching the pattern " + * "$restart_filename_pattern") + end + restart_filename = restart_filename_glob[1] else restart_filename = restart end From e251f8843cea86ee4743c84459a838a1ac2ed80e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 20 Sep 2023 12:08:44 +0100 Subject: [PATCH 15/23] Update README.md to remove reference to running without input.toml file Using an input file is preferred now, as changing options by editing defaults in moment_kinetics_input.toml forces Julia to redo precompilation of the project. --- README.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/README.md b/README.md index 731ce3061..b4fc35b77 100644 --- a/README.md +++ b/README.md @@ -32,13 +32,10 @@ The full documentation is online at [https://mabarnes.github.io/moment_kinetics] ``` this significantly decreases the load time but prevents code changes from taking effect when `moment_kinetics.so` is used without repeating the precompilation (to use this option, add an option `-Jmoment_kinetics.so` when starting julia). 4) To run julia with optimization, type - ``` - $ julia -O3 --project run_moment_kinetics.jl - ``` - Default input options are specified in `moment_kinetics_input.jl`. The defaults can be modified for a particular run by setting options in a TOML file, for example `input.toml`, which can be passed as an argument ``` $ julia -O3 --project run_moment_kinetics.jl input.toml ``` + Options are specified in a TOML file, e.g. `input.toml` here. The defaults are specified in `moment_kinetics_input.jl`. * To run in parallel, just put `mpirun -np ` in front of the call you would normally use, with `` the number of processes to use. * It may be more convenient when running `moment_kinetics` more than once to work from the Julia REPL, e.g. ``` From 7ba85583f3fc11b8c40ffe8f8093e1fda6a9c48b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 20 Sep 2023 12:19:12 +0100 Subject: [PATCH 16/23] Documentation for restart functionality --- README.md | 32 +++++++++++++++++++++++++++----- 1 file changed, 27 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index b4fc35b77..ec27bc97a 100644 --- a/README.md +++ b/README.md @@ -47,27 +47,49 @@ The full documentation is online at [https://mabarnes.github.io/moment_kinetics] ``` julia> run_moment_kinetics("input.toml") ``` -5) To make plots and calculate frequencies/growth rates, run +5) To restart a simulation using `input.toml` from the last time point in the existing run directory, + ``` + $ julia -O3 --project run_moment_kinetics --restart input.toml + ``` + or to restart from a specific output file - either from the same run or (if the settings are compatible) a different one - here `runs/example/example.dfns.h5` + ``` + $ julia -O3 --project run_moment_kinetics input.toml runs/example/example.dfns.h5 + ``` + The output file must include distribution functions. When not using parallel I/O there will be multiple output files from different MPI ranks - any one of these can be passed. + * To do the same from the Julia REPL + ``` + $ julia -O3 --project + julia> run_moment_kinetics("input.toml", restart=true) + ``` + or + ``` + julia> run_moment_kinetics("input.toml", restart="runs/example/example.dfns.h5") + ``` + * When calling the `run_moment_kinetics()` function you can also choose a particular time index to restart from, e.g. + ``` + julia> run_moment_kinetics("input.toml", restart="runs/example/example.dfns.h5", restart_time_index=42) + ``` +6) To make plots and calculate frequencies/growth rates, run ``` $ julia --project run_post_processing.jl runs/ ``` passing the directory to process as a command line argument. Input options for post-processing can be specified in post_processing_input.jl. -6) Parameter scans (see [Running parameter scans](#running-parameter-scans)) or performance tests can be performed by running +7) Parameter scans (see [Running parameter scans](#running-parameter-scans)) or performance tests can be performed by running ``` $ julia -O3 --project driver.jl ``` If running a scan, it can be parallelised by passing the number of processors as an argument. Scan options are set in `scan_inputs.jl`. -7) Post processing can be done for several directories at once using +8) Post processing can be done for several directories at once using ``` $ julia --project post_processing_driver.jl runs/ runs/ ... ``` passing the directories to process as command line arguments. Optionally pass a number as the first argument to parallelise post processing of different directories. Input options for post-processing can be specified in `post_processing_input.jl`. -8) In the course of development, it is sometimes helpful to upgrade the Julia veriosn. Upgrading the version of Julia or upgrading packages may require a fresh installation of `moment_kinetics`. To make a fresh install with the latest package versions it is necessary to remove (or rename) the `Manifest.jl` file in the main directory, and generate a new `Manifest.jl` with step 1) above. It can sometimes be necessary to remove or rename the `.julia/` folder in your root directory for this step to be successful. +9) In the course of development, it is sometimes helpful to upgrade the Julia veriosn. Upgrading the version of Julia or upgrading packages may require a fresh installation of `moment_kinetics`. To make a fresh install with the latest package versions it is necessary to remove (or rename) the `Manifest.jl` file in the main directory, and generate a new `Manifest.jl` with step 1) above. It can sometimes be necessary to remove or rename the `.julia/` folder in your root directory for this step to be successful. -9) One may have to set an environment variable to avoid error messages from the Qt library. If you execute the command +10) One may have to set an environment variable to avoid error messages from the Qt library. If you execute the command ``` $ julia --project run_post_processing.jl runs/your_run_dir/ From 00498013fb1f2f165f30cd2fff5b00bce99cdfe0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 20 Sep 2023 14:27:09 +0100 Subject: [PATCH 17/23] Mention HDF5 in description of "restartfile" command line option Also move "--restart" and "--restart-time-index" so they are not under the `# Options for tests` heading in the source code. --- src/command_line_options.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/command_line_options.jl b/src/command_line_options.jl index 8d3404f00..430fbb8da 100644 --- a/src/command_line_options.jl +++ b/src/command_line_options.jl @@ -16,7 +16,7 @@ const s = ArgParseSettings() arg_type = String default = nothing "restartfile" - help = "Name of NetCDF file to restart from" + help = "Name of output file (HDF5 or NetCDF) to restart from" arg_type = String default = nothing "--debug", "-d" @@ -24,10 +24,6 @@ const s = ArgParseSettings() "integer values activate more checks (and increase run time)" arg_type = Int default = 0 - # Options for tests - "--long" - help = "Include more tests, increasing test run time." - action = :store_true "--restart" help = "Restart from latest output file in run directory (ignored if " * "`restartfile` is passed)" @@ -36,6 +32,10 @@ const s = ArgParseSettings() help = "Time index in output file to restart from, defaults to final time point" arg_type = Int default = -1 + # Options for tests + "--long" + help = "Include more tests, increasing test run time." + action = :store_true "--verbose", "-v" help = "Print verbose output from tests." action = :store_true From e2a2224db44f46fc037fc9c6f0eb0e88ac334c6c Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 25 Sep 2023 09:23:35 +0100 Subject: [PATCH 18/23] Add "sqrt" grid to fundamental theorem of calculus tests with large error tolerance = 1.0e-2. --- test/calculus_tests.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl index eab6e882e..2adf15a07 100644 --- a/test/calculus_tests.jl +++ b/test/calculus_tests.jl @@ -14,7 +14,7 @@ function runtests() println("calculus tests") @testset "fundamental theorem of calculus" begin @testset "$discretization $ngrid $nelement" for - discretization ∈ ("finite_difference", "chebyshev_pseudospectral"), + (discretization, element_spacing_option, etol) ∈ (("finite_difference", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "sqrt", 1.0e-2)), ngrid ∈ (5,6,7,8,9,10), nelement ∈ (1, 2, 3, 4, 5) if discretization == "finite_difference" && (ngrid - 1) * nelement % 2 == 1 @@ -26,7 +26,6 @@ function runtests() continue end - etol = 1.0e-15 # define inputs needed for the test L = 6.0 bc = "periodic" @@ -39,7 +38,7 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - element_spacing_option = "uniform" # dummy value + #element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, discretization, fd_option, bc, adv_input, comm, @@ -676,7 +675,7 @@ function runtests() end @testset "Chebyshev pseudospectral derivatives (4 argument), Neumann" verbose=false begin - @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), + @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"), nelement ∈ (1:5), ngrid ∈ (3:33) # define inputs needed for the test @@ -691,7 +690,7 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - element_spacing_option = "uniform" + #element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, From a61ff3776c361fba5cd70d7a298f0dfd79610163 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 25 Sep 2023 10:15:14 +0100 Subject: [PATCH 19/23] Include sqrt grid in calculus tests as suggested by @johnomotani. --- test/calculus_tests.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl index 2adf15a07..bdcd10453 100644 --- a/test/calculus_tests.jl +++ b/test/calculus_tests.jl @@ -690,8 +690,7 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - #element_spacing_option = "uniform" - input = grid_input("coord", ngrid, nelement, + input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, element_spacing_option) @@ -727,7 +726,7 @@ function runtests() end @testset "Chebyshev pseudospectral derivatives upwinding (5 argument), Neumann" verbose=false begin - @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), + @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"), nelement ∈ (1:5), ngrid ∈ (3:33) # define inputs needed for the test @@ -742,8 +741,7 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - element_spacing_option = "uniform" - input = grid_input("coord", ngrid, nelement, + input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, "chebyshev_pseudospectral", fd_option, bc, adv_input, comm, element_spacing_option) From e24ef9c789642c3a0f5c9f435dbed99440f30f72 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 25 Sep 2023 10:15:40 +0100 Subject: [PATCH 20/23] Include sqrt grid in interpolation tests as suggested by @johnomotani. --- test/interpolation_tests.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 33cd0fbaa..8885224cc 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -27,8 +27,9 @@ adv_input = advection_input("default", 1.0, 0.0, 0.0) function runtests() @testset "interpolation" verbose=use_verbose begin @testset "$discretization, $ntest, $nelement, $zlim" for - (discretization, rtol) ∈ - (("finite_difference", 1.e-5), ("chebyshev_pseudospectral", 1.e-8)), + (discretization, element_spacing_option, rtol) ∈ + (("finite_difference", "uniform", 1.e-5), ("chebyshev_pseudospectral", "uniform", 1.e-8), + ("chebyshev_pseudospectral", "sqrt", 1.e-8)), ntest ∈ (3, 14), nelement ∈ (2, 8), zlim ∈ (L/2.0, L/5.0) # create the 'input' struct containing input info needed to create a coordinate @@ -36,7 +37,7 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - element_spacing_option = "uniform" + #element_spacing_option = "uniform" input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, discretization, fd_option, bc, adv_input, comm, From 05b7bfecc2acac68871fe62e8f7dd0038b23aae4 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 25 Sep 2023 10:17:02 +0100 Subject: [PATCH 21/23] Include sqrt grid in wall boundary condition tests as suggested by @johnomotani. Note that the minimum number of elements for a nontrivial test are 5 and 6 -- the "sqrt" element spacing option has different behaviour for odd and even number of elements, so both should be tested. --- test/wall_bc_tests.jl | 60 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl index 9b51063e2..2bcda81b0 100644 --- a/test/wall_bc_tests.jl +++ b/test/wall_bc_tests.jl @@ -82,6 +82,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1, "z_nelement" => 1, "z_bc" => "wall", "z_discretization" => "finite_difference", + "z_element_spacing_option" => "uniform", "vpa_ngrid" => 400, "vpa_nelement" => 1, "vpa_L" => 8.0, @@ -98,6 +99,32 @@ test_input_chebyshev = merge(test_input_finite_difference, "z_discretization" => "chebyshev_pseudospectral", "z_ngrid" => 9, "z_nelement" => 2, + "z_element_spacing_option" => "uniform", + "vpa_discretization" => "chebyshev_pseudospectral", + "vpa_ngrid" => 17, + "vpa_nelement" => 10, + "vz_discretization" => "chebyshev_pseudospectral", + "vz_ngrid" => 17, + "vz_nelement" => 10)) + +test_input_chebyshev_sqrt_grid_odd = merge(test_input_finite_difference, + Dict("run_name" => "chebyshev_pseudospectral", + "z_discretization" => "chebyshev_pseudospectral", + "z_ngrid" => 9, + "z_nelement" => 5, # minimum nontrival nelement (odd) + "z_element_spacing_option" => "sqrt", + "vpa_discretization" => "chebyshev_pseudospectral", + "vpa_ngrid" => 17, + "vpa_nelement" => 10, + "vz_discretization" => "chebyshev_pseudospectral", + "vz_ngrid" => 17, + "vz_nelement" => 10)) +test_input_chebyshev_sqrt_grid_even = merge(test_input_finite_difference, + Dict("run_name" => "chebyshev_pseudospectral", + "z_discretization" => "chebyshev_pseudospectral", + "z_ngrid" => 9, + "z_nelement" => 6, # minimum nontrival nelement (even) + "z_element_spacing_option" => "sqrt", "vpa_discretization" => "chebyshev_pseudospectral", "vpa_ngrid" => 17, "vpa_nelement" => 10, @@ -125,7 +152,7 @@ function run_test(test_input, expected_phi, tolerance; args...) # update the default inputs # Convert keyword arguments to a unique name - name = test_input["run_name"] + name = test_input["run_name"] * ", with element spacing: " * test_input["z_element_spacing_option"] if length(args) > 0 name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...) @@ -196,12 +223,15 @@ function run_test(test_input, expected_phi, tolerance; args...) input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], test_input["z_nelement"], nrank_per_block, irank, 1.0, test_input["z_discretization"], "", test_input["z_bc"], - adv_input, comm, element_spacing_option) + adv_input, comm, test_input["z_element_spacing_option"]) z, z_spectral = define_coordinate(input) # Cross comparison of all discretizations to same benchmark - phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral) - @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15) + if test_input["z_element_spacing_option"] == "uniform" + # only support this test for uniform element spacing + phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral) + @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15) + end end end @@ -214,12 +244,32 @@ function runtests() run_test(test_input_finite_difference, nothing, 2.e-3) end - @testset "Chebyshev" begin + @testset "Chebyshev uniform" begin run_test(test_input_chebyshev, [-1.1689445031600718, -0.7479504438063098, -0.6947559936893813, -0.6917252442591313, -0.7180152498764835, -0.9980114095597415], 2.e-3) end + + @testset "Chebyshev sqrt grid odd" begin + run_test(test_input_chebyshev_sqrt_grid_odd, + [-1.2047298885671576, -0.9431378294506091, -0.8084332392927167, + -0.7812620422650213, -0.7233303514000929, -0.7003878610612269, + -0.69572751349158, -0.6933148921301019, -0.6992503992521327, + -0.7115787972775218, -0.7596015032228407, -0.795776514029509, + -0.876303297135126, -1.1471244425913258], + 2.e-3) + end + @testset "Chebyshev sqrt grid even" begin + run_test(test_input_chebyshev_sqrt_grid_even, + [-1.213617049279473, -1.0054529928344382, -0.871444761913497, + -0.836017699317097, -0.7552110924643832, -0.7264644073096705, + -0.7149147366621806, -0.6950077192395091, -0.6923364889119271, + -0.6950077192395089, -0.7149147366621814, -0.7264644073096692, + -0.7552110924643836, -0.8360176993170979, -0.8714447619134948, + -1.0054529928344376, -1.2136170492794727], + 2.e-3) + end end end From 178178a5f6351968665974db5c4f7f0bdb144c52 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 25 Sep 2023 11:48:20 +0100 Subject: [PATCH 22/23] Comment why sqrt spacing skips cross-compare test [skip ci] --- test/calculus_tests.jl | 1 - test/wall_bc_tests.jl | 4 +++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl index bdcd10453..ae198003c 100644 --- a/test/calculus_tests.jl +++ b/test/calculus_tests.jl @@ -38,7 +38,6 @@ function runtests() nrank_per_block = 0 # dummy value irank = 0 # dummy value comm = MPI.COMM_NULL # dummy value - #element_spacing_option = "uniform" # dummy value input = grid_input("coord", ngrid, nelement, nelement_local, nrank_per_block, irank, L, discretization, fd_option, bc, adv_input, comm, diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl index 2bcda81b0..c0b563393 100644 --- a/test/wall_bc_tests.jl +++ b/test/wall_bc_tests.jl @@ -228,7 +228,9 @@ function run_test(test_input, expected_phi, tolerance; args...) # Cross comparison of all discretizations to same benchmark if test_input["z_element_spacing_option"] == "uniform" - # only support this test for uniform element spacing + # Only support this test for uniform element spacing. + # phi is better resolved by "sqrt" spacing grid, so disagrees with benchmark data from + # simulation with uniform element spacing. phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral) @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15) end From 3018b2b4331ffca190f28d11ec1aedd1c6dca202 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 25 Sep 2023 11:49:58 +0100 Subject: [PATCH 23/23] Increase CI timeout for 'test' workflow Extra tests added mean the CI was timing out frequently. --- .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 a55d73248..de2c2690e 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -10,7 +10,7 @@ jobs: matrix: os: [ubuntu-latest, macOS-latest] fail-fast: false - timeout-minutes: 30 + timeout-minutes: 40 steps: - uses: actions/checkout@v2