diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index c7f34640a56cc026a9952f777ab4e3f67a3531dd..bd2a2b569ab662eeca8dca1a51e2401ff89f4de5 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -212,10 +212,10 @@ while(CONTINUE)
if W_NAPJ
[Napj, Ts5D, ~] = load_5D_data(filename, 'moments');
tic
- Nipj_ = cat(6,Nipj_,Napj(1,:,:,:,:,:,:)); clear Nipj
+ Nipj_ = cat(6,Nipj_,squeeze(Napj(1,:,:,:,:,:,:))); clear Nipj
toc
if KIN_E
- Nepj_ = cat(6,Nepj_,Napj(2,:,:,:,:,:,:)); clear Nepj
+ Nepj_ = cat(6,Nepj_,squeeze(Napj(2,:,:,:,:,:,:))); clear Nepj
end
end
if W_SAPJ
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 9523e62256c9a29f78fcdaf7e9032310f7e018b9..d4f64d9b41259483c7d5cc646417602e13e1ac3c 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -2,48 +2,54 @@ function [] = profiler(data)
%% load profiling
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
% outfilename = ['/misc/HeLaZ_outputs',filename(3:end)];
-outfilename = data.outfilenames{end};
-CPUTIME = double(h5readatt(outfilename,'/data/input','cpu_time'));
-DT_SIM = h5readatt(outfilename,'/data/input/basic','dt');
+CPUTI=[]; DTSIM=[]; RHSTC=[]; POITC=[]; SAPTC=[]; COLTC=[];
+GRATC=[]; NADTC=[]; ADVTC=[]; GHOTC=[]; CLOTC=[]; CHKTC=[];
+DIATC=[]; STETC=[]; TS0TC=[];
+for i = 1:numel(data.outfilenames)
+ outfilename = data.outfilenames{i};
+ CPUTI = [ CPUTI; double(h5readatt(outfilename,'/data/input','cpu_time'))];
+ DTSIM = [ DTSIM; h5readatt(outfilename,'/data/input/basic','dt')];
+ RHSTC = [ RHSTC; h5read(outfilename,'/profiler/Tc_rhs')];
+ POITC = [ POITC; h5read(outfilename,'/profiler/Tc_poisson')];
+ SAPTC = [ SAPTC; h5read(outfilename,'/profiler/Tc_Sapj')];
+ COLTC = [ COLTC; h5read(outfilename,'/profiler/Tc_coll')];
+ GRATC = [ GRATC; h5read(outfilename,'/profiler/Tc_grad')];
+ NADTC = [ NADTC; h5read(outfilename,'/profiler/Tc_nadiab')];
+ ADVTC = [ ADVTC; h5read(outfilename,'/profiler/Tc_adv_field')];
+ GHOTC = [ GHOTC; h5read(outfilename,'/profiler/Tc_ghost')];
+ CLOTC = [ CLOTC; h5read(outfilename,'/profiler/Tc_clos')];
+ CHKTC = [ CHKTC; h5read(outfilename,'/profiler/Tc_checkfield')];
+ DIATC = [ DIATC; h5read(outfilename,'/profiler/Tc_diag')];
+ STETC = [ STETC; h5read(outfilename,'/profiler/Tc_step')];
+ TS0TC = [ TS0TC; h5read(outfilename,'/profiler/time')];
+end
-rhs_Tc = h5read(outfilename,'/profiler/Tc_rhs');
-poisson_Tc = h5read(outfilename,'/profiler/Tc_poisson');
-Sapj_Tc = h5read(outfilename,'/profiler/Tc_Sapj');
-coll_Tc = h5read(outfilename,'/profiler/Tc_coll');
-process_Tc = h5read(outfilename,'/profiler/Tc_process');
-adv_field_Tc = h5read(outfilename,'/profiler/Tc_adv_field');
-ghost_Tc = h5read(outfilename,'/profiler/Tc_ghost');
-clos_Tc = h5read(outfilename,'/profiler/Tc_clos');
-checkfield_Tc= h5read(outfilename,'/profiler/Tc_checkfield');
-diag_Tc = h5read(outfilename,'/profiler/Tc_diag');
-step_Tc = h5read(outfilename,'/profiler/Tc_step');
-Ts0D = h5read(outfilename,'/profiler/time');
-
-N_T = 11;
+N_T = 12;
-missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ...
- -coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc-process_Tc;
-total_Tc = step_Tc;
+missing_Tc = STETC - RHSTC - ADVTC - GHOTC - CLOTC ...
+ -COLTC - POITC - SAPTC - CHKTC - DIATC - GRATC - NADTC;
+total_Tc = STETC;
-TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(ghost_Tc);...
- diff(clos_Tc); diff(coll_Tc); diff(poisson_Tc); diff(Sapj_Tc); ...
- diff(checkfield_Tc); diff(diag_Tc); diff(process_Tc); diff(missing_Tc)];
+TIME_PER_FCT = [diff(RHSTC); diff(ADVTC); diff(GHOTC);...
+ diff(CLOTC); diff(COLTC); diff(POITC); diff(SAPTC); ...
+ diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(missing_Tc)];
TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]);
TIME_PER_STEP = sum(TIME_PER_FCT,2);
-TIME_PER_CPU = trapz(Ts0D(2:end),TIME_PER_STEP);
+TIME_PER_CPU = trapz(TS0TC(2:end),TIME_PER_STEP);
-rhs_Ta = mean(diff(rhs_Tc));
-adv_field_Ta = mean(diff(adv_field_Tc));
-ghost_Ta = mean(diff(ghost_Tc));
-clos_Ta = mean(diff(clos_Tc));
-coll_Ta = mean(diff(coll_Tc));
-poisson_Ta = mean(diff(poisson_Tc));
-Sapj_Ta = mean(diff(Sapj_Tc));
-checkfield_Ta = mean(diff(checkfield_Tc));
-process_Ta = mean(diff(process_Tc));
-diag_Ta = mean(diff(diag_Tc));
+rhs_Ta = mean(diff(RHSTC));
+adv_field_Ta = mean(diff(ADVTC));
+ghost_Ta = mean(diff(GHOTC));
+clos_Ta = mean(diff(CLOTC));
+coll_Ta = mean(diff(COLTC));
+poisson_Ta = mean(diff(POITC));
+Sapj_Ta = mean(diff(SAPTC));
+checkfield_Ta = mean(diff(CHKTC));
+grad_Ta = mean(diff(GRATC));
+nadiab_Ta = mean(diff(NADTC));
+diag_Ta = mean(diff(DIATC));
miss_Ta = mean(diff(missing_Tc));
total_Ta = mean(diff(total_Tc));
names = {...
@@ -51,33 +57,51 @@ names = {...
'Advf';
'Ghst';
'Clos';
- 'Coll';
+ 'Capj';
'Pois';
'Sapj';
'Chck';
'Diag';
- 'Proc';
+ 'Grad';
+ 'napj';
'Miss';
};
-Ts_A = [rhs_Tc(end) adv_field_Tc(end) ghost_Tc(end) clos_Tc(end) coll_Tc(end)...
- poisson_Tc(end) Sapj_Tc(end) checkfield_Tc(end) diag_Tc(end) process_Tc(end) missing_Tc(end)];
-NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM;
+Ts_A = [RHSTC(end) ADVTC(end) GHOTC(end) CLOTC(end) COLTC(end)...
+ POITC(end) SAPTC(end) CHKTC(end) DIATC(end) GRATC(end) NADTC(end) missing_Tc(end)];
+NSTEP_PER_SAMP= mean(diff(TS0TC))/DTSIM;
%% Plots
if 1
%% Area plot
fig = figure;
% colors = rand(N_T,3);
-colors = lines(N_T);
-p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none');
+% colors = lines(N_T);
+colors = distinguishable_colors(N_T);
+x_ = TS0TC(2:end);
+y_ = TIME_PER_FCT;
+xx_= zeros(2*numel(x_),1);
+yy_= zeros(2*numel(x_),numel(names));
+dx = (x_(2)-x_(1));
+xx_(1) = x_(1)-dx/2; xx_(2) = x_(1)+dx/2;
+yy_(1,:) = y_(1,:)/dx;
+yy_(2,:) = y_(2,:)/dx;
+for i = 2:numel(x_)
+ dx = x_(i) - x_(i-1);
+ xx_(2*i-1) = x_(i)-dx/2;
+ xx_(2*i ) = x_(i)+dx/2;
+ yy_(2*i-1,:) = y_(i,:)/dx;
+ yy_(2*i ,:) = y_(i,:)/dx;
+end
+p1 = area(xx_,yy_,'LineStyle','none');
for i = 1:N_T; p1(i).FaceColor = colors(i,:);
LEGEND{i} = sprintf('%s t=%1.1e[s] (%0.1f %s)',names{i},Ts_A(i),Ts_A(i)/total_Tc(end)*100,'\%');
end;
-legend(LEGEND);
+legend(LEGEND,'Location','bestoutside');
% legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
-xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
-xlim([Ts0D(2),Ts0D(end)]);
-title(sprintf('Proc. 1, total sim. time ~%.0f [h]',CPUTIME/3600))
+xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('N Steps Comp. Time [s]')
+xlim([TS0TC(2),TS0TC(end)]);
+title(sprintf('Gyacomo 2, CPU time: ~%.0f [h] ~%.0f [min] ~%.0f [s]',...
+ CPUTI/3600,CPUTI/60-60*floor(CPUTI/3600),CPUTI-60*floor(CPUTI/60-60*floor(CPUTI/3600))))
hold on
FIGNAME = 'profiler';
% save_figure
diff --git a/src/advance_field_mod.F90 b/src/advance_field_mod.F90
index 1818f8fb8e40bc4ac60ff58af2eb91fe807d1cd6..173b1966d2fce3aedb159b21c962dd0aa5994679 100644
--- a/src/advance_field_mod.F90
+++ b/src/advance_field_mod.F90
@@ -13,7 +13,7 @@ CONTAINS
SUBROUTINE advance_moments
- USE basic, ONLY: t0_adv_field,t1_adv_field,tc_adv_field, dt
+ USE basic, ONLY: dt
USE grid, ONLY:local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
ngp, ngj, ngz, dmax, parray, jarray, dmax
USE model, ONLY: CLOS
@@ -22,7 +22,6 @@ CONTAINS
USE time_integration, ONLY: updatetlevel, A_E, b_E, ntimelevel
IMPLICIT NONE
INTEGER :: ia, ip, ij, ikx,iky,iz, istage, ipi, iji, izi
- CALL cpu_time(t0_adv_field)
SELECT CASE (updatetlevel)
CASE(1)
DO istage=1,ntimelevel
@@ -67,8 +66,5 @@ CONTAINS
END DO
END DO
END SELECT
- ! Execution time end
- CALL cpu_time(t1_adv_field)
- tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field)
END SUBROUTINE advance_moments
END MODULE advance_field_routine
diff --git a/src/auxval.F90 b/src/auxval.F90
index 547b697ab60c84f3ed99ba4f1445241f9cfca170..f9722ce41e5fe0bd6e7acd6333f09ddeac642975 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -52,15 +52,15 @@ subroutine auxval
WRITE(*,'(A9,I3,A10,I3,A10,I3,A9,I3)')&
'my_id = ', my_id, ', rank_p = ', rank_p, ', rank_ky = ', rank_ky,', rank_z = ', rank_z
WRITE(*,'(A22,I3,A11,I3)')&
- ' ips = ', ips , ', ipe = ', ipe
+ ' local_np = ', local_np , ', offset = ', local_np_offset
WRITE(*,'(A22,I3,A11,I3)')&
- ' ijs = ', ijs , ', ije = ', ije
+ ' local_nj = ', local_nj , ', offset = ', local_nj_offset
WRITE(*,'(A22,I3,A11,I3)')&
- ' ikxs = ', ikxs , ', ikxe = ', ikxe
+ ' local_nkx = ', local_nkx , ', offset = ', local_nkx_offset
WRITE(*,'(A22,I3,A11,I3)')&
- ' ikys = ', ikys , ', ikye = ', ikye
+ ' local_nky = ', local_nky , ', offset = ', local_nky_offset
WRITE(*,'(A22,I3,A11,I3)')&
- ' izs = ', izs , ', ize = ', ize
+ ' local_nz = ', local_nz , ', offset = ', local_nz_offset
IF (my_id .NE. num_procs-1) WRITE (*,*) ''
IF (my_id .EQ. num_procs-1) WRITE(*,*) '------------------------------------------'
ENDIF
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index e8953b546d24df99f32e0f6472540ea4cc398742..cc1c14f0aa7f563b438efe8d520cb43c7af1a417 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -31,19 +31,20 @@ MODULE basic
INTEGER, PUBLIC, PROTECTED :: lu_stop = 91 ! stop file, see subroutine TESEND
! To measure computation time
- real(dp), PUBLIC :: start, finish
- real(dp), PUBLIC :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield,&
- t0_step, t0_clos, t0_ghost, t0_coll, t0_process
- real(dp), PUBLIC :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield,&
- t1_step, t1_clos, t1_ghost, t1_coll, t1_process
- real(dp), PUBLIC :: tc_rhs, tc_adv_field, tc_poisson, tc_Sapj, tc_diag, tc_checkfield,&
- tc_step, tc_clos, tc_ghost, tc_coll, tc_process
+ type :: chrono
+ real(dp) :: tstart !start of the chrono
+ real(dp) :: tstop !stop
+ real(dp) :: ttot !cumulative time
+ end type chrono
+
+ type(chrono), PUBLIC, PROTECTED :: chrono_runt, chrono_mrhs, chrono_advf, chrono_pois, chrono_sapj,&
+ chrono_diag, chrono_chck, chrono_step, chrono_clos, chrono_ghst, chrono_coll, chrono_napj, chrono_grad
LOGICAL, PUBLIC, PROTECTED :: GATHERV_OUTPUT = .true.
PUBLIC :: allocate_array, basic_outputinputs,basic_data,&
speak, str, increase_step, increase_cstep, increase_time, display_h_min_s,&
- set_basic_cp, daytim
+ set_basic_cp, daytim, start_chrono, stop_chrono
INTERFACE allocate_array
MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3, &
@@ -72,18 +73,18 @@ CONTAINS
READ(lu_in,basic)
- !Init cumulative timers
- tc_rhs = 0.
- tc_poisson = 0.
- tc_Sapj = 0.
- tc_coll = 0.
- tc_process = 0.
- tc_adv_field = 0.
- tc_ghost = 0.
- tc_clos = 0.
- tc_checkfield = 0.
- tc_diag = 0.
- tc_step = 0.
+ !Init chronometers
+ chrono_mrhs%ttot = 0
+ chrono_pois%ttot = 0
+ chrono_sapj%ttot = 0
+ chrono_napj%ttot = 0
+ chrono_grad%ttot = 0
+ chrono_advf%ttot = 0
+ chrono_ghst%ttot = 0
+ chrono_clos%ttot = 0
+ chrono_chck%ttot = 0
+ chrono_diag%ttot = 0
+ chrono_step%ttot = 0
END SUBROUTINE basic_data
@@ -109,7 +110,7 @@ CONTAINS
CALL attach(fid, TRIM(str), "nrun", nrun)
CALL attach(fid, TRIM(str), "cpu_time", -1)
END SUBROUTINE basic_outputinputs
-
+ !! Increments private attributes
SUBROUTINE increase_step
IMPLICIT NONE
step = step + 1
@@ -130,6 +131,18 @@ CONTAINS
time = time_cp
jobnum = jobnum_cp+1
END SUBROUTINE
+ !! Chrono handling
+ SUBROUTINE start_chrono(timer)
+ IMPLICIT NONE
+ type(chrono) :: timer
+ CALL cpu_time(timer%tstart)
+ END SUBROUTINE
+ SUBROUTINE stop_chrono(timer)
+ IMPLICIT NONE
+ type(chrono) :: timer
+ CALL cpu_time(timer%tstop)
+ timer%ttot = timer%ttot + (timer%tstop-timer%tstart)
+ END SUBROUTINE
!================================================================================
! routine to speak in the terminal
SUBROUTINE speak(message)
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 613ec92ee6a064eb41fcaa641206c68ca8f31cd3..a6c75d03c36e99e1316a039e4af4816d5a66a460 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -16,7 +16,7 @@ MODULE calculus
(/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, odd to even
REAL(dp), dimension(-1:2) :: iz_e2o = &
(/ -1._dp/16._dp, 9._dp/16._dp, 9._dp/16._dp, -1._dp/16._dp /) ! grid interpolation, 4th order, even to odd
- PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4
+ PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4, grad_z_5D, grad_z4_5D
CONTAINS
@@ -80,6 +80,30 @@ CONTAINS
END SUBROUTINE grad_z_e2o
END SUBROUTINE grad_z
+SUBROUTINE grad_z_5D(local_nz,ngz,inv_deltaz,f,ddzf)
+ implicit none
+ ! Compute the periodic boundary condition 4 points centered finite differences
+ ! formula among staggered grid or not.
+ ! not staggered : the derivative results must be on the same grid as the field
+ ! staggered : the derivative is computed from a grid to the other
+ INTEGER, INTENT(IN) :: local_nz, ngz
+ REAL(dp), INTENT(IN) :: inv_deltaz
+ COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
+ COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddzf
+ INTEGER :: iz
+ IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
+ DO iz = 1,local_nz
+ ddzf(:,:,:,:,:,iz) = inv_deltaz*&
+ (dz_usu(-2)*f(:,:,:,:,:,iz ) + dz_usu(-1)*f(:,:,:,:,:,iz+1) &
+ +dz_usu( 0)*f(:,:,:,:,:,iz+2) &
+ +dz_usu( 1)*f(:,:,:,:,:,iz+3) + dz_usu( 2)*f(:,:,:,:,:,iz+4))
+ ENDDO
+ ELSE
+ ddzf = 0._dp
+ ENDIF
+END SUBROUTINE grad_z_5D
+
+
SUBROUTINE grad_z2(local_nz,ngz,inv_deltaz,f,ddz2f)
! Compute the second order fourth derivative for periodic boundary condition
implicit none
@@ -100,6 +124,25 @@ SUBROUTINE grad_z2(local_nz,ngz,inv_deltaz,f,ddz2f)
ddz2f = ddz2f * inv_deltaz**2
END SUBROUTINE grad_z2
+SUBROUTINE grad_z4_5D(local_nz,ngz,inv_deltaz,f,ddz4f)
+ ! Compute the second order fourth derivative for periodic boundary condition
+ implicit none
+ INTEGER, INTENT(IN) :: local_nz, ngz
+ REAL(dp), INTENT(IN) :: inv_deltaz
+ COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(IN) :: f
+ COMPLEX(dp),dimension(:,:,:,:,:,:), INTENT(OUT) :: ddz4f
+ INTEGER :: iz
+ IF(ngz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
+ DO iz = 1,local_nz
+ ddz4f(:,:,:,:,:,iz) = inv_deltaz**4*&
+ (dz4_usu(-2)*f(:,:,:,:,:,iz ) + dz4_usu(-1)*f(:,:,:,:,:,iz+1) &
+ +dz4_usu( 0)*f(:,:,:,:,:,iz+2)&
+ +dz4_usu( 1)*f(:,:,:,:,:,iz+3) + dz4_usu( 2)*f(:,:,:,:,:,iz+4))
+ ENDDO
+ ELSE
+ ddz4f = 0._dp
+ ENDIF
+END SUBROUTINE grad_z4_5D
SUBROUTINE grad_z4(local_nz,ngz,inv_deltaz,f,ddz4f)
! Compute the second order fourth derivative for periodic boundary condition
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 4cff3e6e0a7d2026a3330af0e0ba412622c499f5..cd6e231a604768e6f9667b6c82a8cd11ba0f59a2 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -8,17 +8,14 @@ CONTAINS
! Positive Oob indices are approximated with a model
SUBROUTINE apply_closure_model
- USE basic, ONLY: t0_clos, t1_clos, tc_clos
USE prec_const, ONLY: dp
USE model, ONLY: CLOS
- USE grid, ONLY: local_na, local_nky, local_nkx, local_nz,ngz,&
- local_nj,ngj, jarray,&
+ USE grid, ONLY: local_nj,ngj, jarray,&
local_np,ngp, parray, dmax
USE fields, ONLY: moments
USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
- INTEGER :: iz,ikx,iky,ij,ip,ia
- CALL cpu_time(t0_clos)
+ INTEGER ::ij,ip,ia
IF (CLOS .EQ. 0) THEN
! zero truncation, An+1=0 for n+1>nmax only
CALL ghosts_upper_truncation
@@ -28,125 +25,62 @@ SUBROUTINE apply_closure_model
! only Napj s.t. p+2j <= 3 are evolved
! -> (p,j) allowed are (0,0),(1,0),(0,1),(2,0),(1,1),(3,0)
! =>> Dmax = min(Pmax,2*Jmax+1)
- ! z: DO iz = 1,local_nz+ngz
- ! kx:DO ikx= 1,local_nkx
- ! ky:DO iky= 1,local_nky
j: DO ij = 1,local_nj+ngj
p: DO ip = 1,local_np+ngp
IF ( parray(ip)+2*jarray(ij) .GT. dmax) THEN
- ! a:DO ia = 1,local_na
- ! moments(ia,ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
moments(ia,ip,ij,:,:,:,updatetlevel) = 0._dp
- ! ENDDO a
ENDIF
ENDDO p
ENDDO j
- ! ENDDO ky
- ! ENDDO kx
- ! ENDDO z
ELSE
ERROR STOP '>> ERROR << Closure scheme not found '
ENDIF
CALL ghosts_lower_truncation
- CALL cpu_time(t1_clos)
- tc_clos = tc_clos + (t1_clos - t0_clos)
END SUBROUTINE apply_closure_model
! Positive Oob indices are approximated with a model
SUBROUTINE ghosts_upper_truncation
USE prec_const, ONLY: dp
- USE grid, ONLY: local_na, local_np,ngp,Pmax,&
- local_nj,ngj,Jmax,&
- local_nky,local_nkx,&
- local_nz,ngz,&
- local_pmax, local_jmax
+ USE grid, ONLY: local_np,ngp,local_pmax, pmax,&
+ local_nj,ngj,local_jmax, jmax
USE fields, ONLY: moments
USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
- INTEGER :: iz,ikx,iky,ip,ij,ia,ig
+ INTEGER ::ig
! zero truncation, An+1=0 for n+1>nmax
! applies only for the processes that evolve the highest moment degree
IF(local_jmax .GE. Jmax) THEN
- ! DO iz = 1,local_nz+ngz
- ! DO ikx= 1,local_nkx
- ! DO iky= 1,local_nky
DO ig = 1,ngj/2
- ! DO ip = 1,local_np+ngp
- ! DO ia = 1,local_na
- ! moments(ia,ip,local_nj+ngj/2+ig,iky,ikx,iz,updatetlevel) = 0._dp
moments(:,:,local_nj+ngj/2+ig,:,:,:,updatetlevel) = 0._dp
- ! ENDDO
- ! ENDDO
ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
ENDIF
! applies only for the process that has largest p index
IF(local_pmax .GE. Pmax) THEN
- ! DO iz = 1,local_nz+ngz
- ! DO ikx = 1,local_nkx
- ! DO iky = 1,local_nky
- ! DO ij = 1,local_nj+ngj
DO ig = 1,ngp/2
- ! DO ia = 1,local_na
- ! moments(ia,local_np+ngp/2+ig,ij,iky,ikx,iz,updatetlevel) = 0._dp
moments(:,local_np+ngp/2+ig,:,:,:,:,updatetlevel) = 0._dp
- ! ENDDO
ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
ENDIF
END SUBROUTINE ghosts_upper_truncation
! Negative OoB indices are 0
SUBROUTINE ghosts_lower_truncation
USE prec_const, ONLY: dp
- USE grid, ONLY: local_na,local_np,ngp,&
- local_nj,ngj,&
- local_nky,local_nkx,&
- local_nz,ngz,&
- local_pmin, local_jmin
+ USE grid, ONLY: ngp,ngj,local_pmin,local_jmin
USE fields, ONLY: moments
USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
- INTEGER :: iz,ikx,iky,ip,ia,ij,ig
+ INTEGER :: ig
! zero truncation, An=0 for n<0
IF(local_jmin .EQ. 0) THEN
- ! DO iz = 1,local_nz+ngz
- ! DO ikx = 1,local_nkx
- ! DO iky = 1,local_nky
DO ig = 1,ngj/2
- ! DO ip = 1,local_np+ngp
- ! DO ia = 1,local_na
- ! set to zero the first ghosts cells
- ! moments(ia,ip,ig,iky,ikx,iz,updatetlevel) = 0._dp
moments(:,:,ig,:,:,:,updatetlevel) = 0._dp
- ! ENDDO
- ! ENDDO
ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
ENDIF
! applies only for the process that has lowest p index
IF(local_pmin .EQ. 0) THEN
- ! DO iz = 1,local_nz+ngz
- ! DO ikx = 1,local_nkx
- ! DO iky = 1,local_nky
- ! DO ij = 1,local_nj+ngj
DO ig = 1,ngp/2
- ! DO ia = 1,local_na
- ! moments(ia,ig,ij,iky,ikx,iz,updatetlevel) = 0._dp
moments(:,ig,:,:,:,:,updatetlevel) = 0._dp
ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
ENDIF
END SUBROUTINE ghosts_lower_truncation
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index b4ed42cbd375090d1fd5bd0b899d9edc79decb79..a223da4d13c924a9ff015b7067860a82962d6bb7 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -83,14 +83,10 @@ CONTAINS
END SUBROUTINE coll_outputinputs
SUBROUTINE compute_Capj
- USE basic, ONLY: tc_coll, t0_coll, t1_coll
USE array, ONLY: Capj
USE model, ONLY: nu
USE cosolver_interface, ONLY: compute_cosolver_coll
IMPLICIT NONE
- ! Execution time start
- CALL cpu_time(t0_coll)
-
IF (nu .NE. 0) THEN
SELECT CASE(collision_model)
CASE ('LB')
@@ -107,10 +103,6 @@ CONTAINS
ELSE
Capj = 0._dp
ENDIF
-
- ! Execution time end
- CALL cpu_time(t1_coll)
- tc_coll = tc_coll + (t1_coll - t0_coll)
END SUBROUTINE compute_Capj
!******************************************************************************!
diff --git a/src/control.F90 b/src/control.F90
index d4a965bcdad38429c1d643c7d0aabfe88c7c86c9..b70961cd79b794af3ad77b8b3aed15ef75eadede 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -1,15 +1,17 @@
SUBROUTINE control
! Control the run
- use basic, ONLY: str,daytim,speak,basic_data,start,t0_step,t1_step,tc_step,&
- nlend,step,increase_step,increase_time,increase_cstep
+ use basic, ONLY: str,daytim,speak,basic_data,&
+ nlend,step,increase_step,increase_time,increase_cstep,&
+ chrono_runt,chrono_step, chrono_diag, start_chrono, stop_chrono
use prec_const, ONLY: dp, stdout
USE parallel, ONLY: ppinit
USE mpi
IMPLICIT NONE
REAL(dp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
- CALL cpu_time(start)
+ ! start the chronometer for the total runtime
+ CALL start_chrono(chrono_runt)
!________________________________________________________________________________
! 1. Prologue
! 1.1 Initialize the parallel environment
@@ -60,34 +62,43 @@ SUBROUTINE control
!________________________________________________________________________________
! 2. Main loop
DO
- CALL cpu_time(t0_step) ! Measuring time
+ CALL start_chrono(chrono_step) ! Measuring time per step
+ ! Test if the stopping requirements are met (update nlend)
CALL tesend
IF( nlend ) EXIT ! exit do loop
+ ! Increment steps and csteps (private in basic module)
CALL increase_step
CALL increase_cstep
+ ! Do a full RK step (e.g. 4 substeps for ERK4)
CALL stepon
+ ! Increment time (private in basic module)
CALL increase_time
- CALL diagnose(step)
+ ! Periodic diagnostics
+ CALL start_chrono(chrono_diag)
+ CALL diagnose(step)
+ CALL stop_chrono(chrono_diag)
- CALL cpu_time(t1_step);
- tc_step = tc_step + (t1_step - t0_step)
+ CALL stop_chrono(chrono_step)
END DO
CALL speak('...time integration done')
!________________________________________________________________________________
! 9. Epilogue
-
+ ! Stop total run chronometer (must be done before the last diagnostic)
+ CALL stop_chrono(chrono_runt)
+ ! last diagnostic
CALL diagnose(-1)
+ ! end the run
CALL endrun
-
+ ! display final time
CALL daytim('Done at ')
-
+ ! close mpi environement
CALL ppexit
END SUBROUTINE control
diff --git a/src/cosolver_interface_mod.F90 b/src/cosolver_interface_mod.F90
index 019de6f86a75ae4c6ad0fba930439de1d2c7c2e0..f5b765e9deba5880d4433b34f5d20e476889ca49 100644
--- a/src/cosolver_interface_mod.F90
+++ b/src/cosolver_interface_mod.F90
@@ -69,9 +69,8 @@ CONTAINS
USE grid, ONLY: &
total_na, &
local_np, parray,parray_full, ngp,&
- total_nj, jarray,jarray_full, dmax, ngj, bar, ngz
+ total_nj, jarray,jarray_full, ngj, bar, ngz
USE array, ONLY: nuCself, Cab_F, nadiab_moments
- USE model, ONLY: CLOS
USE species, ONLY: nu_ab
IMPLICIT NONE
INTEGER, INTENT(IN) :: ia,ip,ij,iky,ikx,iz,ikx_C,iky_C,iz_C
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 95dd1774baa4591f4b57874f4640dac1cc22995e..78243c2a8b56d925c7b9d56bb78c0ec943fa7338 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -1,12 +1,11 @@
SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
- USE basic, ONLY: t0_diag,t1_diag,tc_diag, lu_in, finish, start, cstep, dt, time, tmax, display_h_min_s
+ USE basic, ONLY: lu_in, chrono_runt, cstep, dt, time, tmax, display_h_min_s
USE diagnostics_par, ONLY: input_fname
USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
- CALL cpu_time(t0_diag) ! Measuring time
!! Basic diagnose loop for reading input file, displaying advancement and ending
IF ((kstep .EQ. 0)) THEN
INQUIRE(unit=lu_in, name=input_fname)
@@ -14,9 +13,8 @@ SUBROUTINE diagnose(kstep)
ENDIF
!! End diag
IF (kstep .EQ. -1) THEN
- CALL cpu_time(finish)
- ! Display computational time cost
- CALL display_h_min_s(finish-start)
+ ! Display total run time
+ CALL display_h_min_s(chrono_runt%ttot)
! Show last state transport values
IF (my_id .EQ. 0) &
WRITE(*,"(A,G10.2,A8,G10.2,A)") 'Final transport values : | Gxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|'
@@ -27,7 +25,6 @@ SUBROUTINE diagnose(kstep)
IF ((kstep .GE. 0) .AND. (MOD(cstep, INT(1.0/dt)) == 0) .AND. (my_id .EQ. 0)) THEN
WRITE(*,"(A,F6.0,A1,F6.0,A8,G10.2,A8,G10.2,A)")'|t/tmax = ', time,"/",tmax,'| Gxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'|'
ENDIF
- CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
END SUBROUTINE diagnose
SUBROUTINE init_outfile(comm,file0,file,fid)
@@ -91,9 +88,8 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
END SUBROUTINE init_outfile
SUBROUTINE diagnose_full(kstep)
- USE basic, ONLY: speak,&
- cstep,iframe0d,iframe3d,iframe5d,&
- start,finish,crashed
+ USE basic, ONLY: speak,chrono_runt,&
+ cstep,iframe0d,iframe3d,iframe5d,crashed
USE grid, ONLY: &
local_nj,local_nky,local_nkx,local_nz,ngj,ngz,&
parray_full,pmax,jarray_full,jmax,&
@@ -122,7 +118,8 @@ SUBROUTINE diagnose_full(kstep)
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_coll", "cumulative collision computation time")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_process", "cumulative process computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_grad", "cumulative grad computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_nadiab", "cumulative nadiab moments computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field", "cumulative adv. fields computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_ghost", "cumulative communication time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_clos", "cumulative closure computation time")
@@ -252,7 +249,7 @@ SUBROUTINE diagnose_full(kstep)
!_____________________________________________________________________________
! 3. Final diagnostics
ELSEIF (kstep .EQ. -1) THEN
- CALL attach(fidres, "/data/input","cpu_time",finish-start)
+ CALL attach(fidres, "/data/input","cpu_time",chrono_runt%ttot)
! make a checkpoint at last timestep if not crashed
IF(.NOT. crashed) THEN
IF(my_id .EQ. 0) write(*,*) 'Saving last state'
@@ -277,18 +274,19 @@ SUBROUTINE diagnose_0d
CHARACTER :: letter_a
INTEGER :: ia
! Time measurement data
- 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_clos", tc_clos,ionode=0)
- CALL append(fidres, "/profiler/Tc_ghost", tc_ghost,ionode=0)
- CALL append(fidres, "/profiler/Tc_coll", tc_coll,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_checkfield",tc_checkfield,ionode=0)
- CALL append(fidres, "/profiler/Tc_diag", tc_diag,ionode=0)
- CALL append(fidres, "/profiler/Tc_process", tc_process,ionode=0)
- CALL append(fidres, "/profiler/Tc_step", tc_step,ionode=0)
- CALL append(fidres, "/profiler/time", time,ionode=0)
+ CALL append(fidres, "/profiler/Tc_rhs", chrono_mrhs%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_adv_field", chrono_advf%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_clos", chrono_clos%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_ghost", chrono_ghst%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_coll", chrono_coll%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_poisson", chrono_pois%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_Sapj", chrono_sapj%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_checkfield",chrono_chck%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_diag", chrono_diag%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_grad", chrono_grad%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_nadiab", chrono_napj%ttot,ionode=0)
+ CALL append(fidres, "/profiler/Tc_step", chrono_step%ttot,ionode=0)
+ CALL append(fidres, "/profiler/time", time,ionode=0)
! Processing data
CALL append(fidres, "/data/var0d/time", time,ionode=0)
CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index bd83c1196b6e42d323efc545c79da807ffe2bbd1..c578ad1b654e5fc030403ccdab0ef50a1497e83c 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -12,33 +12,25 @@ CONTAINS
SUBROUTINE update_ghosts_moments
USE grid, ONLY: total_nz
USE parallel, ONLY: num_procs_p
- USE basic, ONLY: t0_ghost,t1_ghost,tc_ghost
IMPLICIT NONE
- CALL cpu_time(t0_ghost)
IF (num_procs_p .GT. 1) THEN ! Do it only if we share the p
CALL update_ghosts_p_mom
ENDIF
IF(total_nz .GT. 1) THEN
CALL update_ghosts_z_mom
ENDIF
- CALL cpu_time(t1_ghost)
- tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
END SUBROUTINE update_ghosts_moments
SUBROUTINE update_ghosts_EM
USE model, ONLY : beta
USE grid, ONLY: total_nz
USE fields, ONLY: phi, psi
- USE basic, ONLY: t0_ghost,t1_ghost,tc_ghost
IMPLICIT NONE
- CALL cpu_time(t0_ghost)
IF(total_nz .GT. 1) THEN
CALL update_ghosts_z_3D(phi)
IF(beta .GT. 0._dp) &
CALL update_ghosts_z_3D(psi)
ENDIF
- CALL cpu_time(t1_ghost)
- tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
END SUBROUTINE update_ghosts_EM
@@ -63,23 +55,7 @@ SUBROUTINE update_ghosts_p_mom
USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
ngp,ngj,ngz
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(local_np+ngp) :: slice_p
- !! SUPER SLOW
- ! INTEGER :: ia,iz,ij,iy,ix
- ! DO iz = 1+ngz/2,local_nz+ngz/2
- ! DO ix = 1,local_nkx
- ! DO iy = 1,local_nky
- ! DO ij = 1+ngj/2,local_nj+ngj/2
- ! DO ia = 1,local_na
- ! slice_p = moments(ia,:,ij,iy,ix,iz,updatetlevel)
- ! CALL exchange_ghosts_1D(slice_p,nbr_L,nbr_R,local_np,ngp)
- ! moments(ia,:,ij,iy,ix,iz,updatetlevel) = slice_p
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- INTEGER :: ierr, first, last, ig, count
+ INTEGER :: ierr, first, last, count
first = 1 + ngp/2
last = local_np + ngp/2
! count = local_na*(local_nj+ngj)*local_nky*local_nkx*(local_nz+ngz) ! Number of elements to send
@@ -97,11 +73,11 @@ SUBROUTINE update_ghosts_p_mom
! ENDDO
count = (ngp/2)*local_na*(local_nj+ngj)*local_nky*local_nkx*(local_nz+ngz) ! Number of elements to send
CALL mpi_sendrecv(moments(:,(last-(ngp/2-1)):(last),:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 14, &
- moments(:,(first-ngp/2):(first-1) ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 14, &
- comm0, status, ierr)
+ moments(:,(first-ngp/2):(first-1),:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 14, &
+ comm0, status, ierr)
CALL mpi_sendrecv(moments(:,(first):(first+(ngp/2-1)),:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 16, &
- moments(:,(last+1):(last+ngp/2),:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 16, &
- comm0, status, ierr)
+ moments(:,(last+1):(last+ngp/2) ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 16, &
+ comm0, status, ierr)
END SUBROUTINE update_ghosts_p_mom
!Communicate z+1, z+2 moments to left neighboor and z-1, z-2 moments to right one
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 367f3ea4d1e718c45c616748e129afd2691c8866..95b9c588021094fdad894aa1e66e13430a355816 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -199,7 +199,7 @@ CONTAINS
USE parallel, ONLY: num_procs_p, rank_p
IMPLICIT NONE
LOGICAL, INTENT(IN) :: EM
- INTEGER :: ip, ig
+ INTEGER :: ip
! If no parallel dim (Nz=1) and no EM effects (beta=0), the moment hierarchy
!! is separable between odds and even P and since the energy is injected in
!! P=0 and P=2 for density/temperature gradients there is no need of
@@ -263,7 +263,7 @@ CONTAINS
SUBROUTINE set_jgrid
USE prec_const
IMPLICIT NONE
- INTEGER :: ij, ig
+ INTEGER :: ij
! Total number of J degrees
total_nj = jmax+1
local_jmax = jmax
@@ -320,7 +320,7 @@ CONTAINS
ikys = local_nky_ptr_offset + 1
ikye = ikys + local_nky_ptr - 1
local_nky = ikye - ikys + 1
- local_nky_offset = ikys - 1
+ local_nky_offset = local_nky_ptr_offset
ALLOCATE(kyarray(local_nky))
local_kymax = 0._dp
! Creating a grid ordered as dk*(0 1 2 3)
@@ -397,7 +397,7 @@ CONTAINS
ikxe = total_nkx
local_nkx_ptr = ikxe - ikxs + 1
local_nkx = ikxe - ikxs + 1
- local_nky_offset = ikxs - 1
+ local_nkx_offset = ikxs - 1
ALLOCATE(kxarray(local_nkx))
ALLOCATE(kxarray_full(total_nkx))
IF (Nx .EQ. 1) THEN
@@ -607,17 +607,20 @@ CONTAINS
CHARACTER(len=256) :: str
WRITE(str,'(a)') '/data/input/grid'
CALL creatd(fid, 0,(/0/),TRIM(str),'Grid Input')
- CALL attach(fid, TRIM(str), "pmax", pmax)
- CALL attach(fid, TRIM(str), "jmax", jmax)
- CALL attach(fid, TRIM(str), "Nx", Nx)
- CALL attach(fid, TRIM(str), "Lx", Lx)
- CALL attach(fid, TRIM(str), "Nexc", Nexc)
- CALL attach(fid, TRIM(str), "Ny", Ny)
- CALL attach(fid, TRIM(str), "Ly", Ly)
- CALL attach(fid, TRIM(str), "Nz", Nz)
- CALL attach(fid, TRIM(str), "total_nkx", total_nkx)
- CALL attach(fid, TRIM(str), "Nky", Nky)
- CALL attach(fid, TRIM(str), "SG", SG)
+ CALL attach(fid, TRIM(str), "pmax", pmax)
+ CALL attach(fid, TRIM(str),"deltap", deltap)
+ CALL attach(fid, TRIM(str), "jmax", jmax)
+ CALL attach(fid, TRIM(str), "Nkx", Nkx)
+ CALL attach(fid, TRIM(str), "Nx", Nx)
+ CALL attach(fid, TRIM(str), "Lx", Lx)
+ CALL attach(fid, TRIM(str), "Nexc", Nexc)
+ CALL attach(fid, TRIM(str), "Ny", Ny)
+ CALL attach(fid, TRIM(str), "Nky", Nky)
+ CALL attach(fid, TRIM(str), "Ly", Ly)
+ CALL attach(fid, TRIM(str), "Nz", Nz)
+ CALL attach(fid, TRIM(str), "total_nkx", total_nkx)
+ CALL attach(fid, TRIM(str), "Nky", Nky)
+ CALL attach(fid, TRIM(str), "SG", SG)
END SUBROUTINE grid_outputinputs
FUNCTION bar(p_,j_)
diff --git a/src/inital.F90 b/src/inital.F90
index 54262f736bf5a6dafe6601571e32d7002fbce89a..d9813ba333fb7086a509ce8bf3a08735307a18b4 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -139,11 +139,11 @@ SUBROUTINE init_moments
IF (LINEARITY .NE. 'linear') THEN
DO ikx=1,total_nkx
DO iky=1,local_nky
- DO iz=1,local_nz + ngz
+ DO iz=1,local_nz
izi = iz+ngz/2
- DO ip=1,local_np + ngp
+ DO ip=1,local_np
ipi = ip+ngp/2
- DO ij=1,local_nj + ngj
+ DO ij=1,local_nj
iji = ij+ngj/2
moments(ia,ipi,iji,iky,ikx,izi, :) = moments(ia, ipi,iji,iky,ikx,izi, :)*AA_x(ikx)*AA_y(iky)
ENDDO
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 65b68f1bdd1971febdb52275188af03c43175ed3..815064048b14b04b6c848e04a466674326786177 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -15,7 +15,7 @@ SUBROUTINE compute_moments_eq_rhs
USE prec_const
USE collision
USE time_integration
- USE species, ONLY: dpdx
+ ! USE species, ONLY: dpdx
USE geometry, ONLY: gradz_coeff, dlnBdz, Ckxky!, Gamma_phipar
USE calculus, ONLY: interp_z, grad_z, grad_z2
! USE species, ONLY: dpdx
@@ -31,8 +31,6 @@ SUBROUTINE compute_moments_eq_rhs
COMPLEX(dp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
COMPLEX(dp) :: i_kx,i_ky
COMPLEX(dp) :: Napj, RHS, phikykxz, psikykxz
- ! Measuring execution time
- CALL cpu_time(t0_rhs)
! Spatial loops
z:DO iz = 1,local_nz
izi = iz + ngz/2
@@ -157,24 +155,22 @@ SUBROUTINE compute_moments_eq_rhs
END DO y
END DO x
END DO z
- ! Execution time end
- CALL cpu_time(t1_rhs)
- tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
- ! print*,'sumabs moments i', SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
- ! print*,'moment rhs i 221', moments_rhs(1,1,1,2,2,1,updatetlevel)
- ! print*,'sum real nadiabe', SUM(REAL(nadiab_moments(2,:,:,:,:,:)))
- ! print*,'sumreal momrhs i', SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumimag momrhs i', SUM(IMAG(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumreal phi ', SUM(REAL(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'sumimag phi ', SUM(IMAG(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'phi(2,2,1) ', phi(2,2,1+ngz/2)
- ! print*,'sumreal ddznipj ', SUM(REAL(ddz_napj(1,:,:,:,:,:)))
- ! print*,'sumimag ddznipj ', SUM(IMAG(ddz_napj(1,:,:,:,:,:)))
- ! print*,' ddznipj ',(ddz_napj(1,2+ngp/2,2+ngj/2,2,2,1))
- ! print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
- ! print*,'---'
- ! IF(updatetlevel .EQ. 4) stop
+ ! print*,'sumabs moments i', SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
+ ! print*,'moment rhs i 221', moments_rhs(1,1,1,2,2,1,updatetlevel)
+ ! ! print*,'sum real nadiabe', SUM(REAL(nadiab_moments(2,:,:,:,:,:)))
+ ! print*,'sumreal momrhs i', SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
+ ! print*,'sumimag momrhs i', SUM(IMAG(moments_rhs(1,:,:,:,:,:,:)))
+ ! print*,'sumreal phi ', SUM(REAL(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
+ ! print*,'sumimag phi ', SUM(IMAG(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
+ ! print*,'phi(2,2,1) ', phi(2,2,1+ngz/2)
+ ! print*,'sumreal ddznipj ', SUM(REAL(ddz_napj(1,:,:,:,:,:)))
+ ! print*,'sumimag ddznipj ', SUM(IMAG(ddz_napj(1,:,:,:,:,:)))
+ ! print*,' ddznipj ',(ddz_napj(1,2+ngp/2,2+ngj/2,2,2,1))
+ ! print*,' ddzNDnipj ',SUM(REAL(ddzND_Napj(1,:,:,:,:,:)))
+ ! print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
+ ! print*,'---'
+ ! IF(updatetlevel .EQ. 4) stop
! stop
END SUBROUTINE compute_moments_eq_rhs
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 44f48f8452d1f14492c669045aada5a2b3bdc198..37b3347a7f0456bc9c75f1af9c56eca4d4d3a942 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -1,7 +1,6 @@
MODULE nonlinear
USE array, ONLY : dnjs, Sapj, kernel
USE initial_par, ONLY : ACT_ON_MODES
- USE basic, ONLY : t0_Sapj, t1_Sapj, tc_Sapj
USE fourier, ONLY : bracket_sum_r, bracket_sum_c, planf, planb, poisson_bracket_and_sum
USE fields, ONLY : phi, psi, moments
USE grid, ONLY: local_na, &
@@ -49,9 +48,6 @@ SUBROUTINE compute_Sapj
!! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
!! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps)
!! where # denotes the convolution.
- ! Execution time start
- CALL cpu_time(t0_Sapj)
-
SELECT CASE(LINEARITY)
CASE ('nonlinear')
CALL compute_nonlinear
@@ -60,9 +56,6 @@ SUBROUTINE compute_Sapj
CASE DEFAULT
ERROR STOP '>> ERROR << Linearity not recognized '
END SELECT
- ! Execution time END
- CALL cpu_time(t1_Sapj)
- tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
END SUBROUTINE compute_Sapj
SUBROUTINE compute_nonlinear
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index e37be6cad785028aee2aebddf34240d4f1785b40..9b59750931b72993d33eb68c0186997b396198ca 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -134,7 +134,7 @@ SUBROUTINE evaluate_poisson_op
USE grid, ONLY : local_na, local_nkx, local_nky, local_nz,&
kxarray, kyarray, local_nj, ngj, ngz, ieven
USE species, ONLY : q2_tau
- USE model, ONLY : ADIAB_E, tau_e
+ USE model, ONLY : ADIAB_E
USE prec_const, ONLY: dp
IMPLICIT NONE
REAL(dp) :: pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
@@ -162,7 +162,7 @@ ELSE
ENDDO a
operator_ion = pol_tot
IF(ADIAB_E) THEN ! Adiabatic electron model
- pol_tot = pol_tot + 1._dp/tau_e - 1._dp
+ pol_tot = pol_tot + 1._dp
ENDIF
operator = pol_tot
inv_poisson_op(iky, ikx, iz) = 1._dp/pol_tot
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 47ede3ef7cc09f32f267f0ffc8f3ac71b56603e2..46532eaad54a0fa1daaaad4b6487b487724f6e32 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -1,5 +1,5 @@
MODULE processing
- USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3
+ USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3, onetwelfth, twothird
USE grid, ONLY: &
local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp,nzgrid, &
parray,pmax,ip0, iodd, ieven,&
@@ -14,10 +14,9 @@ MODULE processing
dens, upar, uper, Tpar, Tper, temp
USE geometry, ONLY: Jacobian, iInt_Jacobian
USE time_integration, ONLY: updatetlevel
- USE calculus, ONLY: simpson_rule_z, grad_z, grad_z2, grad_z4, interp_z
+ USE calculus, ONLY: simpson_rule_z, grad_z, grad_z_5D, grad_z2, grad_z4, grad_z4_5D, interp_z
USE model, ONLY: EM, CLOS, beta, HDz_h
USE species, ONLY: tau,q_tau,q_o_sqrt_tau_sigma,sqrt_tau_o_sigma
- USE basic, ONLY: t0_process, t1_process, tc_process
USE parallel, ONLY: num_procs_ky, rank_ky, comm_ky
USE mpi
implicit none
@@ -26,7 +25,7 @@ MODULE processing
REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:), PROTECTED :: hflux_x
INTEGER :: ierr
PUBLIC :: init_process
- PUBLIC :: compute_nadiab_moments_z_gradients_and_interp
+ PUBLIC :: compute_nadiab_moments, compute_gradients_z, compute_interp_z
PUBLIC :: compute_density, compute_upar, compute_uperp
PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
PUBLIC :: compute_radial_transport
@@ -41,8 +40,165 @@ CONTAINS
ALLOCATE( gflux_x(local_na))
ALLOCATE( hflux_x(local_na))
END SUBROUTINE init_process
+!------------------------------ HIGH FREQUENCY ROUTINES -------------
+! The following routines (nadiab computing, interp and grads) must be
+! as fast as possible since they are called every RK substep.
+ ! evaluate the non-adiabatique ion moments
+ !
+ ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
+ !
+ SUBROUTINE compute_nadiab_moments
+ IMPLICIT NONE
+ INTEGER :: ia,ip,ij,iky,ikx,iz, j_int, p_int
+ !non adiab moments
+ DO iz=1,local_nz+ngz
+ DO ikx=1,local_nkx
+ DO iky=1,local_nky
+ DO ij=1,local_nj+ngj
+ DO ip=1,local_np+ngp
+ DO ia = 1,local_na
+ IF(parray(ip) .EQ. 0) THEN
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
+ + q_tau(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*phi(iky,ikx,iz)
+ ELSEIF( (parray(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
+ - q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*psi(iky,ikx,iz)
+ ELSE
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ !! Ensure to kill the moments too high if closue option is set to 1
+ IF(CLOS .EQ. 1) THEN
+ DO iz=1,local_nz+ngz
+ DO ikx=1,local_nkx
+ DO iky=1,local_nky
+ DO ij=1,local_nj+ngj
+ j_int = jarray(ij)
+ DO ip=1,local_np+ngp
+ p_int = parray(ip)
+ DO ia = 1,local_na
+ IF(p_int+2*j_int .GT. dmax) &
+ nadiab_moments(ia,ip,ij,iky,ikx,iz) = 0._dp
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ END SUBROUTINE compute_nadiab_moments
+
+ ! z grid gradients
+ ! SUBROUTINE compute_gradients_z
+ ! IMPLICIT NONE
+ ! INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz,izi
+ ! COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
+ ! COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ ! ! Compute z first derivative
+ ! DO iz=1,local_nz+ngz
+ ! izi = iz+ngz/2
+ ! DO ikx=1,local_nkx
+ ! DO iky=1,local_nky
+ ! DO ij=1,local_nj+ngj
+ ! DO ip=1,local_np+ngp
+ ! DO ia = 1,local_na
+ ! ddz_napj(ia,ip,ij,iky,ikx,iz) = inv_deltaz *(&
+ ! +onetwelfth*nadiab_moments(ia,ip,ij,iky,ikx,izi-2)&
+ ! -twothird*nadiab_moments(ia,ip,ij,iky,ikx,izi-1)&
+ ! +twothird*nadiab_moments(ia,ip,ij,iky,ikx,izi+1)&
+ ! -onetwelfth*nadiab_moments(ia,ip,ij,iky,ikx,izi-2)&
+ ! )
+ ! ddzND_Napj(ia,ip,ij,iky,ikx,iz) = inv_deltaz**4 *(&
+ ! +1._dp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
+ ! -4._dp*moments(ia,ip,ij,iky,ikx,izi-1,updatetlevel)&
+ ! +6._dp*moments(ia,ip,ij,iky,ikx,izi ,updatetlevel)&
+ ! -4._dp*moments(ia,ip,ij,iky,ikx,izi+1,updatetlevel)&
+ ! +1._dp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
+ ! )
+ ! ENDDO
+ ! ENDDO
+ ! ENDDO
+ ! ENDDO
+ ! ENDDO
+ ! ENDDO
+ ! END SUBROUTINE compute_gradients_z
+
+ ! ! z grid gradients
+ SUBROUTINE compute_gradients_z
+ IMPLICIT NONE
+ INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
+ COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
+ COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj+ngj
+ DO ip = 1,local_np+ngp
+ DO ia = 1,local_na
+ IF(nzgrid .GT. 1) THEN
+ p_int = parray(ip+ngp/2)
+ eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
+ ELSE
+ eo = 0
+ ENDIF
+ ! Compute z first derivative
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
+ CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
+ ddz_napj(ia,ip,ij,iky,ikx,:) = f_out(:)
+ ! Parallel numerical diffusion
+ IF (HDz_h) THEN
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
+ ELSE
+ f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
+ ENDIF
+ CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
+ ! get output
+ DO iz = 1,local_nz
+ ddzND_Napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ END SUBROUTINE compute_gradients_z
-! 1D diagnostic to compute the average radial particle transport <n_a v_ExB_x>_xyz
+ ! z grid interpolation
+ SUBROUTINE compute_interp_z
+ IMPLICIT NONE
+ INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
+ COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
+ COMPLEX(dp), DIMENSION(local_nz) :: f_out
+ IF(nzgrid .GT. 1) THEN
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj+ngj
+ DO ip = 1,local_np+ngp
+ DO ia = 1,local_na
+ ! Compute even odd grids interpolation
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
+ CALL interp_z(eo,local_nz,ngz,f_in,f_out)
+ DO iz = 1,local_nz
+ interp_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSE
+ interp_napj(:,:,:,:,:,1:local_nz) = nadiab_moments(:,:,:,:,:,1:local_nz)
+ ENDIF
+ END SUBROUTINE compute_interp_z
+
+ !--------------------- LOW FREQUENCY PROCESSING ROUTINES -------------------!
+ ! The following routines are called by the diagnose routine every nsave3D steps
+ ! (does not need to be optimized)
+ ! 1D diagnostic to compute the average radial particle transport <n_a v_ExB_x>_xyz
SUBROUTINE compute_radial_transport
IMPLICIT NONE
COMPLEX(dp) :: pflux_local, gflux_local, integral
@@ -57,27 +213,27 @@ CONTAINS
! Electrostatic part
IF(CONTAINSp0) THEN
DO iz = 1,local_nz
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- integrant(iz) = integrant(iz) &
- +Jacobian(izi,ieven)*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
- *imagu*kyarray(iky)*CONJG(phi(iky,ikx,izi))
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ integrant(iz) = integrant(iz) &
+ +Jacobian(izi,ieven)*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*CONJG(phi(iky,ikx,izi))
+ ENDDO
+ ENDDO
ENDDO
ENDIF
! Electromagnetic part
IF( EM .AND. CONTAINSp1 ) THEN
DO iz = 1,local_nz ! we take interior points only
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- integrant(iz) = integrant(iz)&
- -Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ij0,iky,ikx,izi,updatetlevel)&
- *imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ integrant(iz) = integrant(iz)&
+ -Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ij0,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))
+ ENDDO
+ ENDDO
ENDDO
ENDIF
! Integrate over z
@@ -90,33 +246,33 @@ CONTAINS
! Electrostatic part
IF(CONTAINSp0) THEN
DO iz = 1,local_nz ! we take interior points only
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1, local_nj
- ini = in + ngj/2 !interior index for ghosted arrays
- integrant(iz) = integrant(iz)+ &
- Jacobian(izi,ieven)*moments(ia,ip0,ini,iky,ikx,izi,updatetlevel)&
- *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(phi(iky,ikx,izi))
- ENDDO
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1, local_nj
+ ini = in + ngj/2 !interior index for ghosted arrays
+ integrant(iz) = integrant(iz)+ &
+ Jacobian(izi,ieven)*moments(ia,ip0,ini,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(phi(iky,ikx,izi))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
ENDIF
! Electromagnetic part
IF( EM .AND. CONTAINSp1 ) THEN
DO iz = 1,local_nz ! we take interior points only
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1, local_nj
- ini = in + ngj/2 !interior index for ghosted arrays
- integrant(iz) = integrant(iz) - &
- Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ini,iky,ikx,izi,updatetlevel)&
- *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,iodd)*CONJG(psi(iky,ikx,izi))
- ENDDO
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1, local_nj
+ ini = in + ngj/2 !interior index for ghosted arrays
+ integrant(iz) = integrant(iz) - &
+ Jacobian(izi,iodd)*sqrt_tau_o_sigma(ia)*moments(ia,ip1,ini,iky,ikx,izi,updatetlevel)&
+ *imagu*kyarray(iky)*kernel(ia,ini,iky,ikx,izi,iodd)*CONJG(psi(iky,ikx,izi))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
ENDIF
! Integrate over z
@@ -150,7 +306,7 @@ CONTAINS
ENDDO
END SUBROUTINE compute_radial_transport
-! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
+! 1D diagnostic to compute the average radial particle heatflux <T_i v_ExB_x>_xyz
SUBROUTINE compute_radial_heatflux
IMPLICIT NONE
COMPLEX(dp) :: hflux_local, integral
@@ -164,22 +320,22 @@ CONTAINS
IF(CONTAINSp0 .AND. CONTAINSp2) THEN
! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Na00 * phi
DO iz = 1,local_nz ! we take interior points only
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1, local_nj
- ini = in+ngj/2 !interior index for ghosted arrays
- n_dp = jarray(ini)
- integrant(iz) = integrant(iz) &
- -Jacobian(izi,ieven)*tau(ia)*imagu*kyarray(iky)*phi(iky,ikx,izi)&
- *kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(&
- 0.5_dp*SQRT2*moments(ia,ip2,ini ,iky,ikx,izi,updatetlevel)&
- +(2._dp*n_dp + 1.5_dp)*moments(ia,ip0,ini ,iky,ikx,izi,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip0,ini+1,iky,ikx,izi,updatetlevel)&
- -n_dp*moments(ia,ip0,ini-1,iky,ikx,izi,updatetlevel))
- ENDDO
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1, local_nj
+ ini = in+ngj/2 !interior index for ghosted arrays
+ n_dp = jarray(ini)
+ integrant(iz) = integrant(iz) &
+ -Jacobian(izi,ieven)*tau(ia)*imagu*kyarray(iky)*phi(iky,ikx,izi)&
+ *kernel(ia,ini,iky,ikx,izi,ieven)*CONJG(&
+ 0.5_dp*SQRT2*moments(ia,ip2,ini ,iky,ikx,izi,updatetlevel)&
+ +(2._dp*n_dp + 1.5_dp)*moments(ia,ip0,ini ,iky,ikx,izi,updatetlevel)&
+ -(n_dp+1._dp)*moments(ia,ip0,ini+1,iky,ikx,izi,updatetlevel)&
+ -n_dp*moments(ia,ip0,ini-1,iky,ikx,izi,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
ELSEIF(CONTAINSp0) THEN
ERROR STOP "Degrees p=0 and p=2 should be owned by the same process"
@@ -187,23 +343,23 @@ CONTAINS
IF(EM .AND. CONTAINSp1 .AND. CONTAINSp3) THEN
!!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
DO iz = 1,local_nz
- izi = iz + ngz/2 !interior index for ghosted arrays
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO in = 1, local_nj
- ini = in + ngj/2 !interior index for ghosted arrays
- n_dp = jarray(ini)
- integrant(iz) = integrant(iz) &
- +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
- *kernel(ia,ini,iky,ikx,izi,iodd)*(&
- 0.5_dp*SQRT2*SQRT3*moments(ia,ip3,ini ,iky,ikx,izi,updatetlevel)&
- +1.5_dp*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
- +(2._dp*n_dp+1._dp)*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
- -(n_dp+1._dp)*moments(ia,ip1,ini+1,iky,ikx,izi,updatetlevel)&
- -n_dp*moments(ia,ip1,ini-1,iky,ikx,izi,updatetlevel))
- ENDDO
- ENDDO
- ENDDO
+ izi = iz + ngz/2 !interior index for ghosted arrays
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO in = 1, local_nj
+ ini = in + ngj/2 !interior index for ghosted arrays
+ n_dp = jarray(ini)
+ integrant(iz) = integrant(iz) &
+ +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
+ *kernel(ia,ini,iky,ikx,izi,iodd)*(&
+ 0.5_dp*SQRT2*SQRT3*moments(ia,ip3,ini ,iky,ikx,izi,updatetlevel)&
+ +1.5_dp*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
+ +(2._dp*n_dp+1._dp)*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
+ -(n_dp+1._dp)*moments(ia,ip1,ini+1,iky,ikx,izi,updatetlevel)&
+ -n_dp*moments(ia,ip1,ini-1,iky,ikx,izi,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
ENDIF
! Add polarisation contribution
@@ -235,105 +391,6 @@ CONTAINS
ENDDO
END SUBROUTINE compute_radial_heatflux
- SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
- ! evaluate the non-adiabatique ion moments
- !
- ! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
- !
- IMPLICIT NONE
- INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz
- COMPLEX(dp), DIMENSION(local_nz+ngz) :: f_in
- COMPLEX(dp), DIMENSION(local_nz) :: f_out
- CALL cpu_time(t0_process)
-
- !non adiab moments
- DO iz=1,local_nz+ngz
- DO ikx=1,local_nkx
- DO iky=1,local_nky
- DO ij=1,local_nj+ngj
- DO ip=1,local_np+ngp
- DO ia = 1,local_na
- IF(parray(ip) .EQ. 0) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- + q_tau(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*phi(iky,ikx,iz)
- ELSEIF( (parray(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- - q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*psi(iky,ikx,iz)
- ELSE
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- !! Ensure to kill the moments too high if closue option is set to 1
- IF(CLOS .EQ. 1) THEN
- DO ij=1,local_nj+ngj
- j_int = jarray(ij)
- DO ip=1,local_np+ngp
- p_int = parray(ip)
- DO ia = 1,local_na
- IF(p_int+2*j_int .GT. dmax) &
- nadiab_moments(ia,ip,ij,:,:,:) = 0._dp
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- !------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO ij = 1,local_nj+ngj
- DO ip = 1,local_np+ngp
- DO ia = 1,local_na
- IF(nzgrid .GT. 1) THEN
- p_int = parray(ip+ngp/2)
- eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
- ELSE
- eo = 0
- ENDIF
- ! Compute z first derivative
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
- ! get output
- DO iz = 1,local_nz
- ddz_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
- ENDDO
- ! Parallel numerical diffusion
- IF (HDz_h) THEN
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- ELSE
- f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
- ENDIF
- CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
- ! get output
- DO iz = 1,local_nz
- ddzND_Napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
- ENDDO
- ! Compute even odd grids interpolation
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- CALL interp_z(eo,local_nz,ngz,f_in,f_out)
- DO iz = 1,local_nz
- interp_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ! Phi parallel gradient (not implemented fully, should be negligible)
- ! DO ikx = 1,local_nkx
- ! DO iky = 1,local_nky
- ! CALL grad_z(0,phi(iky,ikx,1:local_nz+ngz), ddz_phi(iky,ikx,1:local_nz))
- ! ENDDO
- ! ENDDO
-
- ! Execution time end
- CALL cpu_time(t1_process)
- tc_process = tc_process + (t1_process - t0_process)
- END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
-
SUBROUTINE compute_Napjz_spectrum
USE fields, ONLY : moments
USE array, ONLY : Napjz
@@ -347,17 +404,17 @@ CONTAINS
! build local sum
local_sum = 0._dp
DO iz = 1,local_nz
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO ij = 1,local_nj
- DO ip = 1,local_np
- local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
- (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
- * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
- ENDDO
- ENDDO
- ENDDO
- ENDDO
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj
+ DO ip = 1,local_np
+ local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
+ (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
+ * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
ENDDO
! sum reduction
buffer = local_sum
@@ -386,25 +443,25 @@ CONTAINS
!!!!! FLUID MOMENTS COMPUTATIONS !!!!!
! Compute the 2D particle density for electron and ions (sum over Laguerre)
SUBROUTINE compute_density
- IMPLICIT NONE
- COMPLEX(dp) :: dens_
- INTEGER :: ia, iz, iky, ikx, ij
- DO ia=1,local_na
- IF ( CONTAINSp0 ) THEN
- ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- dens_ = 0._dp
- DO ij = 1, local_nj
- dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
- ENDDO
- dens(ia,iky,ikx,iz) = dens_
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ IMPLICIT NONE
+ COMPLEX(dp) :: dens_
+ INTEGER :: ia, iz, iky, ikx, ij
+ DO ia=1,local_na
+ IF ( CONTAINSp0 ) THEN
+ ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ dens_ = 0._dp
+ DO ij = 1, local_nj
+ dens_ = dens_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) * moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
ENDDO
+ dens(ia,iky,ikx,iz) = dens_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
END SUBROUTINE compute_density
! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
@@ -413,21 +470,21 @@ CONTAINS
COMPLEX(dp) :: uperp_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
- IF ( CONTAINSp0 ) THEN
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- uperp_ = 0._dp
- DO ij = 1, local_nj
- uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
- 0.5_dp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
- ENDDO
- uper(ia,iky,ikx,iz) = uperp_
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ IF ( CONTAINSp0 ) THEN
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ uperp_ = 0._dp
+ DO ij = 1, local_nj
+ uperp_ = uperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven) *&
+ 0.5_dp*(moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ uper(ia,iky,ikx,iz) = uperp_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
ENDDO
END SUBROUTINE compute_uperp
@@ -437,19 +494,19 @@ CONTAINS
INTEGER :: ia, iz, iky, ikx, ij
COMPLEX(dp) :: upar_
DO ia=1,local_na
- IF ( CONTAINSp1 ) THEN
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- upar_ = 0._dp
- DO ij = 1, local_nj
- upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
- ENDDO
- upar(ia,iky,ikx,iz) = upar_
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ IF ( CONTAINSp1 ) THEN
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ upar_ = 0._dp
+ DO ij = 1, local_nj
+ upar_ = upar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,iodd)*moments(ia,ip1,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ ENDDO
+ upar(ia,iky,ikx,iz) = upar_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
ENDDO
END SUBROUTINE compute_upar
@@ -461,24 +518,24 @@ CONTAINS
COMPLEX(dp) :: Tperp_
INTEGER :: ia, iz, iky, ikx, ij
DO ia=1,local_na
- IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
- ! Loop to compute T = 1/3*(Tpar + 2Tperp)
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- Tperp_ = 0._dp
- DO ij = 1, local_nj
- j_dp = jarray(ij+ngj/2)
- Tperp_ = Tperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
- ((2_dp*j_dp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -j_dp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
- -(j_dp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
- ENDDO
- Tper(ia,iky,ikx,iz) = Tperp_
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ IF ( CONTAINSp0 .AND. CONTAINSp2 ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ Tperp_ = 0._dp
+ DO ij = 1, local_nj
+ j_dp = jarray(ij+ngj/2)
+ Tperp_ = Tperp_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
+ ((2_dp*j_dp+1)*moments(ia,ip0,ij +ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -j_dp*moments(ia,ip0,ij-1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)&
+ -(j_dp+1)*moments(ia,ip0,ij+1+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ Tper(ia,iky,ikx,iz) = Tperp_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
ENDDO
END SUBROUTINE compute_Tperp
@@ -489,25 +546,24 @@ CONTAINS
REAL(dp) :: j_dp
COMPLEX(dp) :: Tpar_
INTEGER :: ia, iz, iky, ikx, ij
-
DO ia=1,local_na
- IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
- ! Loop to compute T = 1/3*(Tpar + 2Tperp)
- DO iz = 1,local_nz
- DO iky = 1,local_nky
- DO ikx = 1,local_nkx
- Tpar_ = 0._dp
- DO ij = 1, local_nj
- j_dp = REAL(ij-1,dp)
- Tpar_ = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
- (SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
- + moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
- ENDDO
- Tpar(ia,iky,ikx,iz) = Tpar_
- ENDDO
- ENDDO
- ENDDO
- ENDIF
+ IF ( CONTAINSp0 .AND. CONTAINSp0 ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO iz = 1,local_nz
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ Tpar_ = 0._dp
+ DO ij = 1, local_nj
+ j_dp = REAL(ij-1,dp)
+ Tpar_ = Tpar_ + kernel(ia,ij+ngj/2,iky,ikx,iz+ngz/2,ieven)*&
+ (SQRT2 * moments(ia,ip2,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel) &
+ + moments(ia,ip0,ij+ngj/2,iky,ikx,iz+ngz/2,updatetlevel))
+ ENDDO
+ Tpar(ia,iky,ikx,iz) = Tpar_
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
ENDDO
END SUBROUTINE compute_Tpar
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index bf1d7fee97fac1503f2c71f98cb9fd9379fca109..4a574d672d42480d6b0854acb7aa7b3a84013b43 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -25,6 +25,7 @@ CONTAINS
INTEGER :: deltap_cp
INTEGER :: pmax_cp, jmax_cp, n0, Nkx_cp, Nky_cp, Nz_cp, Na_cp, Np_cp, Nj_cp
INTEGER :: ia,ip,ij,iky,ikx,iz, iacp,ipcp,ijcp,iycp,ixcp,izcp, ierr
+ INTEGER :: ipi,iji,izi
REAL(dp):: timer_tot_1,timer_tot_2
COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_cp
CALL cpu_time(timer_tot_1)
@@ -47,7 +48,7 @@ CONTAINS
CALL getatt(fidrst,"/data/input/grid" , "jmax", jmax_cp)
Nj_cp = jmax_cp+1
CALL getatt(fidrst,"/data/input/model", "Na", Na_cp)
- CALL getatt(fidrst,"/data/input/basic" , "startframe5d", n0)
+ CALL getatt(fidrst,"/data/input/basic" , "start_iframe5d", n0)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments", n_ ! start with moments/000001
@@ -73,19 +74,22 @@ CONTAINS
moments = 0._dp;
z: DO iz = 1,local_nz
- izcp = iz+local_nz_offset
+ izcp = iz + local_nz_offset
+ izi = iz + ngz/2
x: DO ikx=1,local_nkx
ixcp = ikx+local_nkx_offset
y: DO iky=1,local_nky
- iycp = iky + local_nkx_offset
+ iycp = iky + local_nky_offset
j: DO ij=1,local_nj
ijcp = ij + local_nj_offset
+ iji = ij + ngj/2
p: DO ip=1,local_np
ipcp = ip + local_np_offset
+ ipi = ip + ngp/2
a: DO ia=1,Na_cp
iacp = ia + local_na_offset
- IF((iacp.LE.Na_cp).AND.(ipcp.LE.Np_cp).AND.(ijcp.LE.Nj_cp).AND.(iycp.LE.Nky_cp).AND.(ixcp.LE.Nkx_cp).AND.(izcp.LE.Nz_cp)) &
- moments(ia,ip,ij,iky,ikx,iz,1) = moments_cp(ia,ip,ij,iky,ikx,iz)
+ ! IF((iacp.LE.Na_cp).AND.(ipcp.LE.Np_cp).AND.(ijcp.LE.Nj_cp).AND.(iycp.LE.Nky_cp).AND.(ixcp.LE.Nkx_cp).AND.(izcp.LE.Nz_cp)) &
+ moments(ia,ipi,iji,iky,ikx,izi,1) = moments_cp(iacp,ipcp,ijcp,iycp,ixcp,izcp)
ENDDO a
ENDDO p
ENDDO j
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index a365937d88e511f40d33a9a03b0c9a17c812016c..73464ae2a432bb2971be251e9f78446f3fbd2d42 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -31,8 +31,6 @@ CONTAINS
INTEGER :: in, ia, ikx, iky, iz, izi, ini
COMPLEX(dp) :: fsa_phi, intf_, rhtot ! current flux averaged phi
COMPLEX(dp), DIMENSION(local_nz) :: rho, integrant ! charge density q_a n_a and aux var
- ! Execution time start
- CALL cpu_time(t0_poisson)
rhtot = 0
!! Poisson can be solved only for process containng p=0
IF ( SOLVE_POISSON ) THEN
@@ -60,11 +58,14 @@ CONTAINS
fsa_phi = 0._dp
IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
! Prepare integrant for z-average
- integrant(iz) = Jacobian(iz+ngz/2,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz)
+ DO iz = 1,local_nz
+ izi = iz+ngz/2
+ integrant(iz) = Jacobian(izi,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz)
+ ENDDO
call simpson_rule_z(local_nz,zweights_SR,integrant,intf_) ! get the flux averaged phi
fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF
- rho(iz) = rho(iz) + fsa_phi
+ rho = rho + fsa_phi
ENDIF
!!!!!!!!!!!!!!! adiabatic ions ?
! IF (ADIAB_I) THEN
@@ -81,9 +82,6 @@ CONTAINS
ENDIF
! Transfer phi to all the others process along p
CALL manual_3D_bcast(phi,local_nky,local_nkx,local_nz+ngz)
- ! Execution time end
- CALL cpu_time(t1_poisson)
- tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
! print*, SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
! print*, SUM(REAL(moments(1,:,:,:,:,:,updatetlevel)))
! print*, SUM(IMAG(moments(1,:,:,:,:,:,updatetlevel)))
@@ -109,8 +107,6 @@ CONTAINS
IMPLICIT NONE
COMPLEX(dp) :: j_a ! current density
INTEGER :: in, ia, ikx, iky, iz, ini, izi
- ! Execution time start
- CALL cpu_time(t0_poisson)
!! Ampere can be solved only with beta > 0 and for process containng p=1 moments
IF ( SOLVE_AMPERE ) THEN
z:DO iz = 1,local_nz
@@ -136,9 +132,6 @@ CONTAINS
IF (contains_kx0 .AND. contains_ky0) psi(iky0,ikx0,:) = 0._dp
! Transfer phi to all the others process along p
CALL manual_3D_bcast(psi,local_nky,local_nkx,local_nz+ngz)
- ! Execution time end
- CALL cpu_time(t1_poisson)
- tc_poisson = tc_poisson + (t1_poisson - t0_poisson)
END SUBROUTINE ampere
END SUBROUTINE solve_EM_fields
diff --git a/src/stepon.F90 b/src/stepon.F90
index 1001468b27cc89e97d05f304391f5d1805b9a45c..e145ecd61f2a1f5d23b636316c91a58099382d69 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -1,19 +1,20 @@
SUBROUTINE stepon
- ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
- USE advance_field_routine, ONLY: advance_time_level, advance_moments
- USE basic, ONLY: nlend
- USE closure, ONLY: apply_closure_model
- USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
- use mpi, ONLY: MPI_COMM_WORLD
- USE time_integration, ONLY: ntimelevel
- USE prec_const, ONLY: dp
- IMPLICIT NONE
-
- INTEGER :: num_step, ierr
- LOGICAL :: mlend
+ ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
+ USE advance_field_routine, ONLY: advance_time_level, advance_moments
+ USE basic, ONLY: nlend, chrono_advf, chrono_pois,&
+ chrono_chck, chrono_clos, chrono_ghst, start_chrono, stop_chrono
+ USE closure, ONLY: apply_closure_model
+ USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
+ use mpi, ONLY: MPI_COMM_WORLD
+ USE time_integration, ONLY: ntimelevel
+ USE prec_const, ONLY: dp
+ IMPLICIT NONE
+
+ INTEGER :: num_step, ierr
+ LOGICAL :: mlend
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
- !----- BEFORE: All fields+ghosts are updated for step = n
+ !----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
CALL assemble_RHS
@@ -25,134 +26,156 @@ SUBROUTINE stepon
! Update moments with the hierarchy RHS (step by step)
! N_n+1 = N_n + N_rhs(n)
- CALL advance_moments
+ CALL start_chrono(chrono_advf)
+ CALL advance_moments
+ CALL stop_chrono(chrono_advf)
+
! Closure enforcement of moments
- CALL apply_closure_model
+ CALL start_chrono(chrono_clos)
+ CALL apply_closure_model
+ CALL stop_chrono(chrono_clos)
+
! Exchanges the ghosts values of N_n+1
- CALL update_ghosts_moments
+ CALL start_chrono(chrono_ghst)
+ CALL update_ghosts_moments
+ CALL stop_chrono(chrono_ghst)
! Update electrostatic potential phi_n = phi(N_n+1) and potential vect psi
- CALL solve_EM_fields
- CALL update_ghosts_EM
+ CALL start_chrono(chrono_pois)
+ CALL solve_EM_fields
+ CALL stop_chrono(chrono_pois)
+ ! Update EM ghosts
+ CALL start_chrono(chrono_ghst)
+ CALL update_ghosts_EM
+ CALL stop_chrono(chrono_ghst)
!- Check before next step
- CALL checkfield_all()
+ CALL start_chrono(chrono_chck)
+ CALL checkfield_all()
+ CALL stop_chrono(chrono_chck)
+
IF( nlend ) EXIT ! exit do loop
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
- !----- AFTER: All fields are updated for step = n+1
+ !----- AFTER: All fields are updated for step = n+1
END DO
- CONTAINS
- !!!! Basic structure to simplify stepon
- SUBROUTINE assemble_RHS
- USE moments_eq_rhs, ONLY: compute_moments_eq_rhs
- USE collision, ONLY: compute_Capj
- USE nonlinear, ONLY: compute_Sapj
- USE processing, ONLY: compute_nadiab_moments_z_gradients_and_interp
- IMPLICIT NONE
- ! compute auxiliary non adiabatic moments array, gradients and interp
- CALL compute_nadiab_moments_z_gradients_and_interp
- ! compute nonlinear term ("if linear" is included inside)
- CALL compute_Sapj
- ! compute collision term ("if coll, if nu >0" is included inside)
- CALL compute_Capj
- ! compute the moments equation rhs
- CALL compute_moments_eq_rhs
- END SUBROUTINE assemble_RHS
-
- SUBROUTINE checkfield_all ! Check all the fields for inf or nan
- USE utility,ONLY: is_nan, is_inf
- USE basic, ONLY: t0_checkfield, t1_checkfield, tc_checkfield
- USE fields, ONLY: phi
- USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
- ngp,ngj,ngz
- USE MPI
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: LINEARITY
- IMPLICIT NONE
- LOGICAL :: checkf_
- INTEGER :: ia, ip, ij, iky, ikx, iz
- REAL :: sum_
- ! Execution time start
- CALL cpu_time(t0_checkfield)
-
- IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule
- IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0
-
- mlend=.FALSE.
- IF(.NOT.nlend) THEN
- sum_ = SUM(ABS(phi))
- checkf_ = (is_nan(sum_,'phi') .OR. is_inf(sum_,'phi'))
- mlend = (mlend .or. checkf_)
- CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
- ENDIF
- ! Execution time end
- CALL cpu_time(t1_checkfield)
- tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield)
- END SUBROUTINE checkfield_all
-
- SUBROUTINE anti_aliasing
- USE fields, ONLY: moments
- USE time_integration, ONLY: updatetlevel
- USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
- ngp,ngj,ngz, AA_x, AA_y
- IMPLICIT NONE
- INTEGER :: ia, ip, ij, iky, ikx, iz
- z: DO iz = 1,local_nz+ngz
- kx:DO ikx= 1,local_nkx
- ky:DO iky=1,local_nky
- j: DO ij=1,local_nj+ngj
- p: DO ip=1,local_np+ngp
- a: DO ia=1,local_na
- moments(ia,ip,ij,iky,ikx,iz,updatetlevel) =&
- AA_x(ikx)* AA_y(iky) * moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
- ENDDO a
- ENDDO p
- ENDDO j
- ENDDO ky
- ENDDO kx
- ENDDO z
- END SUBROUTINE anti_aliasing
-
- SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
- USE fields, ONLY: phi, psi, moments
- USE time_integration, ONLY: updatetlevel
- USE grid, ONLY: local_na,local_np,local_nj,total_nkx,local_nz,&
- ngp,ngj,ngz, ikx0,iky0, contains_ky0
- IMPLICIT NONE
- INTEGER :: ia, ip, ij, ikx, iz
- IF ( contains_ky0 ) THEN
- ! moments
- ! z: DO iz = 1,local_nz+ngz
- ! j: DO ij=1,local_nj+ngj
- ! p: DO ip=1,local_np+ngp
- ! a: DO ia=1,local_na
- DO ikx=2,total_nkx/2 !symmetry at ky = 0
- moments(:,:,:,iky0,ikx,:,updatetlevel) = &
- CONJG(moments(:,:,:,iky0,total_nkx+2-ikx,:,updatetlevel))
- END DO
- ! must be real at origin and top right
- moments(:,:,:, iky0,ikx0,:,updatetlevel) = &
- REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),dp)
- ! ENDDO a
- ! ENDDO p
- ! ENDDO j
- ! ENDDO z
- ! Phi
- DO ikx=2,total_nkx/2 !symmetry at ky = 0
+CONTAINS
+!!!! Basic structure to simplify stepon
+ SUBROUTINE assemble_RHS
+ USE basic, ONLY: chrono_mrhs, chrono_sapj, chrono_coll, chrono_grad, chrono_napj, start_chrono, stop_chrono
+ USE moments_eq_rhs, ONLY: compute_moments_eq_rhs
+ USE collision, ONLY: compute_Capj
+ USE nonlinear, ONLY: compute_Sapj
+ USE processing, ONLY: compute_nadiab_moments, compute_interp_z, compute_gradients_z
+ IMPLICIT NONE
+ ! compute auxiliary non adiabatic moments array
+ CALL start_chrono(chrono_napj)
+ CALL compute_nadiab_moments
+ CALL stop_chrono(chrono_napj)
+
+ ! compute gradients and interp
+ CALL start_chrono(chrono_grad)
+ CALL compute_gradients_z
+ CALL compute_interp_z
+ CALL stop_chrono(chrono_grad)
+
+ ! compute nonlinear term ("if linear" is included inside)
+ CALL start_chrono(chrono_sapj)
+ CALL compute_Sapj
+ CALL stop_chrono(chrono_sapj)
+
+ ! compute collision term ("if coll, if nu >0" is included inside)
+ CALL start_chrono(chrono_coll)
+ CALL compute_Capj
+ CALL stop_chrono(chrono_coll)
+
+ ! compute the moments equation rhs
+ CALL start_chrono(chrono_mrhs)
+ CALL compute_moments_eq_rhs
+ CALL stop_chrono(chrono_mrhs)
+ END SUBROUTINE assemble_RHS
+
+ SUBROUTINE checkfield_all ! Check all the fields for inf or nan
+ USE utility,ONLY: is_nan, is_inf
+ USE fields, ONLY: phi
+ USE MPI
+ USE model, ONLY: LINEARITY
+ IMPLICIT NONE
+ LOGICAL :: checkf_
+ REAL :: sum_
+ IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule
+ IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0
+
+ mlend=.FALSE.
+ IF(.NOT.nlend) THEN
+ sum_ = SUM(ABS(phi))
+ checkf_ = (is_nan(sum_,'phi') .OR. is_inf(sum_,'phi'))
+ mlend = (mlend .or. checkf_)
+ CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
+ ENDIF
+ END SUBROUTINE checkfield_all
+
+ SUBROUTINE anti_aliasing
+ USE fields, ONLY: moments
+ USE time_integration, ONLY: updatetlevel
+ USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
+ ngp,ngj,ngz, AA_x, AA_y
+ IMPLICIT NONE
+ INTEGER :: ia, ip, ij, iky, ikx, iz
+ z: DO iz = 1,local_nz+ngz
+ kx:DO ikx= 1,local_nkx
+ ky:DO iky=1,local_nky
+ j: DO ij=1,local_nj+ngj
+ p: DO ip=1,local_np+ngp
+ a: DO ia=1,local_na
+ moments(ia,ip,ij,iky,ikx,iz,updatetlevel) =&
+ AA_x(ikx)* AA_y(iky) * moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
+ ENDDO a
+ ENDDO p
+ ENDDO j
+ ENDDO ky
+ ENDDO kx
+ ENDDO z
+ END SUBROUTINE anti_aliasing
+
+ SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
+ USE fields, ONLY: phi, psi, moments
+ USE time_integration, ONLY: updatetlevel
+ USE grid, ONLY: total_nkx, ikx0,iky0, contains_ky0
+ IMPLICIT NONE
+ INTEGER :: ikx
+ IF ( contains_ky0 ) THEN
+ ! moments
+ ! z: DO iz = 1,local_nz+ngz
+ ! j: DO ij=1,local_nj+ngj
+ ! p: DO ip=1,local_np+ngp
+ ! a: DO ia=1,local_na
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0
+ moments(:,:,:,iky0,ikx,:,updatetlevel) = &
+ CONJG(moments(:,:,:,iky0,total_nkx+2-ikx,:,updatetlevel))
+ END DO
+ ! must be real at origin and top right
+ moments(:,:,:, iky0,ikx0,:,updatetlevel) = &
+ REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),dp)
+ ! ENDDO a
+ ! ENDDO p
+ ! ENDDO j
+ ! ENDDO z
+ ! Phi
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0
phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:)
- END DO
- ! must be real at origin
- phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
- ! Psi
- DO ikx=2,total_nkx/2 !symmetry at ky = 0
+ END DO
+ ! must be real at origin
+ phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
+ ! Psi
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0
psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:)
- END DO
- ! must be real at origin
- psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
- ENDIF
- END SUBROUTINE enforce_symmetry
+ END DO
+ ! must be real at origin
+ psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
+ ENDIF
+ END SUBROUTINE enforce_symmetry
END SUBROUTINE stepon
diff --git a/src/tesend.F90 b/src/tesend.F90
index 14d673443ff1aa5428e1c0258427ce32f6e197c6..5c04ab7cd8fa9e24da6f2001572ed265135e90bc 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -44,7 +44,7 @@ SUBROUTINE tesend
!________________________________________________________________________________
! 4. Test on run time
CALL cpu_time(tnow)
- mlend = (1.1*(tnow-start)) .GT. maxruntime
+ mlend = (1.1*(tnow-chrono_runt%tstart)) .GT. maxruntime
CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
diff --git a/testcases/cyclone_example/fort.90 b/testcases/cyclone_example/fort.90
deleted file mode 100644
index badc2ce7bff982da60b975610e7c6273e3dd50b3..0000000000000000000000000000000000000000
--- a/testcases/cyclone_example/fort.90
+++ /dev/null
@@ -1,92 +0,0 @@
-&BASIC
- nrun = 99999999
- dt = 0.01
- tmax = 25
- maxruntime = 356400
- job2load = -1
-/
-&GRID
- pmax = 4
- jmax = 2
- Nx = 128
- Lx = 120
- Ny = 48
- Ly = 120
- Nz = 16
- SG = .f.
- Nexc = 0
-/
-&GEOMETRY
- geom = 's-alpha'
- q0 = 1.4
- shear = 0.8
- eps = 0.18
- kappa = 1.0
- s_kappa= 0.0
- delta = 0.0
- s_delta= 0.0
- zeta = 0.0
- s_zeta = 0.0
- parallel_bc = 'dirichlet'
- shift_y= 0.0
-/
-&OUTPUT_PAR
- dtsave_0d = 0.1
- dtsave_1d = -1
- dtsave_2d = -1
- dtsave_3d = 1
- dtsave_5d = 10
- write_doubleprecision = .f.
- write_gamma = .t.
- write_hf = .t.
- write_phi = .t.
- write_Na00 = .t.
- write_Napj = .t.
- write_dens = .t.
- write_fvel = .t.
- write_temp = .t.
-/
-&MODEL_PAR
- ! Collisionality
- CLOS = 0
- NL_CLOS = 0
- LINEARITY = 'nonlinear'
- Na = 1 ! number of species
- mu_x = 1.0
- mu_y = 1.0
- N_HD = 4
- mu_z = 2.0
- HYP_V = 'hypcoll'
- mu_p = 0.0
- mu_j = 0.0
- nu = 0.001
- beta = 0.0
- ADIAB_E = .t.
- tau_e = 1.0
-/
-&SPECIES
- ! ions
- name_ = 'ions'
- tau_ = 1.0
- sigma_= 1.0
- q_ = 1.0
- k_N_ = 2.22
- k_T_ = 6.96
-/
-
-&COLLISION_PAR
- collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
- GK_CO = .f.
- INTERSPECIES = .true.
- mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
-/
-&INITIAL_CON
- INIT_OPT = 'blob'
- ACT_ON_MODES = 'donothing'
- init_background = 0.0
- init_noiselvl = 0.005
- iseed = 42
-/
-&TIME_INTEGRATION_PAR
- numerical_scheme = 'RK4'
-/
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index f0e73003d14d54015a59009e6d50063b046eee53..cf0ff061e658b8cd9d1ec36d94c90b2f9982aaeb 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -1,8 +1,8 @@
%% UNCOMMENT FOR TUTORIAL
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
-resdir = '../testcases/zpinch_example/'; %Name of the directory where the results are located
+resdir = '../testcases/cyclone_example/'; %Name of the directory where the results are located
PARTITION ='';
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+JOBNUMMIN = 03; JOBNUMMAX = 03;
%%
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
@@ -80,7 +80,7 @@ options.NAME = '\phi';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -93,27 +93,28 @@ options.RESOLUTION = 256;
create_film(data,options,'.gif')
end
-if 0
+if 1
%% fields snapshots
% Options
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
% options.NAME = 'T_i';
% options.NAME = 'n_i';
% options.NAME = '\phi^{NZ}';
-% options.NAME = 'N_i^{00}';
+options.NAME = 'N_i^{00}';
% options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
options.COMP = 'avg';
-options.TIME = [15 30 50];
+% options.TIME = [49 50 51];
+options.TIME = data.Ts3D(51:55);
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;