From f994633ea802ce3217155bdf8318142b025d9975 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 3 Jul 2020 18:23:15 +0200
Subject: [PATCH] Full coulomb COSOlver files can be read now
---
matlab/MOLI_time_solver.m | 16 +-
matlab/write_fort90.m | 3 +
src/fourier_grid_mod.F90 | 20 ++
src/inital.F90 | 55 ++++-
src/initial_par_mod.F90 | 7 +
src/memory.F90 | 11 +-
src/moments_eq_rhs.F90 | 412 +++++++++++++++++++++----------------
wk/benchmark_time_solver.m | 294 ++++++++++++++++++++++++++
8 files changed, 628 insertions(+), 190 deletions(-)
create mode 100644 wk/benchmark_time_solver.m
diff --git a/matlab/MOLI_time_solver.m b/matlab/MOLI_time_solver.m
index ff15c259..0bd7d118 100644
--- a/matlab/MOLI_time_solver.m
+++ b/matlab/MOLI_time_solver.m
@@ -45,7 +45,7 @@ options.electrons = 1; % 0 -> adiabatic electrons, 1 no adiabatic electron
options.sw = 1; % 1 -> sound waves on, 0 -> put ion parallel velocity row/column to 0
%% Collision Operator Models and COSOlver Input Parameters
-options.collI = -2; % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
+options.collI = MODEL.CO; % collI=-2 -> Dougherty, -1 -> COSOlver, 0 -> Lenard-Bernstein, other -> hyperviscosity exponent
options.collGK = 0; % collDKGK =1 -> GK collision operator, else DK collision operator
options.COSOlver.GKE = 0;
options.COSOlver.GKI = 0;
@@ -54,7 +54,7 @@ options.COSOlver.GKI = 0;
options.COSOlver.eecolls = 1; % 1 -> electron-electron collisions, 0 -> off
options.COSOlver.iicolls = 1; % 1 -> ion-ion collisions, 0 -> off
options.COSOlver.eicolls = 1; % 1 -> electron-ion collisions (e-i) on, 0 -> off
-options.COSOlver.iecolls = 0; % 1 -> ion-electron collisions (i-e) on, 0 -> off
+options.COSOlver.iecolls = 1; % 1 -> ion-electron collisions (i-e) on, 0 -> off
% Collisional Coulomb sum bounds (only if collI = -1, i.e. Coulomb)
options.COSOlver.lmaxx = 10; % upper bound collision operator first sum first species
@@ -70,12 +70,12 @@ options.COSOlver.nimaxxFLR = 0; % upper bound FLR ion collision
% Set electron/ion test and back-reaction model operator
%
% 0 => Coulomb Collisions
-options.COSOlver.ETEST = 4; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
-options.COSOlver.EBACK = 0;
-options.COSOlver.ITEST = 4;
-options.COSOlver.IBACK = 0;
-options.COSOlver.ESELF = 4;
-options.COSOlver.ISELF = 4;
+options.COSOlver.ETEST = 1; % 0 --> Buffer Operator, 1 --> Coulomb, 2 --> Lorentz
+options.COSOlver.EBACK = 1;
+options.COSOlver.ITEST = 1;
+options.COSOlver.IBACK = 1;
+options.COSOlver.ESELF = 1;
+options.COSOlver.ISELF = 1;
options.COSOlver.OVERWRITE = 0; % overwrite collisional matrices even if exists
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index e7b820ef..6eeb6643 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -51,6 +51,9 @@ fprintf(fid,[' initback_moments=', num2str(INITIAL.initback_moments),' ! x 1e-3
fprintf(fid,' ! Noise amplitude\n');
fprintf(fid,[' initnoise_moments=', num2str(INITIAL.initnoise_moments),'\n']);
fprintf(fid,[' iseed=', num2str(INITIAL.iseed),'\n']);
+fprintf(fid,[' selfmat_file=', INITIAL.selfmat_file,'\n']);
+fprintf(fid,[' eimat_file=', INITIAL.eimat_file,'\n']);
+fprintf(fid,[' iemat_file=', INITIAL.iemat_file,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&TIME_INTEGRATION_PAR\n');
fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90
index c8a40656..b2765095 100644
--- a/src/fourier_grid_mod.F90
+++ b/src/fourier_grid_mod.F90
@@ -20,6 +20,7 @@ MODULE fourier_grid
! Indices of s -> start , e-> END
INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e
INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i
+ INTEGER, PUBLIC, PROTECTED :: ns_e,ne_e, ns_i,ne_i !start/end indices for coll. mat.
INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze
! Toroidal direction
@@ -39,6 +40,7 @@ MODULE fourier_grid
! Public Functions
PUBLIC :: set_krgrid, set_kzgrid, set_pj
PUBLIC :: fourier_grid_readinputs, fourier_grid_outputinputs
+ PUBLIC :: bare, bari
CONTAINS
@@ -91,6 +93,10 @@ CONTAINS
ALLOCATE(jarray_i(ijs_i:ije_i))
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
+
+ ns_e = bare(ips_e,ijs_e); ne_e = bare(ipe_e,ije_e)
+ ns_i = bari(ips_i,ijs_i); ne_i = bari(ipe_i,ije_i)
+
END SUBROUTINE set_pj
SUBROUTINE fourier_grid_readinputs
@@ -132,4 +138,18 @@ CONTAINS
CALL attach(fidres, TRIM(str), "kzmax", kzmax)
END SUBROUTINE fourier_grid_outputinputs
+
+ FUNCTION bare(ip,ij)
+ IMPLICIT NONE
+ INTEGER :: bare, ip, ij
+ bare = (jmaxe+1)*ip + ij + 1
+ END FUNCTION
+
+ FUNCTION bari(ip,ij)
+ IMPLICIT NONE
+ INTEGER :: bari, ip, ij
+ bari = bare(pmaxe,jmaxe) + (jmaxi+1)*ip + ij + 1
+ END FUNCTION
+
+
END MODULE fourier_grid
diff --git a/src/inital.F90 b/src/inital.F90
index c60072f6..1bcebfc2 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -27,7 +27,6 @@ SUBROUTINE init_profiles
USE fields
USE initial_par
USE time_integration
-
USE prec_const
IMPLICIT NONE
@@ -47,21 +46,21 @@ SUBROUTINE init_profiles
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
CALL RANDOM_NUMBER(noise)
- !moments_e( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_e( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
END DO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
CALL RANDOM_NUMBER(noise)
- !moments_i( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_i( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
END DO
! Poke initialization on only Ne00 and Ni00
- moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
-
+ !moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ !moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+
END DO
END DO
@@ -69,9 +68,51 @@ SUBROUTINE init_profiles
END SUBROUTINE init_profiles
-SUBROUTINE load_FC_mat
+SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
USE basic
+ USE fourier_grid
+ USE initial_par
USE array
+ use model
+ USE futils, ONLY : openf, getarr, closef
IMPLICIT NONE
+ INTEGER :: fid, ip,ij, ip2,ij2
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ WRITE(*,*) 'Load self electron collision matrix..'
+ ! get the self electron colision matrix
+ CALL openf(selfmat_file,fid, 'r', 'D');
+ CALL getarr(fid, '/Caapj/Ceepj', Ceepj) ! get array (moli format)
+ CALL closef(fid)
+
+ ! get the Test and Back field electron ion collision matrix
+ WRITE(*,*) 'Load test + field electron collision matrix..'
+ CALL openf(eimat_file,fid, 'r', 'D');
+ CALL getarr(fid, '/Ceipj/CeipjT', CeipjT)
+ CALL getarr(fid, '/Ceipj/CeipjF', CeipjF)
+ CALL closef(fid)
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ WRITE(*,*) 'Load self ion collision matrix..'
+ ! get the self electron colision matrix
+ CALL openf(selfmat_file, fid, 'r', 'D');
+ IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
+ CALL getarr(fid, '/Caapj/Ceepj', Ciipj) ! get array (moli format)
+ ELSE
+ CALL getarr(fid, '/Caapj/Ciipj', Ciipj) ! get array (moli format)
+ ENDIF
+ CALL closef(fid)
+ ! get the Test and Back field ion electron collision matrix
+ WRITE(*,*) 'Load test + field ion collision matrix..'
+ CALL openf(iemat_file, fid, 'r', 'D');
+ CALL getarr(fid, '/Ciepj/CiepjT', CiepjT)
+ CALL getarr(fid, '/Ciepj/CiepjF', CiepjF)
+ CALL closef(fid)
+
+
END SUBROUTINE load_FC_mat
\ No newline at end of file
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index 64953baa..61dbfdb7 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -20,6 +20,10 @@ MODULE initial_par
REAL(dp), PUBLIC, PROTECTED :: init_ampli_density=0.1_dp ! Oscillation amplitude
REAL(dp), PUBLIC, PROTECTED :: init_ampli_temp=0.1_dp
+ CHARACTER(len=128), PUBLIC :: selfmat_file ! COSOlver matrix file names
+ CHARACTER(len=128), PUBLIC :: iemat_file ! COSOlver matrix file names
+ CHARACTER(len=128), PUBLIC :: eimat_file ! COSOlver matrix file names
+
PUBLIC :: initial_outputinputs, initial_readinputs
CONTAINS
@@ -35,6 +39,9 @@ CONTAINS
NAMELIST /INITIAL_CON/ initback_moments
NAMELIST /INITIAL_CON/ initnoise_moments
NAMELIST /INITIAL_CON/ iseed
+ NAMELIST /INITIAL_CON/ selfmat_file
+ NAMELIST /INITIAL_CON/ iemat_file
+ NAMELIST /INITIAL_CON/ eimat_file
READ(lu_in,initial_con)
WRITE(*,initial_con)
diff --git a/src/memory.F90 b/src/memory.F90
index 410abfdc..a094ce92 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -22,7 +22,14 @@ SUBROUTINE memory
! Collision matrix
IF (CO .EQ. -1) THEN
- CALL allocate_array( iCe, ips_i,ipe_i, ijs_i,ije_i )
- CALL allocate_array( iCi, ips_i,ipe_i, ijs_i,ije_i )
+ CALL allocate_array( Ceepj, ns_e,ne_e, 1,(pmaxe+1)*(jmaxe+1))
+ CALL allocate_array( CeipjT, ns_e,ne_e, 1,(pmaxe+1)*(jmaxe+1))
+ CALL allocate_array( CeipjF, ns_e,ne_e, 1,(pmaxi+1)*(jmaxi+1))
+
+ CALL allocate_array( Ciipj, ns_i,ne_i, 1,(pmaxi+1)*(jmaxi+1))
+ CALL allocate_array( CiepjT, ns_i,ne_i, 1,(pmaxi+1)*(jmaxi+1))
+ CALL allocate_array( CiepjF, ns_i,ne_i, 1,(pmaxe+1)*(jmaxe+1))
+
+ write(*,*) 'ns_i=',ns_i,', ns_e=',ns_e
ENDIF
END SUBROUTINE memory
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index e4b95996..e2533b66 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -9,80 +9,89 @@ SUBROUTINE moments_eq_rhs
USE prec_const
IMPLICIT NONE
- INTEGER :: ip,ij, ikr,ikz ! loops indices
+ INTEGER :: ip,ij, ikr,ikz, ip2, ij2 ! loops indices
REAL(dp) :: ip_dp, ij_dp
REAL(dp) :: kr, kz, kperp2
REAL(dp) :: taue_qe_etaB, taui_qi_etaB
REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
- REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
- REAL(dp) :: xphij, xphijp1, xphijm1, xColl
+ REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop
+ REAL(dp) :: xphij, xphijp1, xphijm1 ! ESpot. factors depending on the pj loop
+ REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
+ REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
+!write(*,*) '----------------------------------------'
!Precompute species dependant factors
taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
taui_qi_etaB = tau_i/q_i * eta_B
- sigmae2_taue_o2 = sigma_e**2 * tau_e/2.0 ! factor of the Kernel argument
- sigmai2_taui_o2 = sigma_i**2 * tau_i/2.0
-
+ 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
+ nu_e = nu ! electron-ion collision frequency (where already multiplied by 0.53..)
+ nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ.
+ nu_ee = nu_e/SQRT2 ! e-e coll. frequ.
+ nu_ie = nu*sigma_e**2 ! i-e coll. frequ.
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! 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 pmax)
-
- ! N_e^{p+2,j} multiplicator
- IF (ip+2 .LE. pmaxe+1) THEN
- xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
- ELSE
- xNapp2j = 0.
- ENDIF
-
- ! N_e^{p-2,j} multiplicator
- IF (ip-2 .GE. 1) THEN
- xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
- ELSE
- xNapm2j = 0.
- ENDIF
+ ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxe)
+
+ ! N_e^{p+2,j} coeff
+ xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ ! N_e^{p-2,j} coeff
+ xNapm2j = -taue_qe_etaB * SQRT(ip_dp * (ip_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 jmax)
-
- IF (ij_dp .GT. 0) THEN
- factj = factj * ij_dp; ! Recursive factorial
+ ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxe)
+
+ ! N_e^{p,j+1} coeff
+ xNapjp1 = taue_qe_etaB * (ij_dp + 1._dp)
+ ! N_e^{p,j-1} coeff
+ xNapjm1 = taue_qe_etaB * ij_dp
+ ! N_e^{pj} coeff
+ xNapj = -taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+
+ !! Collision operator pj terms
+ xCapj = -nu_e*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
+ ! Dougherty part
+ IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) 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
+ 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
+ xCa20 = 0._dp
+ xCa01 = 0._dp
+ xCa10 = nu_e
+ ELSE
+ xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp
ENDIF
- ! N_e^{p,j+1} multiplicator
- IF (ij+1 .LE. jmaxe+1) THEN
- xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
- ELSE
- xNapjp1 = 0.
- ENDIF
-
- ! N_e^{p,j-1} multiplicator
- IF (ij-1 .GE. 1) THEN
- xNapjm1 = taue_qe_etaB * ij_dp
+ !! 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
+ ELSE IF (ip .EQ. 3) THEN ! kronecker p2
+ xphij = (eta_T/SQRT2 - SQRT2*eta_B)
+ xphijp1 = 0._dp; xphijm1 = 0._dp
ELSE
- xNapjm1 = 0.
+ xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
ENDIF
- ! N_e^{pj} multiplicator
- xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
-
- ! Collision operator (DK Lenard-Bernstein basis)
- xColl = ip_dp + 2.*ij_dp
-
- ! phi multiplicator for different kernel numbers
- IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
- xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
- xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
- xphijm1 = -(eta_T - eta_B)* ij_dp
- ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
- xphij = (eta_T/SQRT2 - SQRT2*eta_B)
- xphijp1 = 0.; xphijm1 = 0.
+ ! Recursive factorial
+ IF (ij_dp .GT. 0) THEN
+ factj = factj * ij_dp
ELSE
- xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
+ factj = 1._dp
ENDIF
! Loop on kspace
@@ -95,128 +104,157 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_e^{p,j}
- TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
+ TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p+2,j}
- IF (ip+2 .LE. pmaxe+1) THEN
- TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
+ IF (ip+2 .LE. pmaxe+1) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
- TNapp2j = 0.
+ TNapp2j = 0._dp
ENDIF
! term propto N_e^{p-2,j}
- IF (ip-2 .GE. 1) THEN
- TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
+ IF (ip-2 .GE. 1) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
- TNapm2j = 0.
+ TNapm2j = 0._dp
ENDIF
! xterm propto N_e^{p,j+1}
- IF (ij+1 .LE. jmaxe+1) THEN
- TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
+ IF (ij+1 .LE. jmaxe+1) THEN ! OoB check
+ TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
- TNapjp1 = 0.
+ TNapjp1 = 0._dp
ENDIF
! term propto N_e^{p,j-1}
- IF (ij-1 .GE. 1) THEN
- TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
+ IF (ij-1 .GE. 1) THEN ! OoB check
+ TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
- TNapjm1 = 0.
+ TNapjm1 = 0._dp
ENDIF
- ! Dougherty Collision terms
- IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxe .GE. 2) ) THEN ! kronecker p0 * j1
- TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
- TColl20 = 0.0; TColl10 = 0.0;
- ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxe .GE. 1)) THEN ! kronecker p2 * j0
- TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
- TColl10 = 0.0; TColl01 = 0.0;
- ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
- TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
- TColl20 = 0.0; TColl01 = 0.0;
- ELSE
- TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
+ !! 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)
+ ELSE
+ TColl20 = 0._dp
+ ENDIF
+ IF ( (jmaxe .GE. 1) ) THEN ! OoB check
+ TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
+ ELSE
+ TColl01 = 0._dp
+ ENDIF
+ IF ( (pmaxe .GE. 1) ) THEN ! OoB check
+ TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
+ ELSE
+ TColl10 = 0._dp
+ ENDIF
+
+ ! Total collisional term
+ TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
+ + TColl20 + TColl01 + TColl10
+
+ ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb (COSOlver matrix) !!!
+
+ TColl = 0._dp ! Initialization
+
+ ploopee: DO ip2 = 1,pmaxe ! sum the electron-self and electron-ion test terms
+ jloopee: DO ij2 = 1,jmaxe
+ 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)))
+ ENDDO jloopee
+ ENDDO ploopee
+
+ ploopei: DO ip2 = 1,pmaxi ! sum the electron-ion field terms
+ jloopei: DO ij2 = 1,jmaxi
+ TColl = TColl - moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
+ *(nu_e * CeipjF(bare(ip-1,ij-1), bari(ip2-1,ij2-1)))
+ END DO jloopei
+ ENDDO ploopei
+
+ ELSE ! Lenhard Bernstein
+ TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
ENDIF
- ! Total collisional term
- TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
- + TColl01 + TColl10 + TColl20)
!! Electrical potential term
- IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker delta_p^0, delta_p^2
+ IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2
kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
- kerneljp1 = kernelj * b_e2 /(ij_dp + 1.)
+ kerneljp1 = kernelj * b_e2 /(ij_dp + 1._dp)
kerneljm1 = kernelj * ij_dp / b_e2
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
- Tphi = 0
+ Tphi = 0._dp
ENDIF
! Sum of all precomputed terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
- imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
-
+ -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi)&
+ + TColl
+
END DO kzloope
END DO krloope
END DO jloope
END DO ploope
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Ions moments RHS !!!!!!!!!
- ploopi : DO ip = ips_i, ipe_i
- ip_dp = REAL(ip-1,dp)
-
- ! x N_i^{p+2,j}
- IF (ip+2 .LE. pmaxi+1) THEN
- xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
- ELSE
- xNapp2j = 0.
- ENDIF
-
- ! x N_i^{p-2,j}
- IF (ip-2 .GE. 1) THEN
- xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
- ELSE
- xNapm2j = 0.
- ENDIF
-
- factj = 1.0 ! Start of the recursive factorial
-
- jloopi : DO ij = ijs_i, ije_i
- ij_dp = REAL(ij-1,dp)
-
- IF (ij_dp .GT. 0) THEN
- factj = factj * ij_dp; ! Recursive factorial
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ 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)
+
+ ! x N_i^{p+2,j} coeff
+ xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ ! x N_i^{p-2,j} coeff
+ xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_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)
+
+ ! x N_i^{p,j+1} coeff
+ xNapjp1 = taui_qi_etaB * (ij_dp + 1._dp)
+ ! x N_i^{p,j-1} coeff
+ xNapjm1 = taui_qi_etaB * ij_dp
+ ! x N_i^{pj} coeff
+ xNapj = -taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+
+ !! Collision operator pj terms
+ xCapj = -nu_i*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
+ ! Dougherty part
+ IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) 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
+ xCa20 = -nu_i * SQRT2 * 2._dp/3._dp
+ xCa01 = -SQRT2 * xCa20
+ xCa10 = 0._dp
+ ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN
+ xCa20 = 0._dp
+ xCa01 = 0._dp
+ xCa10 = nu_i
+ ELSE
+ xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp
ENDIF
- ! x N_i^{p,j+1}
- IF (ij+1 .LE. jmaxi+1) THEN
- xNapjp1 = taui_qi_etaB * (ij_dp + 1.)
- ELSE
- xNapjp1 = 0.
- ENDIF
-
- ! x N_i^{p,j-1}
- IF (ij-1 .GE. 1) THEN
- xNapjm1 = taui_qi_etaB * ij_dp
+ !! 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
+ ELSE IF (ip .EQ. 3) THEN !krokecker p2
+ xphij = (eta_T/SQRT2 - SQRT2*eta_B)
+ xphijp1 = 0._dp; xphijm1 = 0._dp
ELSE
- xNapjm1 = 0.
+ xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
ENDIF
- ! x N_i^{pj}
- xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
-
- ! Collision operator (DK Lenard-Bernstein basis)
- xColl = ip_dp + 2.*ij_dp
-
- ! x phi
- IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
- xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
- xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
- xphijm1 = -(eta_T - eta_B)* ij_dp
- ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
- xphij = (eta_T/SQRT2 - SQRT2*eta_B)
- xphijp1 = 0.; xphijm1 = 0.
+ ! Recursive factorial
+ IF (ij_dp .GT. 0) THEN
+ factj = factj * ij_dp
ELSE
- xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
+ factj = 1._dp
ENDIF
! Loop on kspace
@@ -229,64 +267,92 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_i^{p,j}
- TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
+ TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p+2,j}
- IF (ip+2 .LE. pmaxi+1) THEN
- TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
+ IF (ip+2 .LE. pmaxi+1) THEN ! OoB check
+ TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
- TNapp2j = 0.
+ TNapp2j = 0._dp
ENDIF
! term propto N_i^{p-2,j}
- IF (ip-2 .GE. 1) THEN
- TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
+ IF (ip-2 .GE. 1) THEN ! OoB check
+ TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
- TNapm2j = 0.
+ TNapm2j = 0._dp
ENDIF
! xterm propto N_i^{p,j+1}
- IF (ij+1 .LE. jmaxi+1) THEN
- TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
+ IF (ij+1 .LE. jmaxi+1) THEN ! OoB check
+ TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
- TNapjp1 = 0.
+ TNapjp1 = 0._dp
ENDIF
! term propto N_i^{p,j-1}
- IF (ij-1 .GE. 1) THEN
- TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
+ IF (ij-1 .GE. 1) THEN ! OoB check
+ TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
- TNapjm1 = 0.
+ TNapjm1 = 0._dp
ENDIF
- ! Dougherty Collision terms
- IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxi .GE. 2)) THEN ! kronecker p0 * j1
- TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
- TColl20 = 0.0; TColl10 = 0.0;
- ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxi .GE. 1)) THEN ! kronecker p2 * j0
- TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
- TColl10 = 0.0; TColl01 = 0.0;
- ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
- TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
- TColl20 = 0.0; TColl01 = 0.0;
- ELSE
- TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
- ENDIF
- ! Total collisional term
- TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
- + TColl01 + TColl10 + TColl20)
-
+ !! 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)
+ ELSE
+ TColl20 = 0._dp
+ ENDIF
+ IF ( (jmaxi .GE. 1) ) THEN ! OoB check
+ TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
+ ELSE
+ TColl01 = 0._dp
+ ENDIF
+ IF ( (pmaxi .GE. 1) ) THEN ! OoB check
+ TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
+ ELSE
+ TColl10 = 0._dp
+ ENDIF
+
+ ! Total collisional term
+ TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
+ + TColl20 + TColl01 + TColl10
+
+ ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb (COSOlver matrix) !!!
+
+ TColl = 0._dp ! Initialization
+
+ ploopii: DO ip2 = 1,pmaxi ! sum the electron-self and electron-ion test terms
+ jloopii: DO ij2 = 1,jmaxi
+ 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)))
+ ENDDO jloopii
+ ENDDO ploopii
+
+ ploopie: DO ip2 = 1,pmaxe ! sum the electron-ion field terms
+ jloopie: DO ij2 = 1,jmaxe
+ TColl = TColl - moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
+ *(nu_ie * CiepjF(bari(ip-1,ij-1), bare(ip2-1,ij2-1)))
+ ENDDO jloopie
+ ENDDO ploopie
+
+ ELSE ! Lenhard Bernstein
+ TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
+ ENDIF
+
!! Electrical potential term
- IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! kronecker delta_p^0, delta_p^2
+ IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2
kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
- kerneljp1 = kernelj * b_i2 /(ij_dp + 1.)
+ kerneljp1 = kernelj * b_i2 /(ij_dp + 1._dp)
kerneljm1 = kernelj * ij_dp / b_i2
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
- Tphi = 0
+ Tphi = 0._dp
ENDIF
+
! Sum of all precomputed terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
- imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
-
+ -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi)&
+ + TColl
+
END DO kzloopi
END DO krloopi
diff --git a/wk/benchmark_time_solver.m b/wk/benchmark_time_solver.m
new file mode 100644
index 00000000..a330b5f3
--- /dev/null
+++ b/wk/benchmark_time_solver.m
@@ -0,0 +1,294 @@
+clear all; close all;
+SIMID = 'benchmark_time_solver'; % Name of the simulations
+addpath(genpath('../matlab')) % ... add
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% outputs options
+OUTPUTS.nsave_0d = 0;
+OUTPUTS.nsave_1d = 0;
+OUTPUTS.nsave_2d = 1;
+OUTPUTS.nsave_5d = 1;
+OUTPUTS.write_Ni00 = '.true.';
+OUTPUTS.write_moments = '.true.';
+OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_doubleprecision = '.true.';
+OUTPUTS.resfile0 = '''results''';
+%% Grid parameters
+GRID.pmaxe = 6;
+GRID.jmaxe = 6;
+GRID.pmaxi = 6;
+GRID.jmaxi = 6;
+GRID.nkr = 1;
+GRID.krmin = 0.;
+GRID.krmax = 0.;
+GRID.nkz = 1;
+GRID.kzmin = 1.0;
+GRID.kzmax = 1.0;
+%% Model parameters
+MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.nu = 1.0; % collisionality nu*L_perp/Cs0
+% temperature ratio T_a/T_e
+MODEL.tau_e = 1.0;
+MODEL.tau_i = 1.0;
+% mass ratio sqrt(m_a/m_i)
+MODEL.sigma_e = 0.0233380;
+MODEL.sigma_i = 1.0;
+% charge q_a/e
+MODEL.q_e =-1.0;
+MODEL.q_i = 1.0;
+% gradients L_perp/L_x
+MODEL.eta_n = 0.0; % Density
+MODEL.eta_T = 0.0; % Temperature
+MODEL.eta_B = 0.0; % Magnetic
+% Debye length
+MODEL.lambdaD = 0.0;
+%% Time integration and intialization parameters
+TIME_INTEGRATION.numerical_scheme = '''RK4''';
+BASIC.nrun = 100000;
+BASIC.dt = 0.01;
+BASIC.tmax = 5.0;
+INITIAL.initback_moments = 0.01;
+INITIAL.initnoise_moments = 0.;
+INITIAL.iseed = 42;
+INITIAL.selfmat_file = ...
+ ['''../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(GRID.pmaxe),...
+ '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
+ num2str(GRID.jmaxi),'_pamaxx_10.h5'''];
+INITIAL.eimat_file = ...
+ ['''../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
+ '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
+ num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
+INITIAL.iemat_file = ...
+ ['''../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
+ '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
+ num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Write input file
+INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Run HeLaZ
+nproc = 1;
+EXEC = ' ../bin/helaz ';
+RUN = ['mpirun -np ' num2str(nproc)];
+CMD = [RUN, EXEC, INPUT];
+system(CMD);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Run MOLI with same input parameters
+params.RK4 = 1;
+run ../matlab/MOLI_time_solver
+if params.RK4; MOLIsolvername = 'RK4';
+else; MOLIsolvername = 'ode15i';
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Analysis and basic figures
+default_plots_options
+SAVEFIG = 1;
+filename = 'results_00.h5';
+TITLE = [];
+TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, '];
+TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
+TITLE = [TITLE,'$\eta_T=',num2str(MODEL.eta_T),'$, '];
+TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
+TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
+if MODEL.CO == -1; CONAME = 'FC';
+elseif MODEL.CO == -2; CONAME = 'DC';
+elseif MODEL.CO == 0; CONAME = 'LB'; end;
+
+bare = @(pp,jj) (GRID.jmaxe +1)*pp + jj+1; % electron 1D-index
+bari = @(pp,jj) bare(GRID.pmaxe, GRID.jmaxe)+(GRID.jmaxi +1)*pp + jj+1; % ion 1D-index
+
+%% Nepj
+%
+% moment = 'moments_e';
+%
+% kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
+% kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
+% time = h5read(filename,'/data/var5d/time');
+% Nepj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
+% for it = 1:numel(time)
+% tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
+% Nepj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+% end
+%
+% fig = figure;
+%
+% %HeLaZ results
+% x1 = time;
+% y1 = x1;
+% ic = 1;
+% for ikr = 1:GRID.nkr
+% for ikz = 1:GRID.nkz
+% for ip = 0:1
+% for ij = 0:1
+% for it = 1:numel(time)
+% y1(it) = abs(Nepj(ip+1,ij+1,ikr,ikz,it));
+% end
+% semilogy(x1,y1,'-','DisplayName',...
+% ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
+% 'color', line_colors(ic,:))
+% hold on
+% ic = ic + 1;
+% end
+% end
+% end
+% end
+%
+% %MOLI results
+% x2 = results.time;
+% y2 = x2;
+% ic = 1;
+% for ip = 0:1
+% for ij = 0:1
+% for it = 1:numel(x2)
+% y2(it) = abs(results.Napj(it,bare(ip,ij)));
+% end
+% semilogy(x2,y2,'--','DisplayName',...
+% ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
+% 'color', line_colors(ic,:))
+% hold on
+% ic = ic + 1;
+% end
+% end
+%
+% title(TITLE);
+% grid on
+% legend('show')
+% xlabel('$t$')
+% ylabel(['$|N_e^{pj}|$'])
+% if SAVEFIG; FIGNAME = ['Nepj_kz_',num2str(GRID.kzmin)]; save_figure; end;
+%
+%% Nipj
+
+moment = 'moments_i';
+
+kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
+kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
+time = h5read(filename,'/data/var5d/time');
+Nipj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
+for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
+ Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+end
+
+fig = figure;
+
+x1 = time;
+x2 = results.time;
+ic = 1;
+
+% Real part
+subplot(321)
+for ip = 0:1
+ for ij = 0:1
+ y1 = squeeze(real(Nipj(ip+1,ij+1,1,1,:)));
+ plot(x1,y1,'-','DisplayName',...
+ ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ y2 = squeeze(real(results.Napj(:,bari(ip,ij))));
+ plot(x2,y2,'--','DisplayName',...
+ ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+end
+grid on
+xlabel('$t$')
+ylabel(['Re$(N_i^{pj})$'])
+
+% Im part
+ic = 1;
+subplot(322)
+for ip = 0:1
+ for ij = 0:1
+ y1 = squeeze(imag(Nipj(ip+1,ij+1,1,1,:)));
+ plot(x1,y1,'-','DisplayName',...
+ ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ y2 = squeeze(imag(results.Napj(:,bari(ip,ij))));
+ plot(x2,y2,'--','DisplayName',...
+ ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+end
+grid on
+xlabel('$t$')
+ylabel(['Im$(N_i^{pj})$'])
+%suptitle(TITLE);
+%if SAVEFIG; FIGNAME = ['Nipj_kz_',num2str(GRID.kzmin)]; save_figure; end;
+
+%% phi
+timephi = h5read(filename,'/data/var2d/time');
+kr = h5read(filename,'/data/var2d/phi/coordkr');
+kz = h5read(filename,'/data/var2d/phi/coordkz');
+phiHeLaZ = zeros(numel(timephi),numel(kr),numel(kz));
+for it = 1:numel(timephi)
+ tmp = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
+ phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
+end
+
+timephiMOLI = results.time;
+phiMOLI = zeros(size(timephiMOLI));
+for it = 1:numel(timephiMOLI)
+ phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
+end
+
+%fig = figure;
+%Real
+subplot(323)
+plot(timephi,real(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
+hold on
+plot(timephiMOLI,real(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
+grid on
+xlabel('$t$')
+ylabel('Re$(\phi)$')
+%Imag
+subplot(324)
+plot(timephi,imag(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
+hold on
+plot(timephiMOLI,imag(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
+grid on
+xlabel('$t$')
+ylabel('Im$(\phi)$')
+%if SAVEFIG; FIGNAME = ['phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
+
+%% phi error
+timephi = h5read(filename,'/data/var2d/time');
+kr = h5read(filename,'/data/var2d/phi/coordkr');
+kz = h5read(filename,'/data/var2d/phi/coordkz');
+phiHeLaZ = zeros(numel(timephi),numel(kr),numel(kz));
+for it = 1:numel(timephi)
+ tmp = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
+ phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
+end
+
+timephiMOLI = results.time;
+phiMOLI = zeros(size(timephiMOLI));
+for it = 1:numel(timephiMOLI)
+ phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
+end
+
+ERR1 = abs(real(phiMOLI - phiHeLaZ));
+ERR2 = abs(imag(phiMOLI - phiHeLaZ));
+
+%fig = figure;
+subplot(325);
+plot(timephiMOLI,ERR1,'-','DisplayName','Real')
+hold on
+plot(timephiMOLI,ERR2,'--','DisplayName','Imag')
+%title(TITLE);
+grid on
+xlabel('$t$')
+ylabel('$|e_\phi|$')
+legend('show')
+%if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
+
+suptitle(TITLE);
+if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
--
GitLab