Merge branch 'master' into Krook-collisions

.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

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 <n>` in front of the call you would normally use, with `<n>` 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.
         ```
@@ -50,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/<directory to process>
     ```
     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/<directory1> runs/<directory2> ...
     ```
     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/

machines/README.md

@@ -67,7 +67,7 @@ $ ./submit-restart.sh <path to input file>.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

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)"

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)"

src/chebyshev.jl

@@ -41,9 +41,25 @@ 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)
+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)
     # initialize chebyshev grid defined on [1,-1]
     # with n grid points chosen to facilitate
     # the fast Chebyshev transform (aka the discrete cosine transform)
@@ -52,27 +68,22 @@ function scaled_chebyshev_grid(ngrid, nelement_global, 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
-    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
-        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
         # after first element, increase minimum index for chebyshev_grid to 2
         # 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
 
@@ -88,11 +99,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
@@ -112,7 +121,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 +332,21 @@ 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
             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
         end
     end
     return wgts

src/command_line_options.jl

@@ -16,22 +16,26 @@ 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"
         help = "Set debugging level, default is 0 (no extra debugging). Higher " *
                "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."
+    "--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
         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

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