From 6a0625fc44dfdf65a0732fe5e24c98ac649dad1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20J=2E=20K=C3=BChn?=
 <62713180+mknaranja@users.noreply.github.com>
Date: Fri, 6 Oct 2023 11:18:06 +0200
Subject: [PATCH] 1 Provide automated scripts for benchmarking and visualizing
 results (#2)

- Creation of a script for benchmarking GMGPolar with MUMPs.
- Creation of python utitilies to extract LIKWID outputs and plot basic curves.
- Added some additional documentation.
---
 batch.sh_paper                                |  17 +-
 include/constants.h                           |  19 +-
 performance/run_gmgpolar.sh                   | 137 ++++++++++---
 plot_tools/plot_benchmarks_scaling.py         | 118 +++++++++++
 plot_tools/read_output.py                     | 190 ++++++++++++++++++
 src/main.cpp                                  |   8 +-
 src/polar_multigrid.cpp                       |   2 +-
 .../CartesianR2GyroSonnendruckerCircular.cpp  |   4 +-
 .../CartesianR2GyroSonnendruckerShafranov.cpp |   4 +-
 ...CartesianR2GyroSonnendruckerTriangular.cpp |   4 +-
 .../CartesianR2SonnendruckerCircular.cpp      |   4 +-
 .../CartesianR2SonnendruckerShafranov.cpp     |   4 +-
 .../CartesianR2SonnendruckerTriangular.cpp    |   4 +-
 .../CartesianR6GyroSonnendruckerCircular.cpp  |   8 +-
 .../CartesianR6GyroSonnendruckerShafranov.cpp |   5 +-
 ...CartesianR6GyroSonnendruckerTriangular.cpp |   5 +-
 .../CartesianR6GyroZoniCircular.cpp           |   1 +
 .../CartesianR6GyroZoniShafranov.cpp          |   1 +
 .../CartesianR6GyroZoniShiftedCircular.cpp    |   1 +
 .../CartesianR6GyroZoniShiftedShafranov.cpp   |   1 +
 .../CartesianR6GyroZoniShiftedTriangular.cpp  |   1 +
 .../CartesianR6GyroZoniTriangular.cpp         |   1 +
 src/test_cases/CartesianR6PoissonCircular.cpp |   1 +
 .../CartesianR6PoissonShafranov.cpp           |   1 +
 .../CartesianR6PoissonTriangular.cpp          |   1 +
 .../CartesianR6SonnendruckerCircular.cpp      |   5 +-
 .../CartesianR6SonnendruckerShafranov.cpp     |   5 +-
 .../CartesianR6SonnendruckerTriangular.cpp    |   5 +-
 src/test_cases/CartesianR6ZoniCircular.cpp    |   1 +
 src/test_cases/CartesianR6ZoniShafranov.cpp   |   1 +
 .../CartesianR6ZoniShiftedCircular.cpp        |   1 +
 .../CartesianR6ZoniShiftedShafranov.cpp       |   1 +
 .../CartesianR6ZoniShiftedTriangular.cpp      |   1 +
 src/test_cases/CartesianR6ZoniTriangular.cpp  |   1 +
 .../PolarR6GyroSonnendruckerCircular.cpp      |   5 +-
 .../PolarR6GyroSonnendruckerShafranov.cpp     |   5 +-
 .../PolarR6GyroSonnendruckerTriangular.cpp    |   5 +-
 src/test_cases/PolarR6GyroZoniCircular.cpp    |   1 +
 src/test_cases/PolarR6GyroZoniShafranov.cpp   |   1 +
 .../PolarR6GyroZoniShiftedCircular.cpp        |   1 +
 .../PolarR6GyroZoniShiftedCulham.cpp          |   1 +
 .../PolarR6GyroZoniShiftedShafranov.cpp       |   1 +
 .../PolarR6GyroZoniShiftedTriangular.cpp      |   1 +
 src/test_cases/PolarR6GyroZoniTriangular.cpp  |   1 +
 src/test_cases/PolarR6PoissonCircular.cpp     |   1 +
 src/test_cases/PolarR6PoissonShafranov.cpp    |   1 +
 src/test_cases/PolarR6PoissonTriangular.cpp   |   1 +
 .../PolarR6SonnendruckerCircular.cpp          |   5 +-
 .../PolarR6SonnendruckerShafranov.cpp         |   5 +-
 .../PolarR6SonnendruckerTriangular.cpp        |   5 +-
 src/test_cases/PolarR6ZoniCircular.cpp        |   1 +
 src/test_cases/PolarR6ZoniShafranov.cpp       |   1 +
 src/test_cases/PolarR6ZoniShiftedCircular.cpp |   1 +
 .../PolarR6ZoniShiftedShafranov.cpp           |   1 +
 .../PolarR6ZoniShiftedTriangular.cpp          |   1 +
 src/test_cases/PolarR6ZoniTriangular.cpp      |   1 +
 treat_output.py                               |   2 +-
 57 files changed, 518 insertions(+), 93 deletions(-)
 create mode 100644 plot_tools/plot_benchmarks_scaling.py
 create mode 100644 plot_tools/read_output.py

diff --git a/batch.sh_paper b/batch.sh_paper
index 3fc45d1a..76db471a 100644
--- a/batch.sh_paper
+++ b/batch.sh_paper
@@ -16,15 +16,8 @@
 debug=0
 v1=1
 v2=1
-cycle=1
-compute_rho=0
-level=-1
 maxiter=300
-periodic=1
-theta_aniso=0
-discr=1
 nr_exp=4
-ntheta_exp=4
 res_norm=3
 R0=1e-6
 DirBC_Interior=0
@@ -80,28 +73,30 @@ done
 mkdir -p outputs
 
 echo "prob alpha_coeff beta_coeff fac_ani extrapolation mod_pk"
-# Triangular-Shafranov
-for mod_pk in 2 1
+# 1) Triangular/Czarny 2) Shafranov
+for mod_pk in 2 1 # 2=Triangular/Czarny, 1=Shafranov
 do
 	# Cartesian + beta 0 + ani 0
-	prob=7
+	prob=7 # Solution (23) of Bourne et al.
 		echo $prob $alpha_coeff $beta_coeff $fac_ani $extrapolation $mod_pk
 		for extrapolation in 1
 		do
 			for divideBy2 in 0 1 2 3 4 5 6		#iterate over the different grid sizes
 			do
+			# note that the divideBy2 option here is only used as a dummy for looping. Grids need to be stored beforehand and are loaded here.
 				echo "./${build_dir}/gmgpolar_simulation  -n "$nr_exp" -a "$fac_ani" --mod_pk "$mod_pk" --DirBC_Interior "$DirBC_Interior" --divideBy2 0 -r "$R0"  --smoother "$smoother" --verbose 2 --debug "$debug" --extrapolation "$extrapolation" --optimized 1 --openmp "$openmp" --v1 "$v1" --v2 "$v2" -R "$R" --prob "$prob"  --maxiter "$maxiter" --alpha_coeff "$alpha_coeff" --beta_coeff "$beta_coeff" --res_norm "$res_norm" --f_grid_r radii_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt --f_grid_theta angles_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt --rel_red_conv "$rel_red_conv" 1> outputs/job.out_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt 2> outputs/job.err_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt"
 				./${build_dir}/gmgpolar_simulation  -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 0 -r $R0  --smoother $smoother --verbose 2 --debug $debug --extrapolation $extrapolation --optimized 1 --openmp $openmp --v1 $v1 --v2 $v2 -R $R --prob $prob  --maxiter $maxiter --alpha_coeff $alpha_coeff --beta_coeff $beta_coeff --res_norm $res_norm --f_grid_r "radii_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt" --f_grid_theta "angles_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt" --rel_red_conv $rel_red_conv 1> "outputs/job.out_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt" 2> "outputs/job.err_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt"
 			done
 		done
 
 	# Polar + beta 0-1 + ani 0-1
-	prob=6
+	prob=6 # Solution (22) of Bourne et al.
 		echo $prob $alpha_coeff $beta_coeff $fac_ani $extrapolation $mod_pk
 		for extrapolation in 1
 		do
 			for divideBy2 in 0 1 2 3 4 5 6		#iterate over the different grid sizes
 			do
+			# note that the divideBy2 option here is only used as a dummy for looping. Grids need to be stored beforehand and are loaded here.
 					echo "./${build_dir}/gmgpolar_simulation -n "$nr_exp" -a "$fac_ani" --mod_pk "$mod_pk" --DirBC_Interior "$DirBC_Interior" --divideBy2 0 -r "$R0"  --smoother "$smoother" --verbose 2 --debug "$debug" --extrapolation "$extrapolation" --optimized 1 --openmp "$openmp" --v1 "$v1" --v2 "$v2" -R "$R" --prob "$prob"  --maxiter "$maxiter" --alpha_coeff "$alpha_coeff" --beta_coeff "$beta_coeff" --res_norm "$res_norm" --f_grid_r radii_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt --f_grid_theta angles_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt --rel_red_conv "$rel_red_conv" 1> outputs/job.out_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt 2> outputs/job.err_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt"
 				./${build_dir}/gmgpolar_simulation -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 0 -r $R0  --smoother $smoother --verbose 2 --debug $debug --extrapolation $extrapolation --optimized 1 --openmp $openmp --v1 $v1 --v2 $v2 -R $R --prob $prob  --maxiter $maxiter --alpha_coeff $alpha_coeff --beta_coeff $beta_coeff --res_norm $res_norm --f_grid_r "radii_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt" --f_grid_theta "angles_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt" --rel_red_conv $rel_red_conv 1> "outputs/job.out_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt" 2> "outputs/job.err_"$fac_ani"_"$mod_pk"_"$prob"_"$beta_coeff"_"$extrapolation"_"$divideBy2"_"$rel_red_conv".txt"
 			done
diff --git a/include/constants.h b/include/constants.h
index 25f757cb..3588e4aa 100644
--- a/include/constants.h
+++ b/include/constants.h
@@ -350,12 +350,14 @@ enum stencil
 //const double PI = 3.141592653589793238463;
 const double PI = M_PI;
 
+// See Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
+// used as Param::mod_pk ("modified polar coordinates") in GMGPolar
 enum geometry_type
 {
-    CIRCULAR   = 0,
-    SHAFRANOV  = 1,
-    TRIANGULAR = 2,
-    CULHAM     = 3
+    CIRCULAR   = 0, // simple circular domain
+    SHAFRANOV  = 1, // Fig. 6a
+    TRIANGULAR = 2, // Fig. 6b (also denoted Czarny)
+    CULHAM     = 3 // Fig. 18
 };
 
 enum alpha_val
@@ -366,12 +368,15 @@ enum alpha_val
     POISSON       = 3,
 };
 
