Skip to content
Snippets Groups Projects
Commit 003c315a authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

First MPI version (not working)

parent 8775464b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1050 additions and 259 deletions
.gitignore
+ 1
− 0
View file @ 003c315a
......@@ -31,3 +31,4 @@ wk/fort.90
.directory
checkpoint/
FM/
wk/test*
local/dirs.inc
+ 3
− 1
View file @ 003c315a
......@@ -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)
local/make.inc
+ 6
− 4
View file @ 003c315a
......@@ -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
#
......
matlab/create_gif_surf.m 0 → 100644
+ 56
− 0
View file @ 003c315a
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
matlab/write_fort90.m
+ 1
− 0
View file @ 003c315a
......@@ -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');
......
src/advance_field.F90
+ 0
− 2
View file @ 003c315a
......@@ -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
......
src/auxval.F90
+ 8
− 1
View file @ 003c315a
......@@ -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
src/basic_mod.F90
+ 16
− 2
View file @ 003c315a
......@@ -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)
......
src/compute_Sapj.F90
+ 16
− 8
View file @ 003c315a
......@@ -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
src/control.F90
+ 36
− 28
View file @ 003c315a
......@@ -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
src/diagnose.F90
+ 154
− 73
View file @ 003c315a
......@@ -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
src/endrun.F90
+ 3
− 3
View file @ 003c315a
......@@ -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
......
src/fourier_mod.F90
+ 74
− 10
View file @ 003c315a
......@@ -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
......
src/grid_mod.F90
+ 37
− 12
View file @ 003c315a
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
......
src/inital.F90
+ 70
− 25
View file @ 003c315a
......@@ -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!"
......
src/moments_eq_rhs.F90
+ 111
− 81
View file @ 003c315a
......@@ -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
src/poisson.F90
+ 8
− 0
View file @ 003c315a
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
src/ppinit.F90
+ 13
− 4
View file @ 003c315a
......@@ -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
src/ppsetup.f90 0 → 100644
+ 437
− 0
View file @ 003c315a
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
src/readinputs.F90
+ 0
− 5
View file @ 003c315a
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment