From ed999c3bc9c51e0b90f055ebb6faa14f10e3f252 Mon Sep 17 00:00:00 2001
From: Antoine <antoine.hoffmann@epfl.ch>
Date: Mon, 4 Sep 2023 16:50:58 +0200
Subject: [PATCH] consequent reorganization of the code -removed the parameters
only modules -simplified the names of the input namelists
---
Makefile | 55 +-
fort.90 | 12 +-
matlab/load/compile_results_low_mem.m | 4 +-
matlab/load/load_params.m | 10 +-
matlab/write_fort90.m | 12 +-
src/adaptive_timestep_mod.F90 | 4 +-
src/auxval.F90 | 17 +-
src/closure_mod.F90 | 4 +-
src/collision_mod.F90 | 4 +-
src/control.F90 | 10 +-
src/diagnose.F90 | 600 ---------------
src/diagnostics_mod.F90 | 684 ++++++++++++++++++
src/diagnostics_par_mod.F90 | 92 ---
src/inital.F90 | 630 ----------------
src/initial_mod.F90 | 651 +++++++++++++++++
src/initial_par_mod.F90 | 60 --
src/main.F90 | 11 +-
src/model_mod.F90 | 8 +-
src/nonlinear_mod.F90 | 25 +-
src/readinputs.F90 | 6 +-
src/restarts_mod.F90 | 13 +-
src/time_integration_mod.F90 | 6 +-
testcases/DLRA_zpinch/base_case/fort_00.90 | 12 +-
testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 | 12 +-
testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 | 12 +-
testcases/ITG_zpinch/fort_00.90 | 12 +-
testcases/cyclone_example/fort_00.90 | 12 +-
testcases/cyclone_example/fort_01.90 | 12 +-
testcases/smallest_problem/fort.90 | 10 +-
testcases/smallest_problem/fort_00.90 | 10 +-
testcases/smallest_problem/fort_01.90 | 10 +-
testcases/smallest_problem/fort_11.90 | 10 +-
testcases/zpinch_3D/fort_00.90 | 12 +-
testcases/zpinch_3D/fort_01.90 | 12 +-
testcases/zpinch_3D/fort_02.90 | 12 +-
testcases/zpinch_example/fort_00.90 | 12 +-
wk/benchmark_and_scan_scripts/old/fort_00.90 | 12 +-
wk/paper_2_scripts_and_results/fort_00.90 | 12 +-
38 files changed, 1538 insertions(+), 1564 deletions(-)
delete mode 100644 src/diagnose.F90
create mode 100644 src/diagnostics_mod.F90
delete mode 100644 src/diagnostics_par_mod.F90
delete mode 100644 src/inital.F90
create mode 100644 src/initial_mod.F90
delete mode 100644 src/initial_par_mod.F90
diff --git a/Makefile b/Makefile
index 6087f542..d2e104d2 100644
--- a/Makefile
+++ b/Makefile
@@ -117,11 +117,11 @@ $(OBJDIR)/diagnose.o : src/srcinfo.h
FOBJ=$(OBJDIR)/advance_field_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o \
$(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/circular_mod.o \
$(OBJDIR)/collision_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o \
-$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \
+$(OBJDIR)/diagnostics_mod.o $(OBJDIR)/endrun.o \
$(OBJDIR)/ExB_shear_flow_mod.o \
$(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/geometry_mod.o \
-$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/inital.o \
-$(OBJDIR)/initial_par_mod.o $(OBJDIR)/lag_interp_mod.o $(OBJDIR)/main.o \
+$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o \
+$(OBJDIR)/initial_mod.o $(OBJDIR)/lag_interp_mod.o $(OBJDIR)/main.o \
$(OBJDIR)/memory.o $(OBJDIR)/miller_mod.o $(OBJDIR)/model_mod.o \
$(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/parallel_mod.o \
$(OBJDIR)/ppexit.o $(OBJDIR)/prec_const_mod.o \
@@ -138,7 +138,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
# Modules compilation
$(OBJDIR)/advance_field_mod.o : src/advance_field_mod.F90 \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_mod.o \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o \
$(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/closure_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field_mod.F90 -o $@
@@ -187,25 +187,21 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/control.o : src/control.F90 \
$(OBJDIR)/auxval.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o \
+ $(OBJDIR)/diagnostics_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
$(OBJDIR)/cosolver_interface_mod.o : src/cosolver_interface_mod.F90 \
$(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/species_mod.o\
- $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/closure_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/cosolver_interface_mod.F90 -o $@
- $(OBJDIR)/diagnose.o : src/diagnose.F90 \
+ $(OBJDIR)/diagnostics_mod.o : src/diagnostics_mod.F90 \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/array_mod.o \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o \
- $(OBJDIR)/time_integration_mod.o\
- $(OBJDIR)/parallel_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnose.F90 -o $@
-
- $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 \
- $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_par_mod.F90 -o $@
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_mod.o $(OBJDIR)/model_mod.o \
+ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/parallel_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_mod.F90 -o $@
$(OBJDIR)/endrun.o : src/endrun.F90 \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
@@ -239,17 +235,14 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
- $(OBJDIR)/inital.o : src/inital.F90 \
+ $(OBJDIR)/initial_mod.o : src/initial_mod.F90 \
$(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
- $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o \
+ $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o \
$(OBJDIR)/solve_EM_fields.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o \
- $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/collision_mod.o $(OBJDIR)/closure_mod.o \
+ $(OBJDIR)/nonlinear_mod.o \
$(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
-
- $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_mod.F90 -o $@
$(OBJDIR)/lag_interp_mod.o : src/lag_interp_mod.F90 \
$(OBJDIR)/prec_const_mod.o
@@ -287,7 +280,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o \
- $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/utility_mod.o $(OBJDIR)/species_mod.o $(OBJDIR)/geometry_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics_mod.F90 -o $@
$(OBJDIR)/parallel_mod.o : src/parallel_mod.F90 \
@@ -307,18 +300,19 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/processing_mod.F90 -o $@
$(OBJDIR)/readinputs.o : src/readinputs.F90 \
- $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o \
+ $(OBJDIR)/diagnostics_mod.o $(OBJDIR)/initial_mod.o $(OBJDIR)/model_mod.o \
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
$(OBJDIR)/restarts_mod.o : src/restarts_mod.F90 \
- $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/restarts_mod.F90 -o $@
$(OBJDIR)/solve_EM_fields.o : src/solve_EM_fields.F90 \
$(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o \
- $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/parallel_mod.o
+ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/parallel_mod.o \
+ $(OBJDIR)/processing_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/solve_EM_fields.F90 -o $@
$(OBJDIR)/species_mod.o : src/species_mod.F90 \
@@ -326,7 +320,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/species_mod.F90 -o $@
$(OBJDIR)/stepon.o : src/stepon.F90 \
- $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field_mod.o \
+ $(OBJDIR)/initial_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field_mod.o \
$(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o \
$(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o\
@@ -343,7 +337,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
$(OBJDIR)/utility_mod.o : src/utility_mod.F90 \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
+ $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
$(OBJDIR)/CLA_mod.o : src/CLA_mod.F90 \
diff --git a/fort.90 b/fort.90
index 3638ba08..c536e536 100644
--- a/fort.90
+++ b/fort.90
@@ -35,7 +35,7 @@
shift_y= 0.0
Npol = 1.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -52,7 +52,7 @@
write_temp = .t.
!diag_mode = 'txtonly'
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 1.0
@@ -67,7 +67,7 @@
ADIAB_E = .f.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
@@ -91,19 +91,19 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'phi' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index 54ee00f4..a90409a8 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -39,8 +39,8 @@ while(CONTINUE)
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
[Parray, Jarray, kx, ky, z] = load_grid_data(filename);
- W_GAMMA = strcmp(h5readatt(filename,'/data/input/diag_par','write_gamma'),'y');
- W_HF = strcmp(h5readatt(filename,'/data/input/diag_par','write_hf' ),'y');
+ W_GAMMA = strcmp(h5readatt(filename,'/data/input/diagnostics','write_gamma'),'y');
+ W_HF = strcmp(h5readatt(filename,'/data/input/diagnostics','write_hf' ),'y');
KIN_E = strcmp(h5readatt(filename,'/data/input/model', 'ADIAB_E' ),'n');
BETA = h5readatt(filename,'/data/input/model','beta');
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index ce0f1794..967baf3c 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -51,11 +51,11 @@ function [DATA] = load_params(DATA,filename)
catch
DATA.inputs.ADIAB_I = 0;
end
- DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
- DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
- DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
- DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y';
- DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y';
+ DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diagnostics','write_gamma') == 'y';
+ DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diagnostics','write_phi') == 'y';
+ DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diagnostics','write_Na00') == 'y';
+ DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diagnostics','write_Napj') == 'y';
+ DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diagnostics','write_Sapj') == 'y';
% Species dependent parameters
DATA.inputs.sigma = zeros(1,DATA.inputs.Na);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 0c0dc4e9..f0bda4cd 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -40,7 +40,7 @@ fprintf(fid,[' Npol = ', num2str(GEOM.Npol),'\n']);
fprintf(fid,[' PB_PHASE= ', GEOM.PB_PHASE,'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&OUTPUT_PAR\n');
+fprintf(fid,'&DIAGNOSTICS\n');
fprintf(fid,[' dtsave_0d = ', num2str(OUTPUTS.dtsave_0d),'\n']);
fprintf(fid,[' dtsave_1d = ', num2str(OUTPUTS.dtsave_1d),'\n']);
fprintf(fid,[' dtsave_2d = ', num2str(OUTPUTS.dtsave_2d),'\n']);
@@ -56,7 +56,7 @@ fprintf(fid,[' write_dens = ', OUTPUTS.write_dens,'\n']);
fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&MODEL_PAR\n');
+fprintf(fid,'&MODEL\n');
fprintf(fid,['LINEARITY = ', MODEL.LINEARITY,'\n']);
fprintf(fid,['RM_LD_T_EQ= ', MODEL.RM_LD_T_EQ,'\n']);
fprintf(fid,[' Na = ', num2str(MODEL.Na),'\n']);
@@ -79,7 +79,7 @@ fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
fprintf(fid,[' MHD_PD = ', MODEL.MHD_PD,'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&CLOSURE_PAR\n');
+fprintf(fid,'&CLOSURE\n');
fprintf(fid,[' hierarchy_closure=',CLOSURE.hierarchy_closure,'\n']);
fprintf(fid,[' dmax =',num2str(CLOSURE.dmax),'\n']);
fprintf(fid,[' nonlinear_closure=',CLOSURE.nonlinear_closure,'\n']);
@@ -107,7 +107,7 @@ if(strcmp(MODEL.ADIAB_E,'.false.'))
fprintf(fid,'/\n');
end
-fprintf(fid,'&COLLISION_PAR\n');
+fprintf(fid,'&COLLISION\n');
fprintf(fid,[' collision_model = ', COLL.collision_model,'\n']);
fprintf(fid,[' GK_CO = ', COLL.GK_CO,'\n']);
fprintf(fid,[' INTERSPECIES = ', COLL.INTERSPECIES,'\n']);
@@ -116,14 +116,14 @@ fprintf(fid,[' collision_kcut = ', num2str(COLL.coll_kcut),'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&INITIAL_CON\n');
+fprintf(fid,'&INITIAL\n');
fprintf(fid,[' INIT_OPT = ', INITIAL.INIT_OPT,'\n']);
fprintf(fid,[' init_background = ', num2str(INITIAL.init_background),'\n']);
fprintf(fid,[' init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
fprintf(fid,[' iseed = ', num2str(INITIAL.iseed),'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&TIME_INTEGRATION_PAR\n');
+fprintf(fid,'&TIME_INTEGRATION\n');
fprintf(fid,[' numerical_scheme = ', TIME_INTEGRATION.numerical_scheme,'\n']);
fprintf(fid,'/\n');
diff --git a/src/adaptive_timestep_mod.F90 b/src/adaptive_timestep_mod.F90
index c8581ea6..f663258e 100644
--- a/src/adaptive_timestep_mod.F90
+++ b/src/adaptive_timestep_mod.F90
@@ -31,9 +31,9 @@ MODULE adaptive_timestep
USE prec_const
USE basic, ONLY : lu_in
IMPLICIT NONE
- namelist /TIME_INTEGRATION_PAR/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
+ namelist /TIME_INTEGRATION/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
- READ(lu_in,time_integration_par)
+ READ(lu_in,time_integration)
CALL set_numerical_scheme
END SUBROUTINE adaptive_timestep_readinputs
diff --git a/src/auxval.F90 b/src/auxval.F90
index 6828d2aa..32966633 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -1,14 +1,17 @@
subroutine auxval
! Set auxiliary values, at beginning of simulation
- USE basic
- USE grid
- USE array
- USE model
+ USE basic, ONLY: str, speak
+ USE grid, ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset, total_nj,&
+ local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset, dmax,&
+ local_nz, local_nz_offset, total_nz, init_grids_data, set_grids
+ !USE array
+ USE model, ONLY: Na, EM, LINEARITY, N_HD
USE fourier, ONLY: init_grid_distr_and_plans
- use prec_const
- USE numerics
- USE geometry
+ use MPI, ONLY: MPI_COMM_WORLD
+ USE numerics, ONLY: build_dnjs_table, build_dv4Hp_table, compute_lin_coeff, &
+ evaluate_EM_op, evaluate_kernels
+ USE geometry, ONLY: Npol, shear, eval_magnetic_geometry
USE closure, ONLY: set_closure_model, hierarchy_closure
USE parallel, ONLY: init_parallel_var, my_id, num_procs, &
num_procs_p, num_procs_z, num_procs_ky, rank_p, rank_ky, rank_z
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 8fd9163c..c580905d 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -17,8 +17,8 @@ CONTAINS
SUBROUTINE closure_readinputs
USE basic, ONLY: lu_in
IMPLICIT NONE
- NAMELIST /CLOSURE_PAR/ hierarchy_closure, dmax, nonlinear_closure, nmax
- READ(lu_in,closure_par)
+ NAMELIST /CLOSURE/ hierarchy_closure, dmax, nonlinear_closure, nmax
+ READ(lu_in,closure)
END SUBROUTINE
SUBROUTINE set_closure_model
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 66df93ec..ef4e0d78 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -27,8 +27,8 @@ CONTAINS
! Read the input parameters
USE basic, ONLY: lu_in
IMPLICIT NONE
- NAMELIST /COLLISION_PAR/ collision_model, GK_CO, INTERSPECIES, mat_file, collision_kcut
- READ(lu_in,collision_par)
+ NAMELIST /COLLISION/ collision_model, GK_CO, INTERSPECIES, mat_file, collision_kcut
+ READ(lu_in,collision)
SELECT CASE(collision_model)
! Lenhard bernstein
CASE ('LB','lenhard-bernstein','Lenhard-Bernstein')
diff --git a/src/control.F90 b/src/control.F90
index c8e307da..bb743ee6 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -4,9 +4,11 @@ SUBROUTINE control
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: xp, stdout
- USE parallel, ONLY: ppinit
- USE mpi
+ use prec_const, ONLY: xp, stdout
+ USE parallel, ONLY: ppinit
+ USE initial, ONLY: initialize
+ USE mpi, ONLY: MPI_COMM_WORLD
+ USE diagnostics, ONLY: diagnose
IMPLICIT NONE
REAL(xp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
@@ -41,7 +43,7 @@ SUBROUTINE control
! 1.5 Initial conditions
CALL speak( 'Create initial state...')
- CALL inital
+ CALL initialize
! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
CALL speak('...initial state created')
diff --git a/src/diagnose.F90 b/src/diagnose.F90
deleted file mode 100644
index 8807f076..00000000
--- a/src/diagnose.F90
+++ /dev/null
@@ -1,600 +0,0 @@
-SUBROUTINE diagnose(kstep)
- ! Diagnostics, writing simulation state to disk
- USE basic, ONLY: lu_in, chrono_runt, cstep, dt, time, tmax, display_h_min_s
- USE diagnostics_par, ONLY: input_fname, diag_mode, nsave_0d
- USE processing, ONLY: pflux_x, hflux_x
- USE parallel, ONLY: my_id
- IMPLICIT NONE
- INTEGER, INTENT(in) :: kstep
- !! Basic diagnose loop for reading input file, displaying advancement and ending
- IF ((kstep .EQ. 0)) THEN
- INQUIRE(unit=lu_in, name=input_fname)
- CLOSE(lu_in)
- ENDIF
- !! End diag
- IF (kstep .EQ. -1) THEN
- ! Display total run time
- CALL display_h_min_s(chrono_runt%ttot)
- END IF
- !! Specific diagnostic calls
- SELECT CASE(diag_mode)
- CASE('full')
- CALL diagnose_full(kstep)
- CASE('txtonly')
- CALL diagnose_txtonly(kstep)
- END SELECT
-
- ! Terminal info
- IF ((kstep .GE. 0) .AND. (MOD(cstep, nsave_0d) == 0) .AND. (my_id .EQ. 0)) THEN
- if(SIZE(pflux_x) .GT. 1) THEN
- WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A8,G10.3,A8,G10.3,A)")&
- '|t = ', time,'| Pxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'| Pxe = ',pflux_x(2),'| Qxe = ',hflux_x(2),'|'
- else
- WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A)")&
- '|t = ', time,'| Px = ',pflux_x(1),'| Qx = ',hflux_x(1),'|'
- endif
- ENDIF
-END SUBROUTINE diagnose
-
-SUBROUTINE init_outfile(comm,file0,file,fid)
- USE diagnostics_par, ONLY: write_doubleprecision, diag_par_outputinputs, input_fname
- USE basic, ONLY: speak, jobnum, basic_outputinputs
- USE grid, ONLY: grid_outputinputs
- USE geometry, ONLY: geometry_outputinputs
- USE model, ONLY: model_outputinputs
- USE closure, ONLY: closure_outputinputs
- USE species, ONLY: species_outputinputs
- USE collision, ONLY: coll_outputinputs
- USE initial_par, ONLY: initial_outputinputs
- USE time_integration,ONLY: time_integration_outputinputs
- USE parallel, ONLY: parallel_outputinputs
- USE futils, ONLY: creatf, creatg, creatd, attach, putfile
- IMPLICIT NONE
- !input
- INTEGER, INTENT(IN) :: comm
- CHARACTER(len=256), INTENT(IN) :: file0
- CHARACTER(len=256), INTENT(OUT) :: file
- INTEGER, INTENT(OUT) :: fid
- CHARACTER(len=256) :: str
- INCLUDE 'srcinfo.h'
- ! Writing output filename
- WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0) ,'_',jobnum,'.h5'
- ! 1.1 Initial run
- ! Main output file creation
- IF (write_doubleprecision) THEN
- CALL creatf(file, fid, real_prec='d', mpicomm=comm)
- ELSE
- CALL creatf(file, fid, mpicomm=comm)
- END IF
- CALL speak(TRIM(file)//' created')
- ! basic data group
- CALL creatg(fid, "/data", "data")
- ! File group
- CALL creatg(fid, "/files", "files")
- CALL attach(fid, "/files", "jobnum", jobnum)
- ! Add the code info and parameters to the file
- CALL creatg(fid, "/data/input", "input")
- CALL creatd(fid, 0,(/0/),"/data/input/codeinfo",'Code Information')
- CALL attach(fid, "/data/input/codeinfo", "version", VERSION) !defined in srcinfo.h
- CALL attach(fid, "/data/input/codeinfo", "branch", BRANCH) !defined in srcinfo.h
- CALL attach(fid, "/data/input/codeinfo", "author", AUTHOR) !defined in srcinfo.h
- CALL attach(fid, "/data/input/codeinfo", "execdate", EXECDATE) !defined in srcinfo.h
- CALL attach(fid, "/data/input/codeinfo", "host", HOST) !defined in srcinfo.h
- CALL basic_outputinputs(fid)
- CALL grid_outputinputs(fid)
- CALL geometry_outputinputs(fid)
- CALL diag_par_outputinputs(fid)
- CALL model_outputinputs(fid)
- CALL closure_outputinputs(fid)
- CALL species_outputinputs(fid)
- CALL coll_outputinputs(fid)
- CALL initial_outputinputs(fid)
- CALL time_integration_outputinputs(fid)
- CALL parallel_outputinputs(fid)
- ! Save STDIN (input file) of this run
- IF(jobnum .LE. 99) THEN
- WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
- ELSE
- WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum
- END IF
- CALL putfile(fid, TRIM(str), TRIM(input_fname),ionode=0)
-END SUBROUTINE init_outfile
-
-SUBROUTINE diagnose_full(kstep)
- USE basic, ONLY: speak,chrono_runt,&
- cstep,iframe0d,iframe3d,iframe5d,crashed
- USE grid, ONLY: &
- parray_full,pmax,jarray_full,jmax, kparray, &
- kyarray_full,kxarray_full,zarray_full, ngz, total_nz, local_nz, ieven,&
- local_Nky, total_nky, local_nkx, total_nkx
- USE geometry, ONLY: gxx, gxy, gxz, gyy, gyz, gzz, &
- hatR, hatZ, hatB, dBdx, dBdy, dBdz, Jacobian, gradz_coeff
- USE diagnostics_par
- USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf!, putarrnd ! Routine de merde, jamais l'utiliser
- USE array
- USE model, ONLY: EM
- USE parallel, ONLY: my_id, comm0, gather_z, gather_xyz_real
- USE collision, ONLY: coll_outputinputs
- IMPLICIT NONE
- INTEGER, INTENT(in) :: kstep
- INTEGER, parameter :: BUFSIZE = 2
- REAL(xp), DIMENSION(total_nz) :: Az_full ! full z array for metric output
- REAL(xp), DIMENSION(total_nky,total_nkx,total_nz) :: kp_full
- INTEGER :: rank = 0, ierr
- INTEGER :: dims(1) = (/0/)
- !____________________________________________________________________________
- ! 1. Initial diagnostics
- IF ((kstep .EQ. 0)) THEN
- CALL init_outfile(comm0, resfile0,resfile,fidres)
- ! Profiler time measurement
- CALL creatg(fidres, "/profiler", "performance analysis")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs", "cumulative rhs computation time")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_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_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_ExBshear", "cumulative ExB shear 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")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative diagnostic time")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
- CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
-#ifdef TEST_SVD
- CALL creatd(fidres, 0, (/0/), "/profiler/Tc_CLA", "cumulative total CLA computation time")
-#endif
- ! Grid info
- CALL creatg(fidres, "/data/grid", "Grid data")
- CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0)
- CALL putarr(fidres, "/data/grid/coordky", kyarray_full, "ky*rho_s0", ionode=0)
- CALL putarr(fidres, "/data/grid/coordz", zarray_full, "z/R", ionode=0)
- CALL putarr(fidres, "/data/grid/coordp" , parray_full, "p", ionode=0)
- CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
- CALL gather_xyz_real(kparray(1:local_Nky,1:local_Nkx,1:local_nz,ieven),kp_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
- CALL putarr(fidres, "/data/grid/coordkp" , kp_full, "kp", ionode=0)
- ! Metric info
- CALL creatg(fidres, "/data/metric", "Metric data")
- CALL gather_z(gxx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gxx", Az_full, "gxx", ionode =0)
- CALL gather_z(gxy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gxy", Az_full, "gxy", ionode =0)
- CALL gather_z(gxz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gxz", Az_full, "gxz", ionode =0)
- CALL gather_z(gyy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gyy", Az_full, "gyy", ionode =0)
- CALL gather_z(gyz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gyz", Az_full, "gyz", ionode =0)
- CALL gather_z(gzz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gzz", Az_full, "gzz", ionode =0)
- CALL gather_z(hatR((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/hatR", Az_full, "hatR", ionode =0)
- CALL gather_z(hatZ((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/hatZ", Az_full, "hatZ", ionode =0)
- CALL gather_z(hatB((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/hatB", Az_full, "hatB", ionode =0)
- CALL gather_z(dBdx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/dBdx", Az_full, "dBdx", ionode =0)
- CALL gather_z(dBdy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/dBdy", Az_full, "dBdy", ionode =0)
- CALL gather_z(dBdz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/dBdz", Az_full, "dBdz", ionode =0)
- CALL gather_z(Jacobian((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/Jacobian", Az_full, "Jacobian", ionode =0)
- CALL gather_z(gradz_coeff((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
- CALL putarr(fidres, "/data/metric/gradz_coeff", Az_full, "gradz_coeff", ionode =0)
- ! var0d group (gyro transport)
- IF (nsave_0d .GT. 0) THEN
- CALL creatg(fidres, "/data/var0d", "0d profiles")
- CALL creatd(fidres, rank, dims, "/data/var0d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var0d/cstep", "iteration number")
- IF (write_gamma) THEN
- CALL creatd(fidres, rank, dims, "/data/var0d/gflux_x", "Radial gyro transport")
- CALL creatd(fidres, rank, dims, "/data/var0d/pflux_x", "Radial part transport")
- ENDIF
- IF (write_hf) THEN
- CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x", "Radial part heat flux")
- ENDIF
- IF (cstep==0) THEN
- iframe0d=0
- ENDIF
- CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
- END IF
- ! var2d group (gyro transport)
- IF (nsave_0d .GT. 0) THEN
- CALL creatg(fidres, "/data/var2d", "2d profiles")
- CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
-#ifdef TEST_SVD
- CALL creatg(fidres, "/data/var2d/sv_ky_pj", "singular values of the moment ky/pj cut")
-#endif
- ENDIF
- ! var3d group (phi,psi, fluid moments, Ni00, Napjz)
- IF (nsave_3d .GT. 0) THEN
- CALL creatg(fidres, "/data/var3d", "3d profiles")
- CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
- IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
- IF (write_phi.AND.EM) CALL creatg(fidres, "/data/var3d/psi", "psi")
- IF (write_Na00) THEN
- CALL creatg(fidres, "/data/var3d/Na00", "gyrocenter density ")
- CALL creatg(fidres, "/data/var3d/Napjz", "pj(z) moment spectrum ")
- ENDIF
- IF (write_dens) THEN
- CALL creatg(fidres, "/data/var3d/dens", "density ")
- ENDIF
- IF (write_fvel) THEN
- CALL creatg(fidres, "/data/var3d/upar", "parallel fluid velocity ")
- CALL creatg(fidres, "/data/var3d/uper", "perpendicular fluid velocity ")
- ENDIF
- IF (write_temp) THEN
- CALL creatg(fidres, "/data/var3d/Tper", "perpendicular temperature ")
- CALL creatg(fidres, "/data/var3d/Tpar", "parallel temperature ")
- CALL creatg(fidres, "/data/var3d/temp", "tiotal temperature ")
- ENDIF
- IF (cstep==0) THEN
- iframe3d=0
- ENDIF
- CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
- END IF
- ! var5d group (moments)
- IF (nsave_5d .GT. 0) THEN
- CALL creatg(fidres, "/data/var5d", "5d profiles")
- CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
- IF (write_Napj) THEN
- CALL creatg(fidres, "/data/var5d/moments", "full moments array")
- ENDIF
- IF (cstep==0) THEN
- iframe5d=0
- END IF
- CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
- END IF
- ENDIF
- !_____________________________________________________________________________
- ! 2. Periodic diagnostics
- !
- IF (kstep .GE. 0) THEN
- ! 2.1 0d history arrays
- IF (nsave_0d .GT. 0) THEN
- IF ( MOD(cstep, nsave_0d) == 0 ) THEN
- CALL diagnose_0d
- END IF
- END IF
- ! 2.2 2d profiles
- IF (nsave_2d .GT. 0) THEN
- IF (MOD(cstep, nsave_2d) == 0) THEN
- CALL diagnose_2d
- ENDIF
- ENDIF
- ! 2.3 3d profiles
- IF (nsave_3d .GT. 0) THEN
- IF (MOD(cstep, nsave_3d) == 0) THEN
- CALL diagnose_3d
- ! Looks at the folder if the file check_phi exists and spits a snapshot
- ! of the current electrostatic potential in a basic text file
- CALL spit_snapshot_check
- ENDIF
- ENDIF
- ! 2.4 5d profiles
- IF (nsave_5d .GT. 0) THEN
- IF (MOD(cstep, nsave_5d) == 0) THEN
- CALL diagnose_5d
- END IF
- END IF
- !_____________________________________________________________________________
- ! 3. Final diagnostics
- ELSEIF (kstep .EQ. -1) THEN
- 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'
- IF (nsave_5d .GT. 0) CALL diagnose_5d
- ENDIF
- ! Close all diagnostic files
- CALL mpi_barrier(MPI_COMM_WORLD, ierr)
- CALL closef(fidres)
- END IF
-END SUBROUTINE diagnose_full
-
-!! This routine outputs only txt file 0D data (flux and time)
-SUBROUTINE diagnose_txtonly(kstep)
- USE basic
- USE diagnostics_par
- USE processing, ONLY: pflux_x, hflux_x, compute_radial_transport, compute_radial_heatflux
- USE parallel, ONLY: my_id, comm0
- USE futils, ONLY: creatf, creatg, creatd, attach, putfile, closef
- IMPLICIT NONE
- INTEGER, INTENT(in) :: kstep
- INTEGER, parameter :: BUFSIZE = 2
- INTEGER :: rank = 0
- INTEGER :: dims(1) = (/0/)
- IF (kstep .GE. 0) THEN
- ! output the transport in a txt file
- IF ( MOD(cstep, nsave_0d) == 0 ) THEN
- CALL compute_radial_transport
- CALL compute_radial_heatflux
- IF (my_id .EQ. 0) &
- WRITE(1,*) time, pflux_x(1), hflux_x(1)
- END IF
- !! Save the last state
- ELSEIF (kstep .EQ. -1) THEN
- CALL init_outfile(comm0, resfile0,resfile,fidres)
- CALL creatg(fidres, "/data/var5d", "5d profiles")
- CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
- CALL creatg(fidres, "/data/var5d/moments", "full moments array")
- CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
- CALL diagnose_5d
- CALL closef(fidres)
- ENDIF
-END SUBROUTINE diagnose_txtonly
-
-!!-------------- Auxiliary routines -----------------!!
-SUBROUTINE diagnose_0d
- USE basic
- USE futils, ONLY: append, attach, getatt, creatd
- USE diagnostics_par
- USE prec_const
- USE processing
- USE model, ONLY: Na
- IMPLICIT NONE
- ! Time measurement data
- CALL append(fidres, "/profiler/Tc_rhs", REAL(chrono_mrhs%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_adv_field", REAL(chrono_advf%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_clos", REAL(chrono_clos%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_ghost", REAL(chrono_ghst%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_coll", REAL(chrono_coll%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_poisson", REAL(chrono_pois%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_Sapj", REAL(chrono_sapj%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_checkfield",REAL(chrono_chck%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_diag", REAL(chrono_diag%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_grad", REAL(chrono_grad%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_nadiab", REAL(chrono_napj%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_ExBshear", REAL(chrono_ExBs%ttot,dp),ionode=0)
- CALL append(fidres, "/profiler/Tc_step", REAL(chrono_step%ttot,dp),ionode=0)
-#ifdef TEST_SVD
- CALL append(fidres, "/profiler/Tc_CLA", REAL(chrono_CLA%ttot,dp),ionode=0)
-#endif
- CALL append(fidres, "/profiler/time", REAL(time,dp),ionode=0)
- ! Processing data
- CALL append(fidres, "/data/var0d/time", REAL(time,dp),ionode=0)
- CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var0d/", "frames",iframe0d)
- iframe0d=iframe0d+1
- CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
- ! Ion transport data
- IF (write_gamma) THEN
- CALL compute_radial_transport
- CALL append(fidres, "/data/var0d/gflux_x",REAL(gflux_x,dp),ionode=0)
- CALL append(fidres, "/data/var0d/pflux_x",REAL(pflux_x,dp),ionode=0)
- ENDIF
- IF (write_hf) THEN
- CALL compute_radial_heatflux
- CALL append(fidres, "/data/var0d/hflux_x",REAL(hflux_x,dp),ionode=0)
- ENDIF
-END SUBROUTINE diagnose_0d
-
-SUBROUTINE diagnose_2d
- USE prec_const
- USE basic
- USE diagnostics_par
- USE futils, ONLY: putarr, append
-#ifdef TEST_SVD
- USE CLA, ONLY: Sf
-#endif
- IMPLICIT NONE
- iframe2d=iframe2d+1
- CALL append(fidres,"/data/var2d/time", REAL(time,dp), ionode=0)
- CALL append(fidres,"/data/var2d/cstep",REAL(cstep,dp),ionode=0)
-#ifdef TEST_SVD
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var2d/sv_ky_pj/", iframe2d
- CALL putarr(fidres, dset_name, Sf, ionode=0)
-#endif
-END SUBROUTINE
-
-SUBROUTINE diagnose_3d
- USE basic
- USE futils, ONLY: append, getatt, attach, putarr!, putarrnd
- USE fields, ONLY: phi, psi, moments
- USE array, ONLY: Napjz,dens,upar,uper,Tpar,Tper,temp
- USE grid, ONLY: CONTAINSp0, ip0,ij0, local_na, total_na,&
- total_np, total_nj, total_nky, total_nkx, total_nz, &
- local_np, local_nky, local_nkx, local_nz, &
- ngz
- USE time_integration, ONLY: updatetlevel
- USE diagnostics_par
- USE prec_const
- USE processing, ONLY: compute_fluid_moments, compute_Napjz_spectrum
- USE model, ONLY: EM
- USE parallel, ONLY: manual_3D_bcast
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(local_na,local_nky,local_nkx,local_nz) :: Na00_
- ! add current time, cstep and frame
- CALL append(fidres, "/data/var3d/time", REAL(time,dp),ionode=0)
- CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
- iframe3d=iframe3d+1
- CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
- ! Write current EM fields
- IF (write_phi) CALL write_field3d_kykxz(phi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'phi')
- IF (write_phi.AND.EM) CALL write_field3d_kykxz(psi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'psi')
- IF (write_Na00) THEN
- CALL compute_Napjz_spectrum
- IF (CONTAINSp0) THEN
- ! gyrocenter density
- Na00_ = moments(:,ip0,ij0,:,:,(1+ngz/2):(local_nz+ngz/2),updatetlevel)
- ELSE
- Na00_ = 0._xp
- ENDIF
- CALL write_field3da_kykxz(Na00_, 'Na00')
- ! <<Napj>x>y spectrum
- CALL write_field3da_pjz(Napjz, 'Napjz')
- ENDIF
- !! Fuid moments
- IF (write_dens .OR. write_fvel .OR. write_temp) &
- CALL compute_fluid_moments
- IF (write_dens) THEN
- CALL write_field3da_kykxz(dens, 'dens')
- ENDIF
- IF (write_fvel) THEN
- CALL write_field3da_kykxz(upar, 'upar')
- CALL write_field3da_kykxz(uper, 'uper')
- ENDIF
- IF (write_temp) THEN
- CALL write_field3da_kykxz(Tpar, 'Tpar')
- CALL write_field3da_kykxz(Tper, 'Tper')
- CALL write_field3da_kykxz(temp, 'temp')
- ENDIF
- CONTAINS
- SUBROUTINE write_field3d_kykxz(field, text)
- USE parallel, ONLY : gather_xyz, num_procs
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(:,:,:), INTENT(IN) :: field
- CHARACTER(*), INTENT(IN) :: text
- COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_full
- CHARACTER(50) :: dset_name
- field_full = 0;
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
- IF (num_procs .EQ. 1) THEN ! no data distribution
- CALL putarr(fidres, dset_name, field, ionode=0)
- ELSE
- CALL gather_xyz(field,field_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
- CALL putarr(fidres, dset_name, field_full, ionode=0)
- ENDIF
- END SUBROUTINE write_field3d_kykxz
-
- SUBROUTINE write_field3da_kykxz(field, text)
- USE parallel, ONLY : gather_xyz, num_procs
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
- CHARACTER(*), INTENT(IN) :: text
- COMPLEX(xp), DIMENSION(total_na,total_nky,total_nkx,total_nz) :: field_full
- COMPLEX(xp), DIMENSION(local_nky,total_nkx,local_nz) :: buff_local
- COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: buff_full
- CHARACTER(50) :: dset_name
- INTEGER :: ia
- field_full = 0;
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
- IF (num_procs .EQ. 1) THEN ! no data distribution
- CALL putarr(fidres, dset_name, field, ionode=0)
- ELSE
- DO ia = 1,total_Na
- buff_local = field(ia,:,:,:)
- CALL gather_xyz(buff_local,buff_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
- field_full(ia,:,:,:) = buff_full
- ENDDO
- CALL putarr(fidres, dset_name, field_full, ionode=0)
- ENDIF
- END SUBROUTINE write_field3da_kykxz
-
- SUBROUTINE write_field3da_pjz(field, text)
- USE parallel, ONLY : gather_pjz, num_procs
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
- CHARACTER(*), INTENT(IN) :: text
- COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nz) :: field_full
- CHARACTER(LEN=50) :: dset_name
- field_full = 0;
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
- IF (num_procs .EQ. 1) THEN ! no data distribution
- CALL putarr(fidres, dset_name, field, ionode=0)
- ELSE
- CALL gather_pjz(field,field_full,total_na,local_np,total_np,total_nj,local_nz,total_nz)
- CALL putarr(fidres, dset_name, field_full, ionode=0)
- ENDIF
- CALL attach(fidres, dset_name, "time", time)
- END SUBROUTINE write_field3da_pjz
-
-END SUBROUTINE diagnose_3d
-
-SUBROUTINE diagnose_5d
-
- USE basic, ONLY: time, iframe5d,cstep
- USE futils, ONLY: append, getatt, attach, putarr !,putarrnd
- USE fields, ONLY: moments
- USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
- local_np, local_nj, local_nky, local_nkx, local_nz, &
- ngp, ngj, ngz, total_na
- USE diagnostics_par
- USE prec_const, ONLY: xp,dp
- IMPLICIT NONE
-
- CALL append(fidres, "/data/var5d/time", REAL(time,dp),ionode=0)
- CALL append(fidres, "/data/var5d/cstep", REAL(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
- iframe5d=iframe5d+1
- CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
-
- IF (write_Napj) THEN
- CALL write_field5d(moments, 'moments')
- ENDIF
-
- CONTAINS
-
- SUBROUTINE write_field5d(field, text)
- USE basic, ONLY: jobnum, dt
- USE futils, ONLY: attach, putarr!, putarrnd
- USE parallel, ONLY: gather_pjxyz, num_procs
- USE prec_const, ONLY: xp
- USE time_integration, ONLY: ntimelevel
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_Nky,local_Nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
- CHARACTER(*), INTENT(IN) :: text
- COMPLEX(xp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
- COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
- CHARACTER(LEN=50) :: dset_name
- field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
- 1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),1)
- field_full = 0;
- WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
- IF (num_procs .EQ. 1) THEN
- CALL putarr(fidres, dset_name, field_sub, ionode=0)
- ELSE!IF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
- CALL gather_pjxyz(field_sub,field_full,total_na,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
- CALL putarr(fidres, dset_name, field_full, ionode=0)
- ! ELSE
- ! CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/))
- ENDIF
- CALL attach(fidres, dset_name, 'cstep', cstep)
- CALL attach(fidres, dset_name, 'time', time)
- CALL attach(fidres, dset_name, 'jobnum', jobnum)
- CALL attach(fidres, dset_name, 'dt', dt)
- CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
- END SUBROUTINE write_field5d
-END SUBROUTINE diagnose_5d
-
-SUBROUTINE spit_snapshot_check
- USE fields, ONLY: phi
- USE grid, ONLY: total_nkx,total_nky,total_nz,&
- local_nky,local_nz, ngz
- USE parallel, ONLY: gather_xyz, my_id
- USE basic
- USE prec_const, ONLY: xp
- IMPLICIT NONE
- LOGICAL :: file_exist
- INTEGER :: fid_check, ikx, iky, iz
- CHARACTER(256) :: check_filename
- COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_to_check
- !! Spit a snapshot of PHI if requested (triggered by creating a file named "check_phi")
- INQUIRE(file='check_phi', exist=file_exist)
- IF( file_exist ) THEN
- IF(my_id.EQ. 0) WRITE(*,*) 'Check file found -> gather phi..'
- CALL gather_xyz(phi(:,:,(1+Ngz/2):(local_nz+Ngz/2)), field_to_check,local_nky,total_nky,total_nkx,local_nz,total_nz)
- IF(my_id.EQ. 0) THEN
- WRITE(check_filename,'(a16)') 'check_phi.out'
- OPEN(fid_check, file=check_filename, form='formatted')
- WRITE(*,*) 'Check file found -> output phi ..'
- WRITE(fid_check,*) total_nky, total_nkx, total_nz
- DO iky = 1,total_nky; DO ikx = 1, total_nkx; DO iz = 1,total_nz
- WRITE(fid_check,*) real(field_to_check(iky,ikx,iz)), ',' , imag(field_to_check(iky,ikx,iz))
- ENDDO; ENDDO; ENDDO
- CLOSE(fid_check)
- WRITE(*,*) 'Check file found -> done.'
- ! delete the check_phi flagfile
- OPEN(fid_check, file='check_phi')
- CLOSE(fid_check, status='delete')
- ENDIF
- ENDIF
-END SUBROUTINE spit_snapshot_check
diff --git a/src/diagnostics_mod.F90 b/src/diagnostics_mod.F90
new file mode 100644
index 00000000..db03142c
--- /dev/null
+++ b/src/diagnostics_mod.F90
@@ -0,0 +1,684 @@
+MODULE diagnostics
+ ! Module for diagnostic parameters
+
+ USE prec_const
+ IMPLICIT NONE
+ PRIVATE
+
+ LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision = .FALSE.
+ LOGICAL, PUBLIC, PROTECTED :: write_gamma, write_hf ! output particle transport and heat flux
+ LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00
+ LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj
+ LOGICAL, PUBLIC, PROTECTED :: write_dens, write_fvel, write_temp
+
+ INTEGER, PUBLIC, PROTECTED :: nsave_0d, nsave_1d, nsave_2d, nsave_3d, nsave_5d ! save data every n step
+ REAL, PUBLIC, PROTECTED :: dtsave_0d, dtsave_1d, dtsave_2d, dtsave_3d, dtsave_5d ! save data every dt time unit
+ ! Change diagnostic mode (full/txtonly)
+ CHARACTER(len=256), PUBLIC, PROTECTED :: diag_mode = "full"
+ ! HDF5 file
+ CHARACTER(len=256), PUBLIC :: resfile,resfile0 = "outputs" ! Head of main result file name
+ CHARACTER(len=256), PUBLIC :: momfile,momfile0 = "moments" ! Head of the moment spectrum file (N_a(p,j,z))
+ CHARACTER(len=256), PUBLIC :: mspfile,mspfile0 = "moments_spectrum" ! Head of the moment spectrum file (N_a(p,j,z))
+ CHARACTER(len=256), PUBLIC :: fldfile,fldfile0 = "fields" ! Head of field (phi,A)
+ CHARACTER(len=256), PUBLIC :: ttrfile,ttrfile0 = "time_traces" ! Head of time traces (gamma_x,Q_x)
+ CHARACTER(len=256), PUBLIC :: ggmfile,ggmfile0 = "grid_geometry" ! Head of time traces (gamma_x,Q_x)
+ CHARACTER(len=256), PUBLIC :: prffile,prffile0 = "profiler" ! Head of time traces (gamma_x,Q_x)
+ CHARACTER(len=256), PUBLIC :: input_fname
+ INTEGER, PUBLIC :: fidres,fidmsp,fidfld,fidttr ! FID for output
+ INTEGER, PUBLIC :: fidmom,fidggm, fidprf
+
+ PUBLIC :: diag_readinputs, diag_outputinputs, diagnose, spit_snapshot_check
+
+CONTAINS
+
+ SUBROUTINE diag_readinputs
+ ! Read the input parameters
+
+ USE basic, ONLY : lu_in, dt
+ USE prec_const
+ IMPLICIT NONE
+
+ NAMELIST /DIAGNOSTICS/ dtsave_0d, dtsave_1d, dtsave_2d, dtsave_3d, dtsave_5d
+ NAMELIST /DIAGNOSTICS/ write_doubleprecision, write_gamma, write_hf, write_phi
+ NAMELIST /DIAGNOSTICS/ write_Na00, write_Napj, write_Sapj
+ NAMELIST /DIAGNOSTICS/ write_dens, write_fvel, write_temp
+ NAMELIST /DIAGNOSTICS/ diag_mode
+
+ READ(lu_in,diagnostics)
+
+ ! set nsave variables from dtsave ones (time unit to steps)
+ nsave_0d = CEILING(dtsave_0d/dt)
+ nsave_1d = CEILING(dtsave_1d/dt)
+ nsave_2d = CEILING(dtsave_2d/dt)
+ nsave_3d = CEILING(dtsave_3d/dt)
+ nsave_5d = CEILING(dtsave_5d/dt)
+
+ END SUBROUTINE diag_readinputs
+
+ SUBROUTINE init_outfile(comm,file0,file,fid)
+ USE basic, ONLY: speak, jobnum, basic_outputinputs
+ USE grid, ONLY: grid_outputinputs
+ USE geometry, ONLY: geometry_outputinputs
+ USE model, ONLY: model_outputinputs
+ USE closure, ONLY: closure_outputinputs
+ USE species, ONLY: species_outputinputs
+ USE collision, ONLY: coll_outputinputs
+ USE initial, ONLY: initial_outputinputs
+ USE time_integration,ONLY: time_integration_outputinputs
+ USE parallel, ONLY: parallel_outputinputs
+ USE futils, ONLY: creatf, creatg, creatd, attach, putfile
+ IMPLICIT NONE
+ !input
+ INTEGER, INTENT(IN) :: comm
+ CHARACTER(len=256), INTENT(IN) :: file0
+ CHARACTER(len=256), INTENT(OUT) :: file
+ INTEGER, INTENT(OUT) :: fid
+ CHARACTER(len=256) :: str
+ INCLUDE 'srcinfo.h'
+ ! Writing output filename
+ WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0) ,'_',jobnum,'.h5'
+ ! 1.1 Initial run
+ ! Main output file creation
+ IF (write_doubleprecision) THEN
+ CALL creatf(file, fid, real_prec='d', mpicomm=comm)
+ ELSE
+ CALL creatf(file, fid, mpicomm=comm)
+ END IF
+ CALL speak(TRIM(file)//' created')
+ ! basic data group
+ CALL creatg(fid, "/data", "data")
+ ! File group
+ CALL creatg(fid, "/files", "files")
+ CALL attach(fid, "/files", "jobnum", jobnum)
+ ! Add the code info and parameters to the file
+ CALL creatg(fid, "/data/input", "input")
+ CALL creatd(fid, 0,(/0/),"/data/input/codeinfo",'Code Information')
+ CALL attach(fid, "/data/input/codeinfo", "version", VERSION) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "branch", BRANCH) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "author", AUTHOR) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "execdate", EXECDATE) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "host", HOST) !defined in srcinfo.h
+ CALL basic_outputinputs(fid)
+ CALL grid_outputinputs(fid)
+ CALL geometry_outputinputs(fid)
+ CALL diag_outputinputs(fid)
+ CALL model_outputinputs(fid)
+ CALL closure_outputinputs(fid)
+ CALL species_outputinputs(fid)
+ CALL coll_outputinputs(fid)
+ CALL initial_outputinputs(fid)
+ CALL time_integration_outputinputs(fid)
+ CALL parallel_outputinputs(fid)
+ ! Save STDIN (input file) of this run
+ IF(jobnum .LE. 99) THEN
+ WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
+ ELSE
+ WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum
+ END IF
+ CALL putfile(fid, TRIM(str), TRIM(input_fname),ionode=0)
+ END SUBROUTINE init_outfile
+
+ SUBROUTINE diagnose(kstep)
+ ! Diagnostics, writing simulation state to disk
+ USE basic, ONLY: lu_in, chrono_runt, cstep, dt, time, tmax, display_h_min_s
+ USE processing, ONLY: pflux_x, hflux_x
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: kstep
+ !! Basic diagnose loop for reading input file, displaying advancement and ending
+ IF ((kstep .EQ. 0)) THEN
+ INQUIRE(unit=lu_in, name=input_fname)
+ CLOSE(lu_in)
+ ENDIF
+ !! End diag
+ IF (kstep .EQ. -1) THEN
+ ! Display total run time
+ CALL display_h_min_s(chrono_runt%ttot)
+ END IF
+ !! Specific diagnostic calls
+ SELECT CASE(diag_mode)
+ CASE('full')
+ CALL diagnose_full(kstep)
+ CASE('txtonly')
+ CALL diagnose_txtonly(kstep)
+ END SELECT
+
+ ! Terminal info
+ IF ((kstep .GE. 0) .AND. (MOD(cstep, nsave_0d) == 0) .AND. (my_id .EQ. 0)) THEN
+ if(SIZE(pflux_x) .GT. 1) THEN
+ WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A8,G10.3,A8,G10.3,A)")&
+ '|t = ', time,'| Pxi = ',pflux_x(1),'| Qxi = ',hflux_x(1),'| Pxe = ',pflux_x(2),'| Qxe = ',hflux_x(2),'|'
+ else
+ WRITE(*,"(A,F8.2,A8,G10.3,A8,G10.3,A)")&
+ '|t = ', time,'| Px = ',pflux_x(1),'| Qx = ',hflux_x(1),'|'
+ endif
+ ENDIF
+ END SUBROUTINE diagnose
+
+ SUBROUTINE diagnose_full(kstep)
+ USE basic, ONLY: speak,chrono_runt,&
+ cstep,iframe0d,iframe3d,iframe5d,crashed
+ USE grid, ONLY: &
+ parray_full,pmax,jarray_full,jmax, kparray, &
+ kyarray_full,kxarray_full,zarray_full, ngz, total_nz, local_nz, ieven,&
+ local_Nky, total_nky, local_nkx, total_nkx
+ USE geometry, ONLY: gxx, gxy, gxz, gyy, gyz, gzz, &
+ hatR, hatZ, hatB, dBdx, dBdy, dBdz, Jacobian, gradz_coeff
+ USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf!, putarrnd ! Routine de merde, jamais l'utiliser
+ USE array
+ USE model, ONLY: EM
+ USE parallel, ONLY: my_id, comm0, gather_z, gather_xyz_real
+ USE collision, ONLY: coll_outputinputs
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: kstep
+ INTEGER, parameter :: BUFSIZE = 2
+ REAL(xp), DIMENSION(total_nz) :: Az_full ! full z array for metric output
+ REAL(xp), DIMENSION(total_nky,total_nkx,total_nz) :: kp_full
+ INTEGER :: rank = 0, ierr
+ INTEGER :: dims(1) = (/0/)
+ !____________________________________________________________________________
+ ! 1. Initial diagnostics
+ IF ((kstep .EQ. 0)) THEN
+ CALL init_outfile(comm0, resfile0,resfile,fidres)
+ ! Profiler time measurement
+ CALL creatg(fidres, "/profiler", "performance analysis")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs", "cumulative rhs computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_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_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_ExBshear", "cumulative ExB shear 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")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative diagnostic time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
+#ifdef TEST_SVD
+ CALL creatd(fidres, 0, (/0/), "/profiler/Tc_CLA", "cumulative total CLA computation time")
+#endif
+ ! Grid info
+ CALL creatg(fidres, "/data/grid", "Grid data")
+ CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0)
+ CALL putarr(fidres, "/data/grid/coordky", kyarray_full, "ky*rho_s0", ionode=0)
+ CALL putarr(fidres, "/data/grid/coordz", zarray_full, "z/R", ionode=0)
+ CALL putarr(fidres, "/data/grid/coordp" , parray_full, "p", ionode=0)
+ CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
+ CALL gather_xyz_real(kparray(1:local_Nky,1:local_Nkx,1:local_nz,ieven),kp_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL putarr(fidres, "/data/grid/coordkp" , kp_full, "kp", ionode=0)
+ ! Metric info
+ CALL creatg(fidres, "/data/metric", "Metric data")
+ CALL gather_z(gxx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxx", Az_full, "gxx", ionode =0)
+ CALL gather_z(gxy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxy", Az_full, "gxy", ionode =0)
+ CALL gather_z(gxz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxz", Az_full, "gxz", ionode =0)
+ CALL gather_z(gyy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gyy", Az_full, "gyy", ionode =0)
+ CALL gather_z(gyz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gyz", Az_full, "gyz", ionode =0)
+ CALL gather_z(gzz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gzz", Az_full, "gzz", ionode =0)
+ CALL gather_z(hatR((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatR", Az_full, "hatR", ionode =0)
+ CALL gather_z(hatZ((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatZ", Az_full, "hatZ", ionode =0)
+ CALL gather_z(hatB((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatB", Az_full, "hatB", ionode =0)
+ CALL gather_z(dBdx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdx", Az_full, "dBdx", ionode =0)
+ CALL gather_z(dBdy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdy", Az_full, "dBdy", ionode =0)
+ CALL gather_z(dBdz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdz", Az_full, "dBdz", ionode =0)
+ CALL gather_z(Jacobian((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/Jacobian", Az_full, "Jacobian", ionode =0)
+ CALL gather_z(gradz_coeff((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gradz_coeff", Az_full, "gradz_coeff", ionode =0)
+ ! var0d group (gyro transport)
+ IF (nsave_0d .GT. 0) THEN
+ CALL creatg(fidres, "/data/var0d", "0d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var0d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var0d/cstep", "iteration number")
+ IF (write_gamma) THEN
+ CALL creatd(fidres, rank, dims, "/data/var0d/gflux_x", "Radial gyro transport")
+ CALL creatd(fidres, rank, dims, "/data/var0d/pflux_x", "Radial part transport")
+ ENDIF
+ IF (write_hf) THEN
+ CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x", "Radial part heat flux")
+ ENDIF
+ IF (cstep==0) THEN
+ iframe0d=0
+ ENDIF
+ CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
+ END IF
+ ! var2d group (gyro transport)
+ IF (nsave_0d .GT. 0) THEN
+ CALL creatg(fidres, "/data/var2d", "2d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
+#ifdef TEST_SVD
+ CALL creatg(fidres, "/data/var2d/sv_ky_pj", "singular values of the moment ky/pj cut")
+#endif
+ ENDIF
+ ! var3d group (phi,psi, fluid moments, Ni00, Napjz)
+ IF (nsave_3d .GT. 0) THEN
+ CALL creatg(fidres, "/data/var3d", "3d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
+ IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
+ IF (write_phi.AND.EM) CALL creatg(fidres, "/data/var3d/psi", "psi")
+ IF (write_Na00) THEN
+ CALL creatg(fidres, "/data/var3d/Na00", "gyrocenter density ")
+ CALL creatg(fidres, "/data/var3d/Napjz", "pj(z) moment spectrum ")
+ ENDIF
+ IF (write_dens) THEN
+ CALL creatg(fidres, "/data/var3d/dens", "density ")
+ ENDIF
+ IF (write_fvel) THEN
+ CALL creatg(fidres, "/data/var3d/upar", "parallel fluid velocity ")
+ CALL creatg(fidres, "/data/var3d/uper", "perpendicular fluid velocity ")
+ ENDIF
+ IF (write_temp) THEN
+ CALL creatg(fidres, "/data/var3d/Tper", "perpendicular temperature ")
+ CALL creatg(fidres, "/data/var3d/Tpar", "parallel temperature ")
+ CALL creatg(fidres, "/data/var3d/temp", "tiotal temperature ")
+ ENDIF
+ IF (cstep==0) THEN
+ iframe3d=0
+ ENDIF
+ CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+ END IF
+ ! var5d group (moments)
+ IF (nsave_5d .GT. 0) THEN
+ CALL creatg(fidres, "/data/var5d", "5d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
+ IF (write_Napj) THEN
+ CALL creatg(fidres, "/data/var5d/moments", "full moments array")
+ ENDIF
+ IF (cstep==0) THEN
+ iframe5d=0
+ END IF
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+ END IF
+ ENDIF
+ !_____________________________________________________________________________
+ ! 2. Periodic diagnostics
+ !
+ IF (kstep .GE. 0) THEN
+ ! 2.1 0d history arrays
+ IF (nsave_0d .GT. 0) THEN
+ IF ( MOD(cstep, nsave_0d) == 0 ) THEN
+ CALL diagnose_0d
+ END IF
+ END IF
+ ! 2.2 2d profiles
+ IF (nsave_2d .GT. 0) THEN
+ IF (MOD(cstep, nsave_2d) == 0) THEN
+ CALL diagnose_2d
+ ENDIF
+ ENDIF
+ ! 2.3 3d profiles
+ IF (nsave_3d .GT. 0) THEN
+ IF (MOD(cstep, nsave_3d) == 0) THEN
+ CALL diagnose_3d
+ ! Looks at the folder if the file check_phi exists and spits a snapshot
+ ! of the current electrostatic potential in a basic text file
+ CALL spit_snapshot_check
+ ENDIF
+ ENDIF
+ ! 2.4 5d profiles
+ IF (nsave_5d .GT. 0) THEN
+ IF (MOD(cstep, nsave_5d) == 0) THEN
+ CALL diagnose_5d
+ END IF
+ END IF
+ !_____________________________________________________________________________
+ ! 3. Final diagnostics
+ ELSEIF (kstep .EQ. -1) THEN
+ 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'
+ IF (nsave_5d .GT. 0) CALL diagnose_5d
+ ENDIF
+ ! Close all diagnostic files
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+ CALL closef(fidres)
+ END IF
+ END SUBROUTINE diagnose_full
+
+ !! This routine outputs only txt file 0D data (flux and time)
+ SUBROUTINE diagnose_txtonly(kstep)
+ USE basic
+ USE processing, ONLY: pflux_x, hflux_x, compute_radial_transport, compute_radial_heatflux
+ USE parallel, ONLY: my_id, comm0
+ USE futils, ONLY: creatf, creatg, creatd, attach, putfile, closef
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: kstep
+ INTEGER, parameter :: BUFSIZE = 2
+ INTEGER :: rank = 0
+ INTEGER :: dims(1) = (/0/)
+ IF (kstep .GE. 0) THEN
+ ! output the transport in a txt file
+ IF ( MOD(cstep, nsave_0d) == 0 ) THEN
+ CALL compute_radial_transport
+ CALL compute_radial_heatflux
+ IF (my_id .EQ. 0) &
+ WRITE(1,*) time, pflux_x(1), hflux_x(1)
+ END IF
+ !! Save the last state
+ ELSEIF (kstep .EQ. -1) THEN
+ CALL init_outfile(comm0, resfile0,resfile,fidres)
+ CALL creatg(fidres, "/data/var5d", "5d profiles")
+ CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
+ CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
+ CALL creatg(fidres, "/data/var5d/moments", "full moments array")
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+ CALL diagnose_5d
+ CALL closef(fidres)
+ ENDIF
+ END SUBROUTINE diagnose_txtonly
+
+ !!-------------- Auxiliary routines -----------------!!
+ SUBROUTINE diagnose_0d
+ USE basic
+ USE futils, ONLY: append, attach, getatt, creatd
+ USE prec_const
+ USE processing
+ USE model, ONLY: Na
+ IMPLICIT NONE
+ ! Time measurement data
+ CALL append(fidres, "/profiler/Tc_rhs", REAL(chrono_mrhs%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_adv_field", REAL(chrono_advf%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_clos", REAL(chrono_clos%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_ghost", REAL(chrono_ghst%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_coll", REAL(chrono_coll%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_poisson", REAL(chrono_pois%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_Sapj", REAL(chrono_sapj%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_checkfield",REAL(chrono_chck%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_diag", REAL(chrono_diag%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_grad", REAL(chrono_grad%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_nadiab", REAL(chrono_napj%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_ExBshear", REAL(chrono_ExBs%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_step", REAL(chrono_step%ttot,dp),ionode=0)
+#ifdef TEST_SVD
+ CALL append(fidres, "/profiler/Tc_CLA", REAL(chrono_CLA%ttot,dp),ionode=0)
+#endif
+ CALL append(fidres, "/profiler/time", REAL(time,dp),ionode=0)
+ ! Processing data
+ CALL append(fidres, "/data/var0d/time", REAL(time,dp),ionode=0)
+ CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var0d/", "frames",iframe0d)
+ iframe0d=iframe0d+1
+ CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
+ ! Ion transport data
+ IF (write_gamma) THEN
+ CALL compute_radial_transport
+ CALL append(fidres, "/data/var0d/gflux_x",REAL(gflux_x,dp),ionode=0)
+ CALL append(fidres, "/data/var0d/pflux_x",REAL(pflux_x,dp),ionode=0)
+ ENDIF
+ IF (write_hf) THEN
+ CALL compute_radial_heatflux
+ CALL append(fidres, "/data/var0d/hflux_x",REAL(hflux_x,dp),ionode=0)
+ ENDIF
+ END SUBROUTINE diagnose_0d
+
+ SUBROUTINE diagnose_2d
+ USE prec_const
+ USE basic
+ USE futils, ONLY: putarr, append
+#ifdef TEST_SVD
+ USE CLA, ONLY: Sf
+#endif
+ IMPLICIT NONE
+ CHARACTER(50) :: dset_name
+ iframe2d=iframe2d+1
+ CALL append(fidres,"/data/var2d/time", REAL(time,dp), ionode=0)
+ CALL append(fidres,"/data/var2d/cstep",REAL(cstep,dp),ionode=0)
+#ifdef TEST_SVD
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var2d/sv_ky_pj/", iframe2d
+ CALL putarr(fidres, dset_name, Sf, ionode=0)
+#endif
+ END SUBROUTINE
+
+ SUBROUTINE diagnose_3d
+ USE basic
+ USE futils, ONLY: append, getatt, attach, putarr!, putarrnd
+ USE fields, ONLY: phi, psi, moments
+ USE array, ONLY: Napjz,dens,upar,uper,Tpar,Tper,temp
+ USE grid, ONLY: CONTAINSp0, ip0,ij0, local_na, total_na,&
+ total_np, total_nj, total_nky, total_nkx, total_nz, &
+ local_np, local_nky, local_nkx, local_nz, &
+ ngz
+ USE time_integration, ONLY: updatetlevel
+ USE prec_const
+ USE processing, ONLY: compute_fluid_moments, compute_Napjz_spectrum
+ USE model, ONLY: EM
+ USE parallel, ONLY: manual_3D_bcast
+ IMPLICIT NONE
+ COMPLEX(xp), DIMENSION(local_na,local_nky,local_nkx,local_nz) :: Na00_
+ ! add current time, cstep and frame
+ CALL append(fidres, "/data/var3d/time", REAL(time,dp),ionode=0)
+ CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
+ iframe3d=iframe3d+1
+ CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+ ! Write current EM fields
+ IF (write_phi) CALL write_field3d_kykxz(phi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'phi')
+ IF (write_phi.AND.EM) CALL write_field3d_kykxz(psi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'psi')
+ IF (write_Na00) THEN
+ CALL compute_Napjz_spectrum
+ IF (CONTAINSp0) THEN
+ ! gyrocenter density
+ Na00_ = moments(:,ip0,ij0,:,:,(1+ngz/2):(local_nz+ngz/2),updatetlevel)
+ ELSE
+ Na00_ = 0._xp
+ ENDIF
+ CALL write_field3da_kykxz(Na00_, 'Na00')
+ ! <<Napj>x>y spectrum
+ CALL write_field3da_pjz(Napjz, 'Napjz')
+ ENDIF
+ !! Fuid moments
+ IF (write_dens .OR. write_fvel .OR. write_temp) &
+ CALL compute_fluid_moments
+ IF (write_dens) THEN
+ CALL write_field3da_kykxz(dens, 'dens')
+ ENDIF
+ IF (write_fvel) THEN
+ CALL write_field3da_kykxz(upar, 'upar')
+ CALL write_field3da_kykxz(uper, 'uper')
+ ENDIF
+ IF (write_temp) THEN
+ CALL write_field3da_kykxz(Tpar, 'Tpar')
+ CALL write_field3da_kykxz(Tper, 'Tper')
+ CALL write_field3da_kykxz(temp, 'temp')
+ ENDIF
+ CONTAINS
+ SUBROUTINE write_field3d_kykxz(field, text)
+ USE parallel, ONLY : gather_xyz, num_procs
+ IMPLICIT NONE
+ COMPLEX(xp), DIMENSION(:,:,:), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+ COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_full
+ CHARACTER(50) :: dset_name
+ field_full = 0;
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+ IF (num_procs .EQ. 1) THEN ! no data distribution
+ CALL putarr(fidres, dset_name, field, ionode=0)
+ ELSE
+ CALL gather_xyz(field,field_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL putarr(fidres, dset_name, field_full, ionode=0)
+ ENDIF
+ END SUBROUTINE write_field3d_kykxz
+
+ SUBROUTINE write_field3da_kykxz(field, text)
+ USE parallel, ONLY : gather_xyz, num_procs
+ IMPLICIT NONE
+ COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+ COMPLEX(xp), DIMENSION(total_na,total_nky,total_nkx,total_nz) :: field_full
+ COMPLEX(xp), DIMENSION(local_nky,total_nkx,local_nz) :: buff_local
+ COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: buff_full
+ CHARACTER(50) :: dset_name
+ INTEGER :: ia
+ field_full = 0;
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+ IF (num_procs .EQ. 1) THEN ! no data distribution
+ CALL putarr(fidres, dset_name, field, ionode=0)
+ ELSE
+ DO ia = 1,total_Na
+ buff_local = field(ia,:,:,:)
+ CALL gather_xyz(buff_local,buff_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ field_full(ia,:,:,:) = buff_full
+ ENDDO
+ CALL putarr(fidres, dset_name, field_full, ionode=0)
+ ENDIF
+ END SUBROUTINE write_field3da_kykxz
+
+ SUBROUTINE write_field3da_pjz(field, text)
+ USE parallel, ONLY : gather_pjz, num_procs
+ IMPLICIT NONE
+ COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+ COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nz) :: field_full
+ CHARACTER(LEN=50) :: dset_name
+ field_full = 0;
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+ IF (num_procs .EQ. 1) THEN ! no data distribution
+ CALL putarr(fidres, dset_name, field, ionode=0)
+ ELSE
+ CALL gather_pjz(field,field_full,total_na,local_np,total_np,total_nj,local_nz,total_nz)
+ CALL putarr(fidres, dset_name, field_full, ionode=0)
+ ENDIF
+ CALL attach(fidres, dset_name, "time", time)
+ END SUBROUTINE write_field3da_pjz
+
+ END SUBROUTINE diagnose_3d
+
+ SUBROUTINE diagnose_5d
+
+ USE basic, ONLY: time, iframe5d,cstep
+ USE futils, ONLY: append, getatt, attach, putarr !,putarrnd
+ USE fields, ONLY: moments
+ USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
+ local_np, local_nj, local_nky, local_nkx, local_nz, &
+ ngp, ngj, ngz, total_na
+ USE prec_const, ONLY: xp,dp
+ IMPLICIT NONE
+
+ CALL append(fidres, "/data/var5d/time", REAL(time,dp),ionode=0)
+ CALL append(fidres, "/data/var5d/cstep", REAL(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
+ iframe5d=iframe5d+1
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+
+ IF (write_Napj) THEN
+ CALL write_field5d(moments, 'moments')
+ ENDIF
+
+ CONTAINS
+
+ SUBROUTINE write_field5d(field, text)
+ USE basic, ONLY: jobnum, dt
+ USE futils, ONLY: attach, putarr!, putarrnd
+ USE parallel, ONLY: gather_pjxyz, num_procs
+ USE prec_const, ONLY: xp
+ USE time_integration, ONLY: ntimelevel
+ IMPLICIT NONE
+ COMPLEX(xp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_Nky,local_Nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
+ CHARACTER(*), INTENT(IN) :: text
+ COMPLEX(xp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
+ COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
+ CHARACTER(LEN=50) :: dset_name
+ field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
+ 1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),1)
+ field_full = 0;
+ WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
+ IF (num_procs .EQ. 1) THEN
+ CALL putarr(fidres, dset_name, field_sub, ionode=0)
+ ELSE!IF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
+ CALL gather_pjxyz(field_sub,field_full,total_na,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL putarr(fidres, dset_name, field_full, ionode=0)
+ ! ELSE
+ ! CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/))
+ ENDIF
+ CALL attach(fidres, dset_name, 'cstep', cstep)
+ CALL attach(fidres, dset_name, 'time', time)
+ CALL attach(fidres, dset_name, 'jobnum', jobnum)
+ CALL attach(fidres, dset_name, 'dt', dt)
+ CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
+ END SUBROUTINE write_field5d
+ END SUBROUTINE diagnose_5d
+
+ SUBROUTINE spit_snapshot_check
+ USE fields, ONLY: phi
+ USE grid, ONLY: total_nkx,total_nky,total_nz,&
+ local_nky,local_nz, ngz
+ USE parallel, ONLY: gather_xyz, my_id
+ USE basic
+ USE prec_const, ONLY: xp
+ IMPLICIT NONE
+ LOGICAL :: file_exist
+ INTEGER :: fid_check, ikx, iky, iz
+ CHARACTER(256) :: check_filename
+ COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_to_check
+ !! Spit a snapshot of PHI if requested (triggered by creating a file named "check_phi")
+ INQUIRE(file='check_phi', exist=file_exist)
+ IF( file_exist ) THEN
+ IF(my_id.EQ. 0) WRITE(*,*) 'Check file found -> gather phi..'
+ CALL gather_xyz(phi(:,:,(1+Ngz/2):(local_nz+Ngz/2)), field_to_check,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ IF(my_id.EQ. 0) THEN
+ WRITE(check_filename,'(a16)') 'check_phi.out'
+ OPEN(fid_check, file=check_filename, form='formatted')
+ WRITE(*,*) 'Check file found -> output phi ..'
+ WRITE(fid_check,*) total_nky, total_nkx, total_nz
+ DO iky = 1,total_nky; DO ikx = 1, total_nkx; DO iz = 1,total_nz
+ WRITE(fid_check,*) real(field_to_check(iky,ikx,iz)), ',' , imag(field_to_check(iky,ikx,iz))
+ ENDDO; ENDDO; ENDDO
+ CLOSE(fid_check)
+ WRITE(*,*) 'Check file found -> done.'
+ ! delete the check_phi flagfile
+ OPEN(fid_check, file='check_phi')
+ CLOSE(fid_check, status='delete')
+ ENDIF
+ ENDIF
+ END SUBROUTINE spit_snapshot_check
+
+
+
+ SUBROUTINE diag_outputinputs(fid)
+ !
+ ! Write the input parameters to the results_xx.h5 file
+ !
+ USE prec_const
+ USE futils, ONLY: attach, creatd
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: fid
+ CHARACTER(len=256) :: str
+ WRITE(str,'(a)') '/data/input/diagnostics'
+ CALL creatd(fidres, 0,(/0/),TRIM(str),'Diagnostics Parameters Input')
+ CALL attach(fid, TRIM(str), "write_doubleprecision", write_doubleprecision)
+ CALL attach(fid, TRIM(str), "nsave_0d", nsave_0d)
+ CALL attach(fid, TRIM(str), "nsave_1d", nsave_1d)
+ CALL attach(fid, TRIM(str), "nsave_2d", nsave_2d)
+ CALL attach(fid, TRIM(str), "nsave_3d", nsave_3d)
+ CALL attach(fid, TRIM(str), "nsave_5d", nsave_5d)
+ CALL attach(fid, TRIM(str), "write_gamma", write_gamma)
+ CALL attach(fid, TRIM(str), "write_hf", write_hf)
+ CALL attach(fid, TRIM(str), "write_phi", write_phi)
+ CALL attach(fid, TRIM(str), "write_Na00", write_Na00)
+ CALL attach(fid, TRIM(str), "write_Napj", write_Napj)
+ CALL attach(fid, TRIM(str), "write_Sapj", write_Sapj)
+ CALL attach(fid, TRIM(str), "write_dens", write_dens)
+ CALL attach(fid, TRIM(str), "write_fvel", write_fvel)
+ CALL attach(fid, TRIM(str), "write_temp", write_temp)
+
+ END SUBROUTINE diag_outputinputs
+
+
+END MODULE diagnostics
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
deleted file mode 100644
index 2a8630f8..00000000
--- a/src/diagnostics_par_mod.F90
+++ /dev/null
@@ -1,92 +0,0 @@
-MODULE diagnostics_par
- ! Module for diagnostic parameters
-
- USE prec_const
- IMPLICIT NONE
- PRIVATE
-
- LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision = .FALSE.
- LOGICAL, PUBLIC, PROTECTED :: write_gamma, write_hf ! output particle transport and heat flux
- LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00
- LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj
- LOGICAL, PUBLIC, PROTECTED :: write_dens, write_fvel, write_temp
-
- INTEGER, PUBLIC, PROTECTED :: nsave_0d, nsave_1d, nsave_2d, nsave_3d, nsave_5d ! save data every n step
- REAL, PUBLIC, PROTECTED :: dtsave_0d, dtsave_1d, dtsave_2d, dtsave_3d, dtsave_5d ! save data every dt time unit
- ! Change diagnostic mode (full/txtonly)
- CHARACTER(len=256), PUBLIC, PROTECTED :: diag_mode = "full"
- ! HDF5 file
- CHARACTER(len=256), PUBLIC :: resfile,resfile0 = "outputs" ! Head of main result file name
- CHARACTER(len=256), PUBLIC :: momfile,momfile0 = "moments" ! Head of the moment spectrum file (N_a(p,j,z))
- CHARACTER(len=256), PUBLIC :: mspfile,mspfile0 = "moments_spectrum" ! Head of the moment spectrum file (N_a(p,j,z))
- CHARACTER(len=256), PUBLIC :: fldfile,fldfile0 = "fields" ! Head of field (phi,A)
- CHARACTER(len=256), PUBLIC :: ttrfile,ttrfile0 = "time_traces" ! Head of time traces (gamma_x,Q_x)
- CHARACTER(len=256), PUBLIC :: ggmfile,ggmfile0 = "grid_geometry" ! Head of time traces (gamma_x,Q_x)
- CHARACTER(len=256), PUBLIC :: prffile,prffile0 = "profiler" ! Head of time traces (gamma_x,Q_x)
- CHARACTER(len=256), PUBLIC :: input_fname
- CHARACTER(len=256), PUBLIC :: rstfile ! restart result file
- INTEGER, PUBLIC :: fidres,fidmsp,fidfld,fidttr ! FID for output
- INTEGER, PUBLIC :: fidmom,fidggm, fidprf
- INTEGER, PUBLIC :: fidrst ! FID for restart file
-
- PUBLIC :: diag_par_readinputs, diag_par_outputinputs
-
-CONTAINS
-
-
- SUBROUTINE diag_par_readinputs
- ! Read the input parameters
-
- USE basic, ONLY : lu_in, dt
- USE prec_const
- IMPLICIT NONE
-
- NAMELIST /OUTPUT_PAR/ dtsave_0d, dtsave_1d, dtsave_2d, dtsave_3d, dtsave_5d
- NAMELIST /OUTPUT_PAR/ write_doubleprecision, write_gamma, write_hf, write_phi
- NAMELIST /OUTPUT_PAR/ write_Na00, write_Napj, write_Sapj
- NAMELIST /OUTPUT_PAR/ write_dens, write_fvel, write_temp
- NAMELIST /OUTPUT_PAR/ diag_mode
-
- READ(lu_in,output_par)
-
- ! set nsave variables from dtsave ones (time unit to steps)
- nsave_0d = CEILING(dtsave_0d/dt)
- nsave_1d = CEILING(dtsave_1d/dt)
- nsave_2d = CEILING(dtsave_2d/dt)
- nsave_3d = CEILING(dtsave_3d/dt)
- nsave_5d = CEILING(dtsave_5d/dt)
-
- END SUBROUTINE diag_par_readinputs
-
-
- SUBROUTINE diag_par_outputinputs(fid)
- !
- ! Write the input parameters to the results_xx.h5 file
- !
- USE prec_const
- USE futils, ONLY: attach, creatd
- IMPLICIT NONE
- INTEGER, INTENT(in) :: fid
- CHARACTER(len=256) :: str
- WRITE(str,'(a)') '/data/input/diag_par'
- CALL creatd(fidres, 0,(/0/),TRIM(str),'Diagnostics Parameters Input')
- CALL attach(fid, TRIM(str), "write_doubleprecision", write_doubleprecision)
- CALL attach(fid, TRIM(str), "nsave_0d", nsave_0d)
- CALL attach(fid, TRIM(str), "nsave_1d", nsave_1d)
- CALL attach(fid, TRIM(str), "nsave_2d", nsave_2d)
- CALL attach(fid, TRIM(str), "nsave_3d", nsave_3d)
- CALL attach(fid, TRIM(str), "nsave_5d", nsave_5d)
- CALL attach(fid, TRIM(str), "write_gamma", write_gamma)
- CALL attach(fid, TRIM(str), "write_hf", write_hf)
- CALL attach(fid, TRIM(str), "write_phi", write_phi)
- CALL attach(fid, TRIM(str), "write_Na00", write_Na00)
- CALL attach(fid, TRIM(str), "write_Napj", write_Napj)
- CALL attach(fid, TRIM(str), "write_Sapj", write_Sapj)
- CALL attach(fid, TRIM(str), "write_dens", write_dens)
- CALL attach(fid, TRIM(str), "write_fvel", write_fvel)
- CALL attach(fid, TRIM(str), "write_temp", write_temp)
-
- END SUBROUTINE diag_par_outputinputs
-
-
-END MODULE diagnostics_par
diff --git a/src/inital.F90 b/src/inital.F90
deleted file mode 100644
index d426e5de..00000000
--- a/src/inital.F90
+++ /dev/null
@@ -1,630 +0,0 @@
-!******************************************************************************!
-!!!!!! initialize the moments and load/build coeff tables
-!******************************************************************************!
-SUBROUTINE inital
-
- USE basic, ONLY: speak
- USE initial_par, ONLY: INIT_OPT
- USE time_integration, ONLY: set_updatetlevel
- USE collision, ONLY: init_collision
- USE closure, ONLY: apply_closure_model
- USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
- USE restarts, ONLY: load_moments, job2load
- USE processing, ONLY: compute_fluid_moments, compute_radial_heatflux, compute_radial_transport
- USE model, ONLY: LINEARITY
- USE nonlinear, ONLY: compute_Sapj, nonlinear_init
- IMPLICIT NONE
-
- CALL set_updatetlevel(1)
-
- !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
- ! through loading a previous state
- IF ( job2load .GE. 0 ) THEN
- CALL speak('Load moments')
- CALL load_moments ! get N_0
- CALL apply_closure_model
- CALL update_ghosts_moments
- CALL solve_EM_fields ! compute phi_0=phi(N_0)
- CALL update_ghosts_EM
- ! through initialization
- ELSE
- SELECT CASE (INIT_OPT)
- ! set phi with noise and set moments to 0
- CASE ('phi')
- CALL speak('Init noisy phi')
- CALL init_phi
- CALL update_ghosts_EM
- CASE ('phi_ppj')
- CALL speak('Init noisy phi')
- CALL init_phi_ppj
- CALL update_ghosts_EM
- ! set moments_00 (GC density) with noise and compute phi afterwards
- CASE('mom00')
- CALL speak('Init noisy gyrocenter density')
- CALL init_gyrodens ! init only gyrocenter density
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- ! set moments_00 (GC density) with a single kx0,ky0 mode
- ! kx0 is setup but by the noise value and ky0 by the background value
- CASE('mom00_single_mode')
- CALL speak('Init noisy gyrocenter density')
- CALL init_single_mode ! init single mode in gyrocenter density
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- ! init all moments randomly (unadvised)
- CASE('allmom')
- CALL speak('Init noisy moments')
- CALL init_moments ! init all moments
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- ! init a gaussian blob in gyrodens
- CASE('blob')
- CALL speak('--init a blob')
- CALL initialize_blob
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- ! init moments 00 with a power law similarly to GENE
- CASE('ppj')
- CALL speak('ppj init ~ GENE')
- call init_ppj
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- CASE('ricci')
- CALL speak('Init Ricci')
- CALL init_ricci ! init only gyrocenter density
- CALL update_ghosts_moments
- CALL solve_EM_fields
- CALL update_ghosts_EM
- END SELECT
- ENDIF
- ! closure of j>J, p>P and j<0, p<0 moments
- CALL speak('Apply closure')
- CALL apply_closure_model
- ! ghosts for p parallelization
- CALL speak('Ghosts communication')
- CALL update_ghosts_moments
- CALL update_ghosts_EM
- !! End of phi and moments initialization
-
- ! Init collision operator
- CALL init_collision
-
- !! Preparing auxiliary arrays at initial state
- ! particle density, fluid velocity and temperature (used in diagnose)
- CALL speak('Computing fluid moments and transport')
- CALL compute_fluid_moments
- CALL compute_radial_transport
- CALL compute_radial_heatflux
-
- ! init auxval for nonlinear
- CALL nonlinear_init
- ! compute nonlinear for t=0 diagnostic
- CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
-
-
-END SUBROUTINE inital
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize all the moments randomly
-!******************************************************************************!
-SUBROUTINE init_moments
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- ngp, ngj, ngz, iky0, contains_ky0, AA_x, AA_y
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp
- USE model, ONLY: LINEARITY
- USE parallel, ONLY: my_id
- IMPLICIT NONE
-
- REAL(xp) :: noise
- INTEGER, DIMENSION(12) :: iseedarr
- INTEGER :: ia,ip,ij,ikx,iky,iz, ipi,iji,izi
-
- ! Seed random number generator
- iseedarr(:)=iseed
- CALL RANDOM_SEED(PUT=iseedarr+my_id)
-
- !**** Broad noise initialization *******************************************
- DO ia=1,local_na
- DO ip=1,local_np
- ipi = ip+ngp/2
- DO ij=1,local_nj
- iji = ij+ngj/2
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz
- izi = iz+ngz/2
- CALL RANDOM_NUMBER(noise)
- moments(ia,ipi,iji,iky,ikx,izi,:) = (init_background + init_noiselvl*(noise-0.5_xp))
- END DO
- END DO
- END DO
- IF ( contains_ky0 ) THEN
- DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
- moments(ia,ipi,iji,iky0,ikx,:,:) = moments(ia,ipi,iji,iky0,total_nkx+2-ikx,:,:)
- END DO
- ENDIF
- END DO
- END DO
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz
- izi = iz+ngz/2
- DO ip=1,local_np
- ipi = ip+ngp/2
- 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
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE init_moments
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize the gyrocenter density randomly
-!******************************************************************************!
-SUBROUTINE init_gyrodens
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- ngp, ngj, ngz, iky0, parray, jarray, contains_ky0, AA_x, AA_y
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE parallel, ONLY: my_id
- IMPLICIT NONE
-
- REAL(xp) :: noise
- INTEGER :: ia,ip,ij,ikx,iky,iz
- INTEGER, DIMENSION(12) :: iseedarr
-
- ! Seed random number generator
- iseedarr(:)=iseed
- CALL RANDOM_SEED(PUT=iseedarr+my_id)
- moments = 0._xp
- !**** Broad noise initialization *******************************************
- DO ia=1,local_na
- DO ip=1+ngp/2,local_np+ngp/2
- DO ij=1+ngj/2,local_nj+ngj/2
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1+ngz/2,local_nz+ngz/2
- CALL RANDOM_NUMBER(noise)
- IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_xp))
- ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
- ENDIF
- END DO
- END DO
- END DO
- IF ( contains_ky0 ) THEN
- DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
- moments(ia, ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:,:)
- END DO
- ENDIF
- END DO
- END DO
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE init_gyrodens
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize a single mode in the gyrocenter density
-!******************************************************************************!
-SUBROUTINE init_single_mode
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- ngp, ngj, ngz, iky0, kxarray, kyarray, parray, jarray, &
- contains_ky0, deltakx, deltaky
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE parallel, ONLY: my_id
- IMPLICIT NONE
-
- REAL(xp) :: amplitude
- INTEGER :: ia,ip,ij,ikx,iky,iz, ikxi, ikyi
- INTEGER, DIMENSION(12) :: iseedarr
-
- moments = 0._xp
- amplitude = 1._xp
- ! find the closest modes to be intialized
- ikxi = NINT(init_noiselvl);
- ikyi = NINT(init_background);
- !****single mode initialization *******************************************
- ! DO ia=1,local_na
- ! DO ip=1+ngp/2,local_np+ngp/2
- ! DO ij=1+ngj/2,local_nj+ngj/2
- ! DO iz=1+ngz/2,local_nz+ngz/2
- ! IF ( ((parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0)) &
- ! .AND. ((ikxi .NE. -1) .AND. (ikyi .NE. -1)) ) &
- ! moments(ia,ip,ij,ikyi,ikxi,iz,:) = amplitude
- ! END DO
- ! IF ( contains_ky0 ) THEN
- ! DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
- ! moments(ia, ip,ij,iky0,total_nkx+2-ikx,:,:) = moments(ia, ip,ij,iky0,ikx,:,:)
- ! END DO
- ! ENDIF
- ! END DO
- ! END DO
- ! ENDDO
- moments(:,:,:,1,1,:,:) = amplitude
- moments(:,:,:,2,2,:,:) = amplitude
- moments(:,:,:,3,3,:,:) = amplitude
- moments(:,:,:,4,4,:,:) = amplitude
- moments(:,:,:,5,5,:,:) = amplitude
-END SUBROUTINE init_single_mode
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize a noisy ES potential and cancel the moments
-!******************************************************************************!
-SUBROUTINE init_phi
- USE grid, ONLY: total_nkx, local_nky, local_nz,&
- ngz, iky0, ikx0, contains_ky0, AA_x, AA_y
- USE fields, ONLY: phi, moments
- USE prec_const, ONLY: xp
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE parallel, ONLY: my_id
- IMPLICIT NONE
- REAL(xp) :: noise
- INTEGER, DIMENSION(12) :: iseedarr
- INTEGER :: iky,ikx,iz
- ! Seed random number generator
- iseedarr(:)=iseed
- CALL RANDOM_SEED(PUT=iseedarr+my_id)
- !**** noise initialization *******************************************
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- CALL RANDOM_NUMBER(noise)
- phi(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_xp))!*AA_x(ikx)*AA_y(iky)
- ENDDO
- END DO
- END DO
- !symmetry at ky = 0 to keep real inverse transform
- IF ( contains_ky0 ) THEN
- DO iz=1,local_nz+ngz
- DO ikx=2,total_nkx/2
- phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
- ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
- END DO
- ENDIF
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- phi(iky,ikx,iz) = phi(iky,ikx,iz)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- !**** ensure no previous moments initialization
- moments = 0._xp
- !**** Zonal Flow initialization *******************************************
- ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
- ! IF(INIT_ZF .GT. 0) THEN
- ! IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
- ! IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
- ! DO ia=1,local_na
- ! DO iz = 1,local_nz+ngz
- ! phi(iky0,INIT_ZF+1,iz) = ZF_AMP*(2._xp*PI)**2/deltakx/deltaky/2._xp * COS((iz-1)/Nz*2._xp*PI)
- ! moments(ia,ip0,ij0,iky0,INIT_ZF+1,iz,:) = kxarray(INIT_ZF+1)**2*phi(iky0,INIT_ZF+1,iz)* COS((iz-1)/Nz*2._xp*PI)
- ! ENDDO
- ! ENDDO
- ! ENDIF
- ! ENDIF
-END SUBROUTINE init_phi
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize a ppj ES potential and cancel the moments
-!******************************************************************************!
-SUBROUTINE init_phi_ppj
- USE grid, ONLY: total_nkx, local_nky, local_nz, AA_x, AA_y,&
- ngz, iky0, ikx0, contains_ky0, ieven, kxarray, kyarray, zarray, deltakx
- USE fields, ONLY: phi, moments
- USE prec_const, ONLY: xp
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE geometry, ONLY: Jacobian, iInt_Jacobian
- IMPLICIT NONE
- REAL(xp) :: kx, ky, z, amp
- INTEGER :: ikx, iky, iz
- amp = 1.0_xp
- !**** ppj initialization *******************************************
- DO iky=1,local_nky
- ky = kyarray(iky)
- DO ikx=1,total_nkx
- kx = kxarray(iky,ikx)
- DO iz=1,local_nz+ngz
- z = zarray(iz,ieven)
- IF (ky .NE. 0) THEN
- phi(iky,ikx,iz) = 0._xp
- ELSE
- phi(iky,ikx,iz) = 0.5_xp*amp*(deltakx/(ABS(kx)+deltakx))
- ENDIF
- ! z-dep and noise
- phi(iky,ikx,iz) = phi(iky,ikx,iz) * &
- (Jacobian(iz,ieven)*iInt_Jacobian)**2
- END DO
- END DO
- END DO
- !symmetry at ky = 0 to keep real inverse transform
- IF ( contains_ky0 ) THEN
- DO iz=1,local_nz+ngz
- DO ikx=2,total_nkx/2
- phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
- ENDDO
- phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
- END DO
- ENDIF
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- phi(iky,ikx,iz) = phi(iky,ikx,iz)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- !**** ensure no previous moments initialization
- moments = 0._xp
-END SUBROUTINE init_phi_ppj
-!******************************************************************************!
-
-!******************************************************************************!
-!******************************************************************************!
-!!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
-!******************************************************************************!
-SUBROUTINE initialize_blob
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz, total_nz,&
- AA_x, AA_y, parray, jarray,&
- ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
- IMPLICIT NONE
- REAL(xp) ::kx, ky, z, sigma_x, sigma_y, gain
- INTEGER :: ia,iky,ikx,iz,ip,ij, p, j
- sigma_y = 0.5
- sigma_x = sigma_y
- gain = 1.0
- ! One can increase the gain if we run 3D sim
- IF(total_nz .GT. 1) THEN
- sigma_y = 1.0
- sigma_x = sigma_y
- gain = 10.0
- ENDIF
-
- DO ia=1,local_na
- DO iky=1,local_nky
- ky = kyarray(iky)
- DO iz=1+ngz/2,local_nz+ngz/2
- z = zarray(iz,ieven)
- DO ikx=1,total_nkx
- kx = kxarray(iky,ikx) + z*shear*ky
- DO ip=1+ngp/2,local_np+ngp/2
- p = parray(ip)
- DO ij=1+ngj/2,local_nj+ngj/2
- j = jarray(ij)
- IF( (iky .NE. iky0) .AND. (p .EQ. 0) .AND. (j .EQ. 0)) THEN
- moments(ia,ip,ij,iky,ikx,iz, :) = moments(ia,ip,ij,iky,ikx,iz, :) &
- + gain * exp(-((kx/sigma_x)**2+(ky/sigma_y)**2)) &
- * AA_x(ikx)*AA_y(iky)* &
- (Jacobian(iz,ieven)*iInt_Jacobian)**2!&
- ! * exp(sigmai2_taui_o2*(kx**2+ky**2))
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE initialize_blob
-!******************************************************************************!
-!******************************************************************************!
-!!!!!!! Initialize the gyrocenter in a ppj gene manner (power law)
-!******************************************************************************!
-SUBROUTINE init_ppj
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- AA_x, AA_y, deltakx, deltaky,contains_ky0,&
- ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp, pi
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
- IMPLICIT NONE
- REAL(xp) :: kx, ky, sigma_z, amp, z
- INTEGER :: ia,iky,ikx,iz,ip,ij
- sigma_z = pi/4._xp
- amp = 1.0_xp
- !**** Broad noise initialization *******************************************
- DO ia=1,local_na
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
- DO iky=1,local_nky
- ky = kyarray(iky)
- DO ikx=1,total_nkx
- kx = kxarray(iky,ikx)
- DO iz=1,local_nz+ngz
- z = zarray(iz,ieven)
- IF (kx .EQ. 0) THEN
- IF(ky .EQ. 0) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
- ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp * deltaky/(ABS(ky)+deltaky)
- ENDIF
- ELSE
- IF(ky .GT. 0) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = (deltakx/(ABS(kx)+deltakx))*(deltaky/(ABS(ky)+deltaky))
- ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp*amp*(deltakx/(ABS(kx)+deltakx))
- ENDIF
- ENDIF
- ! z-dep and noise
- moments(ia,ip,ij,iky,ikx,iz,:) = moments(ia,ip,ij,iky,ikx,iz,:) * &
- (Jacobian(iz,ieven)*iInt_Jacobian)**2
- ! divide by kernel_0 to adjust to particle density (n = kernel_0 N00)
- ! moments(ia,ip,ij,iky,ikx,iz,:) = moments(ia,ip,ij,iky,ikx,iz,:)/kernel(ia,ij,iky,ikx,iz,0)
- END DO
- END DO
- END DO
- IF ( contains_ky0 ) THEN
- DO ikx=2,total_nkx/2 !symmetry at kx = 0 for all z
- moments(ia,ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:, :)
- END DO
- ENDIF
- ELSE
- moments(ia,ip,ij,:,:,:,:) = 0._xp
- ENDIF
- END DO
- END DO
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
-END SUBROUTINE init_ppj
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Initialize Ricci density
-!******************************************************************************!
-SUBROUTINE init_ricci
- USE basic, ONLY: maindir
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- local_nkx_offset, local_nky_offset, kxarray, kyarray, &
- ngp, ngj, ngz, iky0, ikx0, parray, jarray,&
- deltakx, deltaky, contains_ky0, contains_kx0,&
- AA_x, AA_y
- USE fields, ONLY: moments
- USE prec_const, ONLY: xp, imagu
- USE initial_par,ONLY: iseed, init_noiselvl, init_background
- USE model, ONLY: LINEARITY
- IMPLICIT NONE
- REAL(xp), DIMENSION(186,94) :: ricci_mat_real, ricci_mat_imag
- REAL(xp) :: scaling
- INTEGER :: ia,ip,ij,ikx,iky,iz, LPFx, LPFy
- CHARACTER(256) :: filename
-
- ! open data file
- ricci_mat_real = 0; ricci_mat_imag = 0
- filename = TRIM(maindir) // '/Gallery/fourier_ricci_real.txt'
- open(unit = 1 , file = filename)
- read(1,*) ricci_mat_real
- close(1)
- filename = TRIM(maindir) // '/Gallery/fourier_ricci_imag.txt'
- open(unit = 2 , file = filename)
- read(2,*) ricci_mat_imag
- close(2)
- scaling = 0.000001_xp
- moments = 0._xp
- !**** initialization *******************************************
- DO ia=1,local_na
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO ip=1+ngp/2,local_np+ngp/2
- DO ij=1+ngj/2,local_nj+ngj/2
- DO iz=1+ngz/2,local_nz+ngz/2
- IF((ikx+local_nkx_offset .LE. 186) .AND. (iky+local_nky_offset .LE. 94)) THEN
- ! IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
- moments(ia,ip,ij,iky,ikx,iz,:) = scaling*(ricci_mat_real(ikx+local_nkx_offset,iky+local_nky_offset)&
- - imagu*ricci_mat_imag(ikx+local_nkx_offset,iky+local_nky_offset))
- ! ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
- ! ENDIF
- ELSE
- moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
- ENDIF
- END DO
- END DO
- END DO
- END DO
- END DO
- ! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (LINEARITY .NE. 'linear') THEN
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- DO iz=1,local_nz+ngz
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDIF
- ENDDO
- !! Play with some filtering
- LPFx = 0
- LPFy = 0
- DO ikx=1,total_nkx
- DO iky=1,local_nky
- IF ( (abs(kxarray(ikx,iky)) .LT. LPFx*deltakx ) .AND. (abs(kyarray(iky)) .LT. LPFy*deltaky ) )&
- moments(:,:,:,iky,ikx,:,:) = 0._xp
- ENDDO
- ENDDO
-END SUBROUTINE init_ricci
-!******************************************************************************!
\ No newline at end of file
diff --git a/src/initial_mod.F90 b/src/initial_mod.F90
new file mode 100644
index 00000000..4491fea1
--- /dev/null
+++ b/src/initial_mod.F90
@@ -0,0 +1,651 @@
+MODULE initial
+ ! Module for initial parameters
+
+ USE prec_const
+ IMPLICIT NONE
+ PRIVATE
+
+ ! Initialization option (phi/mom00/allmom/blob)
+ CHARACTER(len=32), PUBLIC, PROTECTED :: INIT_OPT = 'phi'
+ ! Act on modes artificially (keep/wipe, zonal, non zonal, entropy mode etc.)
+ ! REMOVED FEATURE (only here for retrocompatibility)
+ CHARACTER(len=32), PUBLIC, PROTECTED :: ACT_ON_MODES = 'nothing'
+ ! Initial amplitude (for specific modes or image init.)
+ REAL(xp), PUBLIC, PROTECTED :: init_amp = 1._xp
+ ! Initial background level (in addition to noise init)
+ REAL(xp), PUBLIC, PROTECTED :: init_background = 0._xp
+ ! Initial noise amplitude (for noise init)
+ REAL(xp), PUBLIC, PROTECTED :: init_noiselvl = 1E-6_xp
+ ! Initialization for random number generator
+ INTEGER, PUBLIC, PROTECTED :: iseed=42
+
+ PUBLIC :: initial_outputinputs, initial_readinputs, initialize
+
+CONTAINS
+
+
+ SUBROUTINE initial_readinputs
+ ! Read the input parameters
+
+ USE basic, ONLY : lu_in
+ USE prec_const
+ IMPLICIT NONE
+
+ NAMELIST /INITIAL/ INIT_OPT,ACT_ON_MODES,&
+ init_amp,init_background,init_noiselvl,iseed
+
+ READ(lu_in,initial)
+
+ END SUBROUTINE initial_readinputs
+
+ !******************************************************************************!
+ !!!!!! initialize the moments and load/build coeff tables
+ !******************************************************************************!
+ SUBROUTINE initialize
+
+ USE basic, ONLY: speak
+ USE time_integration, ONLY: set_updatetlevel
+ USE collision, ONLY: init_collision
+ USE closure, ONLY: apply_closure_model
+ USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
+ USE restarts, ONLY: load_moments, job2load
+ USE processing, ONLY: compute_fluid_moments, compute_radial_heatflux, compute_radial_transport
+ USE model, ONLY: LINEARITY
+ USE nonlinear, ONLY: compute_Sapj, nonlinear_init
+ IMPLICIT NONE
+
+ CALL set_updatetlevel(1)
+
+ !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
+ ! through loading a previous state
+ IF ( job2load .GE. 0 ) THEN
+ CALL speak('Load moments')
+ CALL load_moments ! get N_0
+ CALL apply_closure_model
+ CALL update_ghosts_moments
+ CALL solve_EM_fields ! compute phi_0=phi(N_0)
+ CALL update_ghosts_EM
+ ! through initialization
+ ELSE
+ SELECT CASE (INIT_OPT)
+ ! set phi with noise and set moments to 0
+ CASE ('phi')
+ CALL speak('Init noisy phi')
+ CALL init_phi
+ CALL update_ghosts_EM
+ CASE ('phi_ppj')
+ CALL speak('Init noisy phi')
+ CALL init_phi_ppj
+ CALL update_ghosts_EM
+ ! set moments_00 (GC density) with noise and compute phi afterwards
+ CASE('mom00')
+ CALL speak('Init noisy gyrocenter density')
+ CALL init_gyrodens ! init only gyrocenter density
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ ! set moments_00 (GC density) with a single kx0,ky0 mode
+ ! kx0 is setup but by the noise value and ky0 by the background value
+ CASE('mom00_single_mode')
+ CALL speak('Init noisy gyrocenter density')
+ CALL init_single_mode ! init single mode in gyrocenter density
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ ! init all moments randomly (unadvised)
+ CASE('allmom')
+ CALL speak('Init noisy moments')
+ CALL init_moments ! init all moments
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ ! init a gaussian blob in gyrodens
+ CASE('blob')
+ CALL speak('--init a blob')
+ CALL initialize_blob
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ ! init moments 00 with a power law similarly to GENE
+ CASE('ppj')
+ CALL speak('ppj init ~ GENE')
+ call init_ppj
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ CASE('ricci')
+ CALL speak('Init Ricci')
+ CALL init_ricci ! init only gyrocenter density
+ CALL update_ghosts_moments
+ CALL solve_EM_fields
+ CALL update_ghosts_EM
+ END SELECT
+ ENDIF
+ ! closure of j>J, p>P and j<0, p<0 moments
+ CALL speak('Apply closure')
+ CALL apply_closure_model
+ ! ghosts for p parallelization
+ CALL speak('Ghosts communication')
+ CALL update_ghosts_moments
+ CALL update_ghosts_EM
+ !! End of phi and moments initialization
+
+ ! Init collision operator
+ CALL init_collision
+
+ !! Preparing auxiliary arrays at initial state
+ ! particle density, fluid velocity and temperature (used in diagnose)
+ CALL speak('Computing fluid moments and transport')
+ CALL compute_fluid_moments
+ CALL compute_radial_transport
+ CALL compute_radial_heatflux
+
+ ! init auxval for nonlinear
+ CALL nonlinear_init
+ ! compute nonlinear for t=0 diagnostic
+ CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
+
+
+ END SUBROUTINE initialize
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize all the moments randomly
+ !******************************************************************************!
+ SUBROUTINE init_moments
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ ngp, ngj, ngz, iky0, contains_ky0, AA_x, AA_y
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+
+ REAL(xp) :: noise
+ INTEGER, DIMENSION(12) :: iseedarr
+ INTEGER :: ia,ip,ij,ikx,iky,iz, ipi,iji,izi
+
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
+
+ !**** Broad noise initialization *******************************************
+ DO ia=1,local_na
+ DO ip=1,local_np
+ ipi = ip+ngp/2
+ DO ij=1,local_nj
+ iji = ij+ngj/2
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz
+ izi = iz+ngz/2
+ CALL RANDOM_NUMBER(noise)
+ moments(ia,ipi,iji,iky,ikx,izi,:) = (init_background + init_noiselvl*(noise-0.5_xp))
+ END DO
+ END DO
+ END DO
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
+ moments(ia,ipi,iji,iky0,ikx,:,:) = moments(ia,ipi,iji,iky0,total_nkx+2-ikx,:,:)
+ END DO
+ ENDIF
+ END DO
+ END DO
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz
+ izi = iz+ngz/2
+ DO ip=1,local_np
+ ipi = ip+ngp/2
+ 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
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ END SUBROUTINE init_moments
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize the gyrocenter density randomly
+ !******************************************************************************!
+ SUBROUTINE init_gyrodens
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ ngp, ngj, ngz, iky0, parray, jarray, contains_ky0, AA_x, AA_y
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+
+ REAL(xp) :: noise
+ INTEGER :: ia,ip,ij,ikx,iky,iz
+ INTEGER, DIMENSION(12) :: iseedarr
+
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
+ moments = 0._xp
+ !**** Broad noise initialization *******************************************
+ DO ia=1,local_na
+ DO ip=1+ngp/2,local_np+ngp/2
+ DO ij=1+ngj/2,local_nj+ngj/2
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1+ngz/2,local_nz+ngz/2
+ CALL RANDOM_NUMBER(noise)
+ IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
+ moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_xp))
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ENDIF
+ END DO
+ END DO
+ END DO
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
+ moments(ia, ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:,:)
+ END DO
+ ENDIF
+ END DO
+ END DO
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ END SUBROUTINE init_gyrodens
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize a single mode in the gyrocenter density
+ !******************************************************************************!
+ SUBROUTINE init_single_mode
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ ngp, ngj, ngz, iky0, kxarray, kyarray, parray, jarray, &
+ contains_ky0, deltakx, deltaky
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+
+ REAL(xp) :: amplitude
+ INTEGER :: ia,ip,ij,ikx,iky,iz, ikxi, ikyi
+ INTEGER, DIMENSION(12) :: iseedarr
+
+ moments = 0._xp
+ IF (my_id .EQ. 0) THEN
+ moments(:,:,:,2,2,:,:) = init_amp
+ moments(:,:,:,3,3,:,:) = init_amp
+ moments(:,:,:,4,4,:,:) = init_amp
+ moments(:,:,:,5,5,:,:) = init_amp
+ ENDIF
+ END SUBROUTINE init_single_mode
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize a noisy ES potential and cancel the moments
+ !******************************************************************************!
+ SUBROUTINE init_phi
+ USE grid, ONLY: total_nkx, local_nky, local_nz,&
+ ngz, iky0, ikx0, contains_ky0, AA_x, AA_y
+ USE fields, ONLY: phi, moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE parallel, ONLY: my_id
+ IMPLICIT NONE
+ REAL(xp) :: noise
+ INTEGER, DIMENSION(12) :: iseedarr
+ INTEGER :: iky,ikx,iz
+ ! Seed random number generator
+ iseedarr(:)=iseed
+ CALL RANDOM_SEED(PUT=iseedarr+my_id)
+ !**** noise initialization *******************************************
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ CALL RANDOM_NUMBER(noise)
+ phi(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_xp))!*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ END DO
+ END DO
+ !symmetry at ky = 0 to keep real inverse transform
+ IF ( contains_ky0 ) THEN
+ DO iz=1,local_nz+ngz
+ DO ikx=2,total_nkx/2
+ phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
+ ENDDO
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
+ END DO
+ ENDIF
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ phi(iky,ikx,iz) = phi(iky,ikx,iz)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ !**** ensure no previous moments initialization
+ moments = 0._xp
+ END SUBROUTINE init_phi
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize a ppj ES potential and cancel the moments
+ !******************************************************************************!
+ SUBROUTINE init_phi_ppj
+ USE grid, ONLY: total_nkx, local_nky, local_nz, AA_x, AA_y,&
+ ngz, iky0, ikx0, contains_ky0, ieven, kxarray, kyarray, zarray, deltakx
+ USE fields, ONLY: phi, moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE geometry, ONLY: Jacobian, iInt_Jacobian
+ IMPLICIT NONE
+ REAL(xp) :: kx, ky, z, amp
+ INTEGER :: ikx, iky, iz
+ amp = 1.0_xp
+ !**** ppj initialization *******************************************
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
+ DO iz=1,local_nz+ngz
+ z = zarray(iz,ieven)
+ IF (ky .NE. 0) THEN
+ phi(iky,ikx,iz) = 0._xp
+ ELSE
+ phi(iky,ikx,iz) = 0.5_xp*amp*(deltakx/(ABS(kx)+deltakx))
+ ENDIF
+ ! z-dep and noise
+ phi(iky,ikx,iz) = phi(iky,ikx,iz) * &
+ (Jacobian(iz,ieven)*iInt_Jacobian)**2
+ END DO
+ END DO
+ END DO
+ !symmetry at ky = 0 to keep real inverse transform
+ IF ( contains_ky0 ) THEN
+ DO iz=1,local_nz+ngz
+ DO ikx=2,total_nkx/2
+ phi(iky0,ikx,iz) = phi(iky0,total_nkx+2-ikx,iz)
+ ENDDO
+ phi(iky0,ikx0,iz) = REAL(phi(iky0,ikx0,iz),xp) !origin must be real
+ END DO
+ ENDIF
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ phi(iky,ikx,iz) = phi(iky,ikx,iz)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ !**** ensure no previous moments initialization
+ moments = 0._xp
+ END SUBROUTINE init_phi_ppj
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !******************************************************************************!
+ !!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
+ !******************************************************************************!
+ SUBROUTINE initialize_blob
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz, total_nz,&
+ AA_x, AA_y, parray, jarray,&
+ ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp
+ USE model, ONLY: LINEARITY
+ USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
+ IMPLICIT NONE
+ REAL(xp) ::kx, ky, z, sigma_x, sigma_y, gain
+ INTEGER :: ia,iky,ikx,iz,ip,ij, p, j
+ sigma_y = 0.5
+ sigma_x = sigma_y
+ gain = 1.0
+ ! One can increase the gain if we run 3D sim
+ IF(total_nz .GT. 1) THEN
+ sigma_y = 1.0
+ sigma_x = sigma_y
+ gain = 10.0
+ ENDIF
+
+ DO ia=1,local_na
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO iz=1+ngz/2,local_nz+ngz/2
+ z = zarray(iz,ieven)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx) + z*shear*ky
+ DO ip=1+ngp/2,local_np+ngp/2
+ p = parray(ip)
+ DO ij=1+ngj/2,local_nj+ngj/2
+ j = jarray(ij)
+ IF( (iky .NE. iky0) .AND. (p .EQ. 0) .AND. (j .EQ. 0)) THEN
+ moments(ia,ip,ij,iky,ikx,iz, :) = moments(ia,ip,ij,iky,ikx,iz, :) &
+ + gain * exp(-((kx/sigma_x)**2+(ky/sigma_y)**2)) &
+ * AA_x(ikx)*AA_y(iky)* &
+ (Jacobian(iz,ieven)*iInt_Jacobian)**2!&
+ ! * exp(sigmai2_taui_o2*(kx**2+ky**2))
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ END SUBROUTINE initialize_blob
+ !******************************************************************************!
+ !******************************************************************************!
+ !!!!!!! Initialize the gyrocenter in a ppj gene manner (power law)
+ !******************************************************************************!
+ SUBROUTINE init_ppj
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ AA_x, AA_y, deltakx, deltaky,contains_ky0,&
+ ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp, pi
+ USE model, ONLY: LINEARITY
+ USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
+ IMPLICIT NONE
+ REAL(xp) :: kx, ky, sigma_z, z
+ INTEGER :: ia,iky,ikx,iz,ip,ij
+ sigma_z = pi/4._xp
+ !**** Broad noise initialization *******************************************
+ DO ia=1,local_na
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
+ DO iz=1,local_nz+ngz
+ z = zarray(iz,ieven)
+ IF (kx .EQ. 0) THEN
+ IF(ky .EQ. 0) THEN
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp * deltaky/(ABS(ky)+deltaky)
+ ENDIF
+ ELSE
+ IF(ky .GT. 0) THEN
+ moments(ia,ip,ij,iky,ikx,iz,:) = (deltakx/(ABS(kx)+deltakx))*(deltaky/(ABS(ky)+deltaky))
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0.5_xp*(deltakx/(ABS(kx)+deltakx))
+ ENDIF
+ ENDIF
+ ! z-dep and noise
+ moments(ia,ip,ij,iky,ikx,iz,:) = moments(ia,ip,ij,iky,ikx,iz,:) * &
+ (Jacobian(iz,ieven)*iInt_Jacobian)**2
+ ! divide by kernel_0 to adjust to particle density (n = kernel_0 N00)
+ ! moments(ia,ip,ij,iky,ikx,iz,:) = moments(ia,ip,ij,iky,ikx,iz,:)/kernel(ia,ij,iky,ikx,iz,0)
+ END DO
+ END DO
+ END DO
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,total_nkx/2 !symmetry at kx = 0 for all z
+ moments(ia,ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:, :)
+ END DO
+ ENDIF
+ ELSE
+ moments(ia,ip,ij,:,:,:,:) = 0._xp
+ ENDIF
+ END DO
+ END DO
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ END SUBROUTINE init_ppj
+ !******************************************************************************!
+
+ !******************************************************************************!
+ !!!!!!! Initialize Ricci density
+ ! We take the picture of Paolo Ricci and use it as input for the initialization
+ ! of the gyrodensity. The picture is expected to be 186x96 fourier modes stored
+ ! in the gyacomo/Gallery directory, in two 2D matrices (real.txt and imag.txt)
+ ! with only spaces as separators.
+ ! One can scale the amplitude of the modes with the init_background parameter
+ ! (1e-5 approx is needed to calm down the system if it is in nonlinear evolution)
+ !******************************************************************************!
+ SUBROUTINE init_ricci
+ USE basic, ONLY: maindir
+ USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+ local_nkx_offset, local_nky_offset, kxarray, kyarray, &
+ ngp, ngj, ngz, iky0, ikx0, parray, jarray,&
+ deltakx, deltaky, contains_ky0, contains_kx0,&
+ AA_x, AA_y
+ USE fields, ONLY: moments
+ USE prec_const, ONLY: xp, imagu
+ USE model, ONLY: LINEARITY
+ IMPLICIT NONE
+ REAL(xp), DIMENSION(186,94) :: ricci_mat_real, ricci_mat_imag
+ REAL(xp) :: scaling
+ INTEGER :: ia,ip,ij,ikx,iky,iz, LPFx, LPFy
+ CHARACTER(256) :: filename
+
+ ! open data file
+ ricci_mat_real = 0; ricci_mat_imag = 0
+ filename = TRIM(maindir) // '/Gallery/fourier_ricci_real.txt'
+ OPEN(unit = 1 , file = filename)
+ READ(1,*) ricci_mat_real
+ CLOSE(1)
+ filename = TRIM(maindir) // '/Gallery/fourier_ricci_imag.txt'
+ OPEN(unit = 2 , file = filename)
+ READ(2,*) ricci_mat_imag
+ CLOSE(2)
+ scaling = init_amp*1e-5
+ moments = 0._xp
+ !**** initialization *******************************************
+ DO ia=1,local_na
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO ip=1+ngp/2,local_np+ngp/2
+ DO ij=1+ngj/2,local_nj+ngj/2
+ DO iz=1+ngz/2,local_nz+ngz/2
+ IF((ikx+local_nkx_offset .LE. 186) .AND. (iky+local_nky_offset .LE. 94)) THEN
+ IF (.TRUE.) THEN
+ ! IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
+ moments(ia,ip,ij,iky,ikx,iz,:) = scaling*(ricci_mat_real(ikx+local_nkx_offset,iky+local_nky_offset)&
+ - imagu*ricci_mat_imag(ikx+local_nkx_offset,iky+local_nky_offset))
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ENDIF
+ ELSE
+ moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+ ENDIF
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+ ! Putting to zero modes that are not in the 2/3 Orszag rule
+ IF (LINEARITY .NE. 'linear') THEN
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ DO iz=1,local_nz+ngz
+ DO ip=1,local_np+ngp
+ DO ij=1,local_nj+ngj
+ moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ !! Play with some filtering
+ LPFx = 0
+ LPFy = 0
+ DO ikx=1,total_nkx
+ DO iky=1,local_nky
+ IF ( (abs(kxarray(ikx,iky)) .LT. LPFx*deltakx ) .AND. (abs(kyarray(iky)) .LT. LPFy*deltaky ) )&
+ moments(:,:,:,iky,ikx,:,:) = 0._xp
+ ENDDO
+ ENDDO
+ END SUBROUTINE init_ricci
+ !******************************************************************************!
+
+
+ SUBROUTINE initial_outputinputs(fid)
+ ! Write the input parameters to the results_xx.h5 file
+ USE futils, ONLY: attach, creatd
+ IMPLICIT NONE
+ INTEGER, INTENT(in) :: fid
+ CHARACTER(len=256) :: str
+ WRITE(str,'(a)') '/data/input/intial'
+ CALL creatd(fid, 0,(/0/),TRIM(str),'Initial Input')
+ CALL attach(fid, TRIM(str), "INIT_OPT", INIT_OPT)
+ CALL attach(fid, TRIM(str), "init_background", init_background)
+ CALL attach(fid, TRIM(str), "init_noiselvl", init_noiselvl)
+ CALL attach(fid, TRIM(str), "iseed", iseed)
+ END SUBROUTINE initial_outputinputs
+
+END MODULE initial
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
deleted file mode 100644
index 4fd93b0b..00000000
--- a/src/initial_par_mod.F90
+++ /dev/null
@@ -1,60 +0,0 @@
-MODULE initial_par
- ! Module for initial parameters
-
- USE prec_const
- IMPLICIT NONE
- PRIVATE
-
- ! Initialization option (phi/mom00/allmom/blob)
- CHARACTER(len=32), PUBLIC, PROTECTED :: INIT_OPT = 'phi'
- ! Initialization through a zonal flow phi
- INTEGER, PUBLIC, PROTECTED :: INIT_ZF = 0
- REAL(xp), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_xp
- ! Act on modes artificially (keep/wipe, zonal, non zonal, entropy mode etc.)
- CHARACTER(len=32), PUBLIC, PROTECTED :: ACT_ON_MODES = 'nothing'
- ! Initial background level
- REAL(xp), PUBLIC, PROTECTED :: init_background=0._xp
- ! Initial noise amplitude
- REAL(xp), PUBLIC, PROTECTED :: init_noiselvl=1E-6_xp
- ! Initialization for random number generator
- INTEGER, PUBLIC, PROTECTED :: iseed=42
-
- PUBLIC :: initial_outputinputs, initial_readinputs
-
-CONTAINS
-
-
- SUBROUTINE initial_readinputs
- ! Read the input parameters
-
- USE basic, ONLY : lu_in
- USE prec_const
- IMPLICIT NONE
-
- NAMELIST /INITIAL_CON/ INIT_OPT
- NAMELIST /INITIAL_CON/ ACT_ON_MODES
- NAMELIST /INITIAL_CON/ init_background
- NAMELIST /INITIAL_CON/ init_noiselvl
- NAMELIST /INITIAL_CON/ iseed
-
- READ(lu_in,initial_con)
- !WRITE(*,initial_con)
-
- END SUBROUTINE initial_readinputs
-
-
- SUBROUTINE initial_outputinputs(fid)
- ! Write the input parameters to the results_xx.h5 file
- USE futils, ONLY: attach, creatd
- IMPLICIT NONE
- INTEGER, INTENT(in) :: fid
- CHARACTER(len=256) :: str
- WRITE(str,'(a)') '/data/input/intial'
- CALL creatd(fid, 0,(/0/),TRIM(str),'Initial Input')
- CALL attach(fid, TRIM(str), "INIT_OPT", INIT_OPT)
- CALL attach(fid, TRIM(str), "init_background", init_background)
- CALL attach(fid, TRIM(str), "init_noiselvl", init_noiselvl)
- CALL attach(fid, TRIM(str), "iseed", iseed)
- END SUBROUTINE initial_outputinputs
-
-END MODULE initial_par
diff --git a/src/main.F90 b/src/main.F90
index 83f9d66e..13cd8a9a 100644
--- a/src/main.F90
+++ b/src/main.F90
@@ -1,10 +1,13 @@
-!GYACOMO (Gyrokinetic Advanced Collision Moment solver, 2021)
-!Copyright (C) 2022 A.C.D. Hoffmann
+! GYACOMO (Gyrokinetic Advanced Collision Moment solver, 2021)
+! GNU GENERAL PUBLIC LICENSE
+! Version 3, 29 June 2007
+
+! Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+! Everyone is permitted to copy and distribute verbatim copies
+! of this license document, but changing it is not allowed.
PROGRAM main
! Main program
-
- use prec_const
IMPLICIT NONE
CALL control
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index f74659a8..f468a0cd 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -26,6 +26,8 @@ MODULE model
REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
REAL(xp), PUBLIC, PROTECTED :: ExBrate = 0._xp ! ExB background shearing rate (radially constant shear flow)
+ INTEGER, PUBLIC, PROTECTED :: ikxZF = 0 ! Background zonal mode wavenumber (acts in the nonlinear term)
+ REAL(xp), PUBLIC, PROTECTED :: ZFamp = 200._xp ! Amplitude of the background zonal mode
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift
@@ -46,12 +48,12 @@ CONTAINS
USE prec_const, ONLY: xp
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
+ NAMELIST /MODEL/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
Na, ADIAB_E, ADIAB_I, tau_i, &
mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, &
- nu, k_gB, k_cB, lambdaD, beta, ExBrate
+ nu, k_gB, k_cB, lambdaD, beta, ExBrate, ikxZF, ZFamp
- READ(lu_in,model_par)
+ READ(lu_in,model)
IF (ADIAB_E .AND. ADIAB_I) &
ERROR STOP '>> ERROR << cannot have both adiab e and adiab i models'
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index ccfa5f5d..be57b22f 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -1,6 +1,5 @@
MODULE nonlinear
USE array, ONLY : dnjs, Sapj, kernel
- USE initial_par, ONLY : ACT_ON_MODES
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, &
@@ -8,8 +7,8 @@ MODULE nonlinear
local_nj,ngj,jarray,jmax, local_nj_offset, dmax,&
kyarray, AA_y, local_nky_ptr, local_nky_ptr_offset,inv_Ny,&
local_nkx_ptr,kxarray, AA_x, inv_Nx,&
- local_nz,ngz,zarray,nzgrid
- USE model, ONLY : LINEARITY, EM
+ local_nz,ngz,zarray,nzgrid, deltakx, iky0, contains_kx0, contains_ky0
+ USE model, ONLY : LINEARITY, EM, ikxZF, ZFamp
USE closure, ONLY : evolve_mom, nmaxarray
USE prec_const, ONLY : xp
USE species, ONLY : sqrt_tau_o_sigma
@@ -63,6 +62,7 @@ END SUBROUTINE compute_Sapj
SUBROUTINE compute_nonlinear
IMPLICIT NONE
INTEGER :: iz,ij,ip,eo,ia,ikx,iky,izi,ipi,iji,ini,isi
+ INTEGER :: ikxExBp, ikxExBn ! Negative and positive ExB flow indices
z:DO iz = 1,local_nz
izi = iz + ngz/2
j:DO ij = 1,local_nj ! Loop over Laguerre moments
@@ -82,7 +82,20 @@ SUBROUTINE compute_nonlinear
ini = in+ngj/2
!-----------!! ELECTROSTATIC CONTRIBUTION
! First convolution terms
- F_cmpx(:,:) = phi(:,:,izi) * kernel(ia,ini,:,:,izi,eo)
+ DO ikx = 1,local_nkx_ptr
+ DO iky = 1,local_nky_ptr
+ F_cmpx(iky,ikx) = phi(iky,ikx,izi) * kernel(ia,ini,iky,ikx,izi,eo)
+ ENDDO
+ ENDDO
+ ! Test to implement the ExB shearing as a additional zonal mode in the ES potential
+ IF(ikxZF .GT. 1) THEN
+ ikxExBp = ikxZF
+ ikxExBn = local_nkx_ptr - (ikxExBp-2)
+ IF(contains_kx0 .AND. contains_ky0) THEN
+ F_cmpx(iky0,ikxExBp) = F_cmpx(iky0,ikxExBp) + ZFamp * kernel(ia,ini,iky0,ikxExBp,izi,eo)
+ F_cmpx(iky0,ikxExBn) = F_cmpx(iky0,ikxExBn) + ZFamp * kernel(ia,ini,iky0,ikxExBn,izi,eo)
+ ENDIF
+ ENDIF
! Second convolution terms
G_cmpx = 0._xp ! initialization of the sum
smax = MIN( jarray(ini)+jarray(iji), jmax );
@@ -91,7 +104,7 @@ SUBROUTINE compute_nonlinear
G_cmpx(:,:) = G_cmpx(:,:) + &
dnjs(in,ij,is) * moments(ia,ipi,isi,:,:,izi,updatetlevel)
ENDDO s1
- ! this function add its result to bracket_sum_r
+ ! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r)
!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
IF(EM) THEN
@@ -105,7 +118,7 @@ SUBROUTINE compute_nonlinear
dnjs(in,ij,is) * (sqrt_pp1*moments(ia,ipi+1,isi,:,:,izi,updatetlevel)&
+sqrt_p *moments(ia,ipi-1,isi,:,:,izi,updatetlevel))
ENDDO s2
- ! this function add its result to bracket_sum_r
+ ! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r)
ENDIF
ENDDO n
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 02f8b5d5..d0550189 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -2,11 +2,11 @@ SUBROUTINE readinputs
! Additional data specific for a new run
USE grid, ONLY: grid_readinputs
- USE diagnostics_par, ONLY: diag_par_readinputs
+ USE diagnostics, ONLY: diag_readinputs
USE collision, ONLY: collision_readinputs
USE model, ONLY: model_readinputs
USE species, ONLY: species_readinputs
- USE initial_par, ONLY: initial_readinputs
+ USE initial, ONLY: initial_readinputs
USE time_integration, ONLY: time_integration_readinputs
USE geometry, ONLY: geometry_readinputs
USE closure, ONLY: closure_readinputs
@@ -25,7 +25,7 @@ SUBROUTINE readinputs
CALL geometry_readinputs
! Load diagnostic options from input file
- CALL diag_par_readinputs
+ CALL diag_readinputs
! Load model parameters from input file
CALL model_readinputs
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 491fd53b..7621c74d 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -7,12 +7,15 @@ USE grid, ONLY: local_Na,local_Na_offset,local_np,local_np_offset,&
local_nz,local_nz_offset,ngp,ngj,ngz, deltap,&
total_Na,total_Np,total_Nj,total_Nky,total_Nkx,total_Nz,ieven,&
kyarray_full,kxarray_full,zarray, zarray_full, local_zmin, local_zmax ! for z interp
-USE fields
-USE diagnostics_par
-USE time_integration
-USE prec_const, ONLY : xp,dp,sp
+USE fields, ONLY: moments, phi, psi
+!USE time_integration
+USE prec_const, ONLY : xp,dp,sp,PI
+USE MPI, ONLY: MPI_COMM_WORLD
IMPLICIT NONE
+CHARACTER(len=256), PUBLIC :: rstfile ! restart result file
+INTEGER, PUBLIC :: fidrst ! FID for restart file
+
PUBLIC :: load_moments!, write_restart
CONTAINS
@@ -39,7 +42,7 @@ CONTAINS
COMPLEX(xp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_cp
CALL cpu_time(timer_tot_1)
! Checkpoint filename
- WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',job2load,'.h5'
+ WRITE(rstfile,'(a,i2.2,a3)') 'outputs_',job2load,'.h5'
CALL speak("Resume from "//rstfile)
! Open file
IF(xp .EQ. dp) THEN
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 5b8e5654..0ed20e6e 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -44,10 +44,10 @@ CONTAINS
USE prec_const
USE basic, ONLY : lu_in
IMPLICIT NONE
- NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme
- namelist /TIME_INTEGRATION_PAR/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
+ NAMELIST /TIME_INTEGRATION/ numerical_scheme
+ namelist /TIME_INTEGRATION/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
- READ(lu_in,time_integration_par)
+ READ(lu_in,time_integration)
CALL set_numerical_scheme
END SUBROUTINE time_integration_readinputs
diff --git a/testcases/DLRA_zpinch/base_case/fort_00.90 b/testcases/DLRA_zpinch/base_case/fort_00.90
index 07e2cd32..1dfdd69d 100644
--- a/testcases/DLRA_zpinch/base_case/fort_00.90
+++ b/testcases/DLRA_zpinch/base_case/fort_00.90
@@ -30,7 +30,7 @@
parallel_bc = 'shearless'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 0.1
dtsave_1d = -1
dtsave_2d = 0.1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 1.0
@@ -61,7 +61,7 @@
ADIAB_E = .f.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
!hierarchy_closure='max_degree'
dmax = 2
@@ -86,20 +86,20 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
&CLA
diff --git a/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
index 1fa68c8b..dc0af1ff 100644
--- a/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
+++ b/testcases/DLRA_zpinch/nsv_filter_2/fort_00.90
@@ -30,7 +30,7 @@
parallel_bc = 'shearless'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 0.1
dtsave_1d = -1
dtsave_2d = 0.1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 1.0
@@ -61,7 +61,7 @@
ADIAB_E = .f.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
!hierarchy_closure='max_degree'
dmax = 2
@@ -86,20 +86,20 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
&CLA
diff --git a/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90 b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
index 81ae89e5..3bc58a92 100644
--- a/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
+++ b/testcases/DLRA_zpinch/nsv_filter_6/fort_00.90
@@ -30,7 +30,7 @@
parallel_bc = 'shearless'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 0.1
dtsave_1d = -1
dtsave_2d = 0.1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 1.0
@@ -61,7 +61,7 @@
ADIAB_E = .f.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
!hierarchy_closure='max_degree'
dmax = 2
@@ -86,20 +86,20 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
&CLA
diff --git a/testcases/ITG_zpinch/fort_00.90 b/testcases/ITG_zpinch/fort_00.90
index 536bb949..87253b73 100644
--- a/testcases/ITG_zpinch/fort_00.90
+++ b/testcases/ITG_zpinch/fort_00.90
@@ -30,7 +30,7 @@
parallel_bc = 'shearless'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 0.0
@@ -61,7 +61,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
!hierarchy_closure='truncation'
hierarchy_closure='max_degree'
dmax = 2
@@ -86,19 +86,19 @@
k_N_ = 1.6
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/cyclone_example/fort_00.90 b/testcases/cyclone_example/fort_00.90
index 07b169f3..309ceed7 100644
--- a/testcases/cyclone_example/fort_00.90
+++ b/testcases/cyclone_example/fort_00.90
@@ -31,7 +31,7 @@
parallel_bc = 'dirichlet'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -47,7 +47,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 0.0
@@ -62,7 +62,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='truncation'
@@ -78,20 +78,20 @@
k_T_ = 6.96
/
-&COLLISION_PAR
+&COLLISION
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
+&INITIAL
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
!numerical_scheme = 'RK4'
numerical_scheme = 'SSP_RK3'
/
diff --git a/testcases/cyclone_example/fort_01.90 b/testcases/cyclone_example/fort_01.90
index 07b169f3..309ceed7 100644
--- a/testcases/cyclone_example/fort_01.90
+++ b/testcases/cyclone_example/fort_01.90
@@ -31,7 +31,7 @@
parallel_bc = 'dirichlet'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -47,7 +47,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 0.0
@@ -62,7 +62,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='truncation'
@@ -78,20 +78,20 @@
k_T_ = 6.96
/
-&COLLISION_PAR
+&COLLISION
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
+&INITIAL
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
!numerical_scheme = 'RK4'
numerical_scheme = 'SSP_RK3'
/
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index bc7dce3f..7b62621e 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -30,7 +30,7 @@
parallel_bc = 'dirichlet'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 0.5
dtsave_1d = -1
dtsave_2d = -1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
! Collisionality
CLOS = 0
NL_CLOS = -1
@@ -83,19 +83,19 @@
k_T_ = 2.0!6.96
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index d330d8e0..5e3925a0 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -23,7 +23,7 @@
eps = 0.18
parallel_bc = 'dirichlet'
/
-&OUTPUT_PAR
+&DIAGNOSTICS
nsave_0d = 50
nsave_1d = -1
nsave_2d = -1
@@ -39,7 +39,7 @@
write_temp = .t.
job2load = -1
/
-&MODEL_PAR
+&MODEL
! Collisionality
CLOS = 0
NL_CLOS = -1
@@ -65,19 +65,19 @@
K_Ti = 4!2.22
beta = 0.1
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
gyrokin_CO = .t.
interspecies = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/smallest_problem/fort_01.90 b/testcases/smallest_problem/fort_01.90
index ecaa8711..5670c21b 100644
--- a/testcases/smallest_problem/fort_01.90
+++ b/testcases/smallest_problem/fort_01.90
@@ -23,7 +23,7 @@
eps = 0.18
parallel_bc = 'dirichlet'
/
-&OUTPUT_PAR
+&DIAGNOSTICS
nsave_0d = 1
nsave_1d = -1
nsave_2d = -1
@@ -40,7 +40,7 @@
write_temp = .t.
job2load = -1
/
-&MODEL_PAR
+&MODEL
! Collisionality
CLOS = 0
NL_CLOS = -1
@@ -65,19 +65,19 @@
K_Ti = 2.22
beta = 0
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
!mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/smallest_problem/fort_11.90 b/testcases/smallest_problem/fort_11.90
index f6a8807c..c934f3b4 100644
--- a/testcases/smallest_problem/fort_11.90
+++ b/testcases/smallest_problem/fort_11.90
@@ -30,7 +30,7 @@
parallel_bc = 'dirichlet'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 0.5
dtsave_1d = -1
dtsave_2d = -1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
! Collisionality
CLOS = 0
NL_CLOS = -1
@@ -83,19 +83,19 @@
k_T_ = 2.0!6.96
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'allmom'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/zpinch_3D/fort_00.90 b/testcases/zpinch_3D/fort_00.90
index dc01408d..8e575fe4 100644
--- a/testcases/zpinch_3D/fort_00.90
+++ b/testcases/zpinch_3D/fort_00.90
@@ -35,7 +35,7 @@
shift_y= 0.0
Npol = 1.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -52,7 +52,7 @@
write_temp = .t.
!diag_mode = 'txtonly'
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 1.0
@@ -67,7 +67,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
@@ -91,19 +91,19 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/zpinch_3D/fort_01.90 b/testcases/zpinch_3D/fort_01.90
index 84aa2a6f..60d9f8a4 100644
--- a/testcases/zpinch_3D/fort_01.90
+++ b/testcases/zpinch_3D/fort_01.90
@@ -35,7 +35,7 @@
shift_y= 0.0
Npol = 1.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -52,7 +52,7 @@
write_temp = .t.
!diag_mode = 'txtonly'
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 1.0
@@ -67,7 +67,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
@@ -91,19 +91,19 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/zpinch_3D/fort_02.90 b/testcases/zpinch_3D/fort_02.90
index 3073e0e0..8fe3e08e 100644
--- a/testcases/zpinch_3D/fort_02.90
+++ b/testcases/zpinch_3D/fort_02.90
@@ -35,7 +35,7 @@
shift_y= 0.0
Npol = 1.8
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -52,7 +52,7 @@
write_temp = .t.
!diag_mode = 'txtonly'
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 1.0
@@ -67,7 +67,7 @@
ADIAB_E = .t.
tau_e = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax = -1
nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
@@ -91,19 +91,19 @@
k_N_ = 2.0
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/testcases/zpinch_example/fort_00.90 b/testcases/zpinch_example/fort_00.90
index 3a7fcb7c..2278803d 100644
--- a/testcases/zpinch_example/fort_00.90
+++ b/testcases/zpinch_example/fort_00.90
@@ -30,7 +30,7 @@
parallel_bc = 'dirichlet'
shift_y= 0.0
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -46,7 +46,7 @@
write_fvel = .t.
write_temp = .t.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 1.0
@@ -61,7 +61,7 @@
ADIAB_E = .f.
tau_i = 1.0
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
!hierarchy_closure='max_degree'
dmax = 2
@@ -86,19 +86,19 @@
k_N_ = 1.6
k_T_ = 0.4
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .t.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'blob' !(phi,blob)
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
diff --git a/wk/benchmark_and_scan_scripts/old/fort_00.90 b/wk/benchmark_and_scan_scripts/old/fort_00.90
index 70572048..78e75f4f 100644
--- a/wk/benchmark_and_scan_scripts/old/fort_00.90
+++ b/wk/benchmark_and_scan_scripts/old/fort_00.90
@@ -28,7 +28,7 @@
shift_y = 0
Npol = 1
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -43,7 +43,7 @@
write_dens = .true.
write_temp = .true.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'linear'
RM_LD_T_EQ= .false.
Na = 2
@@ -61,7 +61,7 @@ RM_LD_T_EQ= .false.
beta = 0.03
ADIAB_E = .false.
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax =-1
nonlinear_closure='truncation'
@@ -82,19 +82,19 @@ RM_LD_T_EQ= .false.
K_N_ = 3
K_T_ = 4.5
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG'
GK_CO = .false.
INTERSPECIES = .true.
mat_file = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/oiCa/null'
collision_kcut = 1
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'phi'
init_background = 1e-05
init_noiselvl = 0
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
\ No newline at end of file
diff --git a/wk/paper_2_scripts_and_results/fort_00.90 b/wk/paper_2_scripts_and_results/fort_00.90
index fe18d5ae..b08af283 100644
--- a/wk/paper_2_scripts_and_results/fort_00.90
+++ b/wk/paper_2_scripts_and_results/fort_00.90
@@ -32,7 +32,7 @@
Npol = 1
PB_PHASE= .false.
/
-&OUTPUT_PAR
+&DIAGNOSTICS
dtsave_0d = 1
dtsave_1d = -1
dtsave_2d = -1
@@ -47,7 +47,7 @@
write_dens = .false.
write_temp = .true.
/
-&MODEL_PAR
+&MODEL
LINEARITY = 'linear'
RM_LD_T_EQ= .false.
Na = 2
@@ -68,7 +68,7 @@ RM_LD_T_EQ= .false.
tau_i = 1
MHD_PD = .false.
/
-&CLOSURE_PAR
+&CLOSURE
hierarchy_closure='truncation'
dmax =-1
nonlinear_closure='truncation'
@@ -90,19 +90,19 @@ RM_LD_T_EQ= .false.
K_N_ = 2.22
K_T_ = 6.96
/
-&COLLISION_PAR
+&COLLISION
collision_model = 'DG'
GK_CO = .false.
INTERSPECIES = .true.
mat_file = '/home/ahoffman/gyacomo/wk/paper_2_scripts_and_resuliCa/null'
collision_kcut = 1.75
/
-&INITIAL_CON
+&INITIAL
INIT_OPT = 'phi'
init_background = 0
init_noiselvl = 1e-05
iseed = 42
/
-&TIME_INTEGRATION_PAR
+&TIME_INTEGRATION
numerical_scheme = 'RK4'
/
--
GitLab