+// Defines the manufactured solution to compare the computed error against.
+// see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1
+// or Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 enum problem_type
 {
     FLAT           = 1,
     REFINED_RADIUS = 4,
-    CARTESIAN_R2   = 5,
-    POLAR_R6       = 6,
-    CARTESIAN_R6   = 7,
+    CARTESIAN_R2   = 5, // 
+    POLAR_R6       = 6, // Bourne et al., Eq. (22)
+    CARTESIAN_R6   = 7, // Bourne et al., Eq. (23)
 };
 #endif // CONSTANTS_HXX
diff --git a/performance/run_gmgpolar.sh b/performance/run_gmgpolar.sh
index cf9906a5..f1460d45 100644
--- a/performance/run_gmgpolar.sh
+++ b/performance/run_gmgpolar.sh
@@ -1,58 +1,133 @@
 #!/bin/bash
+#SBATCH --job-name=gmgpolar-setup
+#SBATCH --output=slurm-%A-setup.out
+#SBATCH --error=slurm-%A-setup.err
+#SBATCH -N 1
+#SBATCH -n 1
+#SBATCH -c 1
+#SBATCH -t 5
 
 #fixed variables
+# If Origin is chosen, the node can be set as coarse or fine. Default: Coarse.
 origin_NOT_coarse=0	# origin_NOT_coarse
+# Choose anisotropy in angle-direction. Default: Off.
 theta_aniso=0		# theta_aniso
+# Smoother 3 is our default, 13 is used for some testing, should be not used
+# for production. -> TODO: Rename smoother names
+smoother=3		    # smoother (3,13)
 
-# changing variables ?!
-prob=5			# prob
-R=1.3			# R
-kappa_eps=0		# k
-delta_e=0		# d
-discr=3			# discr
-fac_ani=3		# a
-nr_exp=4		# n
+# default variables
+# If origin is not a particular node of the mesh, Dirichlet boundary conditions
+# can be implemented on the most inner circle
+DirBC_Interior=1	#  DirBC_Interior (0/1)
+# Generalized radius of most inner circle. Defines if origin will be a particular node.
+R0=1e-8		        # r (1e-8/1e-5/1e-2)
+# Generalized radius of maximum outer circle.
+R=1.3			    # R
+# Anisotropy in radial direction.
+fac_ani=3		    # a
+# TODO: which nr_exp and divideby2 do we want to consider?
+nr_exp=4		    # n
 
 #changing variables
-mod_pk=0		# mod_pk (0/1)
-R0=0.1			# r (1e-8/1e-5/1e-2)
-DirBC_Interior=0	#  DirBC_Interior (0/1)
-divideBy2=0		# divideBy2 (3/4/5/6)
-smoother=3		# smoother (3,13)
-extrapolation=0		# E
+mod_pk=1		    # mod_pk=1: Shafranov geometry
+prob=7			    # Prob=7: Simulate solution (23) of Bourne et al.
+# TODO: which alpha and beta to simulate? Alpha aligned with anisotropy?
+alpha_coeff=2
+beta_coeff=1
+
+# set to on
+extrapolation=1		# E
+
+debug=0
+v1=1
+v2=1
+maxiter=300
+res_norm=3
+rel_red_conv=1e-11
+
+nodes=1
+ranks=1     # number of MPI Ranks
+cores=128    # set OpenMP Num Threads to maximum number of cores requested
+
+####################################
+## create grids                   ##
+####################################
+create_grid=0
+if [ $create_grid ]
+then
+cd ..
+mkdir -p angles_files/Rmax"$R"/aniso"$fac_ani"/
+mkdir -p radii_files/Rmax"$R"/aniso"$fac_ani"/
+# Costly function as setup as expensive and sequential. Only run once.
+for divideBy2 in 0 1 2 3 4 5 6 7 8     # create different grid sizes
+do
+	## ATTENTION / REMARK: 
+	## Please note that these calls will abort/segfault as creation of grids and computation in one step
+	## is not yet supported by GMGPolar. We will make this functionality available in a future commit.
+	## Please ignore abort/segfault for the calls in this loop.
+	# mod_pk has no effect on the creation of grids as the set of (r,theta) is
+	# the same for all geometries, only the mapping F(r,theta) -> (x,y) changes.
+	./build/gmgpolar_simulation -n $nr_exp -a $fac_ani --mod_pk 0 --DirBC_Interior $DirBC_Interior --divideBy2 $divideBy2 -r $R0  --smoother $smoother --verbose 2 --debug $debug --extrapolation $extrapolation --optimized 1 $ --v1 $v1 --v2 $v2 -R $R --prob $prob  --maxiter $maxiter --alpha_coeff $alpha_coeff --beta_coeff $beta_coeff --res_norm $res_norm --write_radii_angles 1 --f_grid_r "radii_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt" --f_grid_theta "angles_files/Rmax"$R"/aniso"$fac_ani"/divide"$divideBy2".txt"
+done
+fi
 
 echo "#!/bin/bash" > run_gmgpolar_sbatch.sh
 # create a short name for your job
 echo "#SBATCH --job-name=gmgpolar" >> run_gmgpolar_sbatch.sh
 # stdout file %A=job id
