Make a tutorial on HPC simulations

README.md

@@ -66,7 +66,7 @@ and [Y.-T. Huang](https://github.com/ytdHuang).
 - **Quantum State and Operator Manipulation:** Easily handle quantum states and operators with a rich set of tools, with the same functionalities as `QuTiP`.
 - **Dynamical Evolution:** Advanced solvers for time evolution of quantum systems, thanks to the powerful [`DifferentialEquations.jl`](https://github.com/SciML/DifferentialEquations.jl) package.
 - **GPU Computing:** Leverage GPU resources for high-performance computing. Simulate quantum dynamics directly on the GPU with the same syntax as the CPU case.
-- **Distributed Computing:** Distribute the computation over multiple nodes (e.g., a cluster). For example, you can run hundreds of quantum trajectories in parallel on a cluster, with, again, the same syntax as the simple case.
+- **Distributed Computing:** Distribute the computation over multiple nodes (e.g., a cluster). For example, you can run hundreds of quantum trajectories in parallel on a cluster, with, again, the same syntax as the simple case. See [this tutorial](https://qutip.org/QuantumToolbox.jl/stable/tutorials/cluster) for more information.
 - **Easy Extension:** Easily extend the package, taking advantage of the `Julia` language features, like multiple dispatch and metaprogramming.
 
 ## Installation

docs/make.jl

@@ -62,6 +62,7 @@ const PAGES = [
         ],
         "Miscellaneous Tutorials" => [
             "tutorials/logo.md",
+            "tutorials/cluster.md",
         ],
     ],
     "Resources" => [

docs/src/index.md

@@ -39,7 +39,7 @@ features:
   - icon: <img width="64" height="64" src="https://img.icons8.com/?size=100&id=1W4Bkj363ov0&format=png&color=000000" />
     title: Distributed Computing
     details: Distribute the computation over multiple nodes (e.g., a cluster). Simulate hundreds of quantum trajectories in parallel on a cluster, with, again, the same syntax as the simple case.
-    link: /users_guide/time_evolution/mcsolve
+    link: /tutorials/cluster
 ---
 ```
 

docs/src/tutorials/cluster.md

@@ -0,0 +1,256 @@
+# [Intensive parallelization on a Cluster](@id doc-tutor:Intensive-parallelization-on-a-Cluster)
+
+## Introduction
+
+In this tutorial, we will demonstrate how to seamlessly perform intensive parallelization on a cluster using the **QuantumToolbox.jl** package. Indeed, thanks to the [**Distributed.jl**](https://docs.julialang.org/en/v1/manual/distributed-computing/) and [**ClusterManagers.jl**](https://github.com/JuliaParallel/ClusterManagers.jl) packages, it is possible to parallelize on a cluster with minimal effort. The following examples are applied to a cluster with the [SLURM](https://slurm.schedmd.com/documentation.html) workload manager, but the same principles can be applied to other workload managers, as the [**ClusterManagers.jl**](https://github.com/JuliaParallel/ClusterManagers.jl) package is very versatile.
+
+We will consider two examples:
+
+1. **Parallelization of a Monte Carlo quantum trajectories**
+2. **Parallelization of a Master Equation by sweeping over parameters**
+
+### Monte Carlo Quantum Trajectories
+
+Let's consider a 2-dimensional transferse field Ising model with 4x3 spins. The Hamiltonian is given by
+
+```math
+\hat{H} = \frac{J_z}{2} \sum_{\langle i,j \rangle} \hat{\sigma}_i^z \hat{\sigma}_j^z + h_x \sum_i \hat{\sigma}_i^x \, ,
+```
+
+where the sums are over nearest neighbors, and the collapse operators are given by 
+
+```math
+\hat{c}_i = \sqrt{\gamma} \hat{\sigma}_i^- \, .
+```
+
+We start by creating a file named `run.batch` with the following content:
+
+```bash
+#!/bin/bash
+#SBATCH --job-name=tutorial
+#SBATCH --output=output.out
+#SBATCH --account=your_account
+#SBATCH --nodes=10
+#SBATCH --ntasks-per-node=1
+#SBATCH --cpus-per-task=72
+#SBATCH --mem=128GB
+#SBATCH --time=0:10:00
+#SBATCH --qos=parallel
+
+# Set PATH to include the directory of your custom Julia installation
+export PATH=/home/username/.juliaup/bin:$PATH
+
+# Now run Julia
+julia --project script.jl
+```
+
+where we have to replace `your_account` with the name of your account. This script will be used to submit the job to the cluster. Here, we are requesting 10 nodes with 72 threads each (720 parallel jobs). The `--time` flag specifies the maximum time that the job can run. To see all the available options, you can check the [SLURM documentation](https://slurm.schedmd.com/documentation.html). We also export the path to the custom Julia installation, which is necessary to run the script (replace `username` with your username). Finally, we run the script `script.jl` with the command `julia --project script.jl`.
+
+The `script.jl` contains the following content:
+
+```julia
+using Distributed
+using ClusterManagers
+
+const SLURM_NUM_TASKS = parse(Int, ENV["SLURM_NTASKS"])
+const SLURM_CPUS_PER_TASK = parse(Int, ENV["SLURM_CPUS_PER_TASK"])
+
+exeflags = ["--project=.", "-t $SLURM_CPUS_PER_TASK"]
+addprocs(SlurmManager(SLURM_NUM_TASKS); exeflags=exeflags, topology=:master_worker)
+
+
+println("################")
+println("Hello! You have $(nworkers()) workers with $(remotecall_fetch(Threads.nthreads, 2)) threads each.")
+
+println("----------------")
+
+
+println("################")
+
+flush(stdout)
+
+@everywhere begin
+    using QuantumToolbox
+    using OrdinaryDiffEq
+
+    BLAS.set_num_threads(1)
+end
+
+# Define lattice
+
+Nx = 4
+Ny = 3
+latt = Lattice(Nx = Nx, Ny = Ny)
+
+# Define Hamiltonian and collapse operators
+Jx = 0.0
+Jy = 0.0
+Jz = 1.0
+hx = 0.2
+hy = 0.0
+hz = 0.0
+γ = 1
+
+Sx = mapreduce(i->SingleSiteOperator(sigmax(), i, latt), +, 1:latt.N)
+Sy = mapreduce(i->SingleSiteOperator(sigmay(), i, latt), +, 1:latt.N)
+Sz = mapreduce(i->SingleSiteOperator(sigmaz(), i, latt), +, 1:latt.N)
+
+H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1)
+e_ops = [Sx, Sy, Sz]
+
+# Time Evolution
+
+ψ0 = fock(2^latt.N, 0, dims = ntuple(i->2, Val(latt.N)))
+
+tlist = range(0, 10, 100)
+
+sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=5000, ensemble_method=EnsembleSplitThreads())
+
+##
+
+println("FINISH!")
+
+rmprocs(workers())
+```
+
+In this script, we first load the necessary packages for distributed computing on the cluster (`Distributed.jl` and `ClusterManagers.jl`). Thanks to the environment variables (previously defined in the SLURM script), we can define the number of tasks and the number of CPUs per task. Then, we initialize the distributed network with the `addprocs(SlurmManager(SLURM_NUM_TASKS); exeflags=exeflags, topology=:master_worker)` command. We then import the packages with the `@everywhere` macro, meaning to load them in all the workers. Moreover, in order to avoid conflicts between the multithreading of the BLAS library and the native Julia multithreading, we set the number of threads of the BLAS library to 1 with the `BLAS.set_num_threads(1)` command. More information about this can be found [here](https://docs.julialang.org/en/v1/manual/performance-tips/#man-multithreading-linear-algebra).
+
+With the
+
+```julia
+println("Hello! You have $(nworkers()) workers with $(remotecall_fetch(Threads.nthreads, 2)) threads each.")
+```
+
+command, we test that the distributed network is correctly initialized. The `remotecall_fetch(Threads.nthreads, 2)` command returns the number of threads of the worker with ID 2.
+
+We then write the main part of the script, where we define the lattice through the [`Lattice`](@ref) function. We set the parameters and define the Hamiltonian and collapse operators with the [`DissipativeIsing`](@ref) function. We also define the expectation operators `e_ops` and the initial state `ψ0`. Finally, we perform the Monte Carlo quantum trajectories with the [`mcsolve`](@ref) function. The `ensemble_method=EnsembleSplitThreads()` argument is used to parallelize the Monte Carlo quantum trajectories, by splitting the ensemble of trajectories among the workers. For a more detailed explanation of the different ensemble methods, you can check the [official documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) of the [**DifferentialEquations.jl**](https://github.com/SciML/DifferentialEquations.jl/) package. Finally, the `rmprocs(workers())` command is used to remove the workers after the computation is finished.
+
+The output of the script will be printed in the `output.out` file, which contains an output similar to the following:
+
+```
+################
+Hello! You have 10 workers with 72 threads each.
+----------------
+################
+
+Progress: [==============================] 100.0% --- Elapsed Time: 0h 00m 21s (ETA: 0h 00m 00s)
+
+FINISH!
+```
+
+where we can see that the computation **lasted only 21 seconds**.
+
+### Master Equation by Sweeping Over Parameters
+
+In this example, we will consider a driven Jaynes-Cummings model, describing a two-level atom interacting with a driven cavity mode. The Hamiltonian is given by
+
+```math
+\hat{H} = \omega_c \hat{a}^\dagger \hat{a} + \frac{\omega_q}{2} \hat{\sigma}_z + g (\hat{a} \hat{\sigma}_+ + \hat{a}^\dagger \hat{\sigma}_-) + F \cos(\omega_d t) (\hat{a} + \hat{a}^\dagger) \, ,
+```
+
+and the collapse operators are given by
+
+```math
+\hat{c}_1 = \sqrt{\gamma} \hat{a} \, , \quad \hat{c}_2 = \sqrt{\gamma} \hat{\sigma}_- \, .
+```
+
+The SLURM file is the same as before, but the `script.jl` file now contains the following content:
+
+```julia
+using Distributed
+using ClusterManagers
+
+const SLURM_NUM_TASKS = parse(Int, ENV["SLURM_NTASKS"])
+const SLURM_CPUS_PER_TASK = parse(Int, ENV["SLURM_CPUS_PER_TASK"])
+
+exeflags = ["--project=.", "-t $SLURM_CPUS_PER_TASK"]
+addprocs(SlurmManager(SLURM_NUM_TASKS); exeflags=exeflags, topology=:master_worker)
+
+
+println("################")
+println("Hello! You have $(nworkers()) workers with $(remotecall_fetch(Threads.nthreads, 2)) threads each.")
+
+println("----------------")
+
+
+println("################")
+
+flush(stdout)
+
+@everywhere begin
+    using QuantumToolbox
+    using OrdinaryDiffEq
+
+    BLAS.set_num_threads(1)
+end
+
+@everywhere begin
+    const Nc = 20
+    const ωc = 1.0
+    const g = 0.05
+    const γ = 0.01
+    const F = 0.01
+
+    const a = tensor(destroy(Nc), qeye(2))
+
+    const σm = tensor(qeye(Nc), sigmam())
+    const σp = tensor(qeye(Nc), sigmap())
+
+    H(ωq) = ωc * a' * a + ωq * tensor(num(Nc), qeye(2)) + g * (a' * σm + a * σp)
+
+    coef(p, t) = p.F * cos(p.ωd * t) # coefficient for the time-dependent term
+
+    const c_ops = [sqrt(γ) * a, sqrt(γ) * σm]
+    const e_ops = [a' * a]
+end
+
+# Define the ODE problem and the EnsembleProblem generation function
+
+@everywhere begin
+    ωq_list = range(ωc - 3*g, ωc + 3*g, 100)
+    ωd_list = range(ωc - 3*g, ωc + 3*g, 100)
+
+    const iter = collect(Iterators.product(ωq_list, ωd_list))
+
+    function my_prob_func(prob, i, repeat, channel)
+        ωq, ωd = iter[i]
+        H_i = H(ωq)
+        H_d_i = H_i + QobjEvo(a + a', coef) # Hamiltonian with a driving term
+
+        L = liouvillian(H_d_i, c_ops).data # Make the Liouvillian
+
+        put!(channel, true) # Update the progress bar channel
+
+        remake(prob, f=L, p=(F = F, ωd = ωd))
+    end
+end
+
+ωq, ωd = iter[1]
+H0 = H(ωq) + QobjEvo(a + a', coef)
+ψ0 = tensor(fock(Nc, 0), basis(2, 1)) # Ground State
+tlist = range(0, 20 / γ, 1000)
+
+prob = mesolveProblem(H0, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=Val(false), params=(F = F, ωd = ωd))
+
+### Just to print the progress bar
+progr = ProgressBar(length(iter))
+progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
+###
+ens_prob = EnsembleProblem(prob.prob, prob_func=(prob, i, repeat) -> my_prob_func(prob, i, repeat, progr_channel))
+
+
+@sync begin
+    @async while take!(progr_channel)
+        next!(progr)
+    end
+
+    @async begin
+        sol = solve(ens_prob, Tsit5(), EnsembleSplitThreads(), trajectories = length(iter))
+        put!(progr_channel, false)
+    end
+end
+```
+
+We are using the [`mesolveProblem`](@ref) function to define the master equation problem. We added some code to manage the progress bar, which is updated through a `RemoteChannel`. The `prob_func` argument of the `EnsembleProblem` function is used to define the function that generates the problem for each iteration. The `iter` variable contains the product of the `ωq_list` and `ωd_list` lists, which are used to sweep over the parameters. The `sol = solve(ens_prob, Tsit5(), EnsembleDistributed(), trajectories=length(iter))` command is used to solve the problem with the distributed ensemble method. The output of the script will be printed in the `output.out` file, which contains an output similar to the previous example.
+
+