Parallel (multi-) XPU iterative 1D, 2D and 3D diffusion solvers (resolving linear, step-function and nonlinear diffusion coefficient). This software is part of the PTsolvers project.
The aim of this project is to provide iterative solvers assessing the scalability, performance, and robustness of the accelerated pseudo-transient method with application to diffusion processes. The solution strategy characterises as semi-iterative, implementing the second-order convergence acceleration as introduced by, e.g., Frankel (1950).
This repository, together with PseudoTransientStokes.jl, relates to the original research article published in the Geoscientific Model Development journal:
@Article{gmd-15-5757-2022,
AUTHOR = {R\"ass, L. and Utkin, I. and Duretz, T. and Omlin, S. and Podladchikov, Y. Y.},
TITLE = {Assessing the robustness and scalability of the accelerated pseudo-transient method},
JOURNAL = {Geoscientific Model Development},
VOLUME = {15},
YEAR = {2022},
NUMBER = {14},
PAGES = {5757--5786},
URL = {https://gmd.copernicus.org/articles/15/5757/2022/},
DOI = {10.5194/gmd-15-5757-2022}
}- The diffusion equation
- Scripts
- Optimal iteration parameters
- Additional infos
- Questions, comments and discussions
In this study we will use the (non-linear) diffusion (reaction) equation in 1D, 2D and 3D, solving the following equation:
qH = -D ∇ H
dH/dt = -∇ ⋅ qHwhere D stands for the diffusion coefficient, defined as
D=1in the linear case;D=1orD=1e-4in the linear step function case;D=H^3in the nonlinear (power-law) case.
We use the following initial condition in 1D, 2D, 3D, respectively:
The scripts folder contains the various Julia routines to solve the diffusion equation in 1D (diff_1D_*.jl), 2D (diff_2D_*.jl) and 3D (diff_3D_*.jl). The 3D scripts are grouped in a separate diff_3D including shell scripts to automatise multi-XPU execution. All Julia routines depend on:
- ParallelStencil.jl to allow for backend-agnostic parallelisation on multi-threaded CPUs and Nvidia GPUs (AMD support being wip)
- Plots.jl for basic visualisation
- MAT.jl for file I/O (disk save for publication-ready figures post-processing using Matlab visualisation scripts)
All 3D routines, with exception of one, rely additionally on:
- ImplicitGlobalGrid.jl to implement global grid domain decomposition, relying on
- MPI.jl as communication library for distributed memory, point-to-point communication.
To get started, clone or download this repository, launch Julia in project mode julia --project and instantiate or resolve the dependencies from within the REPL in package mode julia> ].
The 1D and 2D scripts can be launched either from with the REPL:
julia> include("diff_<script_name>.jl")
or executed from the shell as:
julia --project --check-bounds=no -O3 diff_<script_name>.jlNote that for optimal performance, scripts should be launched from the shell making sure bound-checking to be deactivated.
The 3D scripts can be launched in distributed configuration using, e.g., MPI. This requires either to use the Julia MPI launcher mpiexecjl or to rely on system provided MPI (more here).
💡 Note: refer to the documentation of your supercomputing centre for instructions to run Julia at scale. Instructions for running on the Piz Daint GPU supercomputer at the Swiss National Supercomputing Centre can be found here.
All scripts parse environment variables defining important simulation parameters, defaulting to heuristics, namely
USE_RETURN=false, USE_GPU=false, DO_VIZ=true, DO_SAVE=false, DO_SAVE_VIZ=falseand
NX, NY, NZdefaulting to 512 in 1D and 2D and 64 in 3D. Running, e.g., a script from the shell using the GPU backend can be achieved as following:
USE_GPU=true julia --project --check-bounds=no -O3 diff_<script_name>.jlThe folder dispersion_analysis contains the analytical derivations for the values of iteration parameters. We provide these derivations for 1D stationary and transient diffusion problems. Only the case of D=const is considered.
The main output of the script is the theoretically predicted value for the non-dimensional parameter Re, which is used in the diffusion solvers. The figure showing the dependency of the residual decay rate on Re is also displayed:
For users' convenience, we provide two versions of each script, one version written in Matlab and the other in Python.
To launch the Matlab version, the working installation of Matlab and Matlab Symbolic Math Toolbox is required.
The second version is implemented using the open-source computer algebra library SymPy as a Jupyter/IPython notebook. The Jupyter notebooks can be viewed directly at GitHub (example). However, in order to view the notebook on a local computer or to make changes to the scripts, the recent Python installation is required. Also, several Python packages need to be installed: SymPy, NumPy, Jupyter, and Matplotlib. The easiest way to install these packages along with their dependencies is to use the Anaconda platform.
After installing Anaconda, open the terminal, cd into the dispersion_analysis folder and create a new conda environment with the following command:
> conda create -n ptsolvers sympy numpy jupyter matplotlib
This command will install the required packages. After the installation completes, activate the environment:
> conda activate ptsolvers
The final step is to launch the Jupyter server:
> jupyter notebook
This command starts a server and opens the browser window with file manager.
The repository implements a reference tests suite, using ReferenceTests.jl, to verify the correctness of the outputed results with resepct to a reference solution.
To discuss technical issues, please post on Julia Discourse in the Julia at Scale topic or in the #gpu or #distributed channels on the Julia Slack (to join, visit https://julialang.org/slack/).
To discuss numerical/domain-science issues, please post on Julia Discourse in the Numerics topic or the Modelling & Simulations topic or whichever other topic fits best your issue.