-echo "#SBATCH --output=slurm-%A-p$prob-fa$fac_ani-r$nr_exp-mpk$mod_pk-s$smoother-e$extrapolation.out" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH --output=slurm-%A-p$prob-r$nr_exp-mpk$mod_pk-s$smoother-e$extrapolation--N$nodes-R$ranks-maxC$cores.out" >> run_gmgpolar_sbatch.sh
 # stderr file
-echo "#SBATCH --error=slurm-%A-p$prob-fa$fac_ani-r$nr_exp-mpk$mod_pk-s$smoother-e$extrapolation.err" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH --error=slurm-%A-p$prob-r$nr_exp-mpk$mod_pk-s$smoother-e$extrapolation--N$nodes-R$ranks-maxC$cores.err" >> run_gmgpolar_sbatch.sh
 
-echo "#SBATCH -N 1" >> run_gmgpolar_sbatch.sh
-echo "#SBATCH -n 1" >> run_gmgpolar_sbatch.sh
-echo "#SBATCH -c 14" >> run_gmgpolar_sbatch.sh
-echo "#SBATCH -t 6000" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH -N $nodes" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH -n $ranks" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH -c $cores" >> run_gmgpolar_sbatch.sh
+# fix to one thread per core
+echo "#SBATCH --threads-per-core=1" >> run_gmgpolar_sbatch.sh
+# fix CPU frequency to 1.8 Mhz
+echo "#SBATCH --cpu-freq=1800000" >> run_gmgpolar_sbatch.sh
+echo "#SBATCH -t 600" >> run_gmgpolar_sbatch.sh
 echo "#SBATCH --exclusive" >> run_gmgpolar_sbatch.sh
  
 # remove potentially loaded and conflicting modules
 echo "module purge" >> run_gmgpolar_sbatch.sh
 
-# gcc10
-echo "module load PrgEnv/gcc10-openmpi" >> run_gmgpolar_sbatch.sh
+# CARO
+#echo "module load rev/23.05" >> run_gmgpolar_sbatch.sh
+# spack install mumps@XXX+metis~mpi
+echo "module load likwid/5.2.2" >> run_gmgpolar_sbatch.sh
+# Local machine
+# echo "module load PrgEnv/gcc10-openmpi" >> run_gmgpolar_sbatch.sh
+
+# echo "cd .." >> run_gmgpolar_sbatch.sh
+# echo "make -j16" >> run_gmgpolar_sbatch.sh
+
+# to be defined for use case (3/4/5/6)
+# Attention: divideBy is used as a dummy variable to access folders as grids are read in
+echo "let divideBy2=4" >> run_gmgpolar_sbatch.sh
 
-echo "cd .." >> run_gmgpolar_sbatch.sh
-echo "make -j16" >> run_gmgpolar_sbatch.sh
+####################################
+## solve system                   ##
+####################################
 
-# FLOPS-DP counter
-echo "for m in {0..1}; do" >> run_gmgpolar_sbatch.sh
-echo "likwid-perfctr -C 0-$m -g FLOPS_DP ./build_gnu/main --matrix_free 1 -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 $divideBy2 -r $R0  --smoother $smoother -E $extrapolation" >> run_gmgpolar_sbatch.sh
+# reduce cores as cores count from 0
+max_threads=$((cores))
+echo "let m=1" >> run_gmgpolar_sbatch.sh
+# FLOPS-DP counter from 1 to cores many threads
+echo "while [ \$m -le $max_threads ]; do" >> run_gmgpolar_sbatch.sh
+echo "let mminus1=m-1" >> run_gmgpolar_sbatch.sh
+echo "srun --cpus-per-task=\$m likwid-perfctr -f -C 0-\$mminus1 -g FLOPS_DP ./build/gmgpolar_simulation --openmp \$m --matrix_free 1 -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 0 -r $R0 --smoother $smoother -E $extrapolation --verbose 2 --debug $debug --optimized 1 --v1 $v1 --v2 $v2 -R $R --prob $prob --maxiter $maxiter --alpha_coeff $alpha_coeff --beta_coeff $beta_coeff --res_norm $res_norm --f_grid_r "radii_files/Rmax"$R"/aniso"$fac_ani"/divide"\$divideBy2".txt" --f_grid_theta "angles_files/Rmax"$R"/aniso"$fac_ani"/divide"\$divideBy2".txt" --rel_red_conv $rel_red_conv" >> run_gmgpolar_sbatch.sh
+echo "let m=m*2" >> run_gmgpolar_sbatch.sh
 echo "done;" >> run_gmgpolar_sbatch.sh
 
-# memory (saturation) benchmarks
-echo "for m in {0..1}; do" >> run_gmgpolar_sbatch.sh
-echo "likwid-perfctr -C 0-$m -g CACHES ./build_gnu/main -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 $divideBy2 -r $R0  --smoother $smoother -E $extrapol" >> run_gmgpolar_sbatch.sh
+# # Memory (saturation) benchmarks
+echo "let m=1" >> run_gmgpolar_sbatch.sh
+echo "while [ \$m -le $max_threads ]; do" >> run_gmgpolar_sbatch.sh
+echo "let mminus1=m-1" >> run_gmgpolar_sbatch.sh
+echo "srun --cpus-per-task=\$m likwid-perfctr -f -C 0-\$mminus1 -g CACHE ./build/gmgpolar_simulation --openmp \$m --matrix_free 1 -n $nr_exp -a $fac_ani --mod_pk $mod_pk --DirBC_Interior $DirBC_Interior --divideBy2 0 -r $R0 --smoother $smoother -E $extrapolation --verbose 2 --debug $debug --optimized 1 --v1 $v1 --v2 $v2 -R $R --prob $prob --maxiter $maxiter --alpha_coeff $alpha_coeff --beta_coeff $beta_coeff --res_norm $res_norm --f_grid_r "radii_files/Rmax"$R"/aniso"$fac_ani"/divide"\$divideBy2".txt" --f_grid_theta "angles_files/Rmax"$R"/aniso"$fac_ani"/divide"\$divideBy2".txt" --rel_red_conv $rel_red_conv" >> run_gmgpolar_sbatch.sh
+echo "let m=m*2" >> run_gmgpolar_sbatch.sh
 echo "done;" >> run_gmgpolar_sbatch.sh
 
 #submit the job
