diff --git a/.gitignore b/.gitignore
index 0b32f91f21a82f38540d5e890279fb4e5a9c0278..c0a54e54b9940dc111b254aab4898b72e38eccb1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,3 +31,4 @@ wk/fort.90
.directory
checkpoint/
FM/
+wk/test*
diff --git a/local/dirs.inc b/local/dirs.inc
index 6ff72f9689ab58aa7d192eb00508a67ee6ea58e2..26854325c8f4153f8382245a36320ceeae0e0479 100644
--- a/local/dirs.inc
+++ b/local/dirs.inc
@@ -6,7 +6,9 @@ OBJDIR = $(PREFIX)/obj
LIBDIR = $(PREFIX)/lib
MODDIR = $(PREFIX)/mod
FMDIR = $(PREFIX)/FM
-FFTW3DIR = /usr/local/fftw3
+#FFTW3DIR = /usr/local/fftw3
+FFTW3DIR = /home/ahoffman/lib/fftw-3.3.8
CHCKPTDIR = $(PREFIX)/checkpoint
+FUTILS = /home/ahoffman/lib/futils_1.4/src
# Naming ideas : HeLaZ, MoNoLiT (Moment Non Linear Torroidal)
diff --git a/local/make.inc b/local/make.inc
index b1943c29340c8b72027fb998df7f7cc2e5dc1e74..4a88e7762f9d7751e5ddac2a8ed88cd90ff64e88 100644
--- a/local/make.inc
+++ b/local/make.inc
@@ -15,7 +15,7 @@ else
endif
#0/1 = turn parallel MG solver off/on
-USEPARMG=1
+USEPARMG=0
#0/1 = turn OpenMP compilation off/on
USEOPENMP=0
@@ -24,7 +24,7 @@ USEOPENMP=0
USEVERIFICATION=0
#Please set precompiler flags here
-FPPFLAGS=-DPARDISO=$(PARDISO) -DMUMPS=$(USEMUMPS) -DPARMG=$(USEPARMG) -DVERIFICATION=$(USEVERIFICATION)
+#FPPFLAGS=-DPARDISO=$(PARDISO) -DMUMPS=$(USEMUMPS) -DPARMG=$(USEPARMG) -DVERIFICATION=$(USEVERIFICATION)
################################################################
#
@@ -83,11 +83,13 @@ endif
################################################################
IDIRS := -I$(SPC_LOCAL)/include/$(OPTLEVEL)
+#IDIRS := -I$(FUTILS)/include/$(OPTLEVEL) -I$(HDF5)/include -I$(ZLIB)/include
#LIBS := -lbsplines -lfutils -lpppack -lpputils2 \
# -lhdf5_fortran -lhdf5 -lz -ldl -lpthread
LIBS := -lfutils -lhdf5_fortran -lhdf5 -lz -ldl -lpthread
LDIRS := -L$(SPC_LOCAL)/lib/$(OPTLEVEL) -L$(HDF5)/lib
+#LDIRS := -L$(FUTILS)/lib -L$(HDF5)/lib -L$(ZLIB)/lib
# Add Multiple-Precision Library
LIBS += -lfm
@@ -131,8 +133,8 @@ endif
ifdef FFTW3DIR
LIBS += -lfftw3
- LDIRS += -L/usr/local/fftw3/lib
- IDIRS += -I/usr/local/fftw3/include
+ LDIRS += -L$(FFTW3DIR)/lib64
+ IDIRS += -I$(FFTW3DIR)/include
endif
#
diff --git a/matlab/create_gif_surf.m b/matlab/create_gif_surf.m
new file mode 100644
index 0000000000000000000000000000000000000000..087b788c87053b3bb3822e41e27b84e4521f3e94
--- /dev/null
+++ b/matlab/create_gif_surf.m
@@ -0,0 +1,56 @@
+title1 = GIFNAME;
+FIGDIR = BASIC.RESDIR;
+GIFNAME = [FIGDIR, GIFNAME,'.gif'];
+
+% Set colormap boundaries
+hmax = max(max(max(FIELD)));
+hmin = min(min(min(FIELD)));
+
+flag = 0;
+if hmax == hmin
+ disp('Warning : h = hmin = hmax = const')
+else
+% Setup figure frame
+fig = figure('Color','white','Position', [100, 100, 400, 400]);
+ surf(X,Y,FIELD(:,:,1)); % to set up
+ colormap gray
+ axis tight manual % this ensures that getframe() returns a consistent size
+ if INTERP
+ shading interp;
+ end
+ in = 1;
+ nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
+ for n = FRAMES % loop over selected frames
+ scale = max(max(abs(FIELD(:,:,n))));
+ pclr = surf(X,Y,FIELD(:,:,n)/scale); % frame plot
+ if INTERP
+ shading interp;
+ end
+ set(pclr, 'edgecolor','none'); axis square;
+% caxis([min(min(FIELD(:,:,n))),max(max(FIELD(:,:,n)))]);
+ xlabel(XNAME); ylabel(YNAME); %colorbar;
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scaling = ',sprintf('%.1e',scale)]);
+ drawnow
+ % Capture the plot as an image
+ frame = getframe(fig);
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,32);
+ % Write to the GIF File
+ if in == 1
+ imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+ in = in + 1;
+ end
+ disp(' ')
+ disp(['Gif saved @ : ',GIFNAME])
+end
+
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index f9f3fa748ffe478a1f00dc77b41102fbc17f0638..87fdab4b6dedeac1995b5dd1ecef4ae3624cbc74 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -20,6 +20,7 @@ fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
fprintf(fid,[' Lz = ', num2str(GRID.Lz),'\n']);
fprintf(fid,[' kpar = ', num2str(GRID.kpar),'\n']);
+fprintf(fid,[' CANCEL_ODD_P = ', num2str(GRID.CANCEL_ODD_P),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 0c27b59362c14d96f5e61747d0a964c0829b925c..ef7d2e463336550ce2f93a37468e35c762dfd510 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -11,9 +11,7 @@ CONTAINS
USE time_integration
use prec_const
IMPLICIT NONE
-
CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
-
END SUBROUTINE advance_time_level
diff --git a/src/auxval.F90 b/src/auxval.F90
index e3cfe31583dd3929960d1b63ab345bb9f8ad28c3..952057927bebc2c2c3128e1b450429d2c2373393 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -9,11 +9,18 @@ subroutine auxval
IMPLICIT NONE
INTEGER :: irows,irowe, irow, icol
- WRITE(*,*) '=== Set auxiliary values ==='
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
CALL set_krgrid
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
CALL set_kzgrid
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
CALL set_pj
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
CALL memory ! Allocate memory for global arrays
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index eb0541ed487ff0a7cf6402dc6be68271c6a3b216..7499fd80b1b12ea76f1d958c78bf613d6930dfa7 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -4,6 +4,8 @@ MODULE basic
use prec_const
IMPLICIT none
+ ! INCLUDE 'fftw3-mpi.f03'
+
INTEGER :: nrun = 1 ! Number of time steps to run
real(dp) :: tmax = 100000.0 ! Maximum simulation time
real(dp) :: dt = 1.0 ! Time step
@@ -16,6 +18,11 @@ MODULE basic
LOGICAL :: nlend = .FALSE. ! Signal end of run
INTEGER :: ierr ! flag for MPI error
+ INTEGER :: my_id ! identification number of current process
+ INTEGER :: num_procs ! number of MPI processes
+ INTEGER :: grp_world ! world communicator
+ INTEGER :: comm_world
+
INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose)
INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose)
INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose)
@@ -25,8 +32,11 @@ MODULE basic
INTEGER :: lu_in = 90 ! File duplicated from STDIN
INTEGER :: lu_job = 91 ! myjob file
- real :: start, finish ! To measure computation time
- real :: start_est, finish_est, time_est
+ ! To measure computation time
+ real :: start, finish
+ real(dp) :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield, t0_step
+ real(dp) :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield, t1_step
+ real(dp) :: tc_rhs, tc_adv_field, tc_poisson, tc_Sapj, tc_diag, tc_checkfield, tc_step
INTERFACE allocate_array
MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4
@@ -48,6 +58,10 @@ CONTAINS
READ(lu_in,basic)
+ !Init cumulative timers
+ tc_rhs = 0.;tc_adv_field = 0.; tc_poisson = 0.
+ tc_Sapj = 0.; tc_diag = 0.; tc_checkfield = 0.
+
END SUBROUTINE basic_data
!================================================================================
SUBROUTINE daytim(str)
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 6618ce0c664cb962d00c51898c7bfbcc10911e1b..43693ae25c512c215ab5d0c5d9d63d29a71c18ca 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -20,6 +20,9 @@ SUBROUTINE compute_Sapj
REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn
REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
+ ! Execution time start
+ CALL cpu_time(t0_Sapj)
+
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
@@ -30,8 +33,8 @@ SUBROUTINE compute_Sapj
nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
- krloope1: DO ikr = 1,Nkr ! Loop over kr
- kzloope1: DO ikz = 1,Nkz ! Loop over kz
+ krloope1: DO ikr = ikrs,ikre ! Loop over kr
+ kzloope1: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
@@ -65,8 +68,8 @@ SUBROUTINE compute_Sapj
Sepj(ip,ij,:,:) = Sepj(ip,ij,:,:) + CONV ! Add it to Sepj
IF ( NON_LIN ) THEN ! Fully non linear term ikz phi * ikr Napj
- krloope2: DO ikr = 1,Nkr ! Loop over kr
- kzloope2: DO ikz = 1,Nkz ! Loop over kz
+ krloope2: DO ikr = ikrs,ikre ! Loop over kr
+ kzloope2: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
@@ -107,8 +110,8 @@ SUBROUTINE compute_Sapj
nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
- krloopi1: DO ikr = 1,Nkr ! Loop over kr
- kzloopi1: DO ikz = 1,Nkz ! Loop over kz
+ krloopi1: DO ikr = ikrs,ikre ! Loop over kr
+ kzloopi1: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
@@ -140,8 +143,8 @@ SUBROUTINE compute_Sapj
CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier
Sipj(ip,ij,:,:) = Sipj(ip,ij,:,:) + CONV ! Add it to Sipj
- krloopi2: DO ikr = 1,Nkr ! Loop over kr
- kzloopi2: DO ikz = 1,Nkz ! Loop over kz
+ krloopi2: DO ikr = ikrs,ikre ! Loop over kr
+ kzloopi2: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
@@ -179,4 +182,9 @@ SUBROUTINE compute_Sapj
ENDDO jloopi
ENDDO ploopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ ! Execution time end
+ CALL cpu_time(t1_Sapj)
+ tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
+
END SUBROUTINE compute_Sapj
diff --git a/src/control.F90 b/src/control.F90
index 92638b0148dac00f484800c7675fbe3458cb4935..3a78b912e722324fbac8a08fbdba987f473a9ac8 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -9,78 +9,86 @@ SUBROUTINE control
!________________________________________________________________________________
! 1. Prologue
! 1.1 Initialize the parallel environment
- WRITE(*,*) 'Initialize MPI...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Initialize MPI...'
CALL ppinit
- WRITE(*,'(a/)') '...MPI initialized'
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...MPI initialized'
CALL daytim('Start at ')
! 1.2 Define data specific to run
- WRITE(*,*) 'Load basic data...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Load basic data...'
CALL basic_data
- WRITE(*,'(a/)') '...basic data loaded.'
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...basic data loaded.'
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
! 1.3 Read input parameters from input file
- WRITE(*,*) 'Read input parameters...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Read input parameters...'
CALL readinputs
- WRITE(*,'(a/)') '...input parameters read'
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...input parameters read'
+
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
! 1.4 Set auxiliary values (allocate arrays, set grid, ...)
- WRITE(*,*) 'Calculate auxval...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Calculate auxval...'
CALL auxval
- WRITE(*,'(a/)') '...auxval calculated'
- !
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...auxval calculated'
+
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
! 1.5 Initial conditions
- WRITE(*,*) 'Create initial state...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Create initial state...'
CALL inital
- WRITE(*,'(a/)') '...initial state created'
-
- ! ! 1.5.1 Computation time estimation
- ! WRITE(*,*) 'Estimation of computation time...'
- ! CALL cpu_time(start_est)
- ! CALL stepon
- ! CALL cpu_time(finish_est)
- ! time_est = 4.*tmax/dt*(finish_est-start_est)
- ! CALL display_h_min_s(time_est)
- ! WRITE(*,*) '... reinitializing'
- ! CALL inital
- ! WRITE(*,'(a/)') '... done'
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial state created'
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
! 1.6 Initial diagnostics
- WRITE(*,*) 'Initial diagnostics...'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Initial diagnostics...'
CALL diagnose(0)
- WRITE(*,'(a/)') '...initial diagnostics done'
- !
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial diagnostics done'
+
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
CALL FLUSH(stdout)
!________________________________________________________________________________
- WRITE(*,*) 'Time integration loop..'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Time integration loop..'
!________________________________________________________________________________
! 2. Main loop
DO
+CALL cpu_time(t0_step) ! Measuring time
+
step = step + 1
cstep = cstep + 1
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
CALL stepon
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
time = time + dt
CALL tesend
IF( nlend ) EXIT ! exit do loop
+CALL cpu_time(t0_diag) ! Measuring time
CALL diagnose(step)
+CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
+CALL cpu_time(t1_step); tc_step = tc_step + (t1_step - t0_step)
END DO
- WRITE(*,'(a/)') '...time integration done'
+ IF (my_id .EQ. 0) WRITE(*,'(a/)') '...time integration done'
!________________________________________________________________________________
! 9. Epilogue
CALL diagnose(-1)
CALL endrun
- CALL daytim('Done at ')
+ IF (my_id .EQ. 0) CALL daytim('Done at ')
+
CALL ppexit
END SUBROUTINE control
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 8c91e5603f893808ecaa14af699926e27eec9183..b0885f3d6d288900cf37739860cdb0d9f5597fc4 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -18,28 +18,33 @@ SUBROUTINE diagnose(kstep)
INTEGER, INTENT(in) :: kstep
INTEGER, parameter :: BUFSIZE = 2
INTEGER :: rank, dims(1) = (/0/)
- CHARACTER(len=256) :: str, fname
+ INTEGER :: cp_counter = 0
+ CHARACTER(len=256) :: str, fname,test_
+
!_____________________________________________________________________________
! 1. Initial diagnostics
IF ((kstep .EQ. 0)) THEN
- ! jobnum is either 0 from initialization or some integer values read from chkrst(0)
- IF(jobnum .LE. 99) THEN
- WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
- ELSE
- WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
- END IF
-
- ! 1.1 Initial run
+ ! Writing output filename
+ WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
+ ! 1.1 Initial run
+ ! Main output file creation
IF (write_doubleprecision) THEN
- CALL creatf(resfile, fidres, real_prec='d')
+ CALL creatf(resfile, fidres, real_prec='d', mpicomm=MPI_COMM_WORLD)
ELSE
- CALL creatf(resfile, fidres)
+ CALL creatf(resfile, fidres, mpicomm=MPI_COMM_WORLD)
END IF
+ IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
+ ! Checkpoint file creation
+ WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',jobnum,'.h5'
+ CALL creatf(rstfile, fidrst, real_prec='d', mpicomm=MPI_COMM_WORLD)
+ CALL creatg(fidrst, '/Basic', 'Basic data')
+ CALL creatg(fidrst, '/Basic/moments_e', 'electron moments')
+ CALL creatg(fidrst, '/Basic/moments_i', 'ion moments')
- WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
+ IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(rstfile), ' created'
CALL flush(6)
@@ -63,6 +68,17 @@ SUBROUTINE diagnose(kstep)
CALL creatg(fidres, "/files", "files")
CALL attach(fidres, "/files", "jobnum", jobnum)
+ ! Profiler time measurement
+ CALL creatg(fidres, "/profiler", "performance analysis")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs", "cumulative rhs computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field", "cumulative adv. fields computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_poisson", "cumulative poisson computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_Sapj", "cumulative Sapj computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
+
! var2d group (electro. pot., Ni00 moment)
rank = 0
CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
@@ -70,18 +86,30 @@ SUBROUTINE diagnose(kstep)
IF (write_Na00) THEN
CALL creatg(fidres, "/data/var2d/Ne00", "Ne00")
- CALL putarr(fidres, "/data/var2d/Ne00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var2d/Ne00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var2d/Ne00/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var2d/Ne00/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var2d/Ne00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
- CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
END IF
IF (write_phi) THEN
CALL creatg(fidres, "/data/var2d/phi", "phi")
- CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
END IF
! var5d group (moments)
@@ -90,37 +118,53 @@ SUBROUTINE diagnose(kstep)
CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
IF (write_moments) THEN
CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
- CALL putarr(fidres, "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e), "j_e",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_e/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_e/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e), "p_e", ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e), "j_e", ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var5d/moments_e/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var5d/moments_e/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var5d/moments_e/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
- CALL putarr(fidres, "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i), "p_i",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i), "j_i",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var5d/moments_i/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i), "p_i", ionode=0)
+ CALL putarr(fidres, "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i), "j_i", ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var5d/moments_i/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
END IF
IF (write_non_lin) THEN
CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
- CALL putarr(fidres, "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e), "j_e",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sepj/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sepj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e), "p_e", ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e), "j_e", ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var5d/Sepj/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var5d/Sepj/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var5d/Sepj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
- CALL putarr(fidres, "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i), "p_i",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i), "j_i",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sipj/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var5d/Sipj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i), "p_i", ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i), "j_i", ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, "/data/var5d/Sipj/coordkr", krarray(ikrs:ikre), "kr*rho_s0", ionode=0)
+ ELSE
+ CALL putarr(fidres, "/data/var5d/Sipj/coordkr", krarray(ikrs:ikre), "kr*rho_s0", pardim=1)
+ ENDIF
+ CALL putarr(fidres, "/data/var5d/Sipj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0", ionode=0)
END IF
! Add input namelist variables as attributes of /data/input, defined in srcinfo.h
- WRITE(*,*) 'VERSION=', VERSION
- WRITE(*,*) 'BRANCH=', BRANCH
- WRITE(*,*) 'AUTHOR=', AUTHOR
- WRITE(*,*) 'HOST=', HOST
+ IF (my_id .EQ. 0) WRITE(*,*) 'VERSION=', VERSION
+ IF (my_id .EQ. 0) WRITE(*,*) 'BRANCH=', BRANCH
+ IF (my_id .EQ. 0) WRITE(*,*) 'AUTHOR=', AUTHOR
+ IF (my_id .EQ. 0) WRITE(*,*) 'HOST=', HOST
WRITE(str,'(a,i2.2)') "/data/input"
@@ -138,6 +182,7 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, TRIM(str), "dt", dt)
CALL attach(fidres, TRIM(str), "tmax", tmax)
CALL attach(fidres, TRIM(str), "nrun", nrun)
+ CALL attach(fidres, TRIM(str), "cpu_time", -1)
CALL grid_outputinputs(fidres, str)
@@ -180,7 +225,7 @@ SUBROUTINE diagnose(kstep)
IF (kstep .GE. 0) THEN
! Terminal info
- IF (MOD(cstep, INT(1.0/dt)) == 0) THEN
+ IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN
WRITE(*,"(F5.0,A,F5.0)") time,"/",tmax
ENDIF
@@ -211,7 +256,8 @@ SUBROUTINE diagnose(kstep)
! 2.5 Backups
nsave_cp = INT(5/dt)
IF (MOD(cstep, nsave_cp) == 0) THEN
- CALL checkpoint_save
+ CALL checkpoint_save(cp_counter)
+ cp_counter = cp_counter + 1
ENDIF
!_____________________________________________________________________________
@@ -222,12 +268,14 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, "/data/input","cpu_time",finish-start)
! Display computational time cost
- CALL display_h_min_s(finish-start)
- ! WRITE(*,*) 'Time estimated :'
- ! CALL display_h_min_s(time_est)
+ IF (my_id .EQ. 0) CALL display_h_min_s(finish-start)
! Close all diagnostic files
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
CALL closef(fidres)
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+ CALL closef(fidrst)
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
END IF
END SUBROUTINE diagnose
@@ -241,12 +289,15 @@ SUBROUTINE diagnose_0d
USE prec_const
IMPLICIT NONE
- WRITE(*,'(a,1x,i7.7,a1,i7.7,20x,a,1pe10.3,10x,a,1pe10.3)') &
- '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time, 'dt =', dt
- WRITE(*,*)
- ! flush stdout of all ranks. Usually ONLY rank 0 should WRITE, but error messages might be written from other ranks as well
- CALL FLUSH(stdout)
+ CALL append(fidres, "/profiler/Tc_rhs", tc_rhs,ionode=0)
+ CALL append(fidres, "/profiler/Tc_adv_field", tc_adv_field,ionode=0)
+ CALL append(fidres, "/profiler/Tc_poisson", tc_poisson,ionode=0)
+ CALL append(fidres, "/profiler/Tc_Sapj", tc_Sapj,ionode=0)
+ CALL append(fidres, "/profiler/Tc_diag", tc_diag,ionode=0)
+ CALL append(fidres, "/profiler/Tc_checkfield",tc_checkfield,ionode=0)
+ CALL append(fidres, "/profiler/Tc_step", tc_step,ionode=0)
+ CALL append(fidres, "/profiler/time", time,ionode=0)
END SUBROUTINE diagnose_0d
@@ -289,7 +340,11 @@ CONTAINS
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d
- CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze),ionode=0)
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze), ionode=0)
+ ELSE
+ CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze), pardim=1)
+ ENDIF
CALL attach(fidres, dset_name, "time", time)
@@ -340,8 +395,11 @@ SUBROUTINE diagnose_5d
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
- CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze),ionode=0)
-
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), ionode=0)
+ ELSE
+ CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), pardim=3)
+ ENDIF
CALL attach(fidres, dset_name, "time", time)
END SUBROUTINE write_field5d_e
@@ -358,42 +416,65 @@ SUBROUTINE diagnose_5d
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
- CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze),ionode=0)
-
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), ionode=0)
+ ELSE
+ CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), pardim=3)
+ ENDIF
CALL attach(fidres, dset_name, "time", time)
END SUBROUTINE write_field5d_i
END SUBROUTINE diagnose_5d
-SUBROUTINE checkpoint_save
+SUBROUTINE checkpoint_save(cp_step)
USE basic
- USE grid
+ USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
USE diagnostics_par
- USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf
+ USE futils, ONLY: putarr,attach
USE model
USE initial_par
USE fields
USE time_integration
IMPLICIT NONE
-
- WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',jobnum,'.h5'
-
- CALL creatf(rstfile, fidrst, real_prec='d')
- CALL creatg(fidrst, '/Basic', 'Basic data')
- CALL attach(fidrst, '/Basic', 'cstep', cstep)
- CALL attach(fidrst, '/Basic', 'time', time)
- CALL attach(fidrst, '/Basic', 'jobnum', jobnum)
- CALL attach(fidrst, '/Basic', 'dt', dt)
- CALL attach(fidrst, '/Basic', 'iframe2d', iframe2d)
- CALL attach(fidrst, '/Basic', 'iframe5d', iframe5d)
+ INTEGER, INTENT(IN) :: cp_step
+ CHARACTER(LEN=50) :: dset_name
! Write state of system to restart file
- CALL putarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,&
- ikrs:ikre,ikzs:ikze,1),ionode=0)
- CALL putarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,&
- ikrs:ikre,ikzs:ikze,1),ionode=0)
- CALL closef(fidrst)
- WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved"
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", cp_step
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,&
+ ikrs:ikre,ikzs:ikze,1), ionode=0)
+ ELSE
+ CALL putarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,&
+ ikrs:ikre,ikzs:ikze,1), pardim=3)
+ ENDIF
+
+ CALL attach(fidrst, dset_name, 'cstep', cstep)
+ CALL attach(fidrst, dset_name, 'time', time)
+ CALL attach(fidrst, dset_name, 'jobnum', jobnum)
+ CALL attach(fidrst, dset_name, 'dt', dt)
+ CALL attach(fidrst, dset_name, 'iframe2d', iframe2d)
+ CALL attach(fidrst, dset_name, 'iframe5d', iframe5d)
+
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", cp_step
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,&
+ ikrs:ikre,ikzs:ikze,1), ionode=0)
+ ELSE
+ CALL putarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,&
+ ikrs:ikre,ikzs:ikze,1), pardim=3)
+ ENDIF
+
+ CALL attach(fidrst, dset_name, 'cstep', cstep)
+ CALL attach(fidrst, dset_name, 'time', time)
+ CALL attach(fidrst, dset_name, 'jobnum', jobnum)
+ CALL attach(fidrst, dset_name, 'dt', dt)
+ CALL attach(fidrst, dset_name, 'iframe2d', iframe2d)
+ CALL attach(fidrst, dset_name, 'iframe5d', iframe5d)
+
+ IF (my_id .EQ. 0) THEN
+ WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" updated"
+ ENDIF
END SUBROUTINE checkpoint_save
diff --git a/src/endrun.F90 b/src/endrun.F90
index 7a0d5a52429302f5c071c42a6d0160aed238df49..ce6e739d862f091cb4710185528d20f33fc03c3d 100644
--- a/src/endrun.F90
+++ b/src/endrun.F90
@@ -4,13 +4,13 @@ SUBROUTINE endrun
USE basic
use prec_const
IMPLICIT NONE
-
-
+
+
IF( nlend ) THEN
!----------------------------------------------------------------------
! 1. Normal end of run
WRITE(*,'(/a)') ' Normal exit'
-
+
!----------------------------------------------------------------------
! 2. Abnormal exit
ELSE
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 7617963113ef6c467c28a8186a5c1dfb45827a01..a3dcca39e15f6e07e37f33d60e644e013d1c87df 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -5,7 +5,8 @@ MODULE fourier
use, intrinsic :: iso_c_binding
implicit none
- INCLUDE 'fftw3.f03'
+ ! INCLUDE 'fftw3.f03'
+ INCLUDE 'fftw3-mpi.f03'
PRIVATE
PUBLIC :: fft_r2cc
@@ -14,6 +15,7 @@ MODULE fourier
PUBLIC :: ifft2_cc2r
PUBLIC :: convolve_2D_F2F
+
CONTAINS
SUBROUTINE fft_r2cc(fx_in, Fk_out)
@@ -29,8 +31,6 @@ MODULE fourier
END SUBROUTINE fft_r2cc
-
-
SUBROUTINE ifft_cc2r(Fk_in, fx_out)
IMPLICIT NONE
@@ -81,11 +81,10 @@ MODULE fourier
END SUBROUTINE ifft2_cc2r
-
!!! Convolution 2D Fourier to Fourier
! - Compute the convolution using the convolution theorem and MKL
SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
-
+ USE basic
USE grid, ONLY : AA_r, AA_z, Lr, Lz
IMPLICIT NONE
COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: F_2D, G_2D ! input fields
@@ -96,25 +95,90 @@ MODULE fourier
REAL :: a_r
! 2D inverse Fourier transform
- CALL ifft2_cc2r(F_2D,ff);
- CALL ifft2_cc2r(G_2D,gg);
+ IF ( num_procs .EQ. 1 ) THEN
+ CALL ifft2_cc2r(F_2D,ff)
+ CALL ifft2_cc2r(G_2D,gg)
+ ELSE
+ CALL ifft2_cc2r_mpi(F_2D,ff)
+ CALL ifft2_cc2r_mpi(G_2D,gg)
+ ENDIF
! Product in physical space
ffgg = ff * gg;
! 2D Fourier tranform
- CALL fft2_r2cc(ffgg,C_2D);
+ IF ( num_procs .EQ. 1 ) THEN
+ CALL fft2_r2cc(ffgg,C_2D)
+ ELSE
+ CALL fft2_r2cc_mpi(ffgg,C_2D)
+ ENDIF
! Anti aliasing (2/3 rule)
- DO ikr = 1,Nkr
+ DO ikr = ikrs,ikre
a_r = AA_r(ikr)
- DO ikz = 1,Nkz
+ DO ikz = ikzs,ikze
C_2D(ikr,ikz) = C_2D(ikr,ikz) * a_r * AA_z(ikz)
ENDDO
ENDDO
END SUBROUTINE convolve_2D_F2F
+!! MPI routines
+SUBROUTINE fft2_r2cc_mpi( ffx_in, FFk_out )
+
+ IMPLICIT NONE
+ REAL(dp), DIMENSION(Nr,Nz), INTENT(IN) :: ffx_in
+ COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out
+ REAL(dp), DIMENSION(Nkr,Nkz) :: tmp_r
+ type(C_PTR) :: plan
+ integer(C_INTPTR_T) :: L
+ integer(C_INTPTR_T) :: M
+
+ ! L = Nr; M = Nz
+ !
+ ! tmp_r = ffx_in
+ ! !!! 2D Forward FFT ________________________!
+ ! plan = fftw_mpi_plan_dft_r2c_2d(L, M, tmp_r, FFk_out, MPI_COMM_WORLD, FFTW_ESTIMATE)
+ !
+ ! if ((.not. c_associated(plan))) then
+ ! write(*,*) "plan creation error!!"
+ ! stop
+ ! end if
+ !
+ ! CALL fftw_mpi_execute_dft_r2c(plan,tmp_r,FFk_out)
+ ! CALL fftw_destroy_plan(plan)
+
+END SUBROUTINE fft2_r2cc_mpi
+
+
+
+SUBROUTINE ifft2_cc2r_mpi( FFk_in, ffx_out )
+
+ IMPLICIT NONE
+ COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: FFk_in
+ REAL(dp), DIMENSION(Nr,Nz), INTENT(OUT) :: ffx_out
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: tmp_c
+ type(C_PTR) :: plan
+ integer(C_INTPTR_T) :: L
+ integer(C_INTPTR_T) :: M
+
+ ! L = Nr; M = Nz
+ !
+ ! tmp_c = FFk_in
+ ! !!! 2D Backward FFT ________________________!
+ ! plan = fftw_mpi_plan_dft_c2r_2d(L, M, tmp_c, ffx_out, MPI_COMM_WORLD, FFTW_ESTIMATE)
+ !
+ ! if ((.not. c_associated(plan))) then
+ ! write(*,*) "plan creation error!!"
+ ! stop
+ ! end if
+ !
+ ! CALL fftw_mpi_execute_dft_c2r(plan,tmp_c,ffx_out)
+ ! CALL fftw_destroy_plan(plan)
+ !
+ ! ffx_out = ffx_out/Nr/Nz
+
+END SUBROUTINE ifft2_cc2r_mpi
! Empty set/free routines to switch easily with MKL DFTI
SUBROUTINE set_descriptors
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 1f9686eef4df8e6ea9f6c8298cf560c326bb5646..12f9adf0fda718ca175a388e44636263fca32bda 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -1,6 +1,8 @@
MODULE grid
! Grid module for spatial discretization
USE prec_const
+ USE basic
+
IMPLICIT NONE
PRIVATE
@@ -19,7 +21,8 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: Nkz = 16 ! Number of total internal grid points in kz
REAL(dp), PUBLIC, PROTECTED :: Lkz = 1._dp ! vertical length of the fourier box
REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component
-
+ LOGICAL, PUBLIC, PROTECTED :: CANCEL_ODD_P = .false. ! To cancel odd Hermite polynomials
+ INTEGER, PUBLIC, PROTECTED :: pskip = 0 ! variable to skip p degrees or not
! For Orszag filter
REAL(dp), PUBLIC, PROTECTED :: two_third_krmax
REAL(dp), PUBLIC, PROTECTED :: two_third_kzmax
@@ -62,12 +65,18 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER :: ip, ij
- ips_e = 1; ipe_e = pmaxe + 1
- ips_i = 1; ipe_i = pmaxi + 1
+ IF (CANCEL_ODD_P) THEN
+ pskip = 1
+ ELSE
+ pskip = 0
+ ENDIF
+
+ ips_e = 1; ipe_e = pmaxe/(1+pskip) + 1
+ ips_i = 1; ipe_i = pmaxi/(1+pskip) + 1
ALLOCATE(parray_e(ips_e:ipe_e))
ALLOCATE(parray_i(ips_i:ipe_i))
- DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO
- DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO
+ DO ip = ips_e,ipe_e; parray_e(ip) = (1+pskip)*(ip-1); END DO
+ DO ip = ips_i,ipe_i; parray_i(ip) = (1+pskip)*(ip-1); END DO
ijs_e = 1; ije_e = jmaxe + 1
ijs_i = 1; ije_i = jmaxi + 1
@@ -76,6 +85,7 @@ CONTAINS
DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
maxj = MAX(jmaxi, jmaxe)
+
END SUBROUTINE set_pj
SUBROUTINE set_krgrid
@@ -84,11 +94,25 @@ CONTAINS
INTEGER :: ikr
Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
- ! Start and END indices of grid
- ikrs = 1
- ikre = Nkr
+ ! ! Start and END indices of grid
+ IF ( Nkr .GT. 1 ) THEN
+ ikrs = my_id * (Nkr-1)/num_procs + 1
+ ikre = (my_id+1) * (Nkr-1)/num_procs
+ ELSE
+ ikrs = 1; ikre = Nkr
+ ENDIF
+
+ IF (my_id .EQ. num_procs-1) THEN
+ ikre = Nkr
+ ENDIF
+ WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre
+
! Grid spacings
- deltakr = 2._dp*PI/Lr
+ IF (Lr .GT. 0) THEN
+ deltakr = 2._dp*PI/Lr
+ ELSE
+ deltakr = 1.0
+ ENDIF
! Discretized kr positions ordered as dk*(0 1 2)
ALLOCATE(krarray(ikrs:ikre))
@@ -120,6 +144,9 @@ CONTAINS
! Start and END indices of grid
ikzs = 1
ikze = Nkz
+ ! ikzs = my_id * Nz/num_procs + 1
+ ! ikze = (my_id+1) * Nz/num_procs
+ WRITE(*,*) 'ID = ',my_id,' ikzs = ', ikzs, ' ikze = ', ikze
! Grid spacings
deltakz = 2._dp*PI/Lz
@@ -147,8 +174,6 @@ CONTAINS
! Put kz0RT to the nearest grid point on kz
ikz0KH = NINT(kr0KH/deltakr)+1
kr0KH = kzarray(ikz0KH)
- WRITE(*,*) 'ikz0KH = ', ikz0KH
- WRITE(*,*) 'kr0KH = ', kr0KH
END SUBROUTINE set_kzgrid
@@ -159,7 +184,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nr, Lr, Nz, Lz, kpar
+ Nr, Lr, Nz, Lz, kpar, CANCEL_ODD_P
READ(lu_in,grid)
END SUBROUTINE grid_readinputs
diff --git a/src/inital.F90 b/src/inital.F90
index f1fa2efe2a97aad6dbabe5867a932c2c3ce8359e..79d006893a1698b1b7420fc1ad6b4936456155a5 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -36,8 +36,9 @@ SUBROUTINE inital
ENDIF
!!!!!! Load the full coulomb collision operator coefficients !!!!!!
IF (CO .EQ. -1) THEN
- ! WRITE(*,*) '=== Load Full Coulomb matrix ==='
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load Full Coulomb matrix ==='
CALL load_FC_mat
+ IF (my_id .EQ. 0) WRITE(*,*) '..done'
ENDIF
END SUBROUTINE inital
@@ -72,6 +73,12 @@ SUBROUTINE init_moments
END DO
END DO
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ CALL RANDOM_NUMBER(noise)
+ moments_e( 1,1,1,ikz, :) = moments_e( 1,1,1,Nkz+2-ikz, :)
+ moments_i( 1,1,1,ikz, :) = moments_i( 1,1,1,Nkz+2-ikz, :)
+ END DO
+
ELSE
!**** Gaussian initialization (Hakim 2017) *********************************
! sigma = 5._dp ! Gaussian sigma
@@ -90,32 +97,40 @@ SUBROUTINE init_moments
!**** Broad noise initialization *******************************************
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
+
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
- moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
- ! * AA_r(ikr) * AA_z(ikz)
+ moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
END DO
- DO ikz=2,Nkz/2 !symmetry at kr = 0
- CALL RANDOM_NUMBER(noise)
- moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
- END DO
+
+ ! IF ( ikrs .EQ. 1 ) THEN
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ CALL RANDOM_NUMBER(noise)
+ moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
+ END DO
+ ! ENDIF
+
END DO
END DO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
+
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
CALL RANDOM_NUMBER(noise)
- moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
- ! * AA_r(ikr) * AA_z(ikz)
+ moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
END DO
- DO ikz=2,Nkz/2 !symmetry at kr = 0
- CALL RANDOM_NUMBER(noise)
- moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
- END DO
+
+ ! IF ( ikrs .EQ. 1 ) THEN
+ DO ikz=2,Nkz/2 !symmetry at kr = 0
+ CALL RANDOM_NUMBER(noise)
+ moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
+ END DO
+ ! ENDIF
+
END DO
END DO
@@ -128,29 +143,59 @@ END SUBROUTINE init_moments
!******************************************************************************!
SUBROUTINE load_cp
USE basic
- USE futils, ONLY: openf, closef, getarr, getatt
+ USE futils, ONLY: openf, closef, getarr, getatt, isgroup, isdataset
USE grid
USE fields
USE diagnostics_par
USE time_integration
IMPLICIT NONE
+ INTEGER :: rank, sz_, n_
+ INTEGER :: dims(1) = (/0/)
+ CHARACTER(LEN=50) :: dset_name
+
WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5'
WRITE(*,'(3x,a)') "Resume from previous run"
CALL openf(rstfile, fidrst)
- CALL getatt(fidrst, '/Basic', 'cstep', cstep)
- CALL getatt(fidrst, '/Basic', 'time', time)
- CALL getatt(fidrst, '/Basic', 'jobnum', jobnum)
- jobnum = jobnum+1
- CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d)
- CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d)
- iframe2d = iframe2d-1; iframe5d = iframe5d-1
-
- ! Read state of system from restart file
- CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0)
- CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0)
+
+ IF (isgroup(fidrst,'/Basic/moments_e')) THEN
+ n_ = 0
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
+ DO WHILE (isdataset(fidrst, dset_name))
+ n_ = n_ + 1
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
+ ENDDO
+ n_ = n_ - 1
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_
+ CALL getatt(fidrst, dset_name, 'cstep', cstep)
+ CALL getatt(fidrst, dset_name, 'time', time)
+ CALL getatt(fidrst, dset_name, 'jobnum', jobnum)
+ jobnum = jobnum+1
+ CALL getatt(fidrst, dset_name, 'iframe2d',iframe2d)
+ CALL getatt(fidrst, dset_name, 'iframe5d',iframe5d)
+ iframe2d = iframe2d-1; iframe5d = iframe5d-1
+
+ ! Read state of system from restart file
+ CALL getarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0)
+ WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_
+ CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0)
+
+ ELSE
+ CALL getatt(fidrst, '/Basic', 'cstep', cstep)
+ CALL getatt(fidrst, '/Basic', 'time', time)
+ CALL getatt(fidrst, '/Basic', 'jobnum', jobnum)
+ jobnum = jobnum+1
+ CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d)
+ CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d)
+ iframe2d = iframe2d-1; iframe5d = iframe5d-1
+
+ ! Read state of system from restart file
+ CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0)
+ CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0)
+ ENDIF
+
CALL closef(fidrst)
WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 35cc155a37291ed0f793d8d643b6095a80921006..cb5591fa0d1986508019da86f9b8eb687b920d59 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -10,8 +10,8 @@ SUBROUTINE moments_eq_rhs
USE utility, ONLY : is_nan
IMPLICIT NONE
- INTEGER :: ip2, ij2 ! loops indices
- REAL(dp) :: ip_dp, ij_dp
+ INTEGER :: ip2, ij2, p_int, j_int, p2_int, j2_int ! loops indices and polynom. degrees
+ REAL(dp) :: p_dp, j_dp
REAL(dp) :: kr, kz, kperp2
REAL(dp) :: taue_qe_etaB, taui_qi_etaB
REAL(dp) :: sqrtTaue_qe, sqrtTaui_qi, qe_sigmae_sqrtTaue, qi_sigmai_sqrtTaui
@@ -25,6 +25,9 @@ SUBROUTINE moments_eq_rhs
COMPLEX(dp) :: test_nan
REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
+ ! Measuring execution time
+ CALL cpu_time(t0_rhs)
+
!Precompute species dependant factors
IF( q_e .NE. 0._dp ) THEN
taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
@@ -48,44 +51,46 @@ SUBROUTINE moments_eq_rhs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Electrons moments RHS !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
- ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxe)
+ ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices
+ p_int= parray_e(ip) ! Hermite polynom. degree
+ p_dp = REAL(p_int,dp) ! REAL of the Hermite degree
! N_e^{p+1,j} coeff
- xNapp1j = sqrtTaue_qe * SQRT(ip_dp + 1)
+ xNapp1j = sqrtTaue_qe * SQRT(p_dp + 1)
! N_e^{p-1,j} coeff
- xNapm1j = sqrtTaue_qe * SQRT(ip_dp)
+ xNapm1j = sqrtTaue_qe * SQRT(p_dp)
! N_e^{p+2,j} coeff
- xNapp2j = taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ xNapp2j = taue_qe_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
! N_e^{p-2,j} coeff
- xNapm2j = taue_qe_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
+ xNapm2j = taue_qe_etaB * SQRT(p_dp * (p_dp - 1._dp))
factj = 1.0 ! Start of the recursive factorial
- jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
- ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxe)
+ jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
+ j_int= jarray_e(ij) ! Laguerre polynom. degree
+ j_dp = REAL(j_int,dp) ! REAL of degree
! N_e^{p,j+1} coeff
- xNapjp1 = -taue_qe_etaB * (ij_dp + 1._dp)
+ xNapjp1 = -taue_qe_etaB * (j_dp + 1._dp)
! N_e^{p,j-1} coeff
- xNapjm1 = -taue_qe_etaB * ij_dp
+ xNapjm1 = -taue_qe_etaB * j_dp
! N_e^{pj} coeff
- xNapj = taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+ xNapj = taue_qe_etaB * 2._dp*(p_dp + j_dp + 1._dp)
!! Collision operator pj terms
- xCapj = -nu_e*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
+ xCapj = -nu_e*(p_dp + 2._dp*j_dp) !DK Lenard-Bernstein basis
! Dougherty part
IF ( CO .EQ. -2) THEN
- IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20
+ IF ((p_int .EQ. 2) .AND. (j_int .EQ. 0)) THEN ! kronecker pj20
xCa20 = nu_e * 2._dp/3._dp
xCa01 = -SQRT2 * xCa20
xCa10 = 0._dp
- ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01
+ ELSEIF ((p_int .EQ. 0) .AND. (j_int .EQ. 1)) THEN ! kronecker pj01
xCa20 = -nu_e * SQRT2 * 2._dp/3._dp
xCa01 = -SQRT2 * xCa20
xCa10 = 0._dp
- ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN ! kronecker pj10
+ ELSEIF ((p_int .EQ. 1) .AND. (j_int .EQ. 0)) THEN ! kronecker pj10
xCa20 = 0._dp
xCa01 = 0._dp
xCa10 = nu_e
@@ -95,15 +100,15 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrostatic potential pj terms
- IF (ip .EQ. 1) THEN ! kronecker p0
- xphij = (eta_n + 2.*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp) )
- xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp)
- xphijm1 = -(eta_T - eta_B)* ij_dp
+ IF (p_int .EQ. 0) THEN ! kronecker p0
+ xphij = (eta_n + 2.*j_dp*eta_T - 2._dp*eta_B*(j_dp+1._dp) )
+ xphijp1 = -(eta_T - eta_B)*(j_dp+1._dp)
+ xphijm1 = -(eta_T - eta_B)* j_dp
xphijpar= 0._dp
- ELSE IF (ip .EQ. 2) THEN ! kronecker p1
+ ELSE IF (p_int .EQ. 1) THEN ! kronecker p1
xphijpar = qe_sigmae_sqrtTaue
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
- ELSE IF (ip .EQ. 3) THEN ! kronecker p2
+ ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp
ELSE
@@ -111,8 +116,8 @@ SUBROUTINE moments_eq_rhs
ENDIF
! Recursive factorial
- IF (ij_dp .GT. 0) THEN
- factj = factj * ij_dp
+ IF (j_dp .GT. 0) THEN
+ factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
@@ -128,27 +133,27 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_e^{p,j}
TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
- ! term propto N_e^{p+1,j}
- IF (ip+1 .LE. pmaxe+1) THEN ! OoB check
+ ! term propto N_e^{p+1,j} and kparallel
+ IF ( (ip+1 .LE. pmaxe+1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
- ! term propto N_e^{p-1,j}
- IF (ip-1 .GE. 1) THEN ! OoB check
+ ! term propto N_e^{p-1,j} and kparallel
+ IF ( (ip-1 .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
! term propto N_e^{p+2,j}
- IF (ip+2 .LE. pmaxe+1) THEN ! OoB check
- TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
+ IF (ip+(2-pskip) .LE. pmaxe+1) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_e(ip+(2-pskip),ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
! term propto N_e^{p-2,j}
- IF (ip-2 .GE. 1) THEN ! OoB check
- TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
+ IF (ip-(2-pskip) .GE. 1) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_e(ip-(2-pskip),ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
@@ -167,8 +172,8 @@ SUBROUTINE moments_eq_rhs
!! Collision
IF (CO .EQ. -2) THEN ! Dougherty Collision terms
- IF ( (pmaxe .GE. 2) ) THEN ! OoB check
- TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
+ IF ( (pmaxe .GE. 2-pskip) ) THEN ! OoB check
+ TColl20 = xCa20 * moments_e(3-pskip,1,ikr,ikz,updatetlevel)
ELSE
TColl20 = 0._dp
ENDIF
@@ -177,7 +182,7 @@ SUBROUTINE moments_eq_rhs
ELSE
TColl01 = 0._dp
ENDIF
- IF ( (pmaxe .GE. 1) ) THEN ! OoB check
+ IF ( (pmaxe .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB + odd number for Hermite degrees check
TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
ELSE
TColl10 = 0._dp
@@ -192,17 +197,25 @@ SUBROUTINE moments_eq_rhs
TColl = 0._dp ! Initialization
ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms
+ p2_int = parray_e(ip2)
jloopee: DO ij2 = 1,jmaxe+1
+ j2_int = jarray_e(ij2)
+
TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
- *( nu_e * CeipjT(bare(ip-1,ij-1), bare(ip2-1,ij2-1)) &
- +nu_ee * Ceepj (bare(ip-1,ij-1), bare(ip2-1,ij2-1)))
+ *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int)) &
+ +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int)))
+
ENDDO jloopee
ENDDO ploopee
ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms
+ p2_int = parray_i(ip2)
jloopei: DO ij2 = 1,jmaxi+1
+ j2_int = jarray_i(ij2)
+
TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
- *(nu_e * CeipjF(bare(ip-1,ij-1), bari(ip2-1,ij2-1)))
+ *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int)))
+
END DO jloopei
ENDDO ploopei
@@ -211,11 +224,11 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrical potential term
- IF ( (ip .LE. 3) ) THEN ! kronecker p0 p1 p2
- kernelj = be_2**(ij-1) * exp(-be_2)/factj
- kerneljp1 = kernelj * be_2 /(ij_dp + 1._dp)
+ IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2
+ kernelj = be_2**j_int * exp(-be_2)/factj
+ kerneljp1 = kernelj * be_2 /(j_dp + 1._dp)
IF ( be_2 .NE. 0 ) THEN
- kerneljm1 = kernelj * ij_dp / be_2
+ kerneljm1 = kernelj * j_dp / be_2
ELSE
kerneljm1 = 0._dp
ENDIF
@@ -243,47 +256,52 @@ SUBROUTINE moments_eq_rhs
END DO jloope
END DO ploope
+ !_____________________________________________________________________________!
+ !_____________________________________________________________________________!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Ions moments RHS !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ploopi : DO ip = ips_i, ipe_i ! This loop is from 1 to pmaxi+1
- ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxi)
+ !_____________________________________________________________________________!
+ ploopi : DO ip = ips_i, ipe_i ! Hermite loop
+ p_int= parray_i(ip) ! Hermite degree
+ p_dp = REAL(p_int,dp) ! REAL of Hermite degree
! N_i^{p+1,j} coeff
- xNapp1j = sqrtTaui_qi * SQRT(ip_dp + 1)
+ xNapp1j = sqrtTaui_qi * SQRT(p_dp + 1)
! N_i^{p-1,j} coeff
- xNapm1j = sqrtTaui_qi * SQRT(ip_dp)
+ xNapm1j = sqrtTaui_qi * SQRT(p_dp)
! x N_i^{p+2,j} coeff
- xNapp2j = taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ xNapp2j = taui_qi_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
! x N_i^{p-2,j} coeff
- xNapm2j = taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
+ xNapm2j = taui_qi_etaB * SQRT(p_dp * (p_dp - 1._dp))
factj = 1._dp ! Start of the recursive factorial
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
- ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxi)
+ j_int= jarray_i(ij) ! REALof Laguerre degree
+ j_dp = REAL(j_int,dp) ! REALof Laguerre degree
! x N_i^{p,j+1} coeff
- xNapjp1 = -taui_qi_etaB * (ij_dp + 1._dp)
+ xNapjp1 = -taui_qi_etaB * (j_dp + 1._dp)
! x N_i^{p,j-1} coeff
- xNapjm1 = -taui_qi_etaB * ij_dp
+ xNapjm1 = -taui_qi_etaB * j_dp
! x N_i^{pj} coeff
- xNapj = taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+ xNapj = taui_qi_etaB * 2._dp*(p_dp + j_dp + 1._dp)
!! Collision operator pj terms
- xCapj = -nu_i*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
+ xCapj = -nu_i*(p_dp + 2._dp*j_dp) !DK Lenard-Bernstein basis
! Dougherty part
IF ( CO .EQ. -2) THEN
- IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20
+ IF ((p_int .EQ. 2) .AND. (j_int .EQ. 0)) THEN ! kronecker pj20
xCa20 = nu_i * 2._dp/3._dp
xCa01 = -SQRT2 * xCa20
xCa10 = 0._dp
- ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01
+ ELSEIF ((p_int .EQ. 0) .AND. (j_int .EQ. 1)) THEN ! kronecker pj01
xCa20 = -nu_i * SQRT2 * 2._dp/3._dp
xCa01 = -SQRT2 * xCa20
xCa10 = 0._dp
- ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN
+ ELSEIF ((p_int .EQ. 1) .AND. (j_int .EQ. 0)) THEN ! kronecker pj10
xCa20 = 0._dp
xCa01 = 0._dp
xCa10 = nu_i
@@ -293,15 +311,15 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrostatic potential pj terms
- IF (ip .EQ. 1) THEN ! krokecker p0
- xphij = (eta_n + 2._dp*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp))
- xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp)
- xphijm1 = -(eta_T - eta_B)* ij_dp
+ IF (p_int .EQ. 0) THEN ! krokecker p0
+ xphij = (eta_n + 2._dp*j_dp*eta_T - 2._dp*eta_B*(j_dp+1._dp))
+ xphijp1 = -(eta_T - eta_B)*(j_dp+1._dp)
+ xphijm1 = -(eta_T - eta_B)* j_dp
xphijpar = 0._dp
- ELSE IF (ip .EQ. 2) THEN ! kronecker p1
+ ELSE IF (p_int .EQ. 1) THEN ! kronecker p1
xphijpar = qi_sigmai_sqrtTaui
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
- ELSE IF (ip .EQ. 3) THEN !krokecker p2
+ ELSE IF (p_int .EQ. 2) THEN !krokecker p2
xphij = (eta_T/SQRT2 - SQRT2*eta_B)
xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar = 0._dp
ELSE
@@ -309,8 +327,8 @@ SUBROUTINE moments_eq_rhs
ENDIF
! Recursive factorial
- IF (ij_dp .GT. 0) THEN
- factj = factj * ij_dp
+ IF (j_dp .GT. 0) THEN
+ factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
@@ -327,26 +345,26 @@ SUBROUTINE moments_eq_rhs
! term propto N_i^{p,j}
TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p+1,j}
- IF (ip+1 .LE. pmaxi+1) THEN ! OoB check
+ IF ( (ip+1 .LE. pmaxi+1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
! term propto N_i^{p-1,j}
- IF (ip-1 .GE. 1) THEN ! OoB check
+ IF ( (ip-1 .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
! term propto N_i^{p+2,j}
- IF (ip+2 .LE. pmaxi+1) THEN ! OoB check
- TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
+ IF (ip+(2-pskip) .LE. pmaxi+1) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_i(ip+(2-pskip),ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
! term propto N_i^{p-2,j}
- IF (ip-2 .GE. 1) THEN ! OoB check
- TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
+ IF (ip-(2-pskip) .GE. 1) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_i(ip-(2-pskip),ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
@@ -365,8 +383,8 @@ SUBROUTINE moments_eq_rhs
!! Collision
IF (CO .EQ. -2) THEN ! Dougherty Collision terms
- IF ( (pmaxi .GE. 2) ) THEN ! OoB check
- TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
+ IF ( (pmaxi .GE. 2-pskip) ) THEN ! OoB check
+ TColl20 = xCa20 * moments_i(3-pskip,1,ikr,ikz,updatetlevel)
ELSE
TColl20 = 0._dp
ENDIF
@@ -375,7 +393,7 @@ SUBROUTINE moments_eq_rhs
ELSE
TColl01 = 0._dp
ENDIF
- IF ( (pmaxi .GE. 1) ) THEN ! OoB check
+ IF ( (pmaxi .GE. 1) .AND. (.NOT. CANCEL_ODD_P) ) THEN ! OoB check
TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
ELSE
TColl10 = 0._dp
@@ -390,17 +408,25 @@ SUBROUTINE moments_eq_rhs
TColl = 0._dp ! Initialization
ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms
+ p2_int = parray_i(ip2)
jloopii: DO ij2 = 1,jmaxi+1
+ j2_int = jarray_i(ij2)
+
TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
- *( nu_ie * CiepjT(bari(ip-1,ij-1), bari(ip2-1,ij2-1)) &
- +nu_i * Ciipj (bari(ip-1,ij-1), bari(ip2-1,ij2-1)))
+ *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int)) &
+ +nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int)))
+
ENDDO jloopii
ENDDO ploopii
ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms
+ p2_int = parray_e(ip2)
jloopie: DO ij2 = 1,jmaxe+1
+ j2_int = jarray_e(ij2)
+
TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
- *(nu_ie * CiepjF(bari(ip-1,ij-1), bare(ip2-1,ij2-1)))
+ *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int)))
+
ENDDO jloopie
ENDDO ploopie
@@ -409,11 +435,11 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrical potential term
- IF ( (ip .LE. 3) ) THEN ! kronecker p0 p1 p2
- kernelj = bi_2**(ij-1) * exp(-bi_2)/factj
- kerneljp1 = kernelj * bi_2 /(ij_dp + 1._dp)
+ IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2
+ kernelj = bi_2**j_int * exp(-bi_2)/factj
+ kerneljp1 = kernelj * bi_2 /(j_dp + 1._dp)
IF ( bi_2 .NE. 0 ) THEN
- kerneljm1 = kernelj * ij_dp / bi_2
+ kerneljm1 = kernelj * j_dp / bi_2
ELSE
kerneljm1 = 0._dp
ENDIF
@@ -441,4 +467,8 @@ SUBROUTINE moments_eq_rhs
END DO jloopi
END DO ploopi
+ ! Execution time end
+ CALL cpu_time(t1_rhs)
+ tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
+
END SUBROUTINE moments_eq_rhs
diff --git a/src/poisson.F90 b/src/poisson.F90
index c0fdbcfbc7daf7b40ce2ddb4ba2cce2ed4c3c324..0c8db4f38e03e027a1ccdbf84f9e87fcc74e2ad2 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -1,6 +1,7 @@
SUBROUTINE poisson
! Solve poisson equation to get phi
+ USE basic, ONLY: t0_poisson, t1_poisson, tc_poisson
USE time_integration, ONLY: updatetlevel
USE array
USE fields
@@ -21,6 +22,9 @@ SUBROUTINE poisson
REAL(dp) :: gammaD
COMPLEX(dp) :: gammaD_phi
+ ! Execution time start
+ CALL cpu_time(t0_poisson)
+
!Precompute species dependant factors
sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
@@ -91,4 +95,8 @@ SUBROUTINE poisson
! Cancel origin singularity
phi(ikr_0,ikz_0) = 0
+! Execution time end
+CALL cpu_time(t1_poisson)
+tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
+
END SUBROUTINE poisson
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index 3734013cb995294b8c959906f813f4b23e22feb9..c80bf8779406baaa02d373f19e2d61368a1af2f2 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -2,14 +2,23 @@ SUBROUTINE ppinit
! Parallel environment
USE basic
-
+ USE model, ONLY : NON_LIN
use prec_const
IMPLICIT NONE
-
+
INTEGER :: version_prov=-1
- ! CALL mpi_init(ierr)
+ CALL MPI_INIT(ierr)
! CALL MPI_INIT_THREAD(MPI_THREAD_SINGLE,version_prov,ierr)
- CALL MPI_INIT_THREAD(MPI_THREAD_FUNNELED,version_prov,ierr)
+ ! CALL MPI_INIT_THREAD(MPI_THREAD_FUNNELED,version_prov,ierr)
+
+ CALL MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr)
+ CALL MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
+ CALL MPI_COMM_GROUP(MPI_COMM_WORLD,grp_world,ierr)
+ ! CALL MPI_COMM_CREATE_GROUP(MPI_COMM_WORLD,grp_world,0,comm_self,ierr)
+
+ ! IF (NON_LIN .AND. (num_procs .GT. 1)) THEN
+ ! CALL fftw_mpi_init
+ ! ENDIF
END SUBROUTINE ppinit
diff --git a/src/ppsetup.f90 b/src/ppsetup.f90
new file mode 100644
index 0000000000000000000000000000000000000000..cb6d06d32f92f2ac405735ad619f43a1cc7fd2fa
--- /dev/null
+++ b/src/ppsetup.f90
@@ -0,0 +1,437 @@
+SUBROUTINE PPSETUP
+ !
+ ! sets up the MPI topology and distribute the work betwen the nprocs
+ !
+ use prec_const
+ use basic
+ USE basic_mpi
+ use model
+ !
+ IMPLICIT NONE
+ !
+ INTEGER :: njobsie,njobsei,njobsaa,i
+ INTEGER :: iproc =0
+ INTEGER, DIMENSION(nprocs) :: procslabels
+ LOGICAL, DIMENSION(:),allocatable :: logical_mask
+ INTEGER, DIMENSION(:), allocatable :: buff_
+ !
+ ! loc. vars.
+ INTEGER :: lbaremax,lbarimax
+ INTEGER, DIMENSION(2) :: pjs,pje
+ INTEGER :: ii,grp_rank
+ !
+ INTEGER :: dims(2)
+ LOGICAL :: is_periodic(2)
+ INTEGER :: Ndim_energy
+ INTEGER :: Nprocs_lbar
+ INTEGER :: shift_
+ LOGICAL :: reorder
+ !
+ IF(nprocs > 1 ) THEN
+ !
+ ! Generate an 1D topology: the collisional rows are shared with procesors
+ !
+ lbaremax = numbe(Pmaxe,Jmaxe)
+ lbarimax = numbi(Pmaxi,Jmaxi)
+ !
+ ! Create MPI group E/I/SELF
+ Ngroups = 0
+ IF(eicolls) Ngroups = Ngroups +1
+ IF(iecolls) Ngroups = Ngroups +1
+ IF(eecolls .or. iicolls) Ngroups = Ngroups +1
+ !
+ ! Assign group ranks
+ !
+ iproc = 0
+ IF(eicolls) THEN
+ ALLOCATE(ranks_e(1:Cenproc))
+ DO ii=1,Cenproc
+ ranks_e(ii) = ii-1
+ ENDDO
+ iproc = Cenproc
+ !! if(me .eq. 0) write(*,*) ranks_e(:)
+ ENDIF
+ !
+ IF(iecolls) THEN
+ ALLOCATE(ranks_i(1:Cinproc))
+ DO ii=1,Cinproc
+ ranks_i(ii) = iproc + ii -1
+ ENDDO
+ iproc = iproc + Cinproc
+ !! if(me .eq. 0) write(*,*) ranks_i(:)
+ ENDIF
+ !
+ IF(eecolls .or. iicolls) THEN
+ ALLOCATE(ranks_self(1:Caanproc))
+ DO ii=1,Caanproc
+ ranks_self(ii) = iproc + ii -1
+ ENDDO
+ !! if(me .eq. 0) write(*,*) ranks_self(:)
+ ENDIF
+ !
+ ! CREATE Groups
+ CALL MPI_COMM_GROUP(MPI_COMM_WORLD,grp_world,ierr)
+ !
+ IF(eicolls) THEN
+ CALL MPI_GROUP_INCL(grp_world,Cenproc,ranks_e,grp_e,ierr)
+ ! Get my rank in MPI group
+ CALL MPI_GROUP_RANK(grp_e,me_e,ierr)
+ CALL MPI_COMM_CREATE_GROUP(MPI_COMM_WORLD,grp_e,0,comm_e,ierr)
+ IF(me_e .ne. MPI_UNDEFINED ) is_e = .true.
+ ENDIF
+ !
+ IF(iecolls) THEN
+ CALL MPI_GROUP_INCL(grp_world,Cinproc,ranks_i,grp_i,ierr)
+ ! Get my rank in MPI group
+ CALL MPI_GROUP_RANK(grp_i,me_i,ierr)
+ CALL MPI_COMM_CREATE_GROUP(MPI_COMM_WORLD,grp_i,0,comm_i,ierr)
+ IF(me_i.ne. MPI_UNDEFINED) is_i = .true.
+ ENDIF
+ !
+ IF(eecolls .or. iicolls) THEN
+ ! Create group
+ CALL MPI_GROUP_INCL(grp_world,Caanproc,ranks_self,grp_self,ierr)
+ ! Get my rank in MPI group
+ CALL MPI_GROUP_RANK(grp_self,me_self,ierr)
+ CALL MPI_COMM_CREATE_GROUP(MPI_COMM_WORLD,grp_self,0,comm_self,ierr)
+ IF(me_self .ne. MPI_UNDEFINED) is_self = .true.
+ !
+ IF( IFCOULOMB .and. IFGK) THEN
+ ! Create Cartesian Topology
+ ! Naanprocs_j_ii sould be a mutliple of Cannproc
+ Nprocs_rows_ii = floor(real(Caanproc)/Nprocs_j_ii)
+ dims = (/Nprocs_rows_ii, Nprocs_j_ii/) !
+ is_periodic = (/.FALSE.,.FALSE./)
+ reorder = .true.
+ !
+ CALL MPI_CART_CREATE(comm_self,2,dims,is_periodic,reorder,comm_cart_self,ierr)
+ CALL MPI_CART_COORDS(comm_cart_self,me_self,2,mycoords_ii,ierr)
+ !
+ !! debug
+ !! print*,Nprocs_rows_ii, dims
+ ENDIF
+ !
+ ENDIF
+ !
+ ! local 1D indices: start form 1
+ ! electrons
+ IF(eicolls .and. is_e ) THEN
+ IF(MPI_chksz .eq. 0 ) THEN
+ lbare_s_l = me_e*floor(real(lbaremax/Cenproc)) + 1
+ IF( me_e .eq. (Cenproc - 1) ) THEN
+ lbare_e_l = numbe(Pmaxe,Jmaxe)
+ ELSE
+ lbare_e_l = lbare_s_l + floor(real(lbaremax/Cenproc))-1
+ ENDIF
+ ELSE
+ ! fill up from the bottom with chunk size = MPI_chksz
+ shift_ = numbe(Pmaxe,Jmaxe) - MPI_chksz*(Cenproc -1)
+ IF( me_e .eq. 0 ) THEN
+ lbare_s_l = 1
+ lbare_e_l = shift_
+ ELSE
+ lbare_s_l = shift_ + (me_e -1)*MPI_chksz +1
+ lbare_e_l = lbare_s_l + MPI_chksz -1
+ ENDIF
+ ! Debug
+ !! write(*,*) me_e, lbare_s_l,lbare_e_l,numbe(Pmaxe,Jmaxe)
+ ENDIF
+ ENDIF
+
+ ! ions
+ IF(iecolls .and. is_i ) THEN
+ IF(MPI_chksz .eq. 0 ) THEN
+ lbari_s_l = me_i*floor(real(lbarimax/Cinproc)) + 1
+ IF( me_i .eq. (Cinproc - 1) ) THEN
+ lbari_e_l = numbi(Pmaxi,Jmaxi)
+ ELSE
+ lbari_e_l = lbari_s_l + floor(real(lbarimax/Cinproc))-1
+ ENDIF
+ ELSE
+ ! fill up from the bottom with chunk size = MPI_chksz
+ shift_ = numbi(Pmaxi,Jmaxi) - MPI_chksz*(Cinproc -1)
+ IF( me_i .eq. 0 ) THEN
+ lbari_s_l = 1
+ lbari_e_l = shift_
+ ELSE
+ lbari_s_l = shift_ + (me_i -1)*MPI_chksz +1
+ lbari_e_l = lbari_s_l + MPI_chksz -1
+ ENDIF
+ ! Debug
+ !! write(*,*) me_e, lbare_s_l,lbare_e_l,numbe(Pmaxe,Jmaxe)
+ ENDIF
+ ENDIF
+
+ ! self electron collisions
+ IF( eecolls .and. is_self ) THEN
+ IF(MPI_chksz .eq. 0) THEN
+ lbaree_s_l = me_self*floor(real(lbaremax/Caanproc)) + 1
+ IF( me_self .eq. (Caanproc - 1) ) THEN
+ lbaree_e_l = numbe(Pmaxe,Jmaxe)
+ ELSE
+ lbaree_e_l = lbaree_s_l + floor(real(lbaremax/Caanproc))-1
+ ENDIF
+ ELSE
+ ! fill up from the bottom with chunk size = MPI_chksz
+ shift_ = numbe(Pmaxe,Jmaxe) - MPI_chksz*(Caanproc -1)
+ IF( me_self .eq. 0 ) THEN
+ lbaree_s_l = 1
+ lbaree_e_l = shift_
+ ELSE
+ lbaree_s_l = shift_ + (me_self -1)*MPI_chksz +1
+ lbaree_e_l = lbaree_s_l + MPI_chksz -1
+ ENDIF
+ ! Debug
+ !! write(*,*) me_e, lbare_s_l,lbare_e_l,numbe(Pmaxe,Jmaxe)
+ ENDIF
+ ENDIF
+ !
+ ! self ion collisions
+ IF( iicolls .and. is_self ) THEN
+ IF(MPI_chksz .eq. 0 ) THEN
+ IF(IFCOULOMB .and. IFGK) THEN
+ ! MPI implementation for GK Coulomb (start at 1)
+ lbarii_s_l = mycoords_ii(1)*floor(real(lbarimax/Nprocs_rows_ii)) +1
+ IF(mycoords_ii(1) .eq. Nprocs_rows_ii-1) THEN
+ lbarii_e_l = numbi(Pmaxi,Jmaxi)
+ ELSE
+ lbarii_e_l = lbarii_s_l + floor(real(lbarimax/Nprocs_rows_ii)) -1
+ ENDIF
+ !
+ ! local energy component (start at 0)
+ jimaxx_s_l = mycoords_ii(2)*floor(real(JEmaxx/Nprocs_j_ii))
+ IF(mycoords_ii(2) .eq. Nprocs_j_ii-1) THEN
+ jimaxx_e_l = JEmaxx
+ ELSE
+ jimaxx_e_l = jimaxx_s_l + floor(real(JEmaxx/Nprocs_j_ii)) -1
+ ENDIF
+
+ !! Debug
+ !! print*,mycoords_ii(:), lbarii_s_l,lbarii_e_l,jimaxx_s_l,jimaxx_e_l
+
+ !
+ ! Local FLR indices: (not parallelized in the FLR dimensions)
+ nimaxxFLR_s_l = 0
+ nimaxxFLR_e_l = nimaxxFLR
+ npimaxxFLR_s_l = 0
+ npimaxxFLR_e_l = npimaxxFLR
+
+ ELSE
+ lbarii_s_l = me_self*floor(real(lbarimax/Caanproc)) + 1
+ IF( me_self .eq. (Caanproc - 1) ) THEN
+ lbarii_e_l = numbi(Pmaxi,Jmaxi)
+ ELSE
+ lbarii_e_l = lbarii_s_l + floor(real(lbarimax/Caanproc))-1
+ ENDIF
+ ! Local FLR indices
+ nimaxxFLR_s_l = 0
+ nimaxxFLR_e_l = nimaxxFLR
+ npimaxxFLR_s_l = 0
+ npimaxxFLR_e_l = npimaxxFLR
+ jimaxx_s_l = 0
+ jimaxx_e_l = JEmaxx
+ ENDIF
+ ELSE
+ ! fill up from the bottom with chunk size = MPI_chksz
+ ! MPI implementation for GK Coulomb (start at 1)
+ !
+ IF(IFCOULOMB .and. IFGK) THEN
+ shift_ = numbi(Pmaxi,Jmaxi) - MPI_chksz*(Nprocs_rows_ii -1)
+ IF( mycoords_ii(1) .eq. 0 ) THEN
+ lbarii_s_l = 1
+ lbarii_e_l = shift_
+ ELSE
+ lbarii_s_l = shift_ + (mycoords_ii(1) -1)*MPI_chksz +1
+ lbarii_e_l = lbarii_s_l + MPI_chksz -1
+ ENDIF
+ !
+ ! local energy component (start at 0)
+ jimaxx_s_l = mycoords_ii(2)*floor(real(JEmaxx/Nprocs_j_ii))
+ IF(mycoords_ii(2) .eq. Nprocs_j_ii-1) THEN
+ jimaxx_e_l = JEmaxx
+ ELSE
+ jimaxx_e_l = jimaxx_s_l + floor(real(JEmaxx/Nprocs_j_ii)) -1
+ ENDIF
+ ! debug
+ !! print*, lbarii_s_l,lbarii_e_l,jimaxx_s_l,jimaxx_e_l
+ !
+ ! Local FLR indices
+ nimaxxFLR_s_l = 0
+ nimaxxFLR_e_l = nimaxxFLR
+ npimaxxFLR_s_l = 0
+ npimaxxFLR_e_l = npimaxxFLR
+
+ ELSE
+ shift_ = numbi(Pmaxi,Jmaxi) - MPI_chksz*(Caanproc -1)
+ IF( me_self .eq. 0 ) THEN
+ lbarii_s_l = 1
+ lbarii_e_l = shift_
+ ELSE
+ lbarii_s_l = shift_ + (me_self -1)*MPI_chksz +1
+ lbarii_e_l = lbarii_s_l + MPI_chksz -1
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDIF
+ !
+ ! Get the local Hermite-Laguerre indices
+ !
+ ! ELECTRON ELECTRON
+ IF(eicolls .and. (is_e .or. me .eq. 0)) THEN
+ !
+ pjs = invnumbe(lbare_s_l)
+ pje = invnumbe(lbare_e_l)
+ !
+ Pmaxe_s_l = pjs(1)
+ Jmaxe_s_l = pjs(2)
+ Pmaxe_e_l = pje(1)
+ Jmaxe_e_l = pje(2)
+
+ CALL ALLOCATE_ARRAY(Je_s_l,Pmaxe_s_l,Pmaxe_e_l)
+ CALL ALLOCATE_ARRAY(Je_e_l,Pmaxe_s_l,Pmaxe_e_l)
+
+ ! local Laguere indices
+ Je_s_l(Pmaxe_s_l) = Jmaxe_s_l
+ Je_e_l(:) = Jmaxe
+ Je_e_l(Pmaxe_e_l) = Jmaxe_e_l
+ !
+ ENDIF
+
+ ! ION - ELECTRON
+ !
+ IF(iecolls .and. is_i ) THEN
+ pjs = invnumbi(lbari_s_l)
+ Pmaxi_s_l = pjs(1)
+ Jmaxi_s_l = pjs(2)
+ pje = invnumbi(lbari_e_l)
+ Pmaxi_e_l = pje(1)
+ Jmaxi_e_l = pje(2)
+
+ CALL ALLOCATE_ARRAY(Ji_s_l,Pmaxi_s_l,Pmaxi_e_l)
+ CALL ALLOCATE_ARRAY(Ji_e_l,Pmaxi_s_l,Pmaxi_e_l)
+
+ ! local Laguere indices
+ Ji_s_l(Pmaxi_s_l) = Jmaxi_s_l
+ Ji_e_l(:) = Jmaxi
+ Ji_e_l(Pmaxi_e_l) = Jmaxi_e_l
+ !
+ ENDIF
+
+ ! Electron-Electron
+ !
+ IF(eecolls .and. is_self) THEN
+
+ pjs = invnumbe(lbaree_s_l)
+ Pmaxee_s_l = pjs(1)
+ Jmaxee_s_l = pjs(2)
+
+ pje = invnumbe(lbaree_e_l)
+ Pmaxee_e_l = pje(1)
+ Jmaxee_e_l = pje(2)
+ !
+ CALL ALLOCATE_ARRAY(Jee_s_l,Pmaxee_s_l,Pmaxee_e_l)
+ CALL ALLOCATE_ARRAY(Jee_e_l,Pmaxee_s_l,Pmaxee_e_l)
+ !
+ ENDIF
+
+ !
+ ! ION - ION
+ IF(iicolls .and. is_self) THEN
+
+ pjs = invnumbi(lbarii_s_l)
+ Pmaxii_s_l = pjs(1)
+ Jmaxii_s_l = pjs(2)
+
+ pje = invnumbi(lbarii_e_l)
+ Pmaxii_e_l = pje(1)
+ Jmaxii_e_l = pje(2)
+ !
+ CALL ALLOCATE_ARRAY(Jii_s_l,Pmaxii_s_l,Pmaxii_e_l)
+ CALL ALLOCATE_ARRAY(Jii_e_l,Pmaxii_s_l,Pmaxii_e_l)
+
+ ! local Laguere indices
+ Jii_s_l(Pmaxii_s_l) = Jmaxii_s_l
+ Jii_e_l(:) = Jmaxi
+ Jii_e_l(Pmaxii_e_l) = Jmaxii_e_l
+ !
+ ENDIF
+ !
+ ELSE ! ... if NPROC = 1
+
+ ! set local indices to global
+ ! electrons
+ Pmaxe_s_l =0
+ Pmaxe_e_l = Pmaxe
+ Jmaxe_s_l =0
+ Jmaxe_e_l = Jmaxe
+
+ CALL ALLOCATE_ARRAY(Je_s_l,Pmaxe_s_l,Pmaxe_e_l)
+ CALL ALLOCATE_ARRAY(Je_e_l,Pmaxe_s_l,Pmaxe_e_l)
+ Je_s_l(:) = 0
+ Je_e_l(:) = Jmaxe
+
+ ! ion
+ Pmaxi_s_l =0
+ Pmaxi_e_l = Pmaxi
+ Jmaxi_s_l =0
+ Jmaxi_e_l = Jmaxi
+ !
+ CALL ALLOCATE_ARRAY(Ji_s_l,Pmaxi_s_l,Pmaxi_e_l)
+ CALL ALLOCATE_ARRAY(Ji_e_l,Pmaxi_s_l,Pmaxi_e_l)
+ Ji_s_l(:) = 0
+ Ji_e_l(:) = Jmaxi
+ !
+ ! self electrons
+ !
+ Pmaxee_s_l = 0
+ Pmaxee_e_l = Pmaxe
+ Jmaxee_s_l = 0
+ Jmaxee_e_l = Jmaxe
+ !
+ CALL ALLOCATE_ARRAY(Jee_s_l,Pmaxee_s_l,Pmaxee_e_l)
+ CALL ALLOCATE_ARRAY(Jee_e_l,Pmaxee_s_l,Pmaxee_e_l)
+ Jee_s_l(:) = 0
+ Jee_e_l(:) = Jmaxe
+ !
+ ! self ions
+ !
+ Pmaxii_s_l = 0
+ Pmaxii_e_l = Pmaxi
+ Jmaxii_s_l = 0
+ Jmaxii_e_l = Jmaxi
+ !
+ CALL ALLOCATE_ARRAY(Jii_s_l,Pmaxii_s_l,Pmaxii_e_l)
+ CALL ALLOCATE_ARRAY(Jii_e_l,Pmaxii_s_l,Pmaxii_e_l)
+ Jii_s_l(:) = 0
+ Jii_e_l(:) = Jmaxi
+ !
+ ! Set local GK indices to global
+ ! Engery component
+ jimaxx_s_l = 0
+ jimaxx_e_l = JEmaxx
+ !
+ ! FLR indices
+ nimaxxFLR_s_l = 0
+ nimaxxFLR_e_l = nimaxxFLR
+ !
+ npimaxxFLR_s_l = 0
+ npimaxxFLR_e_l = nimaxxFLR
+ !
+ ! Set all com to MPI_WORLD_COMM
+ comm_e = MPI_COMM_WORLD
+ comm_i = MPI_COMM_WORLD
+ comm_self = MPI_COMM_WORLD
+ comm_cart_self = MPI_COMM_WORLD
+
+ !
+ me_e = 0
+ me_i= 0
+ me_self =0
+ !
+ is_e = .true.
+ is_i = .true.
+ is_self = .true.
+ !
+ ENDIF
+ !
+ END SUBROUTINE PPSETUP
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 776f900908fefdfddafb490f74fe992d7d4400f9..44bfcf8ecd377b0fd47aeb7541604fbc47e187bd 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -16,22 +16,17 @@ SUBROUTINE readinputs
! Load grid data from input file
CALL grid_readinputs
- WRITE(*,*) 'check'
! Load diagnostic options from input file
CALL output_par_readinputs
- WRITE(*,*) 'check'
! Load model parameters from input file
CALL model_readinputs
- WRITE(*,*) 'check'
! Load initial condition parameters from input file
CALL initial_readinputs
- WRITE(*,*) 'check'
! Load parameters for time integration from input file
CALL time_integration_readinputs
- WRITE(*,*) 'check'
END SUBROUTINE readinputs
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 7a17f8f70cc2262dddc2de0bb9089c94607edeb2..e6f1a84e19462bbd3c431f5511f11cbe78806a08 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='4817fe5-dirty')
+parameter (VERSION='8775464-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Mon Oct 19 15:30:35 CEST 2020')
+parameter (EXECDATE='Mon Nov 9 13:42:43 CET 2020')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 7a17f8f70cc2262dddc2de0bb9089c94607edeb2..e6f1a84e19462bbd3c431f5511f11cbe78806a08 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='4817fe5-dirty')
+parameter (VERSION='8775464-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Mon Oct 19 15:30:35 CEST 2020')
+parameter (EXECDATE='Mon Nov 9 13:42:43 CET 2020')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index f23c5ff4979fb368086185581c957b2e56d92e11..4826a63f5c2660394fb4abcb0089fd12fda3f24c 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -19,7 +19,11 @@ SUBROUTINE stepon
CALL moments_eq_rhs
! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
CALL advance_time_level
+
! Update the moments with the hierarchy RHS (step by step)
+
+ ! Execution time start
+ CALL cpu_time(t0_adv_field)
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
CALL advance_field(moments_e(ip,ij,:,:,:),moments_rhs_e(ip,ij,:,:,:))
@@ -30,6 +34,10 @@ SUBROUTINE stepon
CALL advance_field(moments_i(ip,ij,:,:,:),moments_rhs_i(ip,ij,:,:,:))
ENDDO
ENDDO
+ ! Execution time end
+ CALL cpu_time(t1_adv_field)
+ tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field)
+
! Update electrostatic potential
CALL poisson
! Update nonlinear term
@@ -37,15 +45,17 @@ SUBROUTINE stepon
CALL compute_Sapj
ENDIF
!(The two routines above are called in inital for t=0)
-
- CALL enforce_symetry() ! Enforcing symmetry on kr = 0
-
CALL checkfield_all()
END DO
CONTAINS
SUBROUTINE checkfield_all ! Check all the fields for inf or nan
+ ! Execution time start
+ CALL cpu_time(t0_checkfield)
+
+ CALL enforce_symetry() ! Enforcing symmetry on kr = 0
+
IF(.NOT.nlend) THEN
nlend=nlend .or. checkfield(phi,' phi')
DO ip=ips_e,ipe_e
@@ -60,10 +70,12 @@ SUBROUTINE stepon
ENDDO
ENDDO
ENDIF
+ ! Execution time end
+ CALL cpu_time(t1_checkfield)
+ tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield)
END SUBROUTINE checkfield_all
SUBROUTINE enforce_symetry
-
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO ikz=2,Nkz/2 !symmetry at kr = 0
diff --git a/src/tesend.F90 b/src/tesend.F90
index 121637b681bce53896cbcfff474333b0102a3c17..5d77e4444f0bd479cb4606f351d1480f8644aa6e 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -9,27 +9,27 @@ SUBROUTINE tesend
!________________________________________________________________________________
! 1. Some processors had set nlend
IF( nlend ) THEN
- WRITE(*,'(/a)') 'rhs are NaN/Inf'
- WRITE(*,*) 'Run terminated at cstep=',cstep
+ WRITE(*,'(/a)') 'rhs are NaN/Inf'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Run terminated at cstep=',cstep
RETURN
END IF
!________________________________________________________________________________
! 2. Test on NRUN
- nlend = step .GT. nrun
- IF ( nlend ) THEN
+ nlend = step .GT. nrun
+ IF ( nlend ) THEN
WRITE(*,'(/a)') 'NRUN steps done'
RETURN
END IF
-
+
!________________________________________________________________________________
! 3. Test on TMAX
nlend = time .GT. tmax
- IF ( nlend ) THEN
- WRITE(*,'(/a)') 'TMAX reached'
+ IF ( nlend ) THEN
+ IF (my_id .EQ. 0) WRITE(*,'(/a)') 'TMAX reached'
RETURN
END IF
- !
+ !
!
END SUBROUTINE tesend
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 03e7098bb88cb7752a14d0c01e58abfd5ff24685..287e304e23fdcc6701d29ee3458bf6c5563d3bc3 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -68,10 +68,10 @@ CONTAINS
CASE ('DOPRI5_ADAPT')
CALL DOPRI5_ADAPT
CASE DEFAULT
- WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
+ IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
END SELECT
- WRITE(*,*) " Time integration with ", numerical_scheme
+ IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme
END SUBROUTINE set_numerical_scheme
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 6c739bb4c5052b3fdae5a3523f85c80217f0bf98..71fab610aa2f413372230d25a89d04f554b4ccb9 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,5 +1,6 @@
%% Load results
% JOBNUM = 0; load_results;
+% JOBNUM = 1; load_results;
compile_results
%% Retrieving max polynomial degree and sampling info
@@ -114,34 +115,38 @@ end
gkr0kz_Ni00 = real(g_(ikr0KH,:));
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+default_plots_options
disp('Plots')
FMT = '.fig';
+
+if 0
%% Time evolutions
fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
+set(gcf, 'Position', [100, 100, 900, 800])
subplot(221);
for ip = 1:Npe
for ij = 1:Nje
plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
+ plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$'];
semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on;
end
end
- grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
+ grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
subplot(222)
for ip = 1:Npi
for ij = 1:Nji
plt = @(x) squeeze(x(ip,ij,:));
- plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
+ plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$'];
semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on;
end
end
- grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$');
+ grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$');
subplot(223)
plot(Ts2D,Flux_ri,'-','DisplayName','$\Gamma_{ri}$'); hold on;
plot(Ts2D,Flux_zi,'-','DisplayName','$\Gamma_{zi}$'); hold on;
plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$')
plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$')
- grid on; xlabel('$t$'); ylabel('$\Gamma$'); %legend('show');
+ grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma$'); %legend('show');
if strcmp(OUTPUTS.write_non_lin,'.true.')
subplot(224)
for ip = 1:Npi
@@ -151,72 +156,122 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on;
end
end
- grid on; xlabel('$t$'); ylabel('$S$'); %legend('show');
+ grid on; xlabel('$t c_s/R$'); ylabel('$\sum_{k_r,k_z}|S_i^{pj}|$'); %legend('show');
else
%% Growth rate
subplot(224)
[~,ikr0KH] = min(abs(kr-KR0KH));
- plot(kz(1:ikr0KH)/kr(ikr0KH),gkr0kz_Ni00(1:ikr0KH)/(KR0KH*A0KH),...
+ plot(kz(1:Nz/2),gkr0kz_Ni00(1:Nz/2),...
'DisplayName',['J = ',num2str(JMAXI)]);
- xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$');
+ xlabel('$k_z$'); ylabel('$\gamma_{Ni}$');
xlim([0. 1.0]); %ylim([0. 0.04]);
end
+suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]);
save_figure
+end
%%
-if 1
-%% Show frame in real space
-tf = 1000; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
-fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)];
- subplot(221); plt = @(x) (((x)));
- pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_e^{00}|$');
- subplot(222); plt = @(x) ((x));
- pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat n_i^{00}|$');
- subplot(223); plt = @(x) ((x));
- pclr = pcolor((RR),(ZZ),plt(phi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts2D(it))); legend('$|\hat\phi|$');
-if strcmp(OUTPUTS.write_non_lin,'.true.')
- subplot(224); plt = @(x) fftshift((abs(x)),2);
- pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none'); colorbar;
- xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it5D))); legend('$|S_i^{00}|$');
+if 0
+%% Photomaton : real space
+tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
+fig = figure; FIGNAME = ['photo_real',sprintf('_t=%.0f',Ts2D(it))]; set(gcf, 'Position', [100, 100, 1500, 500])
+ subplot(131); plt = @(x) (((x)));
+ pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$'); legend('$n_i$');
+
+ subplot(132); plt = @(x) ((x));
+ DATA = plt(ni00(:,:,it))-plt(ne00(:,:,it));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$'); legend('$n_i-n_e$'); set(gca,'ytick',[]);
+
+
+ subplot(133); plt = @(x) ((x));
+ DATA = plt(phi(:,:,it));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$'); set(gca,'ytick',[]); legend('$\phi$');
+
+% if strcmp(OUTPUTS.write_non_lin,'.true.')
+% subplot(133); plt = @(x) fftshift((abs(x)),2);
+% pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+% xlabel('$r/\rho_s$'); legend('$|S_i^{00}|$'); set(gca,'ytick',[])
+% end
+suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB), sprintf(', $t c_s/R=%.0f$',Ts2D(it))]);
+save_figure
end
+
+if 1
+%% Photomaton : real space
+% FIELD = ni00; FNAME = 'ni';
+FIELD = phi; FNAME = 'phi';
+tf = 60; [~,it1] = min(abs(Ts2D-tf));
+tf = 65; [~,it2] = min(abs(Ts2D-tf));
+tf = 200; [~,it3] = min(abs(Ts2D-tf));
+tf = 400; [~,it4] = min(abs(Ts2D-tf));
+fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position', [100, 100, 1500, 400])
+ subplot(141); plt = @(x) ((x));
+ DATA = plt(FIELD(:,:,it1));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
+ subplot(142); plt = @(x) ((x));
+ DATA = plt(FIELD(:,:,it2));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
+ subplot(143); plt = @(x) ((x));
+ DATA = plt(FIELD(:,:,it3));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
+ subplot(144); plt = @(x) ((x));
+ DATA = plt(FIELD(:,:,it4));
+ pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+ xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
+ title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
+% suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+% ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']);
save_figure
end
%%
-t0 = 0;
+t0 = 50;
skip_ = 1;
DELAY = 0.01*skip_;
FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% Density ion
-GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D;
-FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
+FIELDNAME = '$n_i$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$';
create_gif
end
if 0
%% Density electron
-GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D;
-FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$';
+FIELDNAME = '$n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$';
+create_gif
+end
+if 0
+%% Density ion - electron
+GIFNAME = ['ni-ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
+FIELD = real(ni00+ne00); X = RR; Y = ZZ; T = Ts2D;
+FIELDNAME = '$n_i-n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$';
create_gif
end
if 0
%% Phi
GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D;
-FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$';
+FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
%% Density ion frequency
GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
-FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
+FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
create_gif
end
if 0
@@ -224,63 +279,42 @@ if 0
GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.';
X = (kz); T = Ts2D; YMIN = -.1; YMAX = 1.1; XMIN = min(kz); XMAX = max(kz);
-FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r$';
+FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r\rho_s$';
create_gif_1D
-% %% PJ ion moment frequency
-% p_ = 0+1; j_ = 3+1;
-% GIFNAME = ['Ni',num2str(p_),num2str(j_),sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-% FIELD =ifftshift(abs(squeeze(Nipj(p_,j_,:,:,:))),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
-% FIELDNAME = ['$N_i^{',num2str(p_-1),num2str(j_-1),'}$']; XNAME = '$k_r$'; YNAME = '$k_z$';
-% create_gif
-% %% non linear term frequ. space
-% GIFNAME = ['Si00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-% FIELD = fftshift((abs(Si00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
-% FIELDNAME = '$S_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
-% create_gif
-% %% Electron mode growth
-% GIFNAME = ['norm_Nepj',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-% FIELD = Ne_norm; X = PE; Y = JE; T = Ts;
-% FIELDNAME = '$\max_k{\gamma_e}$'; XNAME = '$p$'; YNAME = '$j$';
-% create_gif
end
%%
-if 0
-%% Asymmetry error on kr = 0
-up = 2:Nkz/2;
-down = Nkz:-1:Nkz/2+2;
-eps_asym_r = Ts5D;
-eps_asym_i = Ts5D;
-for it = 1:numel(Ts5D)
- eps_asym_r(it) = max(squeeze(abs(real(Ni00(1,up,it))-real(Ni00(1,down,it)))));
- eps_asym_i(it) = max(abs(imag(Ni00(1,up,it))+imag(Ni00(1,down,it))));
+if 1
+%% Space time diagramm (fig 11 Ivanov 2020)
+t0 = 1; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position', [100, 100, 1200, 400])
+ subplot(211)
+ plot(Ts2D(it0:it1),Flux_ri(it0:it1));
+ ylabel('$\Gamma_r$'); grid on
+% title(['$\eta=',num2str(ETAB),'\quad',...
+% '\nu_{',CONAME,'}=',num2str(NU),'$'])
+% legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
+ ylim([-0.01,1.1*max(Flux_ri(it0:it1))]);
+ subplot(212)
+ [TY,TX] = meshgrid(r,Ts2D(it0:it1));
+ pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,it0:it1),2))'); set(pclr, 'edgecolor','none'); %colorbar;
+ xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
+ legend('$\langle\partial_r \phi\rangle_z$')
+save_figure
end
-figure
- it = 2;
- subplot(311)
- plot(kz(up),real(Ni00(1,up,it))); hold on
- plot(kz(up),real(Ni00(1,down,it)),'--')
- ylabel('Real'); title(['$N_i^{00}(t=',num2str(Ts5D(it)),')$'])
- subplot(312)
- plot(kz(up),imag(Ni00(1,up,it))); hold on
- plot(kz(up),-imag(Ni00(1,down,it)),'--')
- ylabel('Imag'); xlabel('$kz,kr=0$')
- subplot(313)
- plot(kz(2:Nkz/2),real(Ni00(1,up,it))-real(Ni00(1,down,it))); hold on
- plot(kz(2:Nkz/2),imag(Ni00(1,up,it))+imag(Ni00(1,down,it)),'--')
- legend('Re','Im');
- ylabel('err'); xlabel('$kz,kr=0$')
-%%
-figure
- semilogy(Ts5D,eps_asym_r); hold on
- semilogy(Ts5D,eps_asym_i,'--');
- title('$\max_{k_z}|N_i^{00}(0,k_z)-N_i^{00}(0,-k_z)|$');
- legend('Re','Im'); grid on; xlabel('$t$')
+if 0
+%% Flow
+skip = 10;
+plt = @(x) x(1:skip:end,1:skip:end);
+tf = 120; [~,it] = min(abs(Ts2D-tf));
+fig = figure; FIGNAME = 'space_time_drphi';
+ quiver(plt(RR),plt(ZZ),plt(dzphi(:,:,it)),plt(-drphi(:,:,it)),0);
+ xlabel('$r$'), ylabel('$z$')
+ title(sprintf('$v_E$, $t=%.0f c_s/R$',Ts2D(it)));
end
-
if 0
%% Growth rate
fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
diff --git a/wk/analysis_stat_2D.m b/wk/analysis_stat_2D.m
index 4add74bfcbc3da1d0f4e52548fa24192e98c3579..07d4dcf9fc8b5e579ec0cacf6e88f8196f155ac1 100644
--- a/wk/analysis_stat_2D.m
+++ b/wk/analysis_stat_2D.m
@@ -1,29 +1,40 @@
if 0
%% ni ne phase space for different time windows
fig = figure;
- t0 = 20; t1 = 20; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
- plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
- plot(plt(ni00),plt(ne00),'.');
+ t0 = 80; t1 = 100; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+ plt = @(x) abs(squeeze(reshape(x(2,:,it0:it1),[],1)));
+ plot(plt(Ni00),plt(Ne00),'.');
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_i$'); ylabel('$n_e$');
end
if 0
%% histograms
+t0 = 200; t1 = 300; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
fig = figure;
- t0 = 100; t1 = 100; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
- plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
+subplot(221)
+ hist(plt(ne00),100); hold on
+ title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
+ xlabel('$n_e$');
+subplot(222)
hist(plt(ni00),100); hold on
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_i$');
+subplot(223)
+ hist(plt(phi),100); hold on
+ title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
+ xlabel('$\phi$');
+FIGNAME = ['histograms_ne_ni_phi'];
+save_figure
end
if 0
%% ni ne phi 3D phase space
fig = figure;
- t0 = 20; t1 = 20; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
- plt = @(x) squeeze(reshape(x(:,:,it0:it1),[],1));
- plot3(plt(ni00),plt(ne00),plt(phi),'.','MarkerSize',0.5);
+ t0 = 50; [~,it0] = min(abs(t0-Ts2D));
+ plt = @(x) squeeze(reshape(x(:,:,it0),[],1));
+ plot3(plt(ni00),plt(ne00),plt(phi),'.','MarkerSize',1.9);
title(['$',num2str(Ts2D(it0)),'\leq t \leq',num2str(Ts2D(it1)),'$']);
xlabel('$n_i$'); ylabel('$n_e$'); zlabel('$\phi$'); grid on;
end
diff --git a/wk/flux_results.m b/wk/flux_results.m
index 7ac27ee1f07440b2f27cae32bc9b8756f3494626..1cefbb8a6ad67388cc1cceb283e656008485ea45 100644
--- a/wk/flux_results.m
+++ b/wk/flux_results.m
@@ -1,63 +1,92 @@
default_plots_options
if 0
%% Compute time average and std of the mean flow
-t0 = 160; t1 = 200; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
+t0 = 180; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
range = it0:it1;
avg = mean(Flux_ri(range))
stdev = std(Flux_ri(range))^(.5)
figure
hist(Flux_ri(range),20); xlabel('$\Gamma$')
end
-%% Handwritten results
+%% Handwritten results for nu = 0.01
% High definition results (256x128)
-Results_256x128.Gamma = [0.62, 0.60, 2.7, 0.08];
-Results_256x128.N = [ 256, 256, 256, 256];
-Results_256x128.L = [ 50, 66, 66, 50];
-Results_256x128.P = [ 2, 2, 2, 3];
-Results_256x128.J = [ 1, 1, 1, 2];
-Results_256x128.NU = [0.01, 0.01, 0.01, 0.01];
-Results_256x128.etaB = [ 0.5, 0.5, 0.4, 0.5];
-Results_256x128.mrkr = [ 'v', 'v', 'v', 'v'];
-Results_256x128.clr = [ 'b', 'b', 'b', 'g'];
+Results_256x128.Gamma = [0.02, 0.03, 0.20, 0.037, 2.7, 2.25, 4, 5e-4, 2e-3, 0.03];
+Results_256x128.L = [ 66, 66, 66, 50, 66, 66, 66, 66, 66, 66];
+Results_256x128.P = [ 2, 3, 4, 5, 2, 3, 4, 2, 3, 4];
+Results_256x128.J = [ 1, 2, 2, 3, 1, 2, 2, 1, 2, 2];
+Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6];
+Results_256x128.mrkr = [ 'v', '>', '^', 'o', 'v', '>', '^', 'v', '>', '^'];
+Results_256x128.clr = [ 'k', 'k', 'k', 'r', 'r', 'r', 'r', 'k', 'k', 'k'];
% Low definition results (128x64)
-Results_128x64.Gamma = [0.29, 0.05, 7e-4, 0.31, 3.7, 2e-3];
-Results_128x64.N = [128, 128, 128, 128 128, 128];
-Results_128x64.L = [ 25, 25, 25, 33 33, 33];
-Results_128x64.P = [ 2, 2, 2, 2, 2, 2];
-Results_128x64.J = [ 1, 1, 1, 1, 1, 1];
-Results_128x64.NU = [0.01, 0.1, 0.01, 0.01, 0.01, 0.01];
-Results_128x64.etaB = [ 0.5, 0.5, 0.67, 0.5 0.4, 0.6];
-Results_128x64.mrkr = [ 's', 's', 's', 's', 's', 's'];
-Results_128x64.clr = [ 'b', 'b', 'b', 'r', 'r', 'r'];
-
-Ricci_Rogers.Gamma = [ 1 1e-2];
-Ricci_Rogers.etaB = [0.5 1.0];
+% Results_128x64.Gamma = [0.29, 0.05, 7e-4, 0.31, 3.7, 2e-3];
+% Results_128x64.L = [ 25, 25, 25, 33 33, 33];
+% Results_128x64.P = [ 2, 2, 2, 2, 2, 2];
+% Results_128x64.J = [ 1, 1, 1, 1, 1, 1];
+% Results_128x64.NU = [0.01, 0.1, 0.01, 0.01, 0.01, 0.01];
+% Results_128x64.etaB = [ 0.5, 0.5, 0.67, 0.5 0.4, 0.6];
+% Results_128x64.mrkr = [ 's', 's', 's', 's', 's', 's'];
+% Results_128x64.clr = [ 'b', 'b', 'b', 'r', 'r', 'r'];
+% Ricci_Rogers.Gamma = [2.5 1 1e-2];
+% Ricci_Rogers.etaB = [0.4 0.5 1.0];
+Ricci_Rogers.Gamma = [10 1e-1];
+Ricci_Rogers.etaB = [0.5 1.0];
if 1
% Fig 3 of Ricci Rogers 2006
fig = figure;
-
+semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
+hold on;
res = Results_256x128;
for i = 1:numel(res.Gamma)
- if res.NU(i) == 0.01
semilogy(res.etaB(i),res.Gamma(i),...
res.mrkr(i),'DisplayName','256x128', 'color', res.clr(i));
- end
hold on;
end
-res = Results_128x64;
+% res = Results_128x64;
+% for i = 1:numel(res.Gamma)
+% if res.NU(i) == 0.01
+% semilogy(res.etaB(i),res.Gamma(i),...
+% res.mrkr(i),'DisplayName','128x64', 'color', res.clr(i));
+% end
+% hold on;
+% end
+ xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
+end
+grid on; title('$\nu = 0.01$')
+legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$','$P=5$, $J=3$')
+FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png'];
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
+
+%% Handwritten results for nu = 0.1
+Results_256x128.Gamma = [0.026,0.026, 1e-2, 1, 1, 1, 2e-2, 1, 0.15, 3e-3];
+Results_256x128.P = [ 2, 3, 4, 2, 3, 4, 2, 3, 4, 4];
+Results_256x128.J = [ 1, 2, 2, 1, 2, 2, 1, 2, 2, 2];
+Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.7];
+Results_256x128.mrkr = [ 'v', '>', '^', 'v', '>', '^', 'v', '>', '^', '^'];
+Results_256x128.clr = [ 'k', 'k', 'k', 'b', 'b', 'b', 'r', 'b', 'r', 'r'];
+
+% Ricci_Rogers.Gamma = [2 1e-1];
+% Ricci_Rogers.etaB = [0.5 1.0];
+Ricci_Rogers.Gamma = [10 1e-1];
+Ricci_Rogers.etaB = [0.5 1.0];
+
+if 1
+% Fig 3 of Ricci Rogers 2006
+fig = figure;
+semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
+hold on;
+res = Results_256x128;
for i = 1:numel(res.Gamma)
- if res.NU(i) == 0.01
semilogy(res.etaB(i),res.Gamma(i),...
- res.mrkr(i),'DisplayName','128x64', 'color', res.clr(i));
- end
+ res.mrkr(i),'DisplayName','256x128', 'color', res.clr(i));
hold on;
end
+
xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
end
-plot(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--ok');
-grid on; title('$\nu = 0.01$')
-
-FIGNAME = [SIMDIR,'flux_study.png'];
+grid on; title('$\nu = 0.1$')
+legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$')
+FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png'];
saveas(fig,FIGNAME);
disp(['Figure saved @ : ',FIGNAME])
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index f684ea3ed002a4a38469b7a7a93836c8a4e43276..d823444bba59b837ec95babcacf2150103d978cf 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,41 +1,42 @@
&BASIC
nrun = 100000000
- dt = 0.01
- tmax = 300
- RESTART = .false.
+ dt = 0.004
+ tmax = 400
+ RESTART = .true.
/
&GRID
- pmaxe =2
- jmaxe = 1
- pmaxi = 2
- jmaxi = 1
+ pmaxe =5
+ jmaxe = 3
+ pmaxi = 5
+ jmaxi = 3
Nr = 256
Lr = 66
Nz = 256
Lz = 66
kpar = 0
+ CANCEL_ODD_P = .false.
/
&OUTPUT_PAR
- nsave_0d = 0
+ nsave_0d = 250
nsave_1d = 0
- nsave_2d = 50
- nsave_5d = 1000
+ nsave_2d = 250
+ nsave_5d = 2500
write_Na00 = .true.
write_moments = .true.
write_phi = .true.
write_non_lin = .true.
write_doubleprecision = .true.
- resfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/outputs'
- rstfile0 = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint'
+ resfile0 = '../results/ZP/256x128_L_66_Pe_5_Je_3_Pi_5_Ji_3_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/outputs'
+ rstfile0 = '../results/ZP/256x128_L_66_Pe_5_Je_3_Pi_5_Ji_3_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/checkpoint'
job2load = 0
/
&MODEL_PAR
! Collisionality
- CO = -2
+ CO = -1
DK = .false.
NON_LIN = .true.
mu = 0.0005
- nu = 0.01
+ nu = 0.1
tau_e = 1
tau_i = 1
sigma_e = 0.023338
@@ -44,7 +45,7 @@
q_i = 1
eta_n = 1
eta_T = 0
- eta_B = 0.4
+ eta_B = 0.5
lambdaD = 0
kr0KH = 0
A0KH = 0
@@ -54,9 +55,9 @@
initback_moments =0
initnoise_moments =1e-05
iseed =42
- selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5'
- eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
- iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+ selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10.h5'
+ eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+ iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10_tau_1.0000_mu_0.0233.h5'
/
&TIME_INTEGRATION_PAR
numerical_scheme='RK4'
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 4338d4894a11fda0f28b1a3942c5a31d7a12f036..0e08540a042cfc7192fd106326a2b93acb186b32 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -7,32 +7,33 @@ default_plots_options
%% PHYSICAL PARAMETERS
NU = 1e-2; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-ETAB = 1.0;
+ETAB = 0.8;
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
MU = 1e-4; % Hyper diffusivity coefficient
LAMBDAD = 0.0;
NOISE0 = 1.0e-5;
%% GRID PARAMETERS
-N = 50; % Frequency gridpoints (Nkr = N/2)
-L = 30; % Size of the squared frequency domain
-PMAXE = 20;
-JMAXE = 10;
-PMAXI = 20;
-JMAXI = 10;
+N = 40; % Frequency gridpoints (Nkr = N/2)
+L = 60; % Size of the squared frequency domain
+PMAXE = 12;
+JMAXE = 6;
+PMAXI = 12;
+JMAXI = 6;
KREQ0 = 1; % put kr = 0
-KPAR = 0.0; % Parellel wave vector component
+CANCEL_ODD_P = 1;% Cancels the odd polynomials degree
%% TIME PARAMETERS
-TMAX = 150; % Maximal time unit
-DT = 5e-3; % Time step
+TMAX = 300; % Maximal time unit
+DT = 3e-2; % Time step
+SPS0D = 0.5; % Sampling per time unit for 2D arrays
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 0; % Sampling per time unit for 5D arrays
+SPS5D = 0.1; % Sampling per time unit for 5D arrays
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 00;
%% OPTIONS
-SIMID = 'linear_study'; % Name of the simulation
+SIMID = 'test_linear_study'; % Name of the simulation
NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
-CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unused
@@ -40,6 +41,7 @@ KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
NO_E = 0; % Remove electrons dynamic
% DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
JOBNUM = 00;
+KPAR = 0.0 * (1-CANCEL_ODD_P); % Parellel wave vector component
setup
@@ -47,14 +49,18 @@ setup
%% PARAMETER SCANS
if 0
%% Parameter scan over PJ
+NU = 1e-2; % Collision frequency
ETAB = 1.0;
-TMAX = 150;
-PA = [2,3,4,5,6,8,10,12,14,16,20];
-JA = [1,1,2,2,3,4,5,6,7,8,10];
-% PA = [40];
-% JA = [20];
+ETAN = 1.0;
+TMAX = 400;
+PA = [2,6,10,14,20,40,80];
+JA = [1,3,5,7,10,20,20];
+% PA = [4];
+% JA = [2];
Nparam = numel(PA);
param_name = 'PJ';
+CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+DT = 5e-3; % Time step
gamma_Ni = zeros(Nparam,N);
@@ -64,16 +70,16 @@ for i = 1:Nparam
JMAXE = JA(i); JMAXI = JA(i);
setup
% Run linear simulation
- run
+ system('./../bin/helaz');
% Load and process results
load_results
- tend = Ts2D(end); tstart = 0.6*tend;
+ tend = Ts2D(end); tstart = 0.9*tend;
for ikz = 1:N
gamma_Ni(i,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(1,ikz,:))),tstart,tend);
end
gamma_Ni(i,:) = real(gamma_Ni(i,:) .* (gamma_Ni(i,:)>=0.0));
% Clean output
- system(['rm -r ',BASIC.RESDIR])
+% system(['rm -r ',BASIC.RESDIR])
end
if 1
@@ -93,16 +99,19 @@ grid on; xlabel('$k_z$'); ylabel('$\gamma(N_i^{00})$'); xlim([0.0,max(kz)]);
title(['$\eta_B=',num2str(ETAB),'$'])
legend('show')
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
+saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
end
end
if 0
%% Parameter scan over eta_B
PMAXE = 20; PMAXI = 20;
JMAXE = 10; JMAXI = 10;
-TMAX = 50;
+TMAX = 200;
eta_B = [0.1, 0.33, 0.5, 0.67, 1.0];
Nparam = numel(eta_B);
param_name = 'etaB';
+CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+NU = 1e-2; % Collision frequency
gamma_Ni = zeros(Nparam,N);
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index cdb5180ee8ab9499328ae9d6ca2f548522e1edcb..5025785f199d7622f3d96476429bd0bc7ad4ce52 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -4,30 +4,32 @@ addpath(genpath('../matlab')) % ... add
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 1e-2; % Collision frequency
+NU = 1e-1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
ETAB = 0.5; % Magnetic gradient
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
-MU = 5e-6; % Hyper diffusivity coefficient
+MU = 5e-4; % Hyper diffusivity coefficient
NOISE0 = 1.0e-5;
%% GRID PARAMETERS
-N = 32; % Frequency gridpoints (Nkr = N/2)
-L = 10; % Size of the squared frequency domain
-PMAXE = 2; % Highest electron Hermite polynomial degree
-JMAXE = 1; % Highest '' Laguerre ''
-PMAXI = 2; % Highest ion Hermite polynomial degree
-JMAXI = 1; % Highest '' Laguerre ''
+N = 256; % Frequency gridpoints (Nkr = N/2)
+L = 66; % Size of the squared frequency domain
+PMAXE = 5; % Highest electron Hermite polynomial degree
+JMAXE = 3; % Highest '' Laguerre ''
+PMAXI = 5; % Highest ion Hermite polynomial degree
+JMAXI = 3; % Highest '' Laguerre ''
+CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
%% TIME PARAMETERS
-TMAX = 30; % Maximal time unit
-DT = 1e-2; % Time step
-SPS2D = 2; % Sampling per time unit for 2D arrays
+TMAX = 400; % Maximal time unit
+DT = 4e-3; % Time step
+SPS0D = 1.0; % Sampling per time unit for 2D arrays
+SPS2D = 1; % Sampling per time unit for 2D arrays
SPS5D = 0.1; % Sampling per time unit for 5D arrays
-RESTART = 0; % To restart from last checkpoint
+RESTART = 1; % To restart from last checkpoint
JOB2LOAD= 0;
%% OPTIONS
-SIMID = 'test_parallel_lin'; % Name of the simulation
-CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+SIMID = 'ZP'; % Name of the simulation
+CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% unused
@@ -37,6 +39,6 @@ NO_E = 0; % Remove electrons dynamic
KREQ0 = 0; % put kr = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-NON_LIN = 0 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
+NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
setup
diff --git a/wk/profiler.m b/wk/profiler.m
new file mode 100644
index 0000000000000000000000000000000000000000..30aa7bfc8e92e29262583d8a19a448de62a1b524
--- /dev/null
+++ b/wk/profiler.m
@@ -0,0 +1,37 @@
+%% load profiling
+filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
+
+CPUTIME = h5readatt(filename,'/data/input','cpu_time');
+DT_SIM = h5readatt(filename,'/data/input','dt');
+
+
+rhs_Tc = h5read(filename,'/profiler/Tc_rhs');
+adv_field_Tc = h5read(filename,'/profiler/Tc_adv_field');
+poisson_Tc = h5read(filename,'/profiler/Tc_poisson');
+Sapj_Tc = h5read(filename,'/profiler/Tc_Sapj');
+diag_Tc = h5read(filename,'/profiler/Tc_diag');
+checkfield_Tc= h5read(filename,'/profiler/Tc_checkfield');
+step_Tc = h5read(filename,'/profiler/Tc_step');
+Ts0D = h5read(filename,'/profiler/time');
+
+missing_Tc = step_Tc - rhs_Tc - adv_field_Tc -...
+ poisson_Tc - Sapj_Tc -diag_Tc -checkfield_Tc;
+total_Tc = step_Tc;
+%% Plots
+Y = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc);...
+ diff(Sapj_Tc); diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
+Y = reshape(Y,[numel(Y)/7,7]);
+fig = figure;
+
+p1 = area(Ts0D(2:end),100*Y./diff(total_Tc),'LineStyle','none');
+legend('Compute RHS','Adv. fields','Poisson','Sapj','Check+Sym','Diag','Missing')
+xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
+ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
+hold on
+yyaxis right
+p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0);
+ylabel('Step Comp. Time [s]')
+ylim([0,1.1*max(diff(total_Tc))])
+set(gca,'ycolor','r')
+FIGNAME = 'profiler';
+save_figure
\ No newline at end of file
diff --git a/wk/setup.m b/wk/setup.m
index 8b13ed6aae7396edc4800d6a229039038dd7b170..5c48f1a93cc3f1f905977138f1375c7f5c2c3797 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -10,6 +10,11 @@ GRID.Lr = L * (1-KREQ0); % r length
GRID.Nz = N; % z ''
GRID.Lz = L; % z ''
GRID.kpar = KPAR;
+if CANCEL_ODD_P
+GRID.CANCEL_ODD_P = '.true.';
+else
+GRID.CANCEL_ODD_P = '.false.';
+end
% Model parameters
MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
if 0; MODEL.DK = '.true.'; else; MODEL.DK = '.false.';end;
@@ -73,7 +78,7 @@ BASIC.dt = DT;
BASIC.tmax = TMAX; %time normalized to 1/omega_pe
% Outputs parameters
if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
-OUTPUTS.nsave_0d = 0;
+OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT);
OUTPUTS.nsave_1d = 0;
OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);