diff --git a/.gitignore b/.gitignore
index 795de19eb478f0e9498da061cf7ccf9e446e682a..f2c5b1fce9c33f26420cf33fba5ef18fec5fd662 100644
--- a/.gitignore
+++ b/.gitignore
@@ -44,3 +44,5 @@ Gallery/
out*
!scripts/*
*.avi
+gyacomo
+gyacomo_dbg
diff --git a/Makefile b/Makefile
index fd9aa229c40184ff3b9b294aee34764211229b87..9a0ccb3b062b7f27b8518c43e9692ebb44a9fd2d 100644
--- a/Makefile
+++ b/Makefile
@@ -67,12 +67,13 @@ $(OBJDIR)/collision_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o \
$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.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)/main.o $(OBJDIR)/memory.o \
-$(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o \
-$(OBJDIR)/parallel_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o \
-$(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o \
-$(OBJDIR)/restarts_mod.o $(OBJDIR)/solve_EM_fields.o $(OBJDIR)/stepon.o \
-$(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+$(OBJDIR)/initial_par_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)/ppinit.o $(OBJDIR)/prec_const_mod.o \
+$(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/restarts_mod.o \
+$(OBJDIR)/solve_EM_fields.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
+$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(EXEC): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
@@ -83,158 +84,196 @@ $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(EDBG): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
- $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o \
- $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/advance_field.o : src/advance_field.F90 \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o \
+ $(OBJDIR)/fields_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
- $(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/array_mod.o : src/array_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@
- $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o \
- $(OBJDIR)/model_mod.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o\
- $(OBJDIR)/parallel_mod.o
+ $(OBJDIR)/auxval.o : src/auxval.F90 \
+ $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
+ $(OBJDIR)/geometry_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o \
+ $(OBJDIR)/parallel_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
- $(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/basic_mod.o : src/basic_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@
- $(OBJDIR)/calculus_mod.o : src/calculus_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/parallel_mod.o
+ $(OBJDIR)/calculus_mod.o : src/calculus_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/parallel_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/calculus_mod.F90 -o $@
- $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o \
+ $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
- $(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/closure_mod.o : src/closure_mod.F90 \
+ $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
- $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o \
+ $(OBJDIR)/collision_mod.o : src/collision_mod.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o \
$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -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)/ppinit.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)/ppinit.o \
$(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
- $(OBJDIR)/diagnose.o : src/diagnose.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)/diagnose.o : src/diagnose.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
+ $(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)/endrun.o : src/endrun.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/endrun.o : src/endrun.F90 \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@
- $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/fields_mod.o : src/fields_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
- $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
- $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o\
- $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/miller_mod.o \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o \
+ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
- $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o\
+ $(OBJDIR)/geometry_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
- $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o\
- $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/grid_mod.o : src/grid_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/model_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)/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)/solve_EM_fields.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/inital.o : src/inital.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)/solve_EM_fields.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o \
+ $(OBJDIR)/grid_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
+ $(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 $@
+ $(OBJDIR)/lag_interp_mod.o : src/lag_interp_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/lag_interp_mod.F90 -o $@
+
$(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
- $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/collision_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/memory.o : src/memory.F90 $ \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/collision_mod.o\
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
$(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
- $(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/miller_mod.o : src/miller_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/lag_interp_mod.o \
+ $(OBJDIR)/prec_const_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/miller_mod.F90 -o $@
+
+ $(OBJDIR)/model_mod.o : src/model_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
- $(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 $(OBJDIR)/array_mod.o \
- $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o \
- $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o \
+ $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@
- $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o\
+ $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o \
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o\
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -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)/numerics_mod.o : src/numerics_mod.F90 \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o \
+ $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics_mod.F90 -o $@
- $(OBJDIR)/parallel_mod.o : src/parallel_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/parallel_mod.o : src/parallel_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/parallel_mod.F90 -o $@
- $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/coeff_mod.o
+ $(OBJDIR)/ppexit.o : src/ppexit.F90 \
+ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@
- $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/ppinit.o : src/ppinit.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o\
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o \
$(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@
$(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@
- $(OBJDIR)/processing_mod.o : src/processing_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/processing_mod.o : src/processing_mod.F90 \
+ $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/basic_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)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/readinputs.o : src/readinputs.F90 \
+ $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_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)/restarts_mod.o : src/restarts_mod.F90 \
+ $(OBJDIR)/diagnostics_par_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)/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
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/solve_EM_fields.F90 -o $@
- $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/advance_field.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 $(OBJDIR)/time_integration_mod.o \
- $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+ $(OBJDIR)/stepon.o : src/stepon.F90 \
+ $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.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\
+ $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
- $(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/tesend.o : src/tesend.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/tesend.F90 -o $@
- $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o \
- $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
- $(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
- $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/utility_mod.o : src/utility_mod.F90 \
+ $(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
diff --git a/README.md b/README.md
index e88d5631e599d0bb4914ec0bd0e61de761ec19ee..55f0f197d413ddf893e341d0371408623a83525e 100644
--- a/README.md
+++ b/README.md
@@ -24,15 +24,16 @@ How to run it
4. You can run a typical CBC to test the compilation using the basic fort.90 parameter file,
just type ./bin/gyacomo
5. It is possible to run it in parallel (MPI) as mpirun -np N ./bin/gyacomo Np Ny Nz
- where N=Np x Ny x Nz is the number of processes and Np Ny Nz are the parallel dimensions in
- Hermite polynomials, binormal direction and parallel direction, respectively
+ where N=Np x Ny x Nz is the number of processes and Np Ny Nz are the parallel dimensions in Hermite polynomials, binormal direction and parallel direction, respectively
6. You can stop your simulation without breaking output file by creating a blank file call "mystop"
in the directory where the simulation is running. (the file will be removed once read)
7. You can obtain various plots and gifs using gyacomo/wk/header_3D_results.m once the simulation is done.
-// Comment : For some collision operators (Sugama and Full Coulomb) you have to run COSOlver from B.J.Frei first in order to generate the required matrices in gyacomo/iCa folder.
+// Comment : For some collision operators (Sugama and Full Coulomb) you have to run COSOlver from B.J.Frei first in order to generate the required matrices in gyacomo/iCa folder. //
# Changelog
4. GYACOMO
+ 4.1 Miller geometry is added and benchmarked for CBC adiabatic electrons
+
4.0 new naming and opening the code with GNU GPLv3 license
3. HeLaZ 3D
diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index d0beda9da004f60fa15035417bb9a6db523e1c74..bb7cea95cc0fcc5aa75ac597aaf4b338b8d63fa7 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -1,4 +1,4 @@
-function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT)
+function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, OPTIONS)
% Remove AA part
if DATA.Nx > 1
@@ -11,8 +11,8 @@ ikynz = logical(ikynz .* (DATA.ky > 0));
phi = DATA.PHI(ikynz,ikxnz,:,:);
t = DATA.Ts3D;
-[~,its] = min(abs(t-TRANGE(1)));
-[~,ite] = min(abs(t-TRANGE(end)));
+[~,its] = min(abs(t-OPTIONS.TRANGE(1)));
+[~,ite] = min(abs(t-OPTIONS.TRANGE(end)));
w_ky = zeros(sum(ikynz),ite-its);
ce = zeros(sum(ikynz),ite-its);
@@ -34,7 +34,7 @@ for it = its+1:ite
end
[kys, Is] = sort(DATA.ky(ikynz));
-linear_gr.trange = t(its:ite);
+linear_gr.OPTIONS.TRANGE = t(its:ite);
linear_gr.g_ky = real(w_ky(Is,:));
linear_gr.w_ky = imag(w_ky(Is,:));
linear_gr.ce = abs(ce);
@@ -44,20 +44,29 @@ linear_gr.avg_g = mean(real(w_ky(Is,:)),2);
linear_gr.std_w = std (imag(w_ky(Is,:)),[],2);
linear_gr.avg_w = mean(imag(w_ky(Is,:)),2);
-if PLOT >0
+if OPTIONS.NPLOTS >0
figure
-if PLOT > 1
+if OPTIONS.NPLOTS > 1
subplot(1,2,1)
end
+ x_ = linear_gr.ky;
+ plt = @(x) x;
+ OVERK = '';
+ if OPTIONS.GOK == 1
+ plt = @(x) x./x_;
+ OVERK = '/$k_y$';
+ elseif OPTIONS.GOK == 2
+ plt = @(x) x.^2./x_.^3;
+ OVERK = '/$k_y$';
+ end
+ plot(x_,plt(linear_gr.g_ky(:,end)),'-o','DisplayName',['$Re(\omega_{k_y})$',OVERK]); hold on;
+ plot(x_,plt(linear_gr.w_ky(:,end)),'--*','DisplayName',['$Im(\omega_{k_y})$',OVERK]); hold on;
+ plot(x_,plt(linear_gr.ce (:,end)),'x','DisplayName',['$\epsilon$',OVERK]); hold on;
- plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
- plot(linear_gr.ky,linear_gr.w_ky(:,end),'--*','DisplayName','$Im(\omega_{k_y})$'); hold on;
- plot(linear_gr.ky,linear_gr.ce (:,end),'x','DisplayName','$\epsilon$'); hold on;
-
- errorbar(linear_gr.ky,linear_gr.avg_g,linear_gr.std_g,...
+ errorbar(linear_gr.ky,plt(linear_gr.avg_g),plt(linear_gr.std_g),...
'-o','DisplayName','$Re(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
- errorbar(linear_gr.ky,linear_gr.avg_w,linear_gr.std_w,...
+ errorbar(linear_gr.ky,plt(linear_gr.avg_w),plt(linear_gr.std_w),...
'--*','DisplayName','$Im(\omega_{k_y})$',...
'LineWidth',1.5); hold on;
@@ -66,10 +75,10 @@ end
legend('show');
title(DATA.param_title);
-if PLOT > 1
- if PLOT == 2
+if OPTIONS.NPLOTS > 1
+ if OPTIONS.NPLOTS == 2
subplot(1,2,2)
- elseif PLOT == 3
+ elseif OPTIONS.NPLOTS == 3
subplot(2,2,2)
end
plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on;
@@ -77,7 +86,7 @@ if PLOT > 1
xlabel('t'); ylabel('$\gamma(t),\omega(t)$'); xlim([DATA.Ts3D(1) DATA.Ts3D(end)]);
end
-if PLOT > 2
+if OPTIONS.NPLOTS > 2
xlabel([]); xticks([]);
subplot(2,2,4)
[~,ikx0] = min(abs(DATA.kx));
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index a6a738ce2f1305a1f800235a7e400e5cea593ac2..4725ede3741e371fce409a6c9c85fb5398ec46c5 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -2,13 +2,13 @@
SHOW_FILM = 1;
field2plot ='phi';
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME = 40; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME = 1000; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin = 100;
-dt_ = 0.01;
+Tfin = 500;
+dt_ = 0.1;
Nstep = ceil(Tfin/dt_);
% Init tracers
-Np = 20; %number of tracers
+Np = 200; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
@@ -181,8 +181,7 @@ while(t_<Tfin && it <= Nstep)
u___ = linterp(u__0,u__1,a__i);
% push the particle
-% q = sign(-u___(3));
- q = -u___(3);
+ q = sign(-u___(3));
% q =1;
x_ = x_ + dt_*u___(1)*q;
y_ = y_ + dt_*u___(2)*q;
@@ -224,19 +223,20 @@ while(t_<Tfin && it <= Nstep)
end
scale = max(max(abs(F2P))); % Scaling to normalize
pclr = pcolor(XX_,YY_,F2P/scale);
- caxis([-1 1]); colormap(bluewhitered); colorbar;
+ caxis([-1 1]); colormap(bluewhitered); %colorbar;
set(pclr, 'edgecolor','none'); hold on; caxis([-1+dimmed,1+dimmed]); shading interp
for ip = 1:Np
ia0 = max(1,it-Na);
plot(Traj_x(ip,ia0:it),Traj_y(ip,ia0:it),'.','Color',color(ip,:)); hold on
end
for ip = 1:Np
- plot(Traj_x(ip, 1),Traj_y(ip, 1),'xk'); hold on
+% plot(Traj_x(ip, 1),Traj_y(ip, 1),'xk'); hold on
plot(Traj_x(ip,it),Traj_y(ip,it),'ok','MarkerFaceColor',color(ip,:)); hold on
end
title(['$t \approx$', sprintf('%.3d',ceil(data.Ts3D(itu_)))]);
axis equal
xlim([xmin xmax]); ylim([ymin ymax]);
+ xlabel('$x/\rho_s$'); ylabel('$y/\rho_s$');
drawnow
% Capture the plot as an image
frame = getframe(fig);
@@ -303,10 +303,14 @@ xlim(U_TIME + [0 Tfin]);
subplot(222);
itf = floor(Nt/2); %fit end time
% x^2 displacement
- plot(time_,mean(xtot.^2,1),'DisplayName','$\langle x.^2\rangle_p$'); hold on
+ x2tot = mean(xtot.^2,1);
+ plot(time_-time_(1),x2tot,'DisplayName','$\langle x.^2\rangle_p$'); hold on
fit = polyfit(time_(1:itf),mean(xtot(:,1:itf).^2,1),1);
- plot(time_,fit(1)*time_+fit(2),'--k'); hold on
- ylabel('$\langle x^2 \rangle_p$');
+ plot(time_-time_(1),fit(1)*time_+fit(2),'--k'); hold on
+ % Put a linear growth from the first point to check if it super/sub
+ % diffusive
+% loglog(time_-time_(1),(time_-time_(1))*x2tot(2)/(time_(2)-time_(1)),'--k'); hold on
+% ylabel('$\langle x^2 \rangle_p$');
% % y^2 displacement
% fit = polyfit(time_(1:itf),mean(ytot(:,1:itf).^2,1),1);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index b9d69af104a76dfbf926fd4f59ebef3161a25934..a036d18452ace4fd9aaca42f567028ae996b8663 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
-tw = [500 1000];
+tw = [00 5000];
fig = gcf;
axObjs = fig.Children;
@@ -19,13 +19,13 @@ end
% n1 = numel(X_);
[~,n0] = min(abs(X_-tw(1)));
[~,n1] = min(abs(X_-tw(2)));
-mvm = @(x) movmean(x,50);
-shift = X_(n0);
+mvm = @(x) movmean(x,1);
+shift = 0;%X_(n0);
% shift = 0;
+skip = 50;
figure;
-plot(X_(n0:end),Y_(n0:end));
-plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
+plot(mvm(X_(n0:skip:n1)-shift),mvm(Y_(n0:skip:n1))); hold on;
% t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
avg= mean(Y_(n0:n1)); dev = std(Y_(n0:n1));
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 0a4fee17d6e512396b661d009476b1d94a54a42f..37b545481408c1d8c53295ba825db9d53acea4f3 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -8,9 +8,24 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
cstart= h5readatt(filename,'/data/input','start_iframe3d');
data = zeros(numel(ky),numel(kx),numel(z),numel(time));
- for it = 1:numel(time)
- tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
- data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ [KX,KY] = meshgrid(kx,ky);
+
+ switch variablename
+ case 'vEx'
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = 1i.*KY.*(tmp.real + 1i * tmp.imaginary);
+ end
+ case 'vEy'
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/phi/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = -1i.*KX.*(tmp.real + 1i * tmp.imaginary);
+ end
+ otherwise
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+ data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
end
end
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index f10777e343c2e707b75c6b9b0c83f4f9a11b4c11..686686f88caa4865949adc63e23a339aa46432fc 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -65,11 +65,11 @@ subplot(224)
%% FCGK
P_ = 4; J_ = 2;
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
+% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';
% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5';
-kp = 2.0;
+kp = 0.0;
kp_a = h5read(mat_file_name,'/coordkperp');
[~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx);
@@ -99,10 +99,11 @@ if 0
mfns = {...
'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
% '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
- '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
+ '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
% '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
@@ -114,7 +115,7 @@ CONAME_A = {...
% 'FC 10 5 NFLR 4',...
% 'FC 10 5 NFLR 12',...
% 'FC 10 5 NFLR 12 k<2', ...
- 'FC 10 5 NFLR 30', ...
+ 'LD 6 3 NFLR 20', ...
% 'FC 4 2 NFLR 6',...
'FC 4 2 NFLR 12', ...
% 'Hacked SG A',...
@@ -151,13 +152,13 @@ end
if 0
%%
kperp= 1.5;
-mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
+mfns = {'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 10 5 NFLR 12'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'LD 6 3 NFLR 12'};
grow = {};
puls = {};
for j_ = 1:numel(mfns)
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 7b95c932b051e6473e71c699d3a461fe9fe8e14f..377c80fcaaf213f3867d70b5d423eaa133fedb9d 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,4 +1,4 @@
-function [ fig ] = plot_metric( data )
+function [ fig ] = plot_metric( data, options )
names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
'gxx','gxy','gyy','gyz','gzz','hatB','hatR','hatZ'};
@@ -8,23 +8,34 @@ for i_ = 1:numel(names)
namae = names{i_};
geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end
-fig = figure;
-subplot(311)
- for i = 1:5
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
- end
- xlim([min(data.z),max(data.z)]); legend('show'); title('MoNoLiT geometry');
+NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
+if NPLOT > 0
+ fig = figure;
+ if options.SHOW_METRICS
+ subplot(311)
+ for i = 1:5
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
-subplot(312)
- for i = 6:10
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
- end
- xlim([min(data.z),max(data.z)]); legend('show');
+ subplot(312)
+ for i = 6:10
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
-subplot(313)
- for i = 11:13
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ subplot(313)
+ for i = 11:13
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
end
- xlim([min(data.z),max(data.z)]); legend('show');
+ if options.SHOW_FLUXSURF
+ R = [geo_arrays(1,:,12),geo_arrays(1,1,12)];
+ Z = [geo_arrays(1,:,13),geo_arrays(1,1,13)];
+ plot(R,Z);
+ axis equal
+ end
+end
end
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index eb284c69a1d2333e5a792135a60e932b57f46ae1..3a7cc924e4f82bd750cabbe307a03b2284cecae0 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -22,7 +22,9 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
end
time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0));
-
+
+ fig = 0;
+if options.NPLOTS > 0
fig = figure;
subplot(211)
plot(time_seg,transp_seg_avg,'-'); hold on;
@@ -40,6 +42,9 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
plot(time_seg,transp_seg_avg,'-'); hold on;
xlabel('Averaging time'); ylabel('$\langle Q_x\rangle_{\tau}$');
legend(['$Q_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))])
-
+end
+Gx_infty_avg = mean(gamma);
+Gx_infty_std = std (gamma);
+disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 966b6fbdcd2d95dffe6887f69ddeea827d51b6a8..3b8be6c7d14205b12ca26c8e280ac8824a6ca5a1 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -18,6 +18,9 @@ GEOM.geom = ['''',GEOMETRY,''''];
GEOM.q0 = Q0; % q factor
GEOM.shear = SHEAR; % shear
GEOM.eps = EPS; % inverse aspect ratio
+GEOM.kappa = KAPPA; % elongation
+GEOM.delta = DELTA; % triangularity
+GEOM.zeta = ZETA; % squareness
% Model parameters
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
@@ -64,8 +67,9 @@ switch CO
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD'
% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
-% COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
- COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
+% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
+% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
+ COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
end
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index d42c78009fa34251d709f4b063e77475e25ba7cb..29304ef0a05df30e763f385b202b102b8142eaa4 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -30,6 +30,9 @@ fprintf(fid,[' geom = ', GEOM.geom,'\n']);
fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']);
fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']);
fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
+fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
+fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
+fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index bfd9e9bcb69e8331bbbb4e392504617d4fddeea6..076321686c2eb6fe668b8700d58820a0adb509f5 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -150,12 +150,14 @@ SUBROUTINE diagnose_full(kstep)
CALL creatg(fidres, "/data/metric", "Metric data")
CALL putarrnd(fidres, "/data/metric/gxx", gxx(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gxy", gxy(izs:ize,0:1), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxz", gxz(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gyy", gyy(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gyz", gyz(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gzz", gzz(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/hatR", hatR(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/hatB", hatB(izs:ize,0:1), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatB_NL", hatB_NL(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradxB", gradxB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradyB", gradyB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1, 1/))
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 665f18ba08fcc5ccb39437f953d76855a32b8618..72b71046f30035b4117cec67bd8d33f35731b812 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -9,17 +9,34 @@ module geometry
use fields
use basic
use calculus, ONLY: simpson_rule_z
-
+ use miller, ONLY: set_miller_parameters, get_miller
implicit none
PRIVATE
- ! Geometry parameters
+ ! Geometry input parameters
CHARACTER(len=16), &
PUBLIC, PROTECTED :: geom
- REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
- REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
- REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
- LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
- ! Geometrical operators
+ REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
+ REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
+ REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
+ REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
+ ! parameters for Miller geometry
+ REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation
+ REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr
+ REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity
+ REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
+ REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness
+ REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr
+
+ ! GENE unused additional parameters for miller_mod
+ REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0 ! meant to redistribute the points in z
+ REAL(dp), PUBLIC, PROTECTED :: major_R = 1 ! major radius
+ REAL(dp), PUBLIC, PROTECTED :: major_Z = 0 ! vertical elevation
+ REAL(dp), PUBLIC, PROTECTED :: dpdx_pm_geom = 0 ! amplitude mag. eq. pressure grad.
+ REAL(dp), PUBLIC, PROTECTED :: C_y = 0 ! defines y coordinate : Cy (q theta - phi)
+ REAL(dp), PUBLIC, PROTECTED :: C_xy = 0 ! defines x coordinate : B = Cxy Vx x Vy
+
+ ! Geometrical auxiliary variables
+ LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
! Curvature
REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
@@ -38,6 +55,7 @@ implicit none
REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition
INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R
+
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
eval_magnetic_geometry, set_ikx_zBC_map
@@ -47,7 +65,8 @@ CONTAINS
SUBROUTINE geometry_readinputs
! Read the input parameters
IMPLICIT NONE
- NAMELIST /GEOMETRY/ geom, q0, shear, eps
+ NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
+ kappa, s_kappa,delta, s_delta, zeta, s_zeta ! For miller
READ(lu_in,geometry)
IF(shear .NE. 0._dp) SHEARED = .true.
@@ -59,6 +78,7 @@ CONTAINS
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
INTEGER :: fid
+ real(dp) :: G1,G2,G3,Cx,Cy
! Allocate arrays
CALL geometry_allocate_mem
@@ -68,9 +88,6 @@ CONTAINS
call eval_1D_geometry
ELSE
SELECT CASE(geom)
- CASE('circular')
- IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
- call eval_circular_geometry
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
@@ -78,8 +95,15 @@ CONTAINS
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
SHEARED = .FALSE.
+ CASE('miller')
+ IF( my_id .eq. 0 ) WRITE(*,*) 'Miller geometry'
+ call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
+ call get_miller(eps,major_R,major_Z,q0,shear,alpha_MHD,edge_opt,&
+ C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
+ gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ,&
+ Ckxky,gradz_coeff)
CASE DEFAULT
- ERROR STOP 'Error stop: geometry not recognized!!'
+ STOP 'geometry not recognized!!'
END SELECT
ENDIF
!
@@ -87,18 +111,41 @@ CONTAINS
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
DO eo = 0,1
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx = ikxs, ikxe
- kx = kxarray(ikx)
- DO iz = izgs,izge
- kparray(iky, ikx, iz, eo) = &
- SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
- ! there is a factor 1/B from the normalization; important to match GENE
+ DO iz = izgs,izge
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx = ikxs, ikxe
+ kx = kxarray(ikx)
+ kparray(iky, ikx, iz, eo) = &
+ SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
+ ! there is a factor 1/B from the normalization; important to match GENE
ENDDO
- ENDDO
- ENDDO
+ ENDDO
+ ENDDO
+ ! Curvature operator (Frei et al. 2022 eq 2.15)
+ DO iz = izgs,izge
+ G1 = gxy(iz,eo)*gxy(iz,eo)-gxx(iz,eo)*gyy(iz,eo)
+ G2 = gxy(iz,eo)*gxz(iz,eo)-gxx(iz,eo)*gyz(iz,eo)
+ G3 = gyy(iz,eo)*gxz(iz,eo)-gxy(iz,eo)*gyz(iz,eo)
+ Cx = (G1*gradyB(iz,eo) + G2*gradzB(iz,eo))/hatB(iz,eo)
+ Cy = (G3*gradzB(iz,eo) - G1*gradxB(iz,eo))/hatB(iz,eo)
+
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz,eo) = 1._dp / jacobian(iz,eo) / hatB(iz,eo)
+ ! Factor in front of the nonlinear term
+ hatB_NL(iz,eo) = Jacobian(iz,eo)&
+ *(gxx(iz,eo)*gyy(iz,eo) - gxy(iz,eo)**2)/hatB(iz,eo)
+ ENDDO
ENDDO
+
+
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
@@ -116,7 +163,7 @@ CONTAINS
SUBROUTINE eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
- REAL(dp) :: z, kx, ky, alpha_MHD, Gx, Gy
+ REAL(dp) :: z, kx, ky, Gx, Gy
alpha_MHD = 0._dp
parity: DO eo = 0,1
@@ -146,8 +193,7 @@ CONTAINS
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp / hatR(iz,eo)
- hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
+ hatB(iz,eo) = 1._dp / hatR(iz,eo)
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
@@ -161,7 +207,7 @@ CONTAINS
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)* hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)/ hatB(iz,eo)
! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
@@ -176,220 +222,59 @@ CONTAINS
!--------------------------------------------------------------------------------
!
-
- SUBROUTINE eval_circular_geometry
- ! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
+ SUBROUTINE eval_zpinch_geometry
implicit none
- REAL(dp) :: chi, kx, ky, Gx, Gy
+ REAL(dp) :: z, kx, ky, alpha_MHD
+ alpha_MHD = 0._dp
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- chi = zarray(iz,eo) ! = chi
+ parity: DO eo = 0,1
+ zloop: DO iz = izgs,izge
+ z = zarray(iz,eo)
- ! metric in x,y,z
+ ! metric
gxx(iz,eo) = 1._dp
- gyy(iz,eo) = 1._dp + (shear*chi)**2 - 2._dp*eps*COS(chi) - 2._dp*shear*chi*eps*SIN(chi)
- gxy(iz,eo) = shear*chi - eps*SIN(chi)
- gxz(iz,eo) = -SIN(chi)
- gyz(iz,eo) = (1._dp - 2._dp*eps*COS(chi) - shear*chi*eps*SIN(chi))/eps
- gzz(iz,eo) = (1._dp - 2._dp*eps*COS(chi))/eps**2
- dxdR(iz,eo)= COS(chi)
- dxdZ(iz,eo)= SIN(chi)
+ gxy(iz,eo) = 0._dp
+ gxz(iz,eo) = 0._dp
+ gyy(iz,eo) = 1._dp
+ gyz(iz,eo) = 0._dp
+ gzz(iz,eo) = 1._dp
+ dxdR(iz,eo)= COS(z)
+ dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
- hatR(iz,eo) = 1._dp + eps*COS(chi)
- hatZ(iz,eo) = 1._dp + eps*SIN(chi)
+ hatR(iz,eo) = 1._dp
+ hatZ(iz,eo) = 1._dp
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = chi
+ phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
- Jacobian(iz,eo) = q0*hatR(iz,eo)
+ Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp / hatR(iz,eo)
- hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
+ hatB (iz,eo) = 1._dp
+ hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
- gradxB(iz,eo) = -COS(chi) ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = eps * SIN(chi) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
+ gradxB(iz,eo) = -1._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
- Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Kx
- Gy = -COS(chi) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Ky
- DO iky = ikys, ikye
+ DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = -ky
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
- ENDDO zloop
- ENDDO parity
-
- END SUBROUTINE eval_circular_geometry
- !
- !--------------------------------------------------------------------------------
- !
-
- SUBROUTINE eval_circular_geometry_GENE
- ! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
- implicit none
- REAL(dp) :: qbar, dxdr_, dpsidr, dqdr, dqbardr, Cy, dCydr_Cy, Cxy, fac2, &
- cost, sint, dchidr, dchidt, sign_Ip, B, dBdt, dBdr, &
- g11, g22, g12, g33, dBdchi, dBdr_c !GENE variables
- REAL(dp) :: X, z, kx, ky, Gamma1, Gamma2, Gamma3, feps
-
- sign_Ip = 1._dp
- feps = SQRT(1._dp-eps**2)
- qbar = q0 * feps
- dxdr_ = 1._dp
- dpsidr = eps/qbar
- dqdr = q0/eps*shear
- dqbardr = feps*q0/eps*(shear - eps**2/feps**2)
- Cy = eps/q0 * sign_Ip
- dCydr_Cy = 0._dp
- Cxy = 1._dp/feps
- fac2 = dCydr_Cy + dqdr/q0
-
-
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- z = zarray(iz,eo)
-
- cost = (COS(z)-eps)/(1._dp-eps*COS(z))
- sint = feps*sin(sign_Ip*z)/(1._dp-eps*COS(z))
- dchidr = -SIN(z)/feps**2
- dchidt = sign_Ip*feps/(1._dp+eps*cost)
- B = SQRT(1._dp+(eps/qbar)**2)/(1._dp+eps*cost)
- dBdt = eps*sint*B/(1._dp+eps*cost)
-
- ! metric in r,chi,Phi
- g11 = 1._dp
- g22 = dchidr**2 + dchidt**2/eps**2
- g12 = dchidr
- g33 = 1._dp/(1._dp+eps*cost)**2
-
- ! magnetic field derivatives in
- dBdchi=dBdt/dchidt
- dBdr_c=(dBdr-g12*dBdt/dchidt)/g11
-
- ! metric in x,y,z
- gxx(iz,eo) = dxdr_**2*g11
- gyy(iz,eo) = (Cy*q0)**2 * (fac2*z)**2*g11 &
- + 2._dp*fac2*z*g12 &
- + Cy**2*g33
- gxy(iz,eo) = dxdr_*Cy*sign_Ip*q0*(fac2*z*g11 + g12)
- gxz(iz,eo) = dxdr_*g12
- gyz(iz,eo) = Cy*q0*sign_Ip*(fac2*z*g12 + g22)
- gzz(iz,eo) = g22
-
- ! Jacobian
- Jacobian(iz,eo) = Cxy*abs(q0)*(1._dp+eps*cost)**2
-
- ! Background equilibrium magnetic field
- hatB (iz,eo) = B
- hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
- gradxB(iz,eo) = dBdr_c/dxdr_ ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = dBdchi! Gene put a factor hatB or 1/hatR in this
-
- ! Cylindrical coordinates derivatives
- dxdR(iz,eo) = cost
- dxdZ(iz,eo) = sint
- hatR(iz,eo) = 1._dp + eps*cost
- hatZ(iz,eo) = 1._dp + eps*sint
-
- ! toroidal coordinates
- Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = X
- Zc (iz,eo) = hatZ(iz,eo)
-
- Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo)
- Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo) ! Kx
- Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo) ! Ky
- ! Curvature operator
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- ! Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - (COS(z) + (shear*z - alpha_MHD*SIN(z))* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- Ckxky(iky, ikx, iz,eo) = (-sint*kx - cost*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- ENDDO
- ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
-
- ENDDO zloop
- ENDDO parity
-
- END SUBROUTINE eval_circular_geometry_GENE
- !
- !--------------------------------------------------------------------------------
- !
-
- SUBROUTINE eval_zpinch_geometry
- ! evaluate s-alpha geometry model
- implicit none
- REAL(dp) :: z, kx, ky, alpha_MHD
- alpha_MHD = 0._dp
-
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- z = zarray(iz,eo)
-
- ! metric
- gxx(iz,eo) = 1._dp
- gxy(iz,eo) = 0._dp
- gxz(iz,eo) = 0._dp
- gyy(iz,eo) = 1._dp
- gyz(iz,eo) = 0._dp
- gzz(iz,eo) = 1._dp
- dxdR(iz,eo)= COS(z)
- dxdZ(iz,eo)= SIN(z)
-
- ! Relative strengh of radius
- hatR(iz,eo) = 1._dp
- hatZ(iz,eo) = 1._dp
-
- ! toroidal coordinates
- Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = z
- Zc (iz,eo) = hatZ(iz,eo)
-
- ! Jacobian
- Jacobian(iz,eo) = 1._dp
-
- ! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp
- hatB_NL(iz,eo) = 1._dp
-
- ! Derivative of the magnetic field strenght
- gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
-
- ! Curvature operator
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- ENDDO
- ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
-
- ENDDO zloop
- ENDDO parity
+ ENDDO zloop
+ ENDDO parity
END SUBROUTINE eval_zpinch_geometry
!
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 719effe7f116cce087e250a70a99aab4db4f553d..9bba37ad772d93d33a3e345fc6e206c4fdcbe578 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -281,8 +281,8 @@ CONTAINS
jmaxe_dp = real(jmaxe,dp)
jmaxi_dp = real(jmaxi,dp)
! j=0 indices
- DO ij = ijs_e,ij0_e; IF(jarray_e(ij) .EQ. 0) ij0_e = ij; END DO
- DO ij = ijs_i,ij0_i; IF(jarray_i(ij) .EQ. 0) ij0_i = ij; END DO
+ DO ij = ijs_e,ije_e; IF(jarray_e(ij) .EQ. 0) ij0_e = ij; END DO
+ DO ij = ijs_i,ije_i; IF(jarray_i(ij) .EQ. 0) ij0_i = ij; END DO
END SUBROUTINE set_jgrid
diff --git a/src/inital.F90 b/src/inital.F90
index 8ccb5f201ced5420af4d7b7c9c9dfeac13c8a060..6d4fe5af802c7b7298e73c65312a7ed81aa61306 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -243,7 +243,7 @@ SUBROUTINE init_gyrodens
END DO
END DO
IF ( contains_ky0 ) THEN
- DO iky=2,Nky/2 !symmetry at ky = 0 for all z
+ DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
moments_i( ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:,:)
END DO
ENDIF
diff --git a/src/lag_interp_mod.F90 b/src/lag_interp_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..8415df5d63783ba680b21431aec8dbbeb3fe55e5
--- /dev/null
+++ b/src/lag_interp_mod.F90
@@ -0,0 +1,256 @@
+!! This source is taken from GENE https://genecode.org/ !!
+!>lagrange_interpolation contains subroutines to perform
+!!a mid-point lagrange interpolation of order 3
+MODULE lagrange_interpolation
+ USE prec_const
+ IMPLICIT NONE
+ PUBLIC:: lag3interp, lag3deriv, lag3interp_2d
+ PUBLIC:: lag3interp_complex
+ PRIVATE
+
+ INTERFACE lag3interp
+ MODULE PROCEDURE lag3interp_scalar, lag3interp_array
+ END INTERFACE
+
+ INTERFACE lag3deriv
+ MODULE PROCEDURE lag3deriv_scalar, lag3deriv_array
+ END INTERFACE
+
+CONTAINS
+
+ !> Third order lagrange interpolation
+ SUBROUTINE lag3interp_scalar(y_in,x_in,n_in,y_out,x_out)
+ INTEGER, INTENT(IN) :: n_in
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), INTENT(IN) :: x_out
+ REAL(dp), INTENT(OUT) :: y_out
+
+ REAL(dp), DIMENSION(1) :: xout_wrap, yout_wrap
+
+ xout_wrap = x_out
+ call lag3interp_array(y_in,x_in,n_in,yout_wrap,xout_wrap,1)
+ y_out = yout_wrap(1)
+
+ END SUBROUTINE lag3interp_scalar
+
+ !> Third order lagrange interpolation
+ subroutine lag3interp_array(y_in,x_in,n_in,y_out,x_out,n_out)
+ INTEGER, INTENT(IN) :: n_in,n_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(dp), DIMENSION(n_out), INTENT(OUT) :: y_out
+
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ INTEGER :: j,jm,j0,j1,j2
+ INTEGER :: jstart,jfirst,jlast,jstep
+
+ IF (x_in(n_in) > x_in(1)) THEN
+ jstart=3
+ jfirst=1
+ jlast=n_out
+ jstep=1
+ ELSE
+ jstart=n_in-2
+ jfirst=n_out
+ jlast=1
+ jstep=-1
+ END IF
+
+ j1=jstart
+ DO j=jfirst,jlast,jstep
+ x=x_out(j)
+ DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
+ j1=j1+jstep
+ END DO
+
+ j2=j1+jstep
+ j0=j1-jstep
+ jm=j1-2*jstep
+
+ !... extrapolate inside or outside
+
+ x2=x_in(j2)
+ x1=x_in(j1)
+ x0=x_in(j0)
+ xm=x_in(jm)
+
+ aintm=(x-x0)*(x-x1)*(x-x2)/((xm-x0)*(xm-x1)*(xm-x2))
+ aint0=(x-xm)*(x-x1)*(x-x2)/((x0-xm)*(x0-x1)*(x0-x2))
+ aint1=(x-xm)*(x-x0)*(x-x2)/((x1-xm)*(x1-x0)*(x1-x2))
+ aint2=(x-xm)*(x-x0)*(x-x1)/((x2-xm)*(x2-x0)*(x2-x1))
+
+ y_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+ +aint1*y_in(j1)+aint2*y_in(j2)
+
+ END DO
+
+ END SUBROUTINE Lag3interp_array
+
+
+ !> Third order lagrange interpolation for complex arrays
+ SUBROUTINE lag3interp_complex(y_in,x_in,n_in,y_out,x_out,n_out)
+ INTEGER, INTENT(IN) :: n_in,n_out
+ COMPLEX, DIMENSION(n_in), INTENT(IN) :: y_in
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: x_in
+ COMPLEX, DIMENSION(n_out), INTENT(OUT) :: y_out
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ INTEGER :: j,jm,j0,j1,j2
+ INTEGER :: jstart,jfirst,jlast,jstep
+
+ IF (x_in(n_in) > x_in(1)) THEN
+ jstart=3
+ jfirst=1
+ jlast=n_out
+ jstep=1
+ ELSE
+ jstart=n_in-2
+ jfirst=n_out
+ jlast=1
+ jstep=-1
+ END IF
+
+ j1=jstart
+ DO j=jfirst,jlast,jstep
+ x=x_out(j)
+ DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
+ j1=j1+jstep
+ END DO
+
+ j2=j1+jstep
+ j0=j1-jstep
+ jm=j1-2*jstep
+
+ !... extrapolate inside or outside
+
+ x2=x_in(j2)
+ x1=x_in(j1)
+ x0=x_in(j0)
+ xm=x_in(jm)
+
+ aintm=(x-x0)*(x-x1)*(x-x2)/((xm-x0)*(xm-x1)*(xm-x2))
+ aint0=(x-xm)*(x-x1)*(x-x2)/((x0-xm)*(x0-x1)*(x0-x2))
+ aint1=(x-xm)*(x-x0)*(x-x2)/((x1-xm)*(x1-x0)*(x1-x2))
+ aint2=(x-xm)*(x-x0)*(x-x1)/((x2-xm)*(x2-x0)*(x2-x1))
+
+ y_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+ +aint1*y_in(j1)+aint2*y_in(j2)
+
+ END DO
+
+ END SUBROUTINE lag3interp_complex
+
+
+ !>2D interpolation
+ !\TODO check whether a "REAL(dp)" 2D interpolation would
+ !! be more appropriate
+ SUBROUTINE lag3interp_2d(y_in,x1_in,n1_in,x2_in,n2_in,&
+ &y_out,x1_out,n1_out,x2_out,n2_out)
+ INTEGER, INTENT(IN) :: n1_in,n2_in,n1_out,n2_out
+ REAL(dp), DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
+ REAL(dp), DIMENSION(n1_in) :: x1_in
+ REAL(dp), DIMENSION(n2_in) :: x2_in
+ REAL(dp), DIMENSION(n1_out), INTENT(IN) :: x1_out
+ REAL(dp), DIMENSION(n2_out), INTENT(IN) :: x2_out
+ REAL(dp), DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
+
+ !local variables
+ REAL(dp), DIMENSION(n2_in) :: y2_in_tmp
+ REAL(dp), DIMENSION(n2_out) :: y2_out_tmp
+ REAL(dp), DIMENSION(n1_in,n2_out) :: y_tmp
+ INTEGER :: i
+
+ DO i=1,n1_in
+ y2_in_tmp = y_in(i,:)
+ call lag3interp(y2_in_tmp,x2_in,n2_in,&
+ y2_out_tmp,x2_out,n2_out)
+ y_tmp(i,:) = y2_out_tmp
+ ENDDO
+
+ DO i=1,n2_out
+ call lag3interp(y_tmp(:,i),x1_in,n1_in,&
+ y_out(:,i),x1_out,n1_out)
+ END DO
+
+ END SUBROUTINE lag3interp_2d
+
+ !> Third order lagrange interpolation
+ SUBROUTINE lag3deriv_scalar(y_in,x_in,n_in,dydx_out,x_out)
+
+ IMPLICIT NONE
+
+ INTEGER, INTENT(IN) :: n_in
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), INTENT(IN) :: x_out
+ REAL(dp), INTENT(OUT) :: dydx_out
+
+ REAL(dp), DIMENSION(1) :: xout_wrap, dydxout_wrap
+
+ xout_wrap = x_out
+ call lag3deriv_array(y_in,x_in,n_in,dydxout_wrap,xout_wrap,1)
+ dydx_out = dydxout_wrap(1)
+
+ END SUBROUTINE lag3deriv_scalar
+
+
+
+!>Returns Derivative based on a 3rd order lagrange interpolation
+ subroutine lag3deriv_array(y_in,x_in,n_in,dydx_out,x_out,n_out)
+ INTEGER :: n_in,n_out
+ REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
+ REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
+ REAL(dp), DIMENSION(n_out), INTENT(OUT) :: dydx_out
+
+ REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
+ INTEGER :: j,jm,j0,j1,j2
+ INTEGER :: jstart,jfirst,jlast,jstep
+
+ IF (x_in(n_in) > x_in(1)) THEN
+ jstart=3
+ jfirst=1
+ jlast=n_out
+ jstep=1
+ ELSE
+ jstart=n_in-2
+ jfirst=n_out
+ jlast=1
+ jstep=-1
+ END IF
+
+ j1=jstart
+ DO j=jfirst,jlast,jstep
+ x=x_out(j)
+ DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
+ j1=j1+jstep
+ END DO
+
+ j2=j1+jstep
+ j0=j1-jstep
+ jm=j1-2*jstep
+
+ !... extrapolate inside or outside
+
+ x2=x_in(j2)
+ x1=x_in(j1)
+ x0=x_in(j0)
+ xm=x_in(jm)
+
+ aintm=((x-x1)*(x-x2)+(x-x0)*(x-x2)+(x-x0)*(x-x1)) &
+ /((xm-x0)*(xm-x1)*(xm-x2))
+ aint0=((x-x1)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x1)) &
+ /((x0-xm)*(x0-x1)*(x0-x2))
+ aint1=((x-x0)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x0)) &
+ /((x1-xm)*(x1-x0)*(x1-x2))
+ aint2=((x-x0)*(x-x1)+(x-xm)*(x-x1)+(x-xm)*(x-x0)) &
+ /((x2-xm)*(x2-x0)*(x2-x1))
+
+ dydx_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+ +aint1*y_in(j1)+aint2*y_in(j2)
+
+ END DO
+
+ end subroutine Lag3deriv_array
+
+
+end module lagrange_interpolation
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..eb3a9e5f64627815b1b753b3381fead035a46a10
--- /dev/null
+++ b/src/miller_mod.F90
@@ -0,0 +1,667 @@
+!! This source has been adapted from GENE https://genecode.org/ !!
+!>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998)
+!>and [J. Candy, PPCF 51, 105009 (2009)]
+MODULE miller
+ USE prec_const
+ USE basic
+ ! use coordinates,only: gcoor, get_dzprimedz
+ USE grid
+ ! use discretization
+ USE lagrange_interpolation
+ ! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol
+ USE model
+
+ implicit none
+ public:: get_miller, set_miller_parameters
+ public:: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
+ public:: thetaShift
+ public:: mMode, nMode
+ public:: thetak, thetad
+ public:: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
+ public:: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
+ public:: cN_m, sN_m, cNdr_m, sNdr_m
+
+ private
+
+ real(dp) :: rho, kappa, delta, s_kappa, s_delta, drR, drZ, zeta, s_zeta
+
+ INTEGER :: mMode, nMode
+ real(dp) :: thetaShift
+ real(dp) :: thetak, thetad
+ real(dp) :: aSurf, Delta2, Delta3, theta2, theta3, Raxis, Zaxis
+ real(dp) :: Deltam, Deltan, s_Deltam, s_Deltan, thetam, thetan
+
+ INTEGER, PARAMETER :: IND_M=32
+ real(dp), DIMENSION(0:IND_M-1) :: cN_m, sN_m, cNdr_m, sNdr_m
+
+CONTAINS
+
+ !>Set defaults for miller parameters
+ subroutine set_miller_parameters(kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_)
+ real(dp), INTENT(IN) :: kappa_,s_kappa_,delta_,s_delta_,zeta_,s_zeta_
+ rho = -1.0
+ kappa = kappa_
+ s_kappa = s_kappa_
+ delta = delta_
+ s_delta = s_delta_
+ zeta = zeta_
+ s_zeta = s_zeta_
+ drR = 0.0
+ drZ = 0.0
+
+ thetak = 0.0
+ thetad = 0.0
+
+ aSurf = 0.54
+ Delta2 = 1.0
+ Delta3 = 1.0
+ theta2 = 0.0
+ theta3 = 0.0
+ Raxis = 1.0
+ Zaxis = 0.0
+
+ mMode = 2
+ nMode = 3
+ Deltam = 1.0
+ Deltan = 1.0
+ s_Deltam = 0.0
+ s_Deltan = 0.0
+ thetam = 0.0
+ thetan = 0.0
+
+ cN_m = 0.0
+ sN_m = 0.0
+ cNdr_m = 0.0
+ sNdr_m = 0.0
+ end subroutine set_miller_parameters
+
+ !>Get Miller metric, magnetic field, jacobian etc.
+ subroutine get_miller(trpeps,major_R,major_Z,q0,shat,amhd,edge_opt,&
+ C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,&
+ Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_,Ckxky_,gradz_coeff_)
+ !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
+ real(dp), INTENT(INOUT) :: major_R ! major radius
+ real(dp), INTENT(INOUT) :: major_Z ! major Z
+ real(dp), INTENT(INOUT) :: q0 ! safetyfactor
+ real(dp), INTENT(INOUT) :: shat ! safetyfactor
+ real(dp), INTENT(INOUT) :: amhd ! alpha mhd
+ real(dp), INTENT(INOUT) :: edge_opt ! alpha mhd
+ real(dp), INTENT(INOUT) :: dpdx_pm_geom ! amplitude mag. eq. pressure grad.
+ real(dp), INTENT(INOUT) :: C_y, C_xy
+
+ real(dp), dimension(izgs:izge,0:1), INTENT(INOUT) :: &
+ gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
+ dBdx_,dBdy_,Bfield_,jacobian_,&
+ dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
+ gradz_coeff_
+ real(dp), dimension(ikys:ikye,ikxs:ikxe,izgs:izge,0:1), INTENT(INOUT) :: Ckxky_
+ ! No parameter in gyacomo yet
+ integer :: ikxgs =1 ! the left ghost gpdisc%pi1gl
+ real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
+ real(dp) :: sign_Bt_CW=1 ! current sign (only normal current)
+ !!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Auxiliary variables for curvature computation
+ real(dp) :: G1,G2,G3,Cx,Cy,ky,kx
+
+ integer:: np, np_s, Npol_ext, Npol_s
+
+ real(dp), dimension(500*(Npol+2)):: R,Z,R_rho,Z_rho,R_theta,Z_theta,R_theta_theta,Z_theta_theta,dlp,Rc,cosu,sinu,Bphi
+ real(dp), dimension(500*(Npol+2)):: drRcirc, drRelong, drRelongTilt, drRtri, drRtriTilt, drZcirc, drZelong, drZelongTilt
+ real(dp), dimension(500*(Npol+2)):: drZtri, drZtriTilt, dtdtRcirc, dtdtRelong, dtdtRelongTilt, dtdtRtri, dtdtRtriTilt
+ real(dp), dimension(500*(Npol+2)):: dtdtZcirc, dtdtZelong, dtdtZelongTilt, dtdtZtri, dtdtZtriTilt, dtRcirc, dtRelong
+ real(dp), dimension(500*(Npol+2)):: dtRelongTilt, dtRtri, dtRtriTilt, dtZcirc, dtZelong, dtZelongTilt, dtZtri, dtZtriTilt
+ real(dp), dimension(500*(Npol+2)):: Rcirc, Relong, RelongTilt, Rtri, RtriTilt, Zcirc, Zelong, ZelongTilt, Ztri, ZtriTilt
+ real(dp), dimension(500*(Npol+2)):: drrShape, drrAng, drxAng, dryAng, dtdtrShape, dtdtrAng, dtdtxAng
+ real(dp), dimension(500*(Npol+2)):: dtdtyAng, dtrShape, dtrAng, dtxAng, dtyAng, rShape, rAng, xAng, yAng
+ real(dp), dimension(500*(Npol+2)):: theta, thAdj, J_r, B, Bp, D0, D1, D2, D3, nu, chi, psi1, nu1
+ real(dp), dimension(500*(Npol+2)):: tmp_reverse, theta_reverse, tmp_arr
+
+ real(dp), dimension(500*(Npol+1)):: theta_s, thAdj_s, chi_s, theta_s_reverse
+ real(dp), dimension(500*(Npol+1)):: R_s, Z_s, R_theta_s, Z_theta_s, Rc_s, cosu_s, sinu_s, Bphi_s, B_s, Bp_s, dlp_s
+ real(dp), dimension(500*(Npol+1)):: dtRcirc_s, dtRelong_s, dtRelongTilt_s, dtRtri_s, dtRtriTilt_s, dtZcirc_s
+ real(dp), dimension(500*(Npol+1)):: dtZelong_s, dtZelongTilt_s, dtZtri_s, dtZtriTilt_s, Rcirc_s, Relong_s, RelongTilt_s
+ real(dp), dimension(500*(Npol+1)):: Rtri_s, RtriTilt_s, Zcirc_s, Zelong_s, ZelongTilt_s, Ztri_s, ZtriTilt_s, dtrShape_s
+ real(dp), dimension(500*(Npol+1)):: dtrAng_s, dtxAng_s, dtyAng_s, rShape_s, rAng_s, xAng_s, yAng_s
+ real(dp), dimension(500*(Npol+1)):: psi1_s, nu1_s, dchidx_s, dB_drho_s, dB_dl_s, dnu_drho_s, dnu_dl_s, grad_nu_s
+ real(dp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
+ real(dp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
+
+ real(dp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
+ real(dp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out
+ real(dp):: d_inv, drPsi, dxPsi, dq_dx, dq_dpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
+ !real(dp) :: Lnorm, Psi0 ! currently module-wide defined anyway
+ real(dp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
+ real(dp):: rho_a, psiN, drpsiN, CN2, CN3, Rcenter, Zcenter, drRcenter, drZcenter
+ logical:: bMaxShift
+ integer:: i, k, iBmax
+
+ Npol_ext = Npol+2
+ Npol_s = Npol+1
+ np = 500*Npol_ext
+ np_s = 500*Npol_s
+
+ rho = trpeps*major_R
+ if (rho.le.0.0) stop 'flux surface radius not defined'
+ trpeps = rho/major_R
+
+ q0 = sign_Ip_CW * sign_Bt_CW * abs(q0)
+
+ R0=major_R
+ B0=1.0*sign_Bt_CW
+ F=R0*B0
+ Z0=major_Z
+ pi = acos(-1.0)
+ mu_0=4.0*pi
+
+ theta=linspace(-pi*Npol_ext,pi*Npol_ext-2._dp*pi*Npol_ext/np,np)
+ d_inv=asin(delta)
+
+ thetaShift = 0.0
+ iBmax = 1
+
+ !flux surface parametrization
+ thAdj = theta + thetaShift
+ if (zeta/=0.0 .or. s_zeta/=0.0) then
+ R = R0 + rho*Cos(thAdj + d_inv*Sin(thAdj))
+ Z = Z0 + kappa*rho*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_rho = drR + Cos(thAdj + d_inv*Sin(thAdj)) - s_delta*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_rho = drZ + kappa*s_zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ + kappa*Sin(thAdj + zeta*Sin(2*thAdj)) + kappa*s_kappa*Sin(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta = -(rho*(1 + d_inv*Cos(thAdj))*Sin(thAdj + d_inv*Sin(thAdj)))
+ Z_theta = kappa*rho*(1 + 2*zeta*Cos(2*thAdj))*Cos(thAdj + zeta*Sin(2*thAdj))
+
+ R_theta_theta = -(rho*(1 + d_inv*Cos(thAdj))**2*Cos(thAdj + d_inv*Sin(thAdj))) &
+ + d_inv*rho*Sin(thAdj)*Sin(thAdj + d_inv*Sin(thAdj))
+ Z_theta_theta = -4*kappa*rho*zeta*Cos(thAdj + zeta*Sin(2*thAdj))*Sin(2*thAdj) &
+ - kappa*rho*(1 + 2*zeta*Cos(2*thAdj))**2*Sin(thAdj + zeta*Sin(2*thAdj))
+ else
+ Rcirc = rho*Cos(thAdj - thetad + thetak)
+ Zcirc = rho*Sin(thAdj - thetad + thetak)
+ Relong = Rcirc
+ Zelong = Zcirc + (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ RelongTilt = Relong*Cos(thetad - thetak) - Zelong*Sin(thetad - thetak)
+ ZelongTilt = Zelong*Cos(thetad - thetak) + Relong*Sin(thetad - thetak)
+ Rtri = RelongTilt - rho*Cos(thAdj) + rho*Cos(thAdj + delta*Sin(thAdj))
+ Ztri = ZelongTilt
+ RtriTilt = Rtri*Cos(thetad) + Ztri*Sin(thetad)
+ ZtriTilt = Ztri*Cos(thetad) - Rtri*Sin(thetad)
+ R = R0 + RtriTilt
+ Z = Z0 + ZtriTilt
+
+ drRcirc = Cos(thAdj - thetad + thetak)
+ drZcirc = Sin(thAdj - thetad + thetak)
+ drRelong = drRcirc
+ drZelong = drZcirc - (1 - kappa - kappa*s_kappa)*Sin(thAdj - thetad + thetak)
+ drRelongTilt = drRelong*Cos(thetad - thetak) - drZelong*Sin(thetad - thetak)
+ drZelongTilt = drZelong*Cos(thetad - thetak) + drRelong*Sin(thetad - thetak)
+ drRtri = drRelongTilt - Cos(thAdj) + Cos(thAdj + delta*Sin(thAdj)) &
+ - s_delta*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ drZtri = drZelongTilt
+ drRtriTilt = drRtri*Cos(thetad) + drZtri*Sin(thetad)
+ drZtriTilt = drZtri*Cos(thetad) - drRtri*Sin(thetad)
+ R_rho = drR + drRtriTilt
+ Z_rho = drZ + drZtriTilt
+
+ dtRcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtZcirc = rho*Cos(thAdj - thetad + thetak)
+ dtRelong = dtRcirc
+ dtZelong = dtZcirc + (-1 + kappa)*rho*Cos(thAdj - thetad + thetak)
+ dtRelongTilt = dtRelong*Cos(thetad - thetak) - dtZelong*Sin(thetad - thetak)
+ dtZelongTilt = dtZelong*Cos(thetad - thetak) + dtRelong*Sin(thetad - thetak)
+ dtRtri = dtRelongTilt + rho*Sin(thAdj) - rho*(1 + delta*Cos(thAdj))*Sin(thAdj + delta*Sin(thAdj))
+ dtZtri = dtZelongTilt
+ dtRtriTilt = dtRtri*Cos(thetad) + dtZtri*Sin(thetad)
+ dtZtriTilt = dtZtri*Cos(thetad) - dtRtri*Sin(thetad)
+ R_theta = dtRtriTilt
+ Z_theta = dtZtriTilt
+
+ dtdtRcirc = -(rho*Cos(thAdj - thetad + thetak))
+ dtdtZcirc = -(rho*Sin(thAdj - thetad + thetak))
+ dtdtRelong = dtdtRcirc
+ dtdtZelong = dtdtZcirc - (-1 + kappa)*rho*Sin(thAdj - thetad + thetak)
+ dtdtRelongTilt = dtdtRelong*Cos(thetad - thetak) - dtdtZelong*Sin(thetad - thetak)
+ dtdtZelongTilt = dtdtZelong*Cos(thetad - thetak) + dtdtRelong*Sin(thetad - thetak)
+ dtdtRtri = dtdtRelongTilt + rho*Cos(thAdj) - rho*(1 + delta*Cos(thAdj))**2*Cos(thAdj + delta*Sin(thAdj)) &
+ + delta*rho*Sin(thAdj)*Sin(thAdj + delta*Sin(thAdj))
+ dtdtZtri = dtdtZelongTilt
+ dtdtRtriTilt = dtdtRtri*Cos(thetad) + dtdtZtri*Sin(thetad)
+ dtdtZtriTilt = dtdtZtri*Cos(thetad) - dtdtRtri*Sin(thetad)
+ R_theta_theta = dtdtRtriTilt
+ Z_theta_theta = dtdtZtriTilt
+ endif
+
+ !dl/dtheta
+ dlp=(R_theta**2+Z_theta**2)**0.5
+
+ !curvature radius
+ Rc=dlp**3*(R_theta*Z_theta_theta-Z_theta*R_theta_theta)**(-1)
+
+ ! some useful quantities (see papers for definition of u)
+ cosu=Z_theta/dlp
+ sinu=-R_theta/dlp
+
+ !Jacobian J_r = (dPsi/dr) J_psi = (dPsi/dr) / [(nabla fz x nabla psi)* nabla theta]
+ ! = R * (dR/drho dZ/dtheta - dR/dtheta dZ/drho) = R dlp / |nabla r|
+ J_r=R*(R_rho*Z_theta-R_theta*Z_rho)
+
+ !From definition of q = 1/(2 pi) int (B nabla fz) / (B nabla theta) dtheta:
+ !dPsi/dr = sign_Bt sign_Ip / (2 pi q) int F / R^2 J_r dtheta
+ ! = F / (2 pi |q|) int J_r/R^2 dtheta
+ tmp_arr=J_r/R**2
+ drPsi=sign_Ip_CW*F/(2.*pi*Npol_ext*q0)*sum(tmp_arr)*2*pi*Npol_ext/np !dlp_int(tmp_arr,1.0)
+
+ !Poloidal field (Bp = Bvec * nabla l)
+ Bp=sign_Ip_CW * drPsi / J_r * dlp
+
+ !toroidal field
+ Bphi=F/R
+
+ !total modulus of Bfield
+ B=sqrt(Bphi**2+Bp**2)
+
+ bMaxShift = .false.
+ ! if (thetaShift==0.0.and.trim(magn_geometry).ne.'miller_general') then
+ if (thetaShift==0.0) then
+ do i = 2,np-1
+ if (B(iBmax)<B(i)) then
+ iBmax = i
+ end if
+ enddo
+ if (iBmax/=1) then
+ bMaxShift = .true.
+ thetaShift = theta(iBmax)-theta(1)
+ end if
+ end if
+
+ !definition of radial coordinate! dx_drho=1 --> x = r
+ dx_drho=1. !drPsi/Psi0*Lnorm*q0
+ if (my_id==0) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho
+ dxPsi = drPsi/dx_drho
+ C_y = dxPsi*sign_Ip_CW
+ C_xy = abs(B0*dxPsi/C_y)
+
+ if (my_id==0) then
+ write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") &
+ "Setting C_xy = ",C_xy,' C_y = ', C_y," C_x' = ", 1./dxPsi
+ write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi
+ write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi
+ if (thetaShift.ne.0.0) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift
+ endif
+
+
+ !--------shear is expected to be defined as rho/q*dq/drho--------!
+ dq_dx=shat*q0/rho/dx_drho
+ dq_dpsi=dq_dx/dxPsi
+ pprime=-amhd/q0**2/R0/(2*mu_0)*B0**2/drPsi
+
+ !neg. dpdx normalized to magnetic pressure for pressure term
+ dpdx_pm_geom=amhd/q0**2/R0/dx_drho
+
+ !first coefficient of psi in varrho expansion
+ psi1 = R*Bp*sign_Ip_CW
+
+ !integrals for ffprime evaluation
+ do i=1,np
+ tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
+ D0(i)=-F*dlp_int_ind(tmp_arr,dlp,i)
+ tmp_arr=B**2*R/psi1**3
+ D1(i)=-dlp_int_ind(tmp_arr,dlp,i)/F
+ tmp_arr=mu_0*R/psi1**3
+ D2(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
+ tmp_arr=1./(R*psi1)
+ D3(i)=-dlp_int_ind(tmp_arr,dlp,i)*F
+ enddo
+ tmp_arr=(2./Rc-2.*cosu/R)/(R*psi1**2)
+ D0_full=-F*dlp_int(tmp_arr,dlp)
+ tmp_arr=B**2*R/psi1**3
+ D1_full=-dlp_int(tmp_arr,dlp)/F
+ tmp_arr=mu_0*R/psi1**3
+ D2_full=-dlp_int(tmp_arr,dlp)*F
+ tmp_arr=1./(R*psi1)
+ D3_full=-dlp_int(tmp_arr,dlp)*F
+ D0_mid=D0(np/2+1)
+ D1_mid=D1(np/2+1)
+ D2_mid=D2(np/2+1)
+ D3_mid=D3(np/2+1)
+
+ ffprime=-(sign_Ip_CW*dq_dpsi*2.*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full
+
+ if (my_id==0) then
+ write(*,'(A,ES12.4)') "ffprime = ", ffprime
+ endif
+ D0=D0-D0_mid
+ D1=D1-D1_mid
+ D2=D2-D2_mid
+ nu=D3-D3_mid
+
+ nu1=psi1*(D0+D1*ffprime+D2*pprime)
+
+ !straight field line angle defined on equidistant theta grid
+ !alpha = fz + nu = - (q chi - fz) => chi = -nu / q
+ chi=-nu/q0
+
+ !correct small scaling error (<0.5%, due to finite integration resolution)
+ chi=chi*(maxval(theta)-minval(theta))/(maxval(chi)-minval(chi))
+
+ !new grid equidistant in straight field line angle
+ chi_s = linspace(-pi*Npol_s,pi*Npol_s-2*pi*Npol_s/np_s,np_s)
+
+ if (sign_Ip_CW.lt.0.0) then !make chi increasing function to not confuse lag3interp
+ tmp_reverse = chi(np:1:-1)
+ theta_reverse = theta(np:1:-1)
+ call lag3interp(theta_reverse,tmp_reverse,np,theta_s,chi_s,np_s)
+ theta_s_reverse = theta_s(np_s:1:-1)
+ else
+ !lag3interp(y_in,x_in,n_in,y_out,x_out,n_out)
+ call lag3interp(theta,chi,np,theta_s,chi_s,np_s)
+ endif
+ dtheta_dchi_s=deriv_fd(theta_s,chi_s,np_s)
+
+ !arrays equidistant in straight field line angle
+ thAdj_s = theta_s + thetaShift
+
+ if (zeta/=0.0 .or. s_zeta/=0.0) then
+ R_s = R0 + rho*Cos(thAdj_s + d_inv*Sin(thAdj_s))
+ Z_s = Z0 + kappa*rho*Sin(thAdj_s + zeta*Sin(2*thAdj_s))
+
+ R_theta_s = -(dtheta_dchi_s*rho*(1 + d_inv*Cos(thAdj_s))*Sin(thAdj_s + d_inv*Sin(thAdj_s)))
+ Z_theta_s = dtheta_dchi_s*kappa*rho*(1 + 2*zeta*Cos(2*thAdj_s))*Cos(thAdj_s + zeta*Sin(2*thAdj_s))
+ else
+ Rcirc_s = rho*Cos(thAdj_s - thetad + thetak)
+ Zcirc_s = rho*Sin(thAdj_s - thetad + thetak)
+ Relong_s = Rcirc_s
+ Zelong_s = Zcirc_s + (-1 + kappa)*rho*Sin(thAdj_s - thetad + thetak)
+ RelongTilt_s = Relong_s*Cos(thetad - thetak) - Zelong_s*Sin(thetad - thetak)
+ ZelongTilt_s = Zelong_s*Cos(thetad - thetak) + Relong_s*Sin(thetad - thetak)
+ Rtri_s = RelongTilt_s - rho*Cos(thAdj_s) + rho*Cos(thAdj_s + delta*Sin(thAdj_s))
+ Ztri_s = ZelongTilt_s
+ RtriTilt_s = Rtri_s*Cos(thetad) + Ztri_s*Sin(thetad)
+ ZtriTilt_s = Ztri_s*Cos(thetad) - Rtri_s*Sin(thetad)
+ R_s = R0 + RtriTilt_s
+ Z_s = Z0 + ZtriTilt_s
+
+ dtRcirc_s = -(rho*Sin(thAdj_s - thetad + thetak))
+ dtZcirc_s = rho*Cos(thAdj_s - thetad + thetak)
+ dtRelong_s = dtRcirc_s
+ dtZelong_s = dtZcirc_s + (-1 + kappa)*rho*Cos(thAdj_s - thetad + thetak)
+ dtRelongTilt_s = dtRelong_s*Cos(thetad - thetak) - dtZelong_s*Sin(thetad - thetak)
+ dtZelongTilt_s = dtZelong_s*Cos(thetad - thetak) + dtRelong_s*Sin(thetad - thetak)
+ dtRtri_s = dtRelongTilt_s + rho*Sin(thAdj_s) &
+ - rho*(1 + delta*Cos(thAdj_s))*Sin(thAdj_s + delta*Sin(thAdj_s))
+ dtZtri_s = dtZelongTilt_s
+ dtRtriTilt_s = dtRtri_s*Cos(thetad) + dtZtri_s*Sin(thetad)
+ dtZtriTilt_s = dtZtri_s*Cos(thetad) - dtRtri_s*Sin(thetad)
+ R_theta_s = dtheta_dchi_s*dtRtriTilt_s
+ Z_theta_s = dtheta_dchi_s*dtZtriTilt_s
+ endif
+ if (sign_Ip_CW.lt.0.0) then
+ call lag3interp(nu1,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ nu1_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(Bp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ Bp_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(dlp,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ dlp_s = tmp_arr_s(np_s:1:-1)
+ call lag3interp(Rc,theta,np,tmp_arr_s,theta_s_reverse,np_s)
+ Rc_s = tmp_arr_s(np_s:1:-1)
+ else
+ call lag3interp(nu1,theta,np,nu1_s,theta_s,np_s)
+ call lag3interp(Bp,theta,np,Bp_s,theta_s,np_s)
+ call lag3interp(dlp,theta,np,dlp_s,theta_s,np_s)
+ call lag3interp(Rc,theta,np,Rc_s,theta_s,np_s)
+ endif
+
+ psi1_s = R_s*Bp_s*sign_Ip_CW
+
+ dBp_dchi_s=deriv_fd(Bp_s,chi_s,np_s)
+
+ Bphi_s=F/R_s
+ B_s=sqrt(Bphi_s**2+Bp_s**2)
+ cosu_s=Z_theta_s/dlp_s/dtheta_dchi_s
+ sinu_s=-R_theta_s/dlp_s/dtheta_dchi_s
+
+ !radial derivative of straight field line angle
+ dchidx_s=-(nu1_s/psi1_s*dxPsi+chi_s*dq_dx)/q0
+
+ !Bfield derivatives in Mercier-Luc coordinates (varrho,l,fz)
+ dB_drho_s=-1./B_s*(F**2/R_s**3*cosu_s+Bp_s**2/Rc_s+mu_0*psi1_s*pprime)
+ dB_dl_s=1./B_s*(Bp_s*dBp_dchi_s/dtheta_dchi_s/dlp_s+F**2/R_s**3*sinu_s)
+
+ dnu_drho_s=nu1_s
+ dnu_dl_s=-F/(R_s*psi1_s)
+ grad_nu_s=sqrt(dnu_drho_s**2+dnu_dl_s**2)
+
+ !contravariant metric coefficients (varrho,l,fz)->(x,y,z)
+ gxx=(psi1_s/dxPsi)**2
+ gxy=-psi1_s/dxPsi*C_y*sign_Ip_CW*nu1_s
+ gxz=-psi1_s/dxPsi*(nu1_s+psi1_s*dq_dpsi*chi_s)/q0
+ gyy=C_y**2*(grad_nu_s**2+1/R_s**2)
+ gyz=sign_Ip_CW*C_y/q0*(grad_nu_s**2+dq_dpsi*nu1_s*psi1_s*chi_s)
+ gzz=1./q0**2*(grad_nu_s**2+2.*dq_dpsi*nu1_s*psi1_s*chi_s+(dq_dpsi*psi1_s*chi_s)**2)
+
+ jacobian=1./sqrt(gxx*gyy*gzz + 2.*gxy*gyz*gxz - gxz**2*gyy - gyz**2*gxx - gzz*gxy**2)
+
+ !covariant metric coefficients
+ g_xx=jacobian**2*(gyy*gzz-gyz**2)
+ g_xy=jacobian**2*(gxz*gyz-gxy*gzz)
+ g_xz=jacobian**2*(gxy*gyz-gxz*gyy)
+ g_yy=jacobian**2*(gxx*gzz-gxz**2)
+ g_yz=jacobian**2*(gxz*gxy-gxx*gyz)
+ g_zz=jacobian**2*(gxx*gyy-gxy**2)
+
+ !Bfield derivatives
+ !dBdx = e_x * nabla B = J (nabla y x nabla z) * nabla B
+ dBdx=jacobian*C_y/(q0*R_s)*(F/(R_s*psi1_s)*dB_drho_s+(nu1_s+dq_dpsi*chi_s*psi1_s)*dB_dl_s)
+ dBdz=1./B_s*(Bp_s*dBp_dchi_s-F**2/R_s**3*R_theta_s)
+
+ !curvature terms (these are just local and will be recalculated in geometry.F90)
+ K_x = (0.-g_yz/g_zz*dBdz)
+ K_y = (dBdx-g_xz/g_zz*dBdz)
+
+ !(R,Z) derivatives for visualization
+ dxdR_s = dx_drho/drPsi*psi1_s*cosu_s
+ dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s
+
+ if (edge_opt==0.0) then
+ !gene z-grid
+ chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/Nz,Nz)
+ else
+ !new parallel coordinate chi_out==zprime
+ !see also tracer_aux.F90
+ if (Npol>1) STOP "ERROR: Npol>1 has not been implemented for edge_opt=\=0.0"
+ do k=izs,ize
+ chi_out(k)=sinh((-pi+k*2.*pi/Nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
+ enddo
+ !transform metrics according to chain rule
+ do k=1,np_s
+ !>dz'/dz conversion for edge_opt as function of z
+ if (edge_opt.gt.0) then
+ dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1))/&
+ sqrt((edge_opt*chi_s(k))**2+1)
+ else
+ dzprimedz = 1.0
+ endif
+ gxz(k)=gxz(k)*dzprimedz
+ gyz(k)=gyz(k)*dzprimedz
+ gzz(k)=gzz(k)*dzprimedz**2
+ jacobian(k)=jacobian(k)/dzprimedz
+ dBdz(k)=dBdz(k)/dzprimedz
+ enddo
+ endif !edge_opt
+
+ !interpolate down to GENE z-grid
+ call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,Nz)
+ call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,Nz)
+ call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,Nz)
+ call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,Nz)
+ call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,Nz)
+ call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,Nz)
+ call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,Nz)
+ call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,Nz)
+ call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,Nz)
+ call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,Nz)
+ call lag3interp(R_s,chi_s,np_s,R_out,chi_out,Nz)
+ call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,Nz)
+ call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,Nz)
+ call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz)
+ ! Fill the geom arrays with the results
+ do eo=0,1
+ gxx_(izs:ize,eo) =gxx_out(izs:ize)
+ gyy_(izs:ize,eo) =gyy_out(izs:ize)
+ gxz_(izs:ize,eo) =gxz_out(izs:ize)
+ gyz_(izs:ize,eo) =gyz_out(izs:ize)
+ dBdx_(izs:ize,eo) =dBdx_out(izs:ize)
+ dBdy_(izs:ize,eo) =0.
+ gxy_(izs:ize,eo) =gxy_out(izs:ize)
+ gzz_(izs:ize,eo) =gzz_out(izs:ize)
+ Bfield_(izs:ize,eo) =Bfield_out(izs:ize)
+ jacobian_(izs:ize,eo) =jacobian_out(izs:ize)
+ dBdz_(izs:ize,eo) =dBdz_out(izs:ize)
+ R_hat_(izs:ize,eo) =R_out(izs:ize)
+ Z_hat_(izs:ize,eo) =Z_out(izs:ize)
+ dxdR_(izs:ize,eo) = dxdR_out(izs:ize)
+ dxdZ_(izs:ize,eo) = dxdZ_out(izs:ize)
+
+ !! UPDATE GHOSTS VALUES (since the miller function in GENE does not)
+ CALL update_ghosts_z(gxx_(:,eo))
+ CALL update_ghosts_z(gyy_(:,eo))
+ CALL update_ghosts_z(gxz_(:,eo))
+ CALL update_ghosts_z(gxy_(:,eo))
+ CALL update_ghosts_z(gzz_(:,eo))
+ CALL update_ghosts_z(Bfield_(:,eo))
+ CALL update_ghosts_z(dBdx_(:,eo))
+ CALL update_ghosts_z(dBdy_(:,eo))
+ CALL update_ghosts_z(dBdz_(:,eo))
+ CALL update_ghosts_z(jacobian_(:,eo))
+ CALL update_ghosts_z(R_hat_(:,eo))
+ CALL update_ghosts_z(Z_hat_(:,eo))
+ CALL update_ghosts_z(dxdR_(:,eo))
+ CALL update_ghosts_z(dxdZ_(:,eo))
+
+ ! Curvature operator (Frei et al. 2022 eq 2.15)
+ DO iz = izgs,izge
+ G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
+ G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
+ G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
+ Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo)
+ Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo)
+
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff_(iz,eo) = 1._dp / jacobian_(iz,eo) / Bfield_(iz,eo)
+ ENDDO
+ ENDDO
+
+ contains
+
+
+ SUBROUTINE update_ghosts_z(fz_)
+ IMPLICIT NONE
+ ! INTEGER, INTENT(IN) :: nztot_
+ REAL(dp), DIMENSION(izgs:izge), INTENT(INOUT) :: fz_
+ REAL(dp), DIMENSION(-2:2) :: buff
+ INTEGER :: status(MPI_STATUS_SIZE), count
+
+ IF(Nz .GT. 1) THEN
+ IF (num_procs_z .GT. 1) THEN
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ count = 1 ! one point to exchange
+ !!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(fz_(ize), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last
+ buff(-1), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv(fz_(ize-1), count, MPI_DOUBLE, nbr_U, 0, & ! Send to Up the last
+ buff(-2), count, MPI_DOUBLE, nbr_D, 0, & ! Receive from Down the first-2
+ comm0, status, ierr)
+
+ !!!!!!!!!!! Send ghost to down neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(fz_(izs), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first
+ buff(+1), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+1
+ comm0, status, ierr)
+
+ CALL mpi_sendrecv(fz_(izs+1), count, MPI_DOUBLE, nbr_D, 0, & ! Send to Down the first
+ buff(+2), count, MPI_DOUBLE, nbr_U, 0, & ! Recieve from Up the last+2
+ comm0, status, ierr)
+ ELSE
+ buff(-1) = fz_(ize )
+ buff(-2) = fz_(ize-1)
+ buff(+1) = fz_(izs )
+ buff(+2) = fz_(izs+1)
+ ENDIF
+ fz_(ize+1) = buff(+1)
+ fz_(ize+2) = buff(+2)
+ fz_(izs-1) = buff(-1)
+ fz_(izs-2) = buff(-2)
+ ENDIF
+ END SUBROUTINE update_ghosts_z
+
+
+ !> Generate an equidistant array from min to max with n points
+ function linspace(min,max,n) result(out)
+ real(dp):: min, max
+ integer:: n
+ real(dp), dimension(n):: out
+
+ do i=1,n
+ out(i)=min+(i-1)*(max-min)/(n-1)
+ enddo
+ end function linspace
+
+ !> Weighted average
+ real(dp) function average(var,weight)
+ real(dp), dimension(np):: var, weight
+ average=sum(var*weight)/sum(weight)
+ end function average
+
+ !> full theta integral with weight function dlp
+ real(dp) function dlp_int(var,dlp)
+ real(dp), dimension(np):: var, dlp
+ dlp_int=sum(var*dlp)*2*pi*Npol_ext/np
+ end function dlp_int
+
+ !> theta integral with weight function dlp, up to index 'ind'
+ real(dp) function dlp_int_ind(var,dlp,ind)
+ real(dp), dimension(np):: var, dlp
+ integer:: ind
+
+ dlp_int_ind=0.
+ if (ind.gt.1) then
+ dlp_int_ind=dlp_int_ind+var(1)*dlp(1)*pi*Npol_ext/np
+ dlp_int_ind=dlp_int_ind+(sum(var(2:ind-1)*dlp(2:ind-1)))*2*pi*Npol_ext/np
+ dlp_int_ind=dlp_int_ind+var(ind)*dlp(ind)*pi*Npol_ext/np
+ endif
+ end function dlp_int_ind
+
+ !> 1st derivative with 2nd order finite differences
+ function deriv_fd(y,x,n) result(out)
+ integer, intent(in) :: n
+ real(dp), dimension(n):: x,y,out,dx
+
+ !call lag3deriv(y,x,n,out,x,n)
+
+ out=0.
+ do i=2,n-1
+ out(i)=out(i)-y(i-1)/2
+ out(i)=out(i)+y(i+1)/2
+ enddo
+ out(1)=y(2)-y(1)
+ out(n)=y(n)-y(n-1)
+ dx=x(2)-x(1)
+ out=out/dx
+
+ end function deriv_fd
+
+
+ end subroutine get_miller
+
+
+END MODULE miller
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 75a2768787ab2ba5112ba0d4849644cbc6f673b1..5d28283677c9cf1b2b4614330a0a974bf29bc2eb 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -52,6 +52,7 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum
REAL(dp), PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i
REAL(dp), PUBLIC, PROTECTED :: sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i
+ REAL(dp), PUBLIC, PROTECTED :: dpdx = 0 ! radial pressure gradient
PUBLIC :: model_readinputs, model_outputinputs
CONTAINS
@@ -75,8 +76,9 @@ CONTAINS
ENDIF
IF(beta .GT. 0) THEN
- EM = .TRUE.
IF(my_id.EQ.0) print*, 'Electromagnetic effects are included'
+ EM = .TRUE.
+ dpdx = -(tau_i*(k_Ni + k_Ti) + tau_e*(k_Ne + k_Te))
ENDIF
taue_qe = tau_e/q_e ! factor of the magnetic moment coupling
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 8123fe782f59671450d45a94b0d9e7e5fd21711c..948d66c054936fb005eb96f3d6cc243eba8baba4 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -33,7 +33,7 @@ SUBROUTINE moments_eq_rhs_e
USE model
USE prec_const
USE collision
- USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
+ USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
@@ -121,7 +121,9 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all RHS terms
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp&
+ -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp&
+ ! Perpendicular pressure effects
+ -i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau Damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
@@ -179,7 +181,7 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
- USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
+ USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
@@ -269,7 +271,9 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of all RHS terms
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
+ -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp &
+ ! Perpendicular pressure effects
+ -i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 8e9a2ce5cb3d0374acfbf7846a192cc37f56011b..615471762690d6bc3b5000abb29dbed5f45eb6fa 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -141,7 +141,7 @@ ENDIF
eo = MODULO(parray_i(ip),2)
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
j_int=jarray_i(ij)
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute for every moments except for closure 1
! Set non linear sum truncation
IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxi
@@ -152,7 +152,7 @@ ENDIF
ENDIF
bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
-!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
+!-----------!! ELECTROSTATIC CONTRIBUTION
! First convolution terms
F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
! Second convolution terms
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 54f99c63ae489810c974c4bf6a0003ebc99634be..2c34792828e2c1564c921dee25d05a124c74b2a9 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -286,8 +286,8 @@ SUBROUTINE compute_radial_ion_heatflux
ENDDO
ENDIF
! Add polarisation contribution
- integrant(izs:ize) = integrant(izs:ize) + tau_i*imagu*ky_&
- *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
+ ! integrant(izgs:izge) = integrant(izgs:izge) + tau_i*imagu*ky_&
+ ! *CONJG(phi(iky,ikx,izgs:izge))*phi(iky,ikx,izgs:izge) * HF_phi_correction_operator(iky,ikx,izgs:izge)
! Integrate over z
integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
call simpson_rule_z(integrant,integral)
@@ -370,8 +370,8 @@ SUBROUTINE compute_radial_electron_heatflux
ENDDO
ENDIF
! Add polarisation contribution
- integrant(izs:ize) = integrant(izs:ize) + tau_e*imagu*ky_&
- *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
+ ! integrant(izs:ize) = integrant(izs:ize) + tau_e*imagu*ky_&
+ ! *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
! Integrate over z
integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
call simpson_rule_z(integrant,integral)
diff --git a/testcases/cyclone_example/cyclone_example.txt b/testcases/cyclone_example/cyclone_example.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b38bfef642ef3e498137ca790fa392b2d085fc33
--- /dev/null
+++ b/testcases/cyclone_example/cyclone_example.txt
@@ -0,0 +1,4 @@
+This is a testcase reproducing the cyclone base case of Dimits 2000
+Adiabatic electrons, s-alpha geometry, gradlnN = 2.22, gradlnT = 6.96
+With a small P,J=4,2 polynomial basis, one should observe the secondary instability (KHI) at t~50R/cs.
+The saturated heat flux should be located around Qx ~ 30, i.e. Qx/QGB ~ 2 which is close to Dimits results.
diff --git a/testcases/cyclone_example/fort.90 b/testcases/cyclone_example/fort.90
new file mode 100644
index 0000000000000000000000000000000000000000..d562465a47ad3fcff4858ddad9ea05cd9b707fec
Binary files /dev/null and b/testcases/cyclone_example/fort.90 differ
diff --git a/testcases/fort.90 b/testcases/fort.90
new file mode 100644
index 0000000000000000000000000000000000000000..abce02cfd7ef427b8a376899d1c28642821ae22c
--- /dev/null
+++ b/testcases/fort.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 100
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 128
+ Ly = 240
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 1000
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 1.0
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/kobayashi_2015_linear_results.m b/testcases/kobayashi_2015_linear_results.m
deleted file mode 100644
index 82f339125f78966b0d059746be9b56d8eaa292d0..0000000000000000000000000000000000000000
--- a/testcases/kobayashi_2015_linear_results.m
+++ /dev/null
@@ -1,19 +0,0 @@
-% results for blue triangles in Kobayashi 2015
-% x_ = [0.06 0.4 0.6];
-% y_ = [0.0 0.1 0.0];
-% plot(x_,y_,'^b','DisplayName','Kobayashi et al. 2015');
-
-
-% results for red squares in Kobayashi 2015
-% x_ = [0.065 0.15 0.35];
-% y_ = [0.0 0.15 0.0];
-% plot(x_,y_,'sr','DisplayName','Kobayashi et al. 2015');
-
-
-% results for green triangles in Kobayashi 2015
-x_ = [0.06 0.5 1.8];
-y_ = [0.0 0.2 0.0];
-plot(x_,y_,'dg','DisplayName','Kobayashi et al. 2015');
-
-
-
diff --git a/testcases/Hallenbert.m b/testcases/matlab_testscripts/Hallenbert.m
similarity index 100%
rename from testcases/Hallenbert.m
rename to testcases/matlab_testscripts/Hallenbert.m
diff --git a/testcases/benchmark_HeLaZ_gene_transport_zpinch.m b/testcases/matlab_testscripts/benchmark_HeLaZ_gene_transport_zpinch.m
similarity index 100%
rename from testcases/benchmark_HeLaZ_gene_transport_zpinch.m
rename to testcases/matlab_testscripts/benchmark_HeLaZ_gene_transport_zpinch.m
diff --git a/testcases/benchmark_HeLaZ_molix.m b/testcases/matlab_testscripts/benchmark_HeLaZ_molix.m
similarity index 100%
rename from testcases/benchmark_HeLaZ_molix.m
rename to testcases/matlab_testscripts/benchmark_HeLaZ_molix.m
diff --git a/testcases/benchmark_linear_1D_entropy_mode.m b/testcases/matlab_testscripts/benchmark_linear_1D_entropy_mode.m
similarity index 100%
rename from testcases/benchmark_linear_1D_entropy_mode.m
rename to testcases/matlab_testscripts/benchmark_linear_1D_entropy_mode.m
diff --git a/testcases/cyclone_test_case.m b/testcases/matlab_testscripts/cyclone_test_case.m
similarity index 100%
rename from testcases/cyclone_test_case.m
rename to testcases/matlab_testscripts/cyclone_test_case.m
diff --git a/testcases/linear_1D_entropy_mode.m b/testcases/matlab_testscripts/linear_1D_entropy_mode.m
similarity index 100%
rename from testcases/linear_1D_entropy_mode.m
rename to testcases/matlab_testscripts/linear_1D_entropy_mode.m
diff --git a/testcases/linear_damping.m b/testcases/matlab_testscripts/linear_damping.m
similarity index 100%
rename from testcases/linear_damping.m
rename to testcases/matlab_testscripts/linear_damping.m
diff --git a/testcases/miller_example_w_salpha/fort_00.90 b/testcases/miller_example_w_salpha/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..2701de9ddf236b5b9e7da968c91d1958455dd2bd
--- /dev/null
+++ b/testcases/miller_example_w_salpha/fort_00.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 100
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 's-alpha'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/miller_example_w_salpha/fort_01.90 b/testcases/miller_example_w_salpha/fort_01.90
new file mode 100644
index 0000000000000000000000000000000000000000..7deb16c3b42c4b80864dd5eaf78efb4d91892f63
--- /dev/null
+++ b/testcases/miller_example_w_salpha/fort_01.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 400
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = 0
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/fort.90 b/testcases/zpinch_example/fort.90
similarity index 76%
rename from fort.90
rename to testcases/zpinch_example/fort.90
index ae9af96cf7e62ca4cdffaf1472989319302529d2..f769f6b8fb6a2e041fc5dde3ce1edaf7c39810f0 100644
--- a/fort.90
+++ b/testcases/zpinch_example/fort.90
@@ -1,7 +1,7 @@
&BASIC
nrun = 100000000
- dt = 0.01
- tmax = 400
+ dt = 0.005
+ tmax = 100
maxruntime = 356400
/
&GRID
@@ -10,23 +10,23 @@
pmaxi = 4
jmaxi = 2
Nx = 128
- Lx = 120
- Ny = 40
+ Lx = 200
+ Ny = 48
Ly = 60
Nz = 1
SG = .f.
/
&GEOMETRY
geom = 'Z-pinch'
- q0 = 1.4
+ q0 = 0
shear = 0
- eps = 0.18
+ eps = 0
/
&OUTPUT_PAR
nsave_0d = 10
nsave_1d = -1
nsave_2d = -1
- nsave_3d = 20
+ nsave_3d = 100
nsave_5d = 1000
write_doubleprecision = .f.
write_gamma = .t.
@@ -42,9 +42,9 @@
&MODEL_PAR
! Collisionality
CLOS = 0
- NL_CLOS = 0
+ NL_CLOS = -1
LINEARITY = 'nonlinear'
- KIN_E = .f.
+ KIN_E = .t.
mu_x = 1.0
mu_y = 1.0
N_HD = 4
@@ -58,22 +58,23 @@
sigma_i = 1
q_e = -1
q_i = 1
- K_n = 2.22
- K_T = 6.96
- K_E = 0
+ K_Ne = 2.22
+ K_Te = 6.96
+ K_Ni = 2.22
+ K_Ti = 6.96
GradB = 1
CurvB = 1
lambdaD = 0
+ beta = 1e-4
/
&COLLISION_PAR
collision_model = 'SG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
gyrokin_CO = .false.
interspecies = .true.
- mat_file = 'iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
-! mat_file = 'iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
- INIT_OPT = 'ppj'
+ INIT_OPT = 'phi'
ACT_ON_MODES = 'donothing'
init_background = 0
init_noiselvl = 0.00005
diff --git a/wk/Zpinch_coll_scan_kN_1.7.m b/wk/Zpinch_coll_scan_kN_1.7.m
index e223e0e4506c7cdc3279cba5d13af17832cbc0ce..ad00ba144127771160fe1857e219080b70a12441 100644
--- a/wk/Zpinch_coll_scan_kN_1.7.m
+++ b/wk/Zpinch_coll_scan_kN_1.7.m
@@ -1,5 +1,5 @@
-%%
if 0
+%%
figure
Kn = 1.7;
@@ -11,27 +11,42 @@ Kn = 1.7;
% errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama DK (4,2)'); hold on
-% SUGAMA GK 4,2
+% SUGAMA GK 4,2 200x32
nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
-Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.09e-1];
+Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.79e-1];
Gstd_a = [3.04e-2 1.42e-2 1.56e-2 1.23e-2 1.20e-2 1.57e-2 1.63e-2 2.06e-2 2.14e-02 1.78e-2];
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2)'); hold on
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 200x32'); hold on
+
+% SUGAMA GK 4,2 256x64
+nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
+Gavg_a = [1.99e-2 3.03e-2 4.44e-2 5.87e-2 7.87e-2 1.07e-1 1.22e-1 1.38e-1 1.63e-1 1.75e-1];
+Gstd_a = [1.55e-2 2.22e-2 1.12e-2 1.54e-2 2.20e-2 2.76e-2 2.44e-2 3.20e-2 3.12e-2 3.13e-2];
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 256x64'); hold on
-% FCGK 4,2
+
+% SUGAMA GK 6,3 200x32
+nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
+Gavg_a = [2.35e-2 3.57e-2 5.09e-2 6.97e-2 8.65e-2 1.01e-1 1.16e-1 1.33e-1 1.49e-1 1.62e-1];
+Gstd_a = [3.76e-2 3.91e-2 2.96e-2 1.57e-2 1.73e-2 1.94e-2 2.21e-2 2.01e-2 1.93e-2 2.24e-2];
+
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (6,3) 200x32'); hold on
+
+
+% FCGK 4,2 200x32
nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 10.0];
Gavg_a = [8.57e-2 1.45e-1 2.25e-1 2.87e-1 3.48e-1 4.06e-1 4.51e-1 3.65e-1*Kn];
Gstd_a = [2.07e-2 2.61e-2 2.40e-2 3.46e-2 4.30e-2 5.00e-2 5.11e-2 0];
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2)'); hold on
-
-% LDGK ii 6,3
-nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00];
-Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2];
-Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2];
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2) 200x32'); hold on
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3)'); hold on
+% % LDGK ii 6,3 200x32
+% nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00];
+% Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2];
+% Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2];
+%
+% errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3) 200x32'); hold on
% Collisionless
plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
@@ -41,7 +56,7 @@ plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
xlim([0 0.1]);
legend('show');
xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
-
+set(gca,'Yscale','log');
end
if 0
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index f45fc16d1330d14370a68a1e684873034f8b5ed1..53f11c4f9f6cb9dffb0aa4291a7d347b1c5b0f66 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -12,9 +12,9 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
-% folder = '/misc/gene_results/HP_fig_2a_mu_1e-2/';
-% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
-% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
+folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_5e-2/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
@@ -28,7 +28,8 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
-folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
+% folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
+% folder = '/misc/gene_results/miller/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -75,9 +76,9 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = 'v_y';
-% options.NAME = 'n_i^{NZ}';
+options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xy';
@@ -130,15 +131,15 @@ end
if 0
%% Time averaged spectrum
-options.TIME = 300:600;
+options.TIME = [4000 8000];
options.NORM =1;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = 'n_i';
-% options.NAME ='\Gamma_x';
+options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
-options.COMPXY = 'zero';
+options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(gene_data,options);
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index eecf7147d006ff92ac3ec5b7974eaa85a4bbb2a9..2a875f63bc57ccab9a3aacc780124248ef6bcd97 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -28,7 +28,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-i_ = 11;
+i_ = 1;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
@@ -36,37 +36,43 @@ options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times dur
options.NCUT = 4; % Number of cuts for averaging and error estimation
options.NMVA = 100; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-% options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-% options.INTERP = 1;
+options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.INTERP = 0;
+options.RESOLUTION = 256;
fig = plot_radial_transport_and_spacetime(data,options);
save_figure(data,fig,'.png')
end
if 0
%% statistical transport averaging
-options.T = [200 400];
+for i_ = 1:2:21
+% i_ = 3;
+disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
+disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
+options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)];
+options.NPLOTS = 0;
fig = statistical_transport_averaging(data,options);
end
-
+end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
-% options.NAME = 'n_i';
+options.NAME = 'n_i';
options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
-% options.COMP = 'avg';
+options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [000:50:7000];
+options.TIME = [000:0.1:7000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -146,11 +152,11 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [300 600];
+options.TIME = [2000 3000];
options.NORM =1;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
-% options.NAME ='\Gamma_x';
+options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
@@ -200,7 +206,9 @@ end
if 0
%% Metric infos
-fig = plot_metric(data);
+options.SHOW_FLUXSURF = 1;
+options.SHOW_METRICS = 0;
+fig = plot_metric(data,options);
end
if 0
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index b180dd526e307d5a0cb5a7105dd72c0f2cb0165f..cdd4074d3972b6cd85e0275ae582dea29129c9dc 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -191,10 +191,34 @@ resdir ='';
% resdir ='Zpinch_rerun/kN_1.7_SGGK_conv_200x32x7x3_nu_0.01';
% resdir ='Zpinch_rerun/kN_1.7_LDGKii_200x32x7x3_nu_scan';
% resdir ='Zpinch_rerun/nu_0.1_LDGKii_200x48x7x4_kN_scan';
-resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
+% resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
+% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan';
% resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan';
-% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
-%%
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+% % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
+% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
+% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
+%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_mu_0.01';
+
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
+resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
+
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_nu_0.1';
+
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_nu_0.1';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_nu_0.1';
+
+
+JOBNUMMIN = 00; JOBNUMMAX = 03;
resdir = ['results/',resdir];
run analysis_gyacomo
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 2a56c2dff5fab0b74549baddaeedeaa651fa11cd..63cb0bf33d953ee2f8636517e5faaf903fc35339 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -1,64 +1,66 @@
% Directory of the code "mypathtoHeLaZ/HeLaZ/"
-helazdir = '/home/ahoffman/HeLaZ/';
+gyacomodir = '/home/ahoffman/gyacomo/';
% Directory of the simulation (from results)
% if 1% Local results
-% outfile ='volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir ='volcokas/64x32x16x5x3_kin_e_npol_1';
%% Dimits
-% outfile ='shearless_cyclone/128x64x16x5x3_Dim_90';
-% outfile ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60';
-% outfile ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70';
-% outfile ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70';
+% resdir ='shearless_cyclone/128x64x16x5x3_Dim_90';
+% resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60';
+% resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70';
+% resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70';
%% AVS
-% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1';
-% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1';
-% outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine';
-% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
-% outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
-% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine';
+% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
%% shearless CBC
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_120';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_080';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_120';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_080';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_100';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_120';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_080';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_100';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_120';
-% outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
+% resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
%% CBC
-% outfile = 'CBC/64x32x16x5x3';
-% outfile = 'CBC/64x128x16x5x3';
-% outfile = 'CBC/128x64x16x5x3';
-% outfile = 'CBC/96x96x16x3x2_Nexc_6';
-% outfile = 'CBC/128x96x16x3x2';
-% outfile = 'CBC/192x96x16x3x2';
-% outfile = 'CBC/192x96x24x13x7';
-% outfile = 'CBC/kT_11_128x64x16x5x3';
-% outfile = 'CBC/kT_9_256x128x16x3x2';
-% outfile = 'CBC/kT_4.5_128x64x16x13x3';
-% outfile = 'CBC/kT_4.5_192x96x24x13x7';
-% outfile = 'CBC/kT_4.5_128x64x16x13x7';
-% outfile = 'CBC/kT_4.5_128x96x24x15x5';
-% outfile = 'CBC/kT_5.3_192x96x24x13x7';
-% outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
+% resdir = 'CBC/64x32x16x5x3';
+% resdir = 'CBC/64x128x16x5x3';
+% resdir = 'CBC/128x64x16x5x3';
+% resdir = 'CBC/96x96x16x3x2_Nexc_6';
+% resdir = 'CBC/128x96x16x3x2';
+% resdir = 'CBC/192x96x16x3x2';
+% resdir = 'CBC/192x96x24x13x7';
+% resdir = 'CBC/kT_11_128x64x16x5x3';
+% resdir = 'CBC/kT_9_256x128x16x3x2';
+% resdir = 'CBC/kT_4.5_128x64x16x13x3';
+% resdir = 'CBC/kT_4.5_192x96x24x13x7';
+% resdir = 'CBC/kT_4.5_128x64x16x13x7';
+% resdir = 'CBC/kT_4.5_128x96x24x15x5';
+% resdir = 'CBC/kT_5.3_192x96x24x13x7';
+% resdir = 'CBC/kT_13_large_box_128x64x16x5x3';
+% resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
-% outfile = 'CBC/kT_scan_128x64x16x5x3';
-% outfile = 'CBC/kT_scan_192x96x16x3x2';
-% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
-% outfile = 'dbg/nexc_dbg';
-outfile = 'CBC/NM_F4_kT_4.5_192x64x24x6x4';
+% resdir = 'CBC/kT_scan_128x64x16x5x3';
+% resdir = 'CBC/kT_scan_192x96x16x3x2';
+% resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% resdir = 'dbg/nexc_dbg';
+% resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4';
-% outfile = 'CBC_Ke_EM/192x96x24x5x3';
-% outfile = 'CBC_Ke_EM/96x48x16x5x3';
-% outfile = 'CBC_Ke_EM/minimal_res';
+% resdir = 'CBC_Ke_EM/192x96x24x5x3';
+% resdir = 'CBC_Ke_EM/96x48x16x5x3';
+% resdir = 'CBC_Ke_EM/minimal_res';
%% KBM
-% outfile = 'NL_KBM/192x64x24x5x3';
+% resdir = 'NL_KBM/192x64x24x5x3';
%% Linear CBC
-% outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
-
+% resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
+% resdir = 'testcases/miller_example';
+resdir = 'Miller/128x256x3x2_CBC_dt_5e-3';
+% resdir = ['results/',resdir];
JOBNUMMIN = 00; JOBNUMMAX = 10;
-run analysis_HeLaZ
+run analysis_gyacomo
diff --git a/wk/lin_EPY.m b/wk/lin_EPY.m
index 41abd8aa00081b725df41e4b2ed803da719774ea..2ddd756a47e22d4e3b808e4f0c8dd019d6da95fd 100644
--- a/wk/lin_EPY.m
+++ b/wk/lin_EPY.m
@@ -4,11 +4,12 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_EPY'; % Name of the simulation
-RUN = 1; % To run or just to load
+% SIMID = 'dbg'; % Name of the simulation
+RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
PROGDIR = '/home/ahoffman/gyacomo/';
-% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
@@ -17,7 +18,7 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 2.2; % ele Density '''
+K_Ne = 1.6; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
K_Ni = K_Ne; % ion Density gradient drive
K_Ti = K_Ni/4; % ion Temperature '''
@@ -25,7 +26,7 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
+P = 2;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
@@ -45,11 +46,11 @@ SHEAR = 0.8; % magnetic shear
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 50; % Maximal time unit
+TMAX = 500; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 2; % Sampling per time unit for 2D arrays
+SPS3D = 1/5; % Sampling per time unit for 2D arrays
SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
@@ -57,7 +58,7 @@ JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'DG';
+CO = 'SG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -94,7 +95,8 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
@@ -108,9 +110,10 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
%% Short analysis
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+lg = compute_fluxtube_growth_rate(data,options);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
@@ -154,13 +157,13 @@ options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
-if 0
+if 1
%% Mode evolution
options.NORMALIZED = 0;
-options.K2PLOT = 1;
+options.K2PLOT = 10;
options.TIME = [0:1000];
options.NMA = 1;
-options.NMODES = 1;
+options.NMODES = 10;
options.iz = 'avg';
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index bdf3ec0590faa766d52ac04e5038abfb85477b4c..3ebe622d264ab21cb5e9b2f822abf10c897aed18 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -4,18 +4,19 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'lin_ITG'; % Name of the simulation
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-% EXECNAME = 'helaz3';
-EXECNAME = 'helaz3_dbg';
+HELAZDIR = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_dbg';
+EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.0; % Collision frequency
+NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 2.22; % ele Density '''
K_Te = 6.96; % ele Temperature '''
@@ -32,24 +33,26 @@ PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 11; % real space x-gridpoints
+NX = 20; % real space x-gridpoints
NY = 2; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 2*pi/0.3; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
+NZ = 32; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
-% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
+KAPPA = 0.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 25; % Maximal time unit
-DT = 5e-3; % Time step
+TMAX = 1; % Maximal time unit
+DT = 3e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 2D arrays
@@ -127,7 +130,7 @@ end
if 1
%% Ballooning plot
options.time_2_plot = [120];
-options.kymodes = [0.9];
+options.kymodes = [0.3];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
diff --git a/wk/save_iFFT.m b/wk/save_iFFT.m
index d429b62274e33200d342a9d95540c2ef490ec0c6..d9b3dc9c1912c1aee5918ec63cab20506d3961e9 100644
--- a/wk/save_iFFT.m
+++ b/wk/save_iFFT.m
@@ -3,7 +3,7 @@ simdir = '/misc/gyacomo_outputs/results/Zpinch_rerun';
resolu = 'UHD_512x256x2x1';
output = 'outputs_01.h5';
filename = [simdir,'/',resolu,'/',output];
-fieldname = 'Ni00';
+fieldname = 'vEy';
%% OUTPUT
outdir = '.';