-sbatch run_gmgpolar_sbatch.sh
+sbatch run_gmgpolar_sbatch.sh
\ No newline at end of file
diff --git a/plot_tools/plot_benchmarks_scaling.py b/plot_tools/plot_benchmarks_scaling.py
new file mode 100644
index 00000000..cafcf8a9
--- /dev/null
+++ b/plot_tools/plot_benchmarks_scaling.py
@@ -0,0 +1,118 @@
+import pandas as pd
+import numpy as np
+import sys
+# Use backend to not plot on UI
+# import matplotlib
+# matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from os.path import join, exists, dirname
+from os import makedirs
+
+### Plots scaling of FLOPS and Caches (saturation) scaling from 0 to n Cores
+### as read from data frame
+
+colors  =   [
+            [1.00,  0.49,  0.06],
+            [0.17,  0.63,  0.18],
+            [0.83,  0.15,  0.16],
+            [0.13, 0.47, 0.69],
+            [0.58,  0.40,  0.74],
+            [0.53,  0.35,  0.27],
+            [0.92,  0.46,  0.77],
+            [0.50,  0.50,  0.50],
+            [0.66,  0.85,  0.06],                
+            [0.06,  0.85,  0.85],
+            [0.85,  0.15,  0.85],
+            [0.75,  0.75,  0.75]];  
+
+def plot_perf_per_core(path_out, fname, df, benchname, saturation_limit=0, colors=colors):  
+    fontsize = 16
+
+    fig = plt.figure(figsize=(10, 10))
+    ax = fig.add_subplot()
+    plt.plot(df['Cores'], df[benchname[0]])
+
+    if benchname[0] == 'CACHES':
+        plt.plot(df['Cores'], saturation_limit * np.ones(len(df['Cores'])), linestyle='dotted', linewidth=3, color=[0, 0, 0]) 
+        ax.text(1, saturation_limit+3, 'Memory bandwith (AXPY) (' + str(saturation_limit) + ' GBytes/s)', fontsize=14)
+        ax.set_ylim(0, 90)
+
+
+    ax.set_title(benchname[1][0], fontsize=fontsize+6)
+    ax.set_ylabel(benchname[1][0], fontsize=fontsize)
+
+    ax.set_xlabel('Number of cores used', fontsize=fontsize)
+
+
+    path_out = join(path_out, 'figures')
+    if not exists(path_out):
+        makedirs(path_out)
+    plt.savefig(join(path_out, fname + '_' + benchname[0].lower()), bbox_inches='tight')
+    plt.close()   
+
+
+def main():
+
+    problem = 5
+    nr_exp = 4
+    mod_pk = 1
+    smoother = 3
+    extrapolation = 1
+
+    nodes = 1
+    ranks = 1
+    maxCores = 14
+
+    ## saturation_limit is node specific and needs to be adapted.
+    saturation_limit = 80
+
+    fname = 'p' + str(problem) + '-r' + str(nr_exp) + '-mpk' + str(mod_pk) + '-s' + str(
+        smoother) + '-e' + str(extrapolation) + '--N' + str(nodes) + '-R' + str(ranks) + '-maxC' + str(maxCores)   
+    path_to_files_rel = '' # relative to plot script
+    path_to_files = join(dirname(__file__), join(path_to_files_rel))    
+
+    df = pd.read_csv(
+        join(path_to_files, fname + '_benchmarks.csv'),
+        dtype={'Problem': int, 'rExp': int, 'ModPK': int,
+                'Extrapolation': int, 'Nodes': int, 'Ranks': int,
+                'Cores': int, 'its': int})
+
+    # Likwid benchmark columns, more benchmarks are in timings.
+    likwid_benchmarks = {'FLOPS_DP': ['Flop performance in Multi-Threading', 'Flops (GFlops/s)'], 'CACHES': [
+        'Memory bandwidth saturation', 'Memory bandwidth (GBytes/s)']}  # benchmark : [plot title, plot y-label]
+    timing_benchmarks = {'Total_execution_time' : ['Total execution time in Multi-Threading', 'Execution time']}
+    # Problem setting columns
+    setting_cols = ['Problem', 'rExp', 'ModPK', 'Extrapolation', 'Nodes', 'Ranks']
+
+    # check content
+    for bench in likwid_benchmarks.items():
+        bench_rows = np.where(df[bench[0]].isnull()!=True)[0]
+        if len(bench_rows) > 0:
+            df_subframe = df.iloc[bench_rows].copy()
+
+            # Check that the different number of threads/cores where only conducted
+            # on one particular setting of the above columns
+            if np.max(df_subframe.loc[:,setting_cols].nunique()) > 1:
+                sys.exit('Error in input data, more than one setting found.')
+
+            # TODO or not TODO: If the same run was done multiple times with different LIKWID benchmarks
+            # it is not clear which line to take.
+            # This could be extended to take the weighted sum or minimum of all corresponding lines.
+            # However, these timings should be identical or in the same region...
+            # Nonetheless, there could be a nicer table format to store the results ;-)
+            cores_used = df_subframe['Cores'].unique()
+            if len(cores_used) != len(df_subframe['Cores']):
+                sys.exit('Error: Multiple times computed with the same number of threads.')
+
+            plot_perf_per_core(path_to_files, fname, df_subframe, bench, saturation_limit=saturation_limit)
+
+            # Plot particular timings from table for first benchmark. Timings should be similar for FLOPS_DP and CACHES.
+            if bench[0] == list(likwid_benchmarks.keys())[0]:
+                for timing_benchmark in timing_benchmarks.items():
+                    plot_perf_per_core(path_to_files, fname, df_subframe, timing_benchmark)
+
+
+
+if __name__ == '__main__':
+    main()
+    
\ No newline at end of file
diff --git a/plot_tools/read_output.py b/plot_tools/read_output.py
new file mode 100644
index 00000000..6fbb502a
--- /dev/null
+++ b/plot_tools/read_output.py
@@ -0,0 +1,190 @@
+import pandas as pd
+import numpy as np
+import os
+import sys
+
+
+def main():
+
+    jobid = '1256' # provide correct slurm job id
+    problem = 5
+    nr_exp = 4
+    mod_pk = 1
+    smoother = 3
+    extrapolation = 1
+
+    nodes = 1
+    ranks = 1
+    maxCores = 14
+    cores_used = [] # will be filled automatically
+
+    path_to_perf_files_rel = '' # relative from read_output call
+    path_to_perf_files = os.path.join(os.path.dirname(__file__), os.path.join(path_to_perf_files_rel))
+
+    cols_problem = ['Problem', 'rExp', 'ModPK', 'Extrapolation']
+
+    cols_cluster = ['Nodes', 'Ranks', 'Cores']
+
+    cols_convergence = ['its', '2-norm of error', 'inf-norm of error']
+
+    cols_time = []
+    cols_time += ['t_setup', 't_build', 't_facto_Ac',
+                  't_build_P', 't_build_Asc', 't_facto_Asc']
+    cols_time += ['t_total_(fine)', 't_smoothing', 't_residual',
+                  't_restriction', 't_Ac', 't_prolongation', 't_fine_residual',
+                  't_error']
+    cols_time += ['t_applyA']
+    cols_time += ['t_coeff', 't_arr_art_att', 't_sol', 't_detDFinv', 't_trafo']
+    cols_time += ['Total_execution_time']
+
+    cols_benchmark = ['FLOPS_DP', 'CACHES']
+
+    cols = cols_problem + cols_cluster + cols_convergence + cols_time + cols_benchmark
+    cols_dtypes = {col: 'int' for col in cols_problem + cols_cluster}
+    cols_dtypes.update(
+        {col: 'double' for col in cols_convergence + cols_time + cols_benchmark})
+
+    filename_input = 'p' + str(problem) + '-r' + str(nr_exp) + '-mpk' + str(mod_pk) + '-s' + str(
+        smoother) + '-e' + str(extrapolation) + '--N' + str(nodes) + '-R' + str(ranks) + '-maxC' + str(maxCores)
+
+    err = 0
+    with open(os.path.join(path_to_perf_files, 'slurm-' + str(jobid) + '-' + filename_input + '.err')) as f:
+        line = f.readline()
+        while line and err == 0:
+            line = f.readline()
+            if ' error' in line:
+                err = 1
+                print('Error in job script.\n')
+                print(line)
+
+    if err == 0:
+
+        with open(os.path.join(path_to_perf_files, 'slurm-' + str(jobid) + '-' + filename_input + '.out')) as f:
+
+            # create empty data frame from column set
+            perf_results_df = pd.DataFrame(columns=cols)
+            # set data types for new data frame as defined above
+            for k, v in cols_dtypes.items():
+                perf_results_df[k] = perf_results_df[k].astype(v)
+
+            lines = f.readlines()
+            i = 0
+            while i < len(lines)-1:
+
+                # search for next program execution
+                if '--> Iteration' in lines[i]:
+                    while '--> Iteration' in lines[i] and '--> Iteration' in lines[i+1]:
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+                    # iteration scheme ended:
+                    # split on empty spaces and take number of iterations in third place
+                    its = int(lines[i].split()[2][:-1])
+                    while '-norm' not in lines[i]:
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+                    norm2 = -1
+                    if '2-norm' in lines[i]:
+                        norm2 = float(lines[i].split()[-1])
+                    else:
+                        sys.exit('Error. 2-Norm information not available.')
+                    norminf = -1
+                    if 'inf-norm' in lines[i+1]:
+                        norminf = float(lines[i+1].split()[-1])
+                    else:
+                        sys.exit('Error. Inf-Norm information not available.')
+
+                    while 't_setup' not in lines[i]:
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+                    # global time print out, split line by line of information and add to dataframe
+                    time_dict = dict()
+                    while 'Total_execution_time' not in lines[i-1]:
+                        timings = lines[i].split()
+                        for j in range(0, len(timings), 2):
+                            # remove last char (':') from string and convert timing to floating point
+                            time_dict[timings[j][:-1]
+                                      ] = float(timings[j+1].replace(',', ''))
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+
+                    # iterate until LIKWID printout
+                    while 'Group 1: ' not in lines[i]:
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+                    benchmark = lines[i].split()[-1]
+                    # number of LIKWID outputs gives number of used cores
+                    cores = str.count(lines[i+2], 'HWThread')
+
+                    cores_used.append(cores)
+
+                    search_term_likwid = ''
+                    if benchmark == 'FLOPS_DP':
+                        search_term_likwid = 'DP [MFLOP/s]'
+                    elif benchmark == 'CACHES':
+                        search_term_likwid = 'Memory bandwidth [MBytes/s]'
+                    else:
+                        sys.exit('Error. Benchmark type not supported.')
+
+                    if cores > 1:
+                        search_term_likwid = search_term_likwid + ' STAT'
+
+                    while search_term_likwid not in lines[i]:
+                        i = i+1
+                        if i > len(lines)-1:
+                            sys.exit(
+                                'End of file reached without finding output.')
+
+                    benchmark_result = float(lines[i].split('|')[2].split()[0])
+
+                    # with the number of used cores determined, check if
+                    # line for this result is already present
+                    row_setting_avail = (
+                        perf_results_df[cols_problem[0]] == problem) & (
+                        perf_results_df[cols_problem[1]] == nr_exp) & (
+                        perf_results_df[cols_problem[2]] == mod_pk) & (
+                        perf_results_df[cols_cluster[0]] == nodes) & (
+                        perf_results_df[cols_cluster[1]] == ranks) & (
+                        perf_results_df[cols_cluster[2]] == cores)
+
+                    if row_setting_avail.values.sum() <= 1:  # add new line
+                        # add to all rows fix identifiers for:
+                        #   problem, exponent in r, modified coordinates, extrapolation
+                        #   number of nodes and ranks of the job
+                        # as well as variable numbers for other columns
+                        perf_results_df.loc[len(perf_results_df.index),
+                                            cols_problem + cols_cluster +
+                                            cols_convergence + [benchmark]] = [
+                            problem, nr_exp, mod_pk, extrapolation,
+                            nodes, ranks, cores, its, norm2, norminf,
+                            benchmark_result]
+                        # add timings
+                        for time_col, time_val in time_dict.items():
+                            perf_results_df.loc[len(
+                                perf_results_df.index)-1, time_col] = time_val
+                    elif row_setting_avail.values.sum() > 1:  # error
+                        sys.exit(
+                            'Error. More than two lines corresponds to criterion.'
+                            ' We can only have one line for each benchmark: FLOPS_DP and CACHES.')
+
+                i += 1
+            # end while
+
+            perf_results_df.to_csv(
+                os.path.join(
+                    path_to_perf_files, filename_input + '_benchmarks.csv'),
+                index=False)
+        # close file
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/main.cpp b/src/main.cpp
index c64100b5..3b1c8b9e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -227,7 +227,7 @@ int main(int argc, char* argv[])
             error = 1;
         }
         // if (gyro::icntl[Param::verbose] > 0)
-        std::cout << "Total execution time: " << TOC
+        std::cout << "Total_execution_time: " << TOC
                   << "\n---------------------------------------------------------------------------\n";
     }
 
@@ -305,7 +305,7 @@ int main(int argc, char* argv[])
             error = 1;
         }
         if (gyro::icntl[Param::verbose] > 0)
-            std::cout << "Total execution time: " << TOC << "\n";
+            std::cout << "Total_execution_time: " << TOC << "\n";
 
         std::cout << "\n\n\n\n\n";
 
@@ -375,7 +375,7 @@ int main(int argc, char* argv[])
                     error = 1;
                 }
                 if (gyro::icntl[Param::verbose] > 0)
-                    std::cout << "Total execution time: " << TOC << "\n";
+                    std::cout << "Total_execution_time: " << TOC << "\n";
             }
             // }
         }
@@ -454,7 +454,7 @@ int main(int argc, char* argv[])
                                 error = 1;
                             }
                             if (gyro::icntl[Param::verbose] > 0)
-                                std::cout << "Total execution time: " << TOC << "\n";
+                                std::cout << "Total_execution_time: " << TOC << "\n";
                         }
                     }
                 }
diff --git a/src/polar_multigrid.cpp b/src/polar_multigrid.cpp
index 5cd6ed05..8c79e02f 100644
--- a/src/polar_multigrid.cpp
+++ b/src/polar_multigrid.cpp
@@ -186,7 +186,7 @@ void gmgpolar::polar_multigrid()
                 std::cout << "\nt_setup: " << t_setup << ", t_build: " << t_build << ", t_facto_Ac: " << t_facto_Ac
                           << ", t_build_P: " << t_build_P << ", t_build_Asc: " << t_build_Asc
                           << ", t_facto_Asc: " << t_facto_Asc << "\n";
-                std::cout << "t_total (fine): " << t_total << ", t_smoothing: " << t_smoothing
+                std::cout << "t_total_(fine): " << t_total << ", t_smoothing: " << t_smoothing
                           << ", t_residual: " << t_residual << ", t_restriction: " << t_restriction
                           << ", t_Ac: " << t_Ac << ", t_prolongation: " << t_prolongation
                           << ", t_fine_residual: " << t_fine_residual << ", t_error: " << t_error << "\n";
diff --git a/src/test_cases/CartesianR2GyroSonnendruckerCircular.cpp b/src/test_cases/CartesianR2GyroSonnendruckerCircular.cpp
index 77183700..7ce4cb41 100644
--- a/src/test_cases/CartesianR2GyroSonnendruckerCircular.cpp
+++ b/src/test_cases/CartesianR2GyroSonnendruckerCircular.cpp
@@ -257,12 +257,12 @@ void CartesianR2GyroSonnendruckerCircular::rho_pole(double r, std::vector<double
 }
 /*........................................*/
 double CartesianR2GyroSonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2GyroSonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR2GyroSonnendruckerShafranov.cpp b/src/test_cases/CartesianR2GyroSonnendruckerShafranov.cpp
index c12504f1..b940c3c7 100644
--- a/src/test_cases/CartesianR2GyroSonnendruckerShafranov.cpp
+++ b/src/test_cases/CartesianR2GyroSonnendruckerShafranov.cpp
@@ -257,12 +257,12 @@ void CartesianR2GyroSonnendruckerShafranov::rho_pole(double r, std::vector<doubl
 }
 /*........................................*/
 double CartesianR2GyroSonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2GyroSonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR2GyroSonnendruckerTriangular.cpp b/src/test_cases/CartesianR2GyroSonnendruckerTriangular.cpp
index a6b543c1..3e744518 100644
--- a/src/test_cases/CartesianR2GyroSonnendruckerTriangular.cpp
+++ b/src/test_cases/CartesianR2GyroSonnendruckerTriangular.cpp
@@ -257,12 +257,12 @@ void CartesianR2GyroSonnendruckerTriangular::rho_pole(double r, std::vector<doub
 }
 /*........................................*/
 double CartesianR2GyroSonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2GyroSonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR2SonnendruckerCircular.cpp b/src/test_cases/CartesianR2SonnendruckerCircular.cpp
index ff54a9c7..050e2693 100644
--- a/src/test_cases/CartesianR2SonnendruckerCircular.cpp
+++ b/src/test_cases/CartesianR2SonnendruckerCircular.cpp
@@ -257,12 +257,12 @@ void CartesianR2SonnendruckerCircular::rho_pole(double r, std::vector<double> co
 }
 /*........................................*/
 double CartesianR2SonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2SonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR2SonnendruckerShafranov.cpp b/src/test_cases/CartesianR2SonnendruckerShafranov.cpp
index b77b28f2..b48e116e 100644
--- a/src/test_cases/CartesianR2SonnendruckerShafranov.cpp
+++ b/src/test_cases/CartesianR2SonnendruckerShafranov.cpp
@@ -257,12 +257,12 @@ void CartesianR2SonnendruckerShafranov::rho_pole(double r, std::vector<double> c
 }
 /*........................................*/
 double CartesianR2SonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2SonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR2SonnendruckerTriangular.cpp b/src/test_cases/CartesianR2SonnendruckerTriangular.cpp
index 19987f49..718408b0 100644
--- a/src/test_cases/CartesianR2SonnendruckerTriangular.cpp
+++ b/src/test_cases/CartesianR2SonnendruckerTriangular.cpp
@@ -257,12 +257,12 @@ void CartesianR2SonnendruckerTriangular::rho_pole(double r, std::vector<double>
 }
 /*........................................*/
 double CartesianR2SonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR2SonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6GyroSonnendruckerCircular.cpp b/src/test_cases/CartesianR6GyroSonnendruckerCircular.cpp
index 9255d1f9..8d9eb3c7 100644
--- a/src/test_cases/CartesianR6GyroSonnendruckerCircular.cpp
+++ b/src/test_cases/CartesianR6GyroSonnendruckerCircular.cpp
@@ -1,8 +1,8 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroSonnendruckerCircular.h"
 #include <stdlib.h>
 #include <math.h>
 
-
 /*........................................*/
 double CartesianR6GyroSonnendruckerCircular::x(double r, double theta, double unused_1, double unused_2, double Rmax) const
 {
@@ -257,12 +257,12 @@ void CartesianR6GyroSonnendruckerCircular::rho_pole(double r, std::vector<double
 }
 /*........................................*/
 double CartesianR6GyroSonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6GyroSonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
@@ -275,7 +275,7 @@ double CartesianR6GyroSonnendruckerCircular::coeffs2(double r, double Rmax) cons
 }
 /*........................................*/
 void CartesianR6GyroSonnendruckerCircular::coeffs2(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ 
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = pow((0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111)), (double)((-1)));
diff --git a/src/test_cases/CartesianR6GyroSonnendruckerShafranov.cpp b/src/test_cases/CartesianR6GyroSonnendruckerShafranov.cpp
index 54a6c5bc..99b56e8b 100644
--- a/src/test_cases/CartesianR6GyroSonnendruckerShafranov.cpp
+++ b/src/test_cases/CartesianR6GyroSonnendruckerShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroSonnendruckerShafranov.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void CartesianR6GyroSonnendruckerShafranov::rho_pole(double r, std::vector<doubl
 }
 /*........................................*/
 double CartesianR6GyroSonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6GyroSonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6GyroSonnendruckerTriangular.cpp b/src/test_cases/CartesianR6GyroSonnendruckerTriangular.cpp
index 00067267..e58741fc 100644
--- a/src/test_cases/CartesianR6GyroSonnendruckerTriangular.cpp
+++ b/src/test_cases/CartesianR6GyroSonnendruckerTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroSonnendruckerTriangular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void CartesianR6GyroSonnendruckerTriangular::rho_pole(double r, std::vector<doub
 }
 /*........................................*/
 double CartesianR6GyroSonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6GyroSonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6GyroZoniCircular.cpp b/src/test_cases/CartesianR6GyroZoniCircular.cpp
index fde42da7..22a3610a 100644
--- a/src/test_cases/CartesianR6GyroZoniCircular.cpp
+++ b/src/test_cases/CartesianR6GyroZoniCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6GyroZoniShafranov.cpp b/src/test_cases/CartesianR6GyroZoniShafranov.cpp
index d7b6de0f..4fc457ef 100644
--- a/src/test_cases/CartesianR6GyroZoniShafranov.cpp
+++ b/src/test_cases/CartesianR6GyroZoniShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6GyroZoniShiftedCircular.cpp b/src/test_cases/CartesianR6GyroZoniShiftedCircular.cpp
index 14e843bf..dd8f07e8 100644
--- a/src/test_cases/CartesianR6GyroZoniShiftedCircular.cpp
+++ b/src/test_cases/CartesianR6GyroZoniShiftedCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniShiftedCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6GyroZoniShiftedShafranov.cpp b/src/test_cases/CartesianR6GyroZoniShiftedShafranov.cpp
index 4d780756..59c024e1 100644
--- a/src/test_cases/CartesianR6GyroZoniShiftedShafranov.cpp
+++ b/src/test_cases/CartesianR6GyroZoniShiftedShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniShiftedShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6GyroZoniShiftedTriangular.cpp b/src/test_cases/CartesianR6GyroZoniShiftedTriangular.cpp
index 816002a2..e13c99cc 100644
--- a/src/test_cases/CartesianR6GyroZoniShiftedTriangular.cpp
+++ b/src/test_cases/CartesianR6GyroZoniShiftedTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniShiftedTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6GyroZoniTriangular.cpp b/src/test_cases/CartesianR6GyroZoniTriangular.cpp
index 11ee7407..760e29f6 100644
--- a/src/test_cases/CartesianR6GyroZoniTriangular.cpp
+++ b/src/test_cases/CartesianR6GyroZoniTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6GyroZoniTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6PoissonCircular.cpp b/src/test_cases/CartesianR6PoissonCircular.cpp
index 5ad5e880..aabe16dd 100644
--- a/src/test_cases/CartesianR6PoissonCircular.cpp
+++ b/src/test_cases/CartesianR6PoissonCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6PoissonCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6PoissonShafranov.cpp b/src/test_cases/CartesianR6PoissonShafranov.cpp
index f6d03db7..7eb01852 100644
--- a/src/test_cases/CartesianR6PoissonShafranov.cpp
+++ b/src/test_cases/CartesianR6PoissonShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6PoissonShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6PoissonTriangular.cpp b/src/test_cases/CartesianR6PoissonTriangular.cpp
index 9f96e861..82cc184d 100644
--- a/src/test_cases/CartesianR6PoissonTriangular.cpp
+++ b/src/test_cases/CartesianR6PoissonTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6PoissonTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6SonnendruckerCircular.cpp b/src/test_cases/CartesianR6SonnendruckerCircular.cpp
index cec53984..290a990c 100644
--- a/src/test_cases/CartesianR6SonnendruckerCircular.cpp
+++ b/src/test_cases/CartesianR6SonnendruckerCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6SonnendruckerCircular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void CartesianR6SonnendruckerCircular::rho_pole(double r, std::vector<double> co
 }
 /*........................................*/
 double CartesianR6SonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6SonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6SonnendruckerShafranov.cpp b/src/test_cases/CartesianR6SonnendruckerShafranov.cpp
index 927f0291..63749cff 100644
--- a/src/test_cases/CartesianR6SonnendruckerShafranov.cpp
+++ b/src/test_cases/CartesianR6SonnendruckerShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6SonnendruckerShafranov.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void CartesianR6SonnendruckerShafranov::rho_pole(double r, std::vector<double> c
 }
 /*........................................*/
 double CartesianR6SonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6SonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6SonnendruckerTriangular.cpp b/src/test_cases/CartesianR6SonnendruckerTriangular.cpp
index fc4994f2..9312a869 100644
--- a/src/test_cases/CartesianR6SonnendruckerTriangular.cpp
+++ b/src/test_cases/CartesianR6SonnendruckerTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6SonnendruckerTriangular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void CartesianR6SonnendruckerTriangular::rho_pole(double r, std::vector<double>
 }
 /*........................................*/
 double CartesianR6SonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void CartesianR6SonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/CartesianR6ZoniCircular.cpp b/src/test_cases/CartesianR6ZoniCircular.cpp
index 841f1112..f7721f45 100644
--- a/src/test_cases/CartesianR6ZoniCircular.cpp
+++ b/src/test_cases/CartesianR6ZoniCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6ZoniShafranov.cpp b/src/test_cases/CartesianR6ZoniShafranov.cpp
index a1d0e707..03fbb0c2 100644
--- a/src/test_cases/CartesianR6ZoniShafranov.cpp
+++ b/src/test_cases/CartesianR6ZoniShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6ZoniShiftedCircular.cpp b/src/test_cases/CartesianR6ZoniShiftedCircular.cpp
index 2ecafe81..ff90dfd1 100644
--- a/src/test_cases/CartesianR6ZoniShiftedCircular.cpp
+++ b/src/test_cases/CartesianR6ZoniShiftedCircular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniShiftedCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6ZoniShiftedShafranov.cpp b/src/test_cases/CartesianR6ZoniShiftedShafranov.cpp
index e8b46d1a..390a3d20 100644
--- a/src/test_cases/CartesianR6ZoniShiftedShafranov.cpp
+++ b/src/test_cases/CartesianR6ZoniShiftedShafranov.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniShiftedShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6ZoniShiftedTriangular.cpp b/src/test_cases/CartesianR6ZoniShiftedTriangular.cpp
index df8725e5..f6eb3ed4 100644
--- a/src/test_cases/CartesianR6ZoniShiftedTriangular.cpp
+++ b/src/test_cases/CartesianR6ZoniShiftedTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniShiftedTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/CartesianR6ZoniTriangular.cpp b/src/test_cases/CartesianR6ZoniTriangular.cpp
index 635e78b6..42aaa013 100644
--- a/src/test_cases/CartesianR6ZoniTriangular.cpp
+++ b/src/test_cases/CartesianR6ZoniTriangular.cpp
@@ -1,3 +1,4 @@
+// CartesianR6 simulates solution (23) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "CartesianR6ZoniTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroSonnendruckerCircular.cpp b/src/test_cases/PolarR6GyroSonnendruckerCircular.cpp
index 12b41cbc..04b296e2 100644
--- a/src/test_cases/PolarR6GyroSonnendruckerCircular.cpp
+++ b/src/test_cases/PolarR6GyroSonnendruckerCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroSonnendruckerCircular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6GyroSonnendruckerCircular::rho_pole(double r, std::vector<double> co
 }
 /*........................................*/
 double PolarR6GyroSonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6GyroSonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6GyroSonnendruckerShafranov.cpp b/src/test_cases/PolarR6GyroSonnendruckerShafranov.cpp
index 7f45111c..45337b32 100644
--- a/src/test_cases/PolarR6GyroSonnendruckerShafranov.cpp
+++ b/src/test_cases/PolarR6GyroSonnendruckerShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroSonnendruckerShafranov.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6GyroSonnendruckerShafranov::rho_pole(double r, std::vector<double> c
 }
 /*........................................*/
 double PolarR6GyroSonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6GyroSonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6GyroSonnendruckerTriangular.cpp b/src/test_cases/PolarR6GyroSonnendruckerTriangular.cpp
index 64a565ee..a2b5fd84 100644
--- a/src/test_cases/PolarR6GyroSonnendruckerTriangular.cpp
+++ b/src/test_cases/PolarR6GyroSonnendruckerTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroSonnendruckerTriangular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6GyroSonnendruckerTriangular::rho_pole(double r, std::vector<double>
 }
 /*........................................*/
 double PolarR6GyroSonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6GyroSonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6GyroZoniCircular.cpp b/src/test_cases/PolarR6GyroZoniCircular.cpp
index cffa89c6..95086adc 100644
--- a/src/test_cases/PolarR6GyroZoniCircular.cpp
+++ b/src/test_cases/PolarR6GyroZoniCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroZoniShafranov.cpp b/src/test_cases/PolarR6GyroZoniShafranov.cpp
index aaa43189..e5fe6c0c 100644
--- a/src/test_cases/PolarR6GyroZoniShafranov.cpp
+++ b/src/test_cases/PolarR6GyroZoniShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroZoniShiftedCircular.cpp b/src/test_cases/PolarR6GyroZoniShiftedCircular.cpp
index bf189091..858a7f5e 100644
--- a/src/test_cases/PolarR6GyroZoniShiftedCircular.cpp
+++ b/src/test_cases/PolarR6GyroZoniShiftedCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniShiftedCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroZoniShiftedCulham.cpp b/src/test_cases/PolarR6GyroZoniShiftedCulham.cpp
index 86ba7c9d..a1b72161 100644
--- a/src/test_cases/PolarR6GyroZoniShiftedCulham.cpp
+++ b/src/test_cases/PolarR6GyroZoniShiftedCulham.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniShiftedCulham.h"
 #include <stdlib.h>
 #include <array>
diff --git a/src/test_cases/PolarR6GyroZoniShiftedShafranov.cpp b/src/test_cases/PolarR6GyroZoniShiftedShafranov.cpp
index 146885d7..99d7520e 100644
--- a/src/test_cases/PolarR6GyroZoniShiftedShafranov.cpp
+++ b/src/test_cases/PolarR6GyroZoniShiftedShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniShiftedShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroZoniShiftedTriangular.cpp b/src/test_cases/PolarR6GyroZoniShiftedTriangular.cpp
index 7e7aecf2..fd0e8603 100644
--- a/src/test_cases/PolarR6GyroZoniShiftedTriangular.cpp
+++ b/src/test_cases/PolarR6GyroZoniShiftedTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniShiftedTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6GyroZoniTriangular.cpp b/src/test_cases/PolarR6GyroZoniTriangular.cpp
index a97f722f..037d654a 100644
--- a/src/test_cases/PolarR6GyroZoniTriangular.cpp
+++ b/src/test_cases/PolarR6GyroZoniTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6GyroZoniTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6PoissonCircular.cpp b/src/test_cases/PolarR6PoissonCircular.cpp
index cd1da441..5f2f5d6d 100644
--- a/src/test_cases/PolarR6PoissonCircular.cpp
+++ b/src/test_cases/PolarR6PoissonCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6PoissonCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6PoissonShafranov.cpp b/src/test_cases/PolarR6PoissonShafranov.cpp
index bf6db6a7..cc8c7bc4 100644
--- a/src/test_cases/PolarR6PoissonShafranov.cpp
+++ b/src/test_cases/PolarR6PoissonShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6PoissonShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6PoissonTriangular.cpp b/src/test_cases/PolarR6PoissonTriangular.cpp
index 8ae65e37..0a42f162 100644
--- a/src/test_cases/PolarR6PoissonTriangular.cpp
+++ b/src/test_cases/PolarR6PoissonTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6PoissonTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6SonnendruckerCircular.cpp b/src/test_cases/PolarR6SonnendruckerCircular.cpp
index 5b9013eb..8bc1339f 100644
--- a/src/test_cases/PolarR6SonnendruckerCircular.cpp
+++ b/src/test_cases/PolarR6SonnendruckerCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6SonnendruckerCircular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6SonnendruckerCircular::rho_pole(double r, std::vector<double> const&
 }
 /*........................................*/
 double PolarR6SonnendruckerCircular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6SonnendruckerCircular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6SonnendruckerShafranov.cpp b/src/test_cases/PolarR6SonnendruckerShafranov.cpp
index b1efe2d1..4a735338 100644
--- a/src/test_cases/PolarR6SonnendruckerShafranov.cpp
+++ b/src/test_cases/PolarR6SonnendruckerShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6SonnendruckerShafranov.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6SonnendruckerShafranov::rho_pole(double r, std::vector<double> const
 }
 /*........................................*/
 double PolarR6SonnendruckerShafranov::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6SonnendruckerShafranov::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6SonnendruckerTriangular.cpp b/src/test_cases/PolarR6SonnendruckerTriangular.cpp
index 2a69ae95..7d10608e 100644
--- a/src/test_cases/PolarR6SonnendruckerTriangular.cpp
+++ b/src/test_cases/PolarR6SonnendruckerTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6SonnendruckerTriangular.h"
 #include <stdlib.h>
 #include <math.h>
@@ -257,12 +258,12 @@ void PolarR6SonnendruckerTriangular::rho_pole(double r, std::vector<double> cons
 }
 /*........................................*/
 double PolarR6SonnendruckerTriangular::coeffs1(double r, double Rmax) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r/Rmax) - 11.1111111111111);
 }
 /*........................................*/
 void PolarR6SonnendruckerTriangular::coeffs1(std::vector<double> const& r, double Rmax, std::vector<double>& sol) const
-{
+{ // see Kuehn et al. https://doi.org/10.1007/s10915-022-01802-1, Eq. (2.3) with Rmax=1.3
     for (std::size_t i=0; i < sol.size(); ++i)
     {
         sol[i] = 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r[i]/Rmax) - 11.1111111111111);
diff --git a/src/test_cases/PolarR6ZoniCircular.cpp b/src/test_cases/PolarR6ZoniCircular.cpp
index 8d0c3f6a..7b3c40ae 100644
--- a/src/test_cases/PolarR6ZoniCircular.cpp
+++ b/src/test_cases/PolarR6ZoniCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6ZoniShafranov.cpp b/src/test_cases/PolarR6ZoniShafranov.cpp
index d4b0c2fc..de87f64e 100644
--- a/src/test_cases/PolarR6ZoniShafranov.cpp
+++ b/src/test_cases/PolarR6ZoniShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6ZoniShiftedCircular.cpp b/src/test_cases/PolarR6ZoniShiftedCircular.cpp
index dcd04554..3cbf6e23 100644
--- a/src/test_cases/PolarR6ZoniShiftedCircular.cpp
+++ b/src/test_cases/PolarR6ZoniShiftedCircular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniShiftedCircular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6ZoniShiftedShafranov.cpp b/src/test_cases/PolarR6ZoniShiftedShafranov.cpp
index b2480620..473d174d 100644
--- a/src/test_cases/PolarR6ZoniShiftedShafranov.cpp
+++ b/src/test_cases/PolarR6ZoniShiftedShafranov.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniShiftedShafranov.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6ZoniShiftedTriangular.cpp b/src/test_cases/PolarR6ZoniShiftedTriangular.cpp
index 328d5a94..437fc5c7 100644
--- a/src/test_cases/PolarR6ZoniShiftedTriangular.cpp
+++ b/src/test_cases/PolarR6ZoniShiftedTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniShiftedTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/src/test_cases/PolarR6ZoniTriangular.cpp b/src/test_cases/PolarR6ZoniTriangular.cpp
index 4b3cf5f6..24de7daf 100644
--- a/src/test_cases/PolarR6ZoniTriangular.cpp
+++ b/src/test_cases/PolarR6ZoniTriangular.cpp
@@ -1,3 +1,4 @@
+// PolarR6 simulates solution (22) of Bourne et al. https://doi.org/10.1016/j.jcp.2023.112249
 #include "PolarR6ZoniTriangular.h"
 #include <stdlib.h>
 #include <math.h>
diff --git a/treat_output.py b/treat_output.py
index a2f9e1a3..8da4d0b8 100644
--- a/treat_output.py
+++ b/treat_output.py
@@ -99,7 +99,7 @@
         latex.write(error_norm_inf)
         latex.write(" & ")
 
-    time = line.find('Total execution time:')
+    time = line.find('Total_execution_time:')
     if (time != -1):
         #print(f'{line}')
         #execution time is found at: time+22, until the end