From e04f06e1e9b9aca0a470508e07185b263c60942a Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 30 Sep 2020 08:28:29 +0200
Subject: [PATCH] Backup before putting non linear solver
---
Makefile | 117 ++++--
local/dirs.inc | 4 +-
local/make.inc | 19 +-
matlab/AssociatedLaguerrePoly.m | 17 +
matlab/MOLI_time_solver.m | 4 +-
matlab/compute_Sapj.m | 37 ++
matlab/conv_thm_2D.m | 14 +
matlab/create_gif.m | 37 +-
matlab/dnjs.m | 6 +-
matlab/kernel.m | 15 +
matlab/load_2D_data.m | 13 +
matlab/load_5D_data.m | 13 +
matlab/save_figure.m | 2 +-
matlab/write_fort90.m | 54 +--
src/array_mod.F90 | 11 +-
src/auxval.F90 | 2 +-
src/basic_mod.F90 | 38 +-
src/coeff_mod.F90 | 95 +++++
src/coeff_mod.f90 | 633 --------------------------------
src/compute_Sapj.F90 | 109 ++++++
src/control.F90 | 23 +-
src/convolution_mod.F90 | 221 +++++++++++
src/diagnose.F90 | 155 ++++----
src/diagnostics_par_mod.F90 | 21 +-
src/fields_mod.F90 | 3 +-
src/fourier_grid_mod.F90 | 22 +-
src/inital.F90 | 177 ++++++---
src/initial_par_mod.F90 | 30 +-
src/memory.F90 | 10 +-
src/model_mod.F90 | 8 +-
src/moments_eq_rhs.F90 | 41 ++-
src/poisson.F90 | 7 +-
src/prec_const_mod.F90 | 2 +-
src/readinputs.F90 | 2 +-
src/stepon.F90 | 31 +-
src/time_integration_mod.F90 | 2 +-
wk/Test_Non_Lin.m | 135 +++----
wk/benchmark_kperp_scan.m | 2 +
wk/benchmark_time_solver.m | 394 ++++++++++----------
wk/compute_Sapj.m | 41 ---
wk/fort.90 | 60 +--
wk/kr_kz_time_evolution.m | 92 ++---
wk/kz_linear.m | 1 +
wk/launcher.m | 8 +-
44 files changed, 1341 insertions(+), 1387 deletions(-)
create mode 100644 matlab/AssociatedLaguerrePoly.m
create mode 100644 matlab/compute_Sapj.m
create mode 100644 matlab/conv_thm_2D.m
create mode 100644 matlab/kernel.m
create mode 100644 matlab/load_2D_data.m
create mode 100644 matlab/load_5D_data.m
create mode 100644 src/coeff_mod.F90
delete mode 100644 src/coeff_mod.f90
create mode 100644 src/compute_Sapj.F90
create mode 100644 src/convolution_mod.F90
delete mode 100644 wk/compute_Sapj.m
diff --git a/Makefile b/Makefile
index a0820e0a..d104f738 100644
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,18 @@
include local/dirs.inc
include local/make.inc
+MKLROOT = /usr/local/intel/composerxe/composer_xe_2015/mkl
EXEC = $(BINDIR)/helaz
F90 = mpif90
+# Add Multiple-Precision Library
+EXTLIBS += -L$(FMDIR)/lib
+EXTINC += -I$(FMDIR)/mod
+
+EXTLIBS += -L$(FFTWDIR)/lib64
+EXTINC += -I$(FFTWDIR)/include
all: dirs src/srcinfo.h $(EXEC)
@@ -16,8 +23,9 @@ dirs:
mkdir -p $(BINDIR)
mkdir -p $(OBJDIR)
mkdir -p $(MODDIR)
+ mkdir -p $(BACKUPDIR)
-src/srcinfo.h:
+src/srcinfo.h:
( cd src/srcinfo; $(MAKE);)
clean: cleanobj cleanmod
@@ -37,65 +45,100 @@ cleanbin:
$(OBJDIR)/diagnose.o : src/srcinfo.h
FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
-$(OBJDIR)/control.o $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \
-$(OBJDIR)/fields_mod.o $(OBJDIR)/inital.o \
-$(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
-$(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/ppexit.o \
-$(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/fourier_grid_mod.o \
-$(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
-$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
-
-$(EXEC): $(FOBJ)
+$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/convolution_mod.o \
+$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
+$(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
+$(OBJDIR)/model_mod.o $(OBJDIR)/mkl_dfti.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \
+$(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \
+$(OBJDIR)/fourier_grid_mod.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 $@
-$(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_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)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/fourier_grid_mod.o
+
+ $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/fourier_grid_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)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.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)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/convolution_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/compute_Sapj.F90 -o $@
+
+ $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.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)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+
+ $(OBJDIR)/convolution_mod.o : src/convolution_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/mkl_dfti.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/convolution_mod.F90 -o $@
+
+ $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_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)/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)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
+
+ $(OBJDIR)/fourier_grid_mod.o : src/fourier_grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_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)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_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)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.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)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/fourier_grid_mod.o
+
+ $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/fourier_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)/mkl_dfti.o : $(MKLROOT)/include/mkl_dfti.f90
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) $(MKLROOT)/include/mkl_dfti.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.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
+
+ $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
-$(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+
+ $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_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/poisson.F90 -o $@
-$(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
+
+ $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_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)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+
+ $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_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
+
+ $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_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)/fourier_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)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
-$(OBJDIR)/fourier_grid_mod.o : src/fourier_grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_grid_mod.F90 -o $@
-$(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+
+ $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_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)/fourier_grid_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
+ $(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
diff --git a/local/dirs.inc b/local/dirs.inc
index 0667deb9..cafb3079 100644
--- a/local/dirs.inc
+++ b/local/dirs.inc
@@ -5,4 +5,6 @@ BINDIR = $(PREFIX)/bin
OBJDIR = $(PREFIX)/obj
LIBDIR = $(PREFIX)/lib
MODDIR = $(PREFIX)/mod
-
+FMDIR = $(PREFIX)/FM
+FFTWDIR = /usr/local/fftw-3.3.5
+BACKUPDIR = $(PREFIX)/backup
diff --git a/local/make.inc b/local/make.inc
index e5a03bdc..5a62151b 100644
--- a/local/make.inc
+++ b/local/make.inc
@@ -43,7 +43,7 @@ ifeq ($(OPTLEVEL), O) #optimized
F90FLAGS += -O3 -xHOST
endif
ifeq ($(COMPTYPE), g) #gnu
- F90FLAGS += -ffree-line-length-0 -O3
+ F90FLAGS += -ffree-line-length-0 -O3
endif
ifeq ($(COMPTYPE), c) #cray
F90FLAGS +=
@@ -56,7 +56,7 @@ ifeq ($(OPTLEVEL), g) #debug
F90FLAGS += -g -traceback -CB
endif
ifeq ($(COMPTYPE), g) #gnu
- F90FLAGS += -ffree-line-length-0 -g -fbacktrace -fcheck=all -pedantic -Wall
+ F90FLAGS += -ffree-line-length-0 -g -fbacktrace -fcheck=all -pedantic -Wall
endif
ifeq ($(COMPTYPE), c) #cray
F90FLAGS += -g -O0
@@ -86,13 +86,16 @@ IDIRS := -I$(SPC_LOCAL)/include/$(OPTLEVEL)
#LIBS := -lbsplines -lfutils -lpppack -lpputils2 \
# -lhdf5_fortran -lhdf5 -lz -ldl -lpthread
-LIBS := -lfutils -lhdf5_fortran -lhdf5 -lz -ldl -lpthread
+LIBS := -lfutils -lhdf5_fortran -lhdf5 -lz -ldl -lpthread
LDIRS := -L$(SPC_LOCAL)/lib/$(OPTLEVEL) -L$(HDF5)/lib
+# Add Multiple-Precision Library
+LIBS += -lfm
+
#
# setup of MKL library
# working well with intel compiler and intelmpi
-# MKL_LIB has to be set in .bashrc
+# MKL_LIB has to be set in .bashrc
#
ifdef MKL_LIB
@@ -121,7 +124,7 @@ ifdef MKL_LIB
endif
endif
- LDIRS += -L$(MKL_LIB)
+ LDIRS += -L$(MKL_LIB)
IDIRS += -I$(MKL_LIB)/../../include
endif
@@ -134,11 +137,11 @@ ifeq ($(USEMUMPS),1)
IDIRS += -I$(MUMPS)/include
LDIRS += -L$(MUMPS)/lib
ifdef PARMETIS
- LIBS += -lparmetis -lmetis
+ LIBS += -lparmetis -lmetis
LDIRS += -L$(PARMETIS)/lib -L$(METIS)/lib
endif
ifdef SCOTCH
- LIBS += -lesmumps -lscotch -lscotcherr
+ LIBS += -lesmumps -lscotch -lscotcherr
LDIRS += -L$(SCOTCH)/lib
endif
endif
@@ -194,5 +197,5 @@ endif
#Flag for finding external libraries in LDIR
EXTLIBS = $(LDIRS) -Wl,--start-group $(LIBS) -Wl,--end-group
-#Flag for finding external include files
+#Flag for finding external include files
EXTINC = $(IDIRS)
diff --git a/matlab/AssociatedLaguerrePoly.m b/matlab/AssociatedLaguerrePoly.m
new file mode 100644
index 00000000..7a5b36ee
--- /dev/null
+++ b/matlab/AssociatedLaguerrePoly.m
@@ -0,0 +1,17 @@
+% AssociatedLaguerrePoly.m by David Terr, Raytheon, 5-11-04
+% Given nonnegative integers n and k with k<=n, compute the associated
+% Laguerre polynomial L{n,k}. Return the result as a vector whose mth
+% element is the coefficient of x^(n+1-m).
+% polyval(AssociatedLaguerrePoly(n,k),x) evaluates L{n,k}(x).
+% Note: This program requires downloading binomial.m first.
+function Lnk = AssociatedLaguerrePoly(n,k)
+if k==0
+ Lnk = LaguerrePoly(n);
+else
+ Lnk = zeros(n+1,1);
+
+ for m=0:n
+ Lnk(n+1-m) = (-1)^m * binomial(k+n,n-m) / factorial(m);
+ end
+
+end
\ No newline at end of file
diff --git a/matlab/MOLI_time_solver.m b/matlab/MOLI_time_solver.m
index b2f25e39..d9b4d806 100644
--- a/matlab/MOLI_time_solver.m
+++ b/matlab/MOLI_time_solver.m
@@ -92,8 +92,8 @@ params.nu = MODEL.nu; % electron/ion collision frequency ... on
params.nuoveromegapemax = inf; % Maximum ratio between electron/ion collision frequency and electron plasma frequency [See Banks et al. (2017)]. Set to inf if not desired !!!
params.mu = MODEL.sigma_e; % sqrt(m_e/m_i)
params.kpar = 0.0; % normalized parallel wave number to the major radius
-params.kperp = GRID.kzmin; % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
-params.kr = GRID.krmin; % Radial component of perpendicular vector
+params.kperp = kz_MOLI; % normalized perpendicular toroidal wave number to the soundLarmor radius. Note: If ions ==0 (e.g. EPW), kperp --> b
+params.kr = kr_MOLI; % Radial component of perpendicular vector
params.alphaD = 0.0; % (k*Debye length)^2
params.Rn = MODEL.eta_n; % Major Radius / Background density gradient length
params.RTe = MODEL.eta_T; % Major Radius * normalized kperp / Background electron temperature gradient length
diff --git a/matlab/compute_Sapj.m b/matlab/compute_Sapj.m
new file mode 100644
index 00000000..0a708aab
--- /dev/null
+++ b/matlab/compute_Sapj.m
@@ -0,0 +1,37 @@
+function [ Sapj ] = compute_Sapj(p, j, Kr, Kz, Napj, specie, phi, MODEL, GRID)
+%COMPUT_SAPJ compute the non linear term for moment pj specie a
+% perform a convolution by product of inverse fourier transform and then
+% put the results back into k space with an FFT
+ Jmax = GRID.jmaxe * (specie=='e')...
+ + GRID.jmaxi * (specie=='i');
+ %padding
+ Pad = 2.0;
+ Sapj = zeros(numel(Kr), numel(Kz));
+ F = zeros(numel(Kr), numel(Kz));
+ G = zeros(numel(Kr), numel(Kz));
+ for n = 0:Jmax % Sum over Laguerre
+
+ for ikr = 1:numel(Kr)
+ for ikz = 1:numel(Kz)
+ kr = Kr(ikr); kz = Kz(ikz);
+ BA = sqrt(kr^2+kz^2)*...
+ (MODEL.sigma_e*sqrt(2) * (specie=='e')...
+ + sqrt(2*MODEL.tau_i) * (specie=='i'));
+
+ %First conv term
+ F(ikr,ikz) = (kz-kr).*phi(ikr,ikz).*kernel(n,BA);
+ %Second conv term
+ G(ikr,ikz) = 0.0;
+ for s = 0:min(n+j,Jmax)
+ G(ikr,ikz) = ...
+ G(ikr,ikz) + dnjs(n,j,s) .* squeeze(Napj(p+1,s+1,ikr,ikz));
+ end
+ G(ikr,ikz) = (kz-kr) .* G(ikr,ikz);
+ end
+ end
+ %Conv theorem
+ Sapj = Sapj + conv_thm_2D(F,G,Pad);
+ end
+end
+
+
diff --git a/matlab/conv_thm_2D.m b/matlab/conv_thm_2D.m
new file mode 100644
index 00000000..6f2b4452
--- /dev/null
+++ b/matlab/conv_thm_2D.m
@@ -0,0 +1,14 @@
+function [ conv ] = conv_thm_2D( F, G, Pad )
+%conv_thm_2D computes the convolution between F and G with a zero pading Pad
+% kspace -> real -> product -> kspace
+
+[Nr, Nz] = size(F);
+
+f = ifft2(F,Nr*Pad, Nz*Pad);
+g = ifft2(G,Nr*Pad, Nz*Pad);
+
+conv_pad = fft2(f.*g); % convolution becomes product
+
+conv = conv_pad(1:Nr,1:Nz); % remove padding
+end
+
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 8073feac..f83f2f2f 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -1,4 +1,4 @@
-function [ ] = create_gif(x, y, t, h, BASIC, GRID, MODEL, delay, GIFNAME)
+function [ ] = create_gif(x, y, t, h, BASIC, GRID, MODEL, delay, GIFNAME, saturate)
title1 = GIFNAME;
FIGDIR = ['../results/', BASIC.SIMID,'/'];
if ~exist(FIGDIR, 'dir')
@@ -28,26 +28,35 @@ else
pcolor(x,y,h(:,:,1)); % to set up
colormap jet
colorbar
- if not(flag)
+ if not(flag) && saturate
caxis([hmin hmax])
end
axis tight manual % this ensures that getframe() returns a consistent size
hold on
-
+
+ n = 0;
+ nbytes = fprintf(2,'frame %d/%d',n,numel(h(1,1,:)));
for n = 1:numel(h(1,1,:)) % time loop
- pcolor(x,y,h(:,:,n)); % frame plot
+ pclr = pcolor(x,y,h(:,:,n)); % frame plot
+ set(pclr, 'edgecolor','none')
title([title1,', $t \approx$', sprintf('%.3d',ceil(t(n)))])
drawnow
- % Capture the plot as an image
- frame = getframe(fig);
- im = frame2im(frame);
- [imind,cm] = rgb2ind(im,256);
- % Write to the GIF File
- if n == 1
- imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
- else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',delay);
- end
+ % Capture the plot as an image
+ frame = getframe(fig);
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,64);
+ % Write to the GIF File
+ if n == 1
+ imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',delay);
+ end
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(h(1,1,:)));
end
disp(['Gif saved @ : ',GIFNAME])
diff --git a/matlab/dnjs.m b/matlab/dnjs.m
index d9f8cc03..d333c1b2 100644
--- a/matlab/dnjs.m
+++ b/matlab/dnjs.m
@@ -4,8 +4,8 @@ function [ result ] = dnjs( n, j, s )
% last element of Coeffs is the larger one
L3 = flip(LaguerrePoly(Coeffs(end)));
% retrive smaller order coeff by taking the firsts (flipped array)
- L2 = L3(1:Coeffs(2));
- L1 = L3(1:Coeffs(1));
+ L2 = L3(1:Coeffs(2)+1);
+ L1 = L3(1:Coeffs(1)+1);
% build a factorial array to compute once every needed factorials
Factar = zeros(n+j+s+1,1);
@@ -21,7 +21,7 @@ function [ result ] = dnjs( n, j, s )
for il2 = 1:numel(L2)
for il1 = 1:numel(L1)
result = result + Factar(il1 + il2 + il3 -2)...
- /sqrt(L1(il1) * L2(il2) * L3(il3));
+ *L1(il1) * L2(il2) * L3(il3);
end
end
end
diff --git a/matlab/kernel.m b/matlab/kernel.m
new file mode 100644
index 00000000..f0530feb
--- /dev/null
+++ b/matlab/kernel.m
@@ -0,0 +1,15 @@
+function kernel_=kernel(n,x)
+% Evaluate the kernel function of order n. The normalized argument of the kernel function
+% is b_a = k_\perp \sigma_a \sqrt{2 \tau_a}.
+%
+% Note: - if n<0, then return 0.
+
+% Retrieve physical parameters
+
+if n >= 0
+ kernel_=(x./2).^(2*n).*exp(-x.^2/4)./factorial(n);
+else
+ kernel_ =0;
+end
+
+end
diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m
new file mode 100644
index 00000000..ebdc1538
--- /dev/null
+++ b/matlab/load_2D_data.m
@@ -0,0 +1,13 @@
+function [ data, kr, kz, time ] = load_2D_data( filename, variablename )
+%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
+ time = h5read(filename,'/data/var2d/time');
+ kr = h5read(filename,['/data/var2d/',variablename,'/coordkr']);
+ kz = h5read(filename,['/data/var2d/',variablename,'/coordkz']);
+ data = zeros(numel(kr),numel(kz),numel(time));
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(it,'%06d')]);
+ data(:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
+
+end
+
diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m
new file mode 100644
index 00000000..6f51cc1c
--- /dev/null
+++ b/matlab/load_5D_data.m
@@ -0,0 +1,13 @@
+function [ data, p, j, kr, kz, time ] = load_5D_data( filename, variablename )
+%LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ
+ time = h5read(filename,'/data/var5d/time');
+ p = h5read(filename,['/data/var5d/',variablename,'/coordp']);
+ j = h5read(filename,['/data/var5d/',variablename,'/coordj']);
+ kr = h5read(filename,['/data/var5d/',variablename,'/coordkr']);
+ kz = h5read(filename,['/data/var5d/',variablename,'/coordkz']);
+ data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time));
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(it,'%06d')]);
+ data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ end
+end
\ No newline at end of file
diff --git a/matlab/save_figure.m b/matlab/save_figure.m
index 1c0661e5..4d569a63 100644
--- a/matlab/save_figure.m
+++ b/matlab/save_figure.m
@@ -1,6 +1,6 @@
%% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME)
% and parameters
-FIGDIR = ['../results/', SIMID,'/'];
+FIGDIR = ['../results/', BASIC.SIMID,'/'];
if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 6eeb6643..b1536986 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -5,19 +5,20 @@ fid = fopen(INPUT,'wt');
fprintf(fid,'&BASIC\n');
fprintf(fid,[' nrun=', num2str(BASIC.nrun),'\n']);
fprintf(fid,[' dt=', num2str(BASIC.dt),'\n']);
-fprintf(fid,[' tmax=', num2str(BASIC.tmax),' ! time normalized to 1/omega_pe\n']);
+fprintf(fid,[' tmax=', num2str(BASIC.tmax),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&GRID\n');
-fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),' ! Electron Hermite moments \n']);
-fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),' ! Electron Laguerre moments \n']);
-fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),' ! Ion Hermite moments \n']);
-fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),' ! Ion Laguerre moments\n']);
+fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),'\n']);
+fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']);
+fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']);
+fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']);
fprintf(fid,[' nkr = ', num2str(GRID.nkr),'\n']);
fprintf(fid,[' krmin = ', num2str(GRID.krmin),'\n']);
fprintf(fid,[' krmax = ', num2str(GRID.krmax),'\n']);
fprintf(fid,[' nkz = ', num2str(GRID.nkz),'\n']);
fprintf(fid,[' kzmin = ', num2str(GRID.kzmin),'\n']);
fprintf(fid,[' kzmax = ', num2str(GRID.kzmax),'\n']);
+fprintf(fid,[' Pad = ', num2str(GRID.Pad),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
@@ -27,33 +28,36 @@ fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
fprintf(fid,[' write_Ni00 = ', OUTPUTS.write_Ni00,'\n']);
fprintf(fid,[' write_moments = ', OUTPUTS.write_moments,'\n']);
fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']);
+fprintf(fid,[' write_non_lin = ', OUTPUTS.write_non_lin,'\n']);
fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']);
+fprintf(fid,[' rstfile0 = ', OUTPUTS.rstfile0,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&MODEL_PAR\n');
fprintf(fid,' ! Collisionality\n');
-fprintf(fid,[' CO = ', num2str(MODEL.CO),' ! Collision operator (-1:Full Coulomb, 0: Dougherty)\n']);
-fprintf(fid,[' nu = ', num2str(MODEL.nu),' ! Normalized collision frequency normalized to plasma frequency\n']);
-fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),' ! T_e/T_e\n']);
-fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),' ! T_i/T_e temperature ratio\n']);
-fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),' ! sqrt(m_e/m_i) mass ratio\n']);
-fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),' ! sqrt(m_i/m_i)\n']);
-fprintf(fid,[' q_e = ', num2str(MODEL.q_e),' ! Electrons charge\n']);
-fprintf(fid,[' q_i = ', num2str(MODEL.q_i),' ! Ions charge\n']);
-fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),' ! L_perp / L_n Density gradient\n']);
-fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),' ! L_perp / L_T Temperature gradient\n']);
-fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),' ! L_perp / L_B Magnetic gradient and curvature\n']);
-fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),' ! Debye length\n']);
+fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']);
+fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']);
+fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']);
+fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']);
+fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
+fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
+fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
+fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
+fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
+fprintf(fid,[' eta_n = ', num2str(MODEL.eta_n),'\n']);
+fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),'\n']);
+fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),'\n']);
+fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&INITIAL_CON\n');
-fprintf(fid,' ! Background value\n');
-fprintf(fid,[' initback_moments=', num2str(INITIAL.initback_moments),' ! x 1e-3\n']);
-fprintf(fid,' ! Noise amplitude\n');
-fprintf(fid,[' initnoise_moments=', num2str(INITIAL.initnoise_moments),'\n']);
-fprintf(fid,[' iseed=', num2str(INITIAL.iseed),'\n']);
-fprintf(fid,[' selfmat_file=', INITIAL.selfmat_file,'\n']);
-fprintf(fid,[' eimat_file=', INITIAL.eimat_file,'\n']);
-fprintf(fid,[' iemat_file=', INITIAL.iemat_file,'\n']);
+fprintf(fid,[' RESTART =', INITIAL.RESTART,'\n']);
+fprintf(fid,[' only_Na00 =', INITIAL.only_Na00,'\n']);
+fprintf(fid,[' initback_moments =', num2str(INITIAL.initback_moments),'\n']);
+fprintf(fid,[' initnoise_moments =', num2str(INITIAL.initnoise_moments),'\n']);
+fprintf(fid,[' iseed =', num2str(INITIAL.iseed),'\n']);
+fprintf(fid,[' selfmat_file =', INITIAL.selfmat_file,'\n']);
+fprintf(fid,[' eimat_file =', INITIAL.eimat_file,'\n']);
+fprintf(fid,[' iemat_file =', INITIAL.iemat_file,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&TIME_INTEGRATION_PAR\n');
fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index b00087f2..8e9681ca 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -1,6 +1,5 @@
MODULE array
- !USE mumps_bsplines
use prec_const
implicit none
@@ -13,6 +12,12 @@ MODULE array
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj, CiepjT
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF
-
-END MODULE array
+ ! dnjs coefficient storage (in, ij, is)
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dnjs
+
+ ! Non linear term array (ip,ij,ikr,ikz)
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Sepj ! electron
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Sipj ! ion
+
+END MODULE array
diff --git a/src/auxval.F90 b/src/auxval.F90
index 8925ac1c..412809bd 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -9,7 +9,7 @@ subroutine auxval
IMPLICIT NONE
INTEGER :: irows,irowe, irow, icol
- WRITE(*,'(a/)') '=== Set auxiliary values ==='
+ WRITE(*,*) '=== Set auxiliary values ==='
call set_krgrid
call set_kzgrid
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 20e71f3a..859e81ac 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -4,16 +4,17 @@ MODULE basic
use prec_const
IMPLICIT none
- INTEGER :: nrun=1 ! Number of time steps to run
+ LOGICAL :: RESTART=.true. ! To restart the simulation or look for saved state
+ INTEGER :: nrun=1 ! Number of time steps to run
real(dp) :: tmax=100000.0 ! Maximum simulation time
real(dp) :: dt=1.0 ! Time step
real(dp) :: time=0 ! Current simulation time (Init from restart file)
-
+
INTEGER :: jobnum=0 ! Job number
INTEGER :: step=0 ! Calculation step of this run
INTEGER :: cstep=0 ! Current step number (Init from restart file)
LOGICAL :: nlend=.FALSE. ! Signal end of run
-
+
INTEGER :: ierr ! flag for MPI error
INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose)
INTEGER :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose)
@@ -24,6 +25,8 @@ MODULE basic
INTEGER :: lu_in = 90 ! File duplicated from STDIN
INTEGER :: lu_job = 91 ! myjob file
+ real :: start, finish ! To measure computation time
+
INTERFACE allocate_array
MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4
MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5
@@ -39,11 +42,11 @@ CONTAINS
use prec_const
IMPLICIT NONE
-
- NAMELIST /BASIC/ nrun, dt, tmax
-
+
+ NAMELIST /BASIC/ nrun, dt, tmax
+
READ(lu_in,basic)
- WRITE(*,basic)
+ !WRITE(*,basic)
END SUBROUTINE basic_data
!================================================================================
@@ -52,7 +55,7 @@ CONTAINS
use prec_const
IMPLICIT NONE
-
+
CHARACTER(len=*) , INTENT(in) :: str
CHARACTER(len=16) :: d, t, dat, time
!________________________________________________________________________________
@@ -78,7 +81,7 @@ CONTAINS
SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2)
IMPLICIT NONE
real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
- INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
+ INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=0.0_dp
END SUBROUTINE allocate_array_dp2
@@ -120,7 +123,7 @@ CONTAINS
SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2)
IMPLICIT NONE
DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
- INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
+ INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=CMPLX(0.0_dp,0.0_dp)
END SUBROUTINE allocate_array_dc2
@@ -162,7 +165,7 @@ CONTAINS
SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2)
IMPLICIT NONE
INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
- INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
+ INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
ALLOCATE(a(is1:ie1,is2:ie2))
a=0
END SUBROUTINE allocate_array_i2
@@ -189,7 +192,7 @@ CONTAINS
INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
a=0
- END SUBROUTINE allocate_array_i5
+ END SUBROUTINE allocate_array_i5
!========================================
@@ -233,15 +236,4 @@ CONTAINS
a=.false.
END SUBROUTINE allocate_array_l5
- RECURSIVE FUNCTION Factorial(n) RESULT(Fact)
- IMPLICIT NONE
- INTEGER :: Fact
- INTEGER, INTENT(IN) :: n
- IF (n == 0) THEN
- Fact = 1
- ELSE
- Fact = n * Factorial(n-1)
- END IF
- END FUNCTION Factorial
-
END MODULE basic
diff --git a/src/coeff_mod.F90 b/src/coeff_mod.F90
new file mode 100644
index 00000000..b1c8e276
--- /dev/null
+++ b/src/coeff_mod.F90
@@ -0,0 +1,95 @@
+MODULE coeff
+
+! this module contains routines to compute normalization coefficients and velocity integrals
+!
+
+USE PREC_CONST
+use BASIC
+USE MODEL
+USE FMZM
+
+PUBLIC
+
+CONTAINS
+
+ !-----------------------------------------------------------
+ !> @brief
+ !> Computes the associated-Laguerre/Laguerre to Laguerre basis transformation coefficient \f$ \overline{d}^{|m|}_{nkf} \f$,
+ !>
+ !> \f[ \overline{d}^{|m|}_{nkf} = \sum_{n_1 = 0}^n \sum_{k_1 =0}^k \sum_{f_1=0}^f L_{kk_1}^{-1/2} L_{nn_1}^{|m|-1/2} L_{f f_1}^{-1/2} (n_1 + k_1 + f_1)!. \f]
+ !> @todo Checked against Mathematica
+ !> @param[in] n
+ !> @param[in] k
+ !> @param[in] s
+ !> @param[in] m
+ !> @param[out] XOUT
+ !-----------------------------------------------------------
+ FUNCTION ALL2L(n,k,s,m) RESULT(XOUT)
+ !
+ ! Computes the associated-Laguerre/Laguerre to Laguerre coefficient, i.e.
+ !
+ ! L_n^m L_k = sum_{s=0}^{n+k}\overline{d}_{nks}^m L_s
+ !
+ !
+ !
+ IMPLICIT NONE
+ !
+ INTEGER, intent(in) :: n,k,s,m ! ... degree and orders of Laguerre polynomials
+ TYPE(FM) :: XOUT ! ... coefficient values
+ TYPE(FM), SAVE :: AUXN,AUXK,AUXS
+ INTEGER :: n1,k1,s1
+ !
+ CALL FM_ENTER_USER_FUNCTION(XOUT)
+ ! Compute coefficients
+ !
+ XOUT = TO_FM('0.0')
+ !
+ DO n1=0,n
+ AUXN =Lpjl(REAL(m,dp)-0.5_dp,REAL(n,dp),REAL(n1,dp))
+ DO k1=0,k
+ AUXK =Lpjl(-0.5_dp,REAL(k,dp),REAL(k1,dp))
+ DO s1=0,s
+ AUXS = Lpjl(-0.5_dp,REAL(s,dp),REAL(s1,dp))
+ XOUT = XOUT + FACTORIAL(TO_FM(n1 + k1 + s1 ))*AUXN*AUXK*AUXS
+ ENDDO
+ ENDDO
+ ENDDO
+ !
+ CALL FM_EXIT_USER_FUNCTION(XOUT)
+ !
+ END FUNCTION ALL2L
+
+ !-----------------------------------------------------------
+ !> @brief
+ !> Computes the associated Laguerre serie coefficient \f$ L_{jl}^p \f$
+ !>
+ !> \f[ L_{jl}^{p} = \frac{(-1)^l (p + j + 1/2)!}{(j - l)!(l +p + 1/2 )! l!}. \f]
+ !> such that
+ !> \f[ L^{p+1/2}_j(x) = \sum_{l=0}^j L^p_{jl} x^l \f]
+ !> @param[in] XOUT
+ !-----------------------------------------------------------
+ FUNCTION Lpjl(p,j,l) RESULT(XOUT)
+ !
+ ! Computes the associated Laguerre serie coefficient,
+ !
+
+ IMPLICIT NONE
+
+ !
+ REAL(dp), intent(in) :: p,j,l
+ TYPE(FM) :: XOUT
+ !
+ CALL FM_ENTER_USER_FUNCTION(XOUT)
+
+ ! compute coeff
+
+ XOUT = TO_FM('0.0')
+
+ XOUT = TO_FM(((-1)**l))*FACTORIAL(TO_FM(2*p + 2*j + 1)/2)/&
+ (FACTORIAL(TO_FM( j -l ))*FACTORIAL(TO_FM(2*l + 2*p + 1)/2)*FACTORIAL(TO_FM(l)))
+
+ CALL FM_EXIT_USER_FUNCTION(XOUT)
+
+ END FUNCTION Lpjl
+
+END MODULE coeff
diff --git a/src/coeff_mod.f90 b/src/coeff_mod.f90
deleted file mode 100644
index 8db9837e..00000000
--- a/src/coeff_mod.f90
+++ /dev/null
@@ -1,633 +0,0 @@
-MODULE coeff
-
-! this module contains routines to compute normalization coefficients and velocity integrals
-!
-
-USE PREC_CONST
-use BASIC
-USE MODEL
-USE FMZM
-USE basis_transformation
-
-PUBLIC
-
-CONTAINS
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the associated-Laguerre/Laguerre/x to Laguerre basis transformation coefficient \f$ d^{|m|}_{nks} \f$,
- !>
- !> \f[ L_n^{|m|}(x) L_k(x) x^{|m|} =\sum_{s=0}^{n+|m| +k} d^{|m|}_{nks} L_s(x),\f]
- !> where
- !> \f[d^{|m|}_{nks} = \int_0^\infty d x e^{-x} L_n^{|m|}(x) L_s(x) L_k(x) x^{|m|}. \f]
- !> yielding the closed analyticalform,
- !> \f[ d_{nks}^{|m|} =\sum_{n_1=0}^{n} \sum_{k_1=0}^k \sum_{s_1=0}^s L_{kk_1}^{-1/2} L_{nn_1}^{|m|-1/2} L_{ss_1}^{-1/2} (n_1 + k_1 + s_1 + |m|)!, \f]
- !
- !> @param[in] n
- !> @param[in] k
- !> @param[in] s
- !> @param[in] m
- !>
- !> @param[out] XOUT Output result
- !-----------------------------------------------------------
- FUNCTION ALLX2L(n,k,s,m) RESULT(XOUT)
- !
- ! Computes the associated-Laguerre/Laguerre/x to Laguerre coefficient, i.e.
- !
- ! L_n^m L_k x^m = \sum_{s=0}^{n+m+k}d_{nks}^m L_s
- !
- !
- !
- IMPLICIT NONE
-
- !
- INTEGER, intent(in) :: n,k,s,m ! ... degree and orders of Laguerre polynomials
- TYPE(FM) :: XOUT ! ... coefficient values
- TYPE(FM),SAVE:: AUXN,AUXK,AUXS
- INTEGER :: n1,k1,s1
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- ! Compute coefficients
- !
- XOUT = TO_FM('0.0')
- !
- DO n1=0,n
- AUXN =Lpjl(REAL(m,dp)-0.5_dp,REAL(n,dp),REAL(n1,dp))
- DO k1=0,k
- AUXK =Lpjl(-0.5_dp,REAL(k,dp),REAL(k1,dp))
- DO s1=0,s
- AUXS = Lpjl(-0.5_dp,REAL(s,dp),REAL(s1,dp))
- XOUT = XOUT + FACTORIAL(TO_FM(n1 + k1 + s1 + m))*AUXN*AUXK*AUXS
- ENDDO
- ENDDO
- ENDDO
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION ALLX2L
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the associated-Laguerre/Laguerre to Laguerre basis transformation coefficient \f$ \overline{d}^{|m|}_{nkf} \f$,
- !>
- !> \f[ \overline{d}^{|m|}_{nkf} = \sum_{n_1 = 0}^n \sum_{k_1 =0}^k \sum_{f_1=0}^f L_{kk_1}^{-1/2} L_{nn_1}^{|m|-1/2} L_{f f_1}^{-1/2} (n_1 + k_1 + f_1)!. \f]
- !> @todo Checked against Mathematica
- !> @param[in] n
- !> @param[in] k
- !> @param[in] s
- !> @param[in] m
- !> @param[out] XOUT
- !-----------------------------------------------------------
- FUNCTION ALL2L(n,k,s,m) RESULT(XOUT)
- !
- ! Computes the associated-Laguerre/Laguerre to Laguerre coefficient, i.e.
- !
- ! L_n^m L_k = sum_{s=0}^{n+k}\overline{d}_{nks}^m L_s
- !
- !
- !
- IMPLICIT NONE
- !
- INTEGER, intent(in) :: n,k,s,m ! ... degree and orders of Laguerre polynomials
- TYPE(FM) :: XOUT ! ... coefficient values
- TYPE(FM), SAVE :: AUXN,AUXK,AUXS
- INTEGER :: n1,k1,s1
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
-
- ! Compute coefficients
- !
- XOUT = TO_FM('0.0')
- !
- DO n1=0,n
- AUXN =Lpjl(REAL(m,dp)-0.5_dp,REAL(n,dp),REAL(n1,dp))
- DO k1=0,k
- AUXK =Lpjl(-0.5_dp,REAL(k,dp),REAL(k1,dp))
- DO s1=0,s
- AUXS = Lpjl(-0.5_dp,REAL(s,dp),REAL(s1,dp))
- XOUT = XOUT + FACTORIAL(TO_FM(n1 + k1 + s1 ))*AUXN*AUXK*AUXS
- ENDDO
- ENDDO
- ENDDO
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION ALL2L
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the associated Laguerre serie coefficient \f$ L_{jl}^p \f$
- !>
- !> \f[ L_{jl}^{p} = \frac{(-1)^l (p + j + 1/2)!}{(j - l)!(l +p + 1/2 )! l!}. \f]
- !> such that
- !> \f[ L^{p+1/2}_j(x) = \sum_{l=0}^j L^p_{jl} x^l \f]
- !> @param[in] XOUT
- !-----------------------------------------------------------
- FUNCTION Lpjl(p,j,l) RESULT(XOUT)
- !
- ! Computes the associated Laguerre serie coefficient,
- !
-
- IMPLICIT NONE
-
- !
- REAL(dp), intent(in) :: p,j,l
- TYPE(FM) :: XOUT
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
-
- ! compute coeff
-
- XOUT = TO_FM('0.0')
-
- XOUT = TO_FM(((-1)**l))*FACTORIAL(TO_FM(2*p + 2*j + 1)/2)/&
- (FACTORIAL(TO_FM( j -l ))*FACTORIAL(TO_FM(2*l + 2*p + 1)/2)*FACTORIAL(TO_FM(l)))
-
- CALL FM_EXIT_USER_FUNCTION(XOUT)
-
- END FUNCTION Lpjl
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the particle species 'a' FLR coefficients
- !>
- !> \f[ \mathcal{A}_{an}^m = \frac{ \sigma_m n!}{(n+|m|)!}\left(\frac{b_a}{2} \right)^{|m|} \mathcal{K}_{n}(b_a) \f]
- !>
- !> @param[in] m Spherical Harmonics order
- !> @param[in] n FLR order
- !> @param[in] ba Thermal Normalized Perpendicular Wavenumber
- !> @param[out] Amn
- !-----------------------------------------------------------
- FUNCTION curlyAmn(m,n,ba) RESULT(XOUT)
- ! compute FLR coeff. for particle a
- ! CHECKED
- IMPLICIT NONE
- !
- INTEGER,INTENT(in) :: m,n
- TYPE(FM),INTENT(in) :: ba
- !
- TYPE(FM):: XOUT
- TYPE(FM), SAVE :: ker
- INTEGER :: sigmam
- !
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM('0.0')
- !
- sigmam = sign(1,m)**m
- !
- ker = kerneln(ba,n)
- !
- IF( m .eq. 0) THEN ! NO FLR EFFECTS
- XOUT = ker
- ELSE
- XOUT = sigmam*(FACTORIAL(TO_FM(n))/FACTORIAL(TO_FM(n +ABS(m))))*((ba/2)**ABS(m))*ker
- ENDIF
- !
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION curlyAmn
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the particle species 'b' FLR coefficients
- !> @param[in] m Spherical Harmonics order
- !> @param[in] n FLR order
- !> @param[in] bb Thermal Normalized Perpendicular Wavenumber
- !> @param[out] Bmn
- !-----------------------------------------------------------
- FUNCTION curlyBmn(m,n,b_) RESULT(XOUT)
- !
- ! compute FLR coeff. for particle b
- ! CHECKED
- IMPLICIT NONE
- !
- INTEGER,INTENT(in) :: m,n
- TYPE(FM),INTENT(in) :: b_
- !
- TYPE(FM):: XOUT
- TYPE(FM),SAVE :: ker
- INTEGER :: sigmam
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM('0.0')
- !
- sigmam = sign(1,m)**m
- !
- ker = kerneln(b_,n)
- !
- IF(m .eq. 0) THEN !
- XOUT = ker
- ELSE
- XOUT = sigmam*(FACTORIAL(TO_FM(n))/FACTORIAL(TO_FM(n +ABS(m))))*((b_/2)**ABS(m))*ker
- ENDIF
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
-
- END FUNCTION curlyBmn
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the nth-order kernel function
- !>
- !> \f[ \mathcal{K}_n(x) = \frac{1}{n!} \left(\frac{x}{2} \right)^{2n } e^{- x^2/4} \f]
- !> @param[in] n kernel order
- !> @param[in] b normalized perpendicular wavevector
- !> @param[out] kerneln
- !-----------------------------------------------------------
- FUNCTION kerneln(b,n) RESULT(XOUT)
- ! compute the nth-order kernel function
- ! CHECKED
- IMPLICIT NONE
- !
- INTEGER,INTENT(in):: n
- TYPE(FM), INTENT(in) :: b
- !
- TYPE(FM) :: XOUT
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
-
- XOUT = TO_FM(0)
- ! CHECKED
-
- IF (n .LT. 0) THEN ! ... not defined return 0
- XOUT = TO_FM('0.0')
- ELSE
- IF( b .eq. TO_FM('0.0') .and. n .eq. 0) THEN ! NO FLR effetcs
- XOUT = TO_FM('1.0')
- ELSE
- XOUT = ((b/2)**n)*EXP(- (b/2)*(b/2))/FACTORIAL(TO_FM(n))
- ENDIF
- ENDIF
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
-
- END FUNCTION kerneln
-
-
-
- FUNCTION normepm(p,m) RESULT(XOUT)
- ! compute the norm of e^{pm} \cdot e^{pm}
- IMPLICIT NONE
- !
- INTEGER,intent(in) :: p,m
- !
- TYPE(FM) :: XOUT
-
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- !
- XOUT = TO_FM(0)
- !
- XOUT = 2**(2*p)*FACTORIAL(TO_FM(p))*FACTORIAL( TO_FM(2*p -1)/2)/FACTORIAL(TO_FM(2*p))/SQRT(PIFM)
-
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
-
- END FUNCTION normepm
-
- FUNCTION Injlkmp(n,j,l,k,m,p) RESULT(XOUT)
- ! Compute speed integrals in Lorentz operator
- ! Note: Only the zeroth-order harmonics, i.e. consider m=0 !!
- ! Not defined if p ==0
- IMPLICIT NONE
- !
- INTEGER,INTENT(in) :: n,j,l,k,m,p
- TYPE(FM) :: XOUT
- TYPE(FM),SAVE :: NX1,NX2,T4inv,d0nkf_,L1,L2
- CHARACTER(len=T4len) :: T4str_
- !
- INTEGER :: f,g,h,h1,j1
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM(0)
-
- floop:DO f=0,n+k
- ! Associated Laguerre basis transformation
- d0nkf_ = ALL2L(n,k,f,0)
- gloop:DO g=0,l+2*f
- IF( g .eq. p) THEN
- hloop:DO h=0,f+l/2
- ! Inverse Basis transformation
- ! _______________________________________________________________________________________
- ! Previous method: very slow at the execution time
- ! T4inv = READ_T4_FROM_CSV(l,f,p,h,1) ! ... (T^-1)^gh_lf = ... T^lf_gh
- !
- !_______________________________________________________________________________________
- ! New method: Much faster than previous method
- ! write(*,*) l,f,g,h
- T4str_ = GET_T4(l,f,g,h) ! ... (T^-1)^gh_lf = ... T^lf_gh
-
- T4inv = SQRT(PIFM)*(2**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*(2*g + 1)/2/FACTORIAL(TO_FM(2*g+2*h +1)/2)*TO_FM(T4str_)
- !________________________________________________________________________________________
- IF( T4inv .ne. TO_FM(0)) THEN
-
- h1loop:DO h1=0,h
- L1= Lpjl(real(p,dp),real(h,dp),real(h1,dp))
- j1loop: DO j1=0,j
- L2 = Lpjl(real(p,dp),real(j,dp),real(j1,dp))
- NX1 = FACTORIAL(TO_FM(p + h1 + j1 -1))
- XOUT = XOUT + L2*L1*T4inv*d0nkf_*NX1*2/TO_FM(2*p + 1)/SQRT(PIFM)
- ENDDO j1loop
- ENDDO h1loop
- ENDIF
- ENDDO hloop
- ENDIF
- ENDDO gloop
- ENDDO floop
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- !
-END FUNCTION Injlkmp
-
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the spherical harmonics normalization coefficient
- !>
- !> \f[ \f]
- !> @param[in] p
- !> @param[in] j
- !> @param[out] Results
- !-----------------------------------------------------------
- FUNCTION sigmapj(p,j) RESULT(XOUT)
- ! compute spherical harmonics normalization coefficient
- !
- IMPLICIT NONE
- !
- INTEGER,INTENT(in):: p
- INTEGER, INTENT(in) :: j
- !
- TYPE(FM) :: XOUT
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
-
- XOUT = TO_FM(0)
- !
- XOUT = FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(2*p+2*j+1)/2)/(TO_FM(2)**p)/FACTORIAL(TO_FM(2*p+1)/2)/FACTORIAL(TO_FM(j))
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
-
- END FUNCTION sigmapj
- !
- ! _________________ LORENTZ COEFFICIENTS_________________________
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the drift kinetic \mathcal{A}_{\parallel ei}^{lk0}
- !>
- !> \f[ \f]
- !> @param[in] p
- !> @param[in] j
- !> @param[out] Results
- !-----------------------------------------------------------
- FUNCTION curlyAlk0pareiDK(l,k) RESULT(XOUT)
- IMPLICIT NONE
- !
- INTEGER,INTENT(in):: l,k
- !
- TYPE(FM) :: XOUT
- !
- ! LOCAL VARIABLES
- INTEGER :: h
- TYPE(FM),SAVE :: NX0,NX1,NX2,T4inv
- CHARACTER(len=T4len) :: T4invstr_
-
-
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM(0)
- !
- NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*TO_FM(2)**l)
- !
- IF(l .ge. 1 .or. k .ge. 1) THEN
- DO h=0, k + l/2
- ! INVERSE BASIS TRANSFORMATION
- T4invstr_ = GET_T4(1,h,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
- T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*h +1)/2)
-
- XOUT = XOUT + T4inv*NX0*FACTORIAL(TO_FM(2*h + 1)/2)/FACTORIAL(TO_FM(h))
- !
- ENDDO
- ENDIF
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION curlyAlk0pareiDK
-
-
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the gyrokinetic \mathcal{A}_{\parallel ei}^{lk0}
- !>
- !> \f[ \f]
- !> @param[in] p
- !> @param[in] j
- !> @param[out] Results
- !-----------------------------------------------------------
- FUNCTION curlyAlk0pareiGK(l,k) RESULT(XOUT)
- IMPLICIT NONE
- !
- INTEGER,INTENT(in):: l,k
- !
- TYPE(FM) :: XOUT
- !
- ! LOCAL VARIABLES
- INTEGER :: h,g,f,n
- TYPE(FM),SAVE :: NX0,T4inv,ker_,d0nkf_
- CHARACTER(len=T4len) :: T4invstr_
-
-
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM(0)
- !
- NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*TO_FM(2)**l)
- !
- nloop: DO n=0, nimaxxFLR
- ker_ = kerneln(bbfm_,n)
- floop: DO f=0, n + k
- d0nkf_ = ALLX2L(n,k,f,0)
- gloop: DO g=0, l+2*f
- IF(g .eq. 1) THEN
- hloop: DO h=0, f + l/2
- ! INVERSE BASIS TRANSFORMATION
- T4invstr_ = GET_T4(1,h,l,f) ! ... (T^-1)^pt_lk = ... T^lk_pt
- T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM(2*1 + 1)/2/FACTORIAL(TO_FM(2*1+2*h +1)/2)
-
- XOUT = XOUT + T4inv*NX0*d0nkf_*ker_*FACTORIAL(TO_FM(2*h + 1)/2)/FACTORIAL(TO_FM(h))
- !
- ENDDO hloop
- ENDIF
- ENDDO gloop
- ENDDO floop
- ENDDO nloop
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION curlyAlk0pareiGK
-
- !-----------------------------------------------------------
- !> @brief
- !> Computes the gyrokinetic \mathcal{A}_{\perp ei}^{lk0}
- !>
- !> \f[ \f]
- !> @param[in] p
- !> @param[in] j
- !> @param[out] Results
- !-----------------------------------------------------------
- FUNCTION curlyAlk0perpeiGK(l,k) RESULT(XOUT)
- IMPLICIT NONE
- !
- INTEGER,INTENT(in):: l,k
- !
- TYPE(FM) :: XOUT
- !
- ! LOCAL VARIABLES
- INTEGER :: h,g,f,n,h1
- TYPE(FM),SAVE :: NX0,T4inv,ker_,L0hh1,d1nkf_
- CHARACTER(len=T4len) :: T4invstr_
-
-
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM(0)
- !
- NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*TO_FM(2)**l)
- !
- nloop: DO n=0, nimaxxFLR
- ker_ = kerneln(bbfm_,n)
- floop: DO f=0, n + k +1
- d1nkf_ = ALLX2L(n,k,f,1)
- hloop: DO h=0, f + l/2
- ! INVERSE BASIS TRANSFORMATION
- T4invstr_ = GET_T4(0,h,l,f) ! ... (T^-1)^pt_lk = ... T^lk_pt
- T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*h +1)/2)
- h1loop: DO h1=1,h
- L0hh1 =Lpjl(real(0,dp),real(h,dp),real(h1,dp))
- XOUT = XOUT + T4inv*NX0*d1nkf_*ker_*bbfm_*L0hh1/2/(n+1)*FACTORIAL(TO_FM(h1-1))
- !
- ENDDO h1loop
- ENDDO hloop
- ENDDO floop
- ENDDO nloop
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION curlyAlk0perpeiGK
- !
- !__________________________________________________________________
- !
- FUNCTION X2prodLegendre(g,h) RESULT(XOUT)
- !
- ! compute \int d x P_g(x) P_h(x) x^2 = d_h^g
- !
- ! Ref : Arfkten G. B. and Weber H. J., Mathematical methods for physicist, p. 806
- !
- IMPLICIT NONE
- !
- INTEGER,INTENT(in) :: g,h
- TYPE(FM) :: XOUT
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = TO_FM(0)
- !
- IF(g .eq. h) THEN
- XOUT = 2*TO_FM((2*TO_FM(h)**2 + 2*h -1))/(2*h-1)/(2*h +1 )/(2*h +3)
- ELSEIF (g .eq. h -2) THEN
- XOUT = TO_FM(2*h*(h-1))/(2*h -3)/(2*h-1)/(2*h+1)
- ELSEIF( g .eq. h +2) THEN
- XOUT = TO_FM(2*(h+1)*(h+2))/(2*h +1)/(2*h + 3)/(2*h +5)
- ENDIF
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION X2prodLegendre
- !
- !__________________________________________________________________
- !
- FUNCTION CJ(j) RESULT(XOUT)
- ! Associated Laguerre Normalization factor
- IMPLICIT NONE
- !
- INTEGER,INTENT(in) :: j
- TYPE(FM) :: XOUT
- !
- CALL FM_ENTER_USER_FUNCTION(XOUT)
- !
- XOUT = FACTORIAL( TO_FM(2*j+3))*3*TO_FM(2)**(3*j +3)*FACTORIAL(TO_FM(j))/FACTORIAL(TO_FM(4*j +6))
- !
- CALL FM_EXIT_USER_FUNCTION(XOUT)
- !
- END FUNCTION
- !____________________________________________________
- !
- FUNCTION BinomialFM( n, k ) result(XOUT )
- ! Compute the binomial coefficient for positive and negative integers
- implicit none
- integer, intent(in) :: n, k
- TYPE(FM) :: XOUT
- type(FM), save :: f1, f2, f3
- ! See http://mathworld.wolfram.com/BinomialCoefficient.html
-
- call FM_ENTER_USER_FUNCTION(XOUT)
-
- XOUT = TO_FM(0)
- if( k>=0 .and. n>=k ) then
- f1 = FACTORIAL(TO_FM(n))
- f2 = FACTORIAL(TO_FM(k))
- f3 = FACTORIAL(TO_FM(n-k))
- XOUT = f1/(f2*f3)
- else if( n<0 ) then
- if( k>=0 ) then
- f1 = FACTORIAL(TO_FM( -n+k-1 ))
- f2 = FACTORIAL( TO_FM(k ))
- f3 = FACTORIAL( TO_FM(-n-1 ))
- XOUT = TO_FM(-1)**k * f1 / ( f2*f3 )
- else if( k<=n ) then
- f1 = FACTORIAL( TO_FM(-k-1 ))
- f2 = FACTORIAL( TO_FM(n-k ))
- f3 = FACTORIAL( -n-1 )
- XOUT = TO_FM(-1)**(n-k)*f1/(f2*f3)
- end if
- end if
- call FM_EXIT_USER_FUNCTION(XOUT)
- END FUNCTION BinomialFM
-!___________________________________________________
-END MODULE coeff
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
new file mode 100644
index 00000000..ca2787e3
--- /dev/null
+++ b/src/compute_Sapj.F90
@@ -0,0 +1,109 @@
+SUBROUTINE compute_Sapj(Sepj, Sipj)
+
+ USE array, ONLY : dnjs
+ USE basic
+ USE convolution
+ USE fields!, ONLY : phi, moments_e, moments_i
+ USE fourier_grid
+ USE model
+ USE prec_const
+ USE time_integration!, ONLY : updatetlevel
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: F_, G_, CONV
+ COMPLEX(dp), DIMENSION(ips_e:ipe_e, ijs_e:ije_e, Nkr, Nkz) :: Sepj
+ COMPLEX(dp), DIMENSION(ips_i:ipe_i, ijs_i:ije_i, Nkr, Nkz) :: Sipj
+ INTEGER :: ip, ij, in, is, ikr, ikz
+ REAL(dp):: kr, kz, kernel, be_2, bi_2, factn
+ REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
+
+ !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
+ sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
+ ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
+ jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
+ Sepj(ip,ij,:,:) = 0._dp
+ factn = 1
+
+ nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
+
+ krloope: DO ikr = 1,Nkr ! Loop over kr
+ kzloope: DO ikz = 1,Nkz ! Loop over kz
+ kr = krarray(ikr)
+ kz = kzarray(ikz)
+ be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
+ kernel = be_2**(in-1)/factn * EXP(-be_2)
+
+ ! First convolution term
+ F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) * kernel
+ ! Second convolution term
+ G_(ikr,ikz) = 0._dp ! initialization of the sum
+ DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
+ G_(ikr,ikz) = G_(ikr,ikz) + &
+ dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
+ ENDDO
+ G_(ikr,ikz) = (kz - kr) * G_(ikr,ikz)
+ ENDDO kzloope
+ ENDDO krloope
+
+ CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve and go back to Fourier space
+ !CALL convolve_2D_F2R( F_, G_, CONV ) ! .. or convolve and keep the results in real space**
+
+ Sepj(ip,ij,:,:) = Sepj(ip,ij,:,:) + CONV ! Add it to Sepj (Real space here)
+
+ IF ( in + 1 .LE. jmaxe+1 ) THEN
+ factn = real(in,dp) * factn ! factorial(n+1)
+ ENDIF
+ ENDDO nloope
+
+ !CALL forward_FFT(Sepj(ip,ij,:,:)) !**then put the sum back to fourier space (spares n FFT)
+
+ ENDDO jloope
+ ENDDO ploope
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
+ sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
+
+ ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
+ jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
+ Sipj(ip,ij,:,:) = 0._dp
+ factn = 1
+
+ nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
+
+ krloopi: DO ikr = 1,Nkr ! Loop over kr
+ kzloopi: DO ikz = 1,Nkz ! Loop over kz
+ kr = krarray(ikr)
+ kz = kzarray(ikz)
+ bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
+ kernel = bi_2**(in-1)/factn * EXP(-bi_2)
+
+ F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) * kernel
+
+ G_(ikr,ikz) = 0._dp ! initialization of the sum
+ DO is = 1, MIN( in+ij-1, jmaxi+1 )
+ G_(ikr,ikz) = G_(ikr,ikz) + &
+ dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
+ ENDDO
+ G_(ikr,ikz) = (kz - kr) * G_(ikr,ikz)
+ ENDDO kzloopi
+ ENDDO krloopi
+
+ CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve and back to Fourier
+ !CALL convolve_2D_F2R( F_, G_, CONV ) ! or convolve and keep the results in real space**
+
+ Sipj(ip,ij,:,:) = Sipj(ip,ij,:,:) + CONV ! Add it to Sipj (Real space here)
+
+ IF ( in + 1 .LE. jmaxi+1 ) THEN
+ factn = real(in,dp) * factn ! factorial(n+1)
+ ENDIF
+
+ ENDDO nloopi
+
+ !CALL forward_FFT(Sipj(ip,ij,:,:)) ! **then put the sum back to fourier space (spares n FFT)
+
+ ENDDO jloopi
+ ENDDO ploopi
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+END SUBROUTINE compute_Sapj
diff --git a/src/control.F90 b/src/control.F90
index 273ce8ad..9218566a 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -2,16 +2,16 @@ SUBROUTINE control
! Control the run
use basic
-
use prec_const
IMPLICIT NONE
+ call cpu_time(start)
!________________________________________________________________________________
! 1. Prologue
! 1.1 Initialize the parallel environment
WRITE(*,*) 'Initialize MPI...'
CALL ppinit
- WRITE(*,*) '...MPI initialized'
+ WRITE(*,'(a/)') '...MPI initialized'
call daytim('Start at ')
@@ -19,46 +19,49 @@ SUBROUTINE control
! 1.2 Define data specific to run
WRITE(*,*) 'Load basic data...'
CALL basic_data
- WRITE(*,*) '...basic data loaded.'
+ WRITE(*,'(a/)') '...basic data loaded.'
! 1.3 Read input parameters from input file
WRITE(*,*) 'Read input parameters...'
CALL readinputs
- WRITE(*,*) '...input parameters read'
+ WRITE(*,'(a/)') '...input parameters read'
! 1.4 Set auxiliary values (allocate arrays, set laplace operator, ...)
WRITE(*,*) 'Calculate auxval...'
CALL auxval
- WRITE(*,*) '...auxval calculated'
+ WRITE(*,'(a/)') '...auxval calculated'
!
! 1.5 Initial conditions
WRITE(*,*) 'Create initial state...'
CALL inital
- WRITE(*,*) '...initial state created'
+ WRITE(*,'(a/)') '...initial state created'
! 1.6 Initial diagnostics
WRITE(*,*) 'Initial diagnostics...'
CALL diagnose(0)
- WRITE(*,*) '...initial diagnostics done'
+ WRITE(*,'(a/)') '...initial diagnostics done'
!
CALL FLUSH(stdout)
+ !________________________________________________________________________________
+ WRITE(*,*) 'Time integration loop..'
!________________________________________________________________________________
! 2. Main loop
DO
step = step + 1
cstep = cstep + 1
+
CALL stepon
time = time + dt
- CALL diagnose(step)
-
CALL tesend
IF( nlend ) EXIT ! exit do loop
- ! CALL write_restart ! if want to write a restart file every so often (in case of crash)
+ CALL diagnose(step)
+
END DO
+ WRITE(*,'(a/)') '...time integration done'
!________________________________________________________________________________
! 9. Epilogue
diff --git a/src/convolution_mod.F90 b/src/convolution_mod.F90
new file mode 100644
index 00000000..b50e9173
--- /dev/null
+++ b/src/convolution_mod.F90
@@ -0,0 +1,221 @@
+MODULE convolution
+ USE prec_const
+ implicit none
+
+ CONTAINS
+
+ !!! Convolution 2D Fourier to Fourier
+ ! - Compute the convolution using the convolution theorem and MKL
+ SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
+
+ USE prec_const
+ USE fourier_grid, ONLY : Nkr, Nkz, Pad
+ USE MKL_DFTI
+
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: F_2D, G_2D, C_2D
+ COMPLEX(dp), DIMENSION(Pad*Nkr*Pad*Nkz) :: F_1D, G_1D, C_1D
+
+ INTEGER :: ix, iy, Mkr, Mkz
+ INTEGER :: Status, L(2), L_pad(2)
+
+ type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle,My_Desc2_Handle,My_Desc3_Handle
+
+ L = (/ Nkr, Nkz /)
+ Mkr = Pad*Nkr; Mkz = Pad*Nkz
+ L_pad = (/ Mkr, Mkz /)
+
+ !WRITE(*,*) 'Reshaping..'
+ DO ix = 1,Mkr
+ DO iy = 1,Mkz
+ IF ( (ix .LE. Nkr) .AND. (iy .LE. Nkz) ) THEN
+ F_1D(Mkr*(iy-1) + ix) = F_2D(ix,iy)
+ G_1D(Mkr*(iy-1) + ix) = G_2D(ix,iy)
+ ELSE
+ F_1D(Mkr*(iy-1) + ix) = 0._dp
+ G_1D(Mkr*(iy-1) + ix) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+
+ Status = DftiCreateDescriptor(My_Desc1_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 2, L_pad)
+ Status = DftiCommitDescriptor(My_Desc1_Handle)
+ !WRITE(*,*) 'Backward FFT on F..'
+ Status = DftiComputeBackward (My_Desc1_Handle, F_1D)
+ Status = DftiFreeDescriptor (My_Desc1_Handle)
+
+ Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 2, L_pad)
+ Status = DftiCommitDescriptor(My_Desc2_Handle)
+ !WRITE(*,*) 'Backward FFT on G..'
+ Status = DftiComputeBackward (My_Desc2_Handle, G_1D)
+ Status = DftiFreeDescriptor (My_Desc2_Handle)
+
+ !CALL backward_FFT( F_2D )
+ !CALL backward_FFT( G_2D )
+
+ !WRITE(*,*) 'C = F/(Mkr *Mkz) x G/(Mkr *Mkz)..'
+ DO ix=1,Mkr*Mkz
+ C_1D(ix) = F_1D(ix)/Mkr/Mkz * G_1D(ix)/Mkr/Mkz
+ ENDDO
+
+ !WRITE(*,*) 'Forward FFT on C..'
+ Status = DftiCreateDescriptor( My_Desc3_Handle, DFTI_DOUBLE,&
+ DFTI_COMPLEX, 2, L_pad)
+ Status = DftiCommitDescriptor(My_Desc3_Handle)
+ Status = DftiComputeForward (My_Desc3_Handle, C_1D)
+ Status = DftiFreeDescriptor (My_Desc3_Handle)
+
+ DO ix=1,Nkr
+ DO iy=1,Nkz
+ C_2D(ix,iy) = C_1D(Mkr*(iy-1) + ix)
+ ENDDO
+ ENDDO
+
+ END SUBROUTINE convolve_2D_F2F
+
+ !!! Convolution 2D Fourier to Real
+ ! - Same as convolve_2D_F2F but does not compute the forward FFT at the end
+ ! - Used to same computation in the compute_Sapj loop
+ ! - Output : C_2D that is defined in the real space
+ SUBROUTINE convolve_2D_F2R( F_2D, G_2D, C_2D )
+
+ USE prec_const
+ USE fourier_grid, ONLY : Nkr, Nkz, Pad
+ USE MKL_DFTI
+
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: F_2D, G_2D, C_2D
+ COMPLEX(dp), DIMENSION(Pad*Nkr*Pad*Nkz) :: F_1D, G_1D, C_1D
+
+ INTEGER :: ix, iy, Mkr, Mkz
+ INTEGER :: Status, L(2), L_pad(2)
+
+ type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle,My_Desc2_Handle
+
+ L = (/ Nkr, Nkz /)
+ Mkr = Pad*Nkr; Mkz = Pad*Nkz
+ L_pad = (/ Mkr, Mkz /)
+
+ !WRITE(*,*) 'Reshaping..'
+ DO ix = 1,Mkr
+ DO iy = 1,Mkz
+ IF ( (ix .LE. Nkr) .AND. (iy .LE. Nkz) ) THEN
+ F_1D(Mkr*(iy-1) + ix) = F_2D(ix,iy)
+ G_1D(Mkr*(iy-1) + ix) = G_2D(ix,iy)
+ ELSE
+ F_1D(Mkr*(iy-1) + ix) = 0._dp
+ G_1D(Mkr*(iy-1) + ix) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+
+ CALL backward_FFT( F_2D )
+ CALL backward_FFT( G_2D )
+
+ ! We adapt the 1D output vector to a 2D form
+ DO ix=1,Nkr
+ DO iy=1,Nkz
+ !C_2D(ix,iy) = F_1D(Mkr*(iy-1) + ix)/Mkr/Mkz * G_1D(Mkr*(iy-1) + ix)/Mkr/Mkz
+ C_2D(ix,iy) = F_2D(ix,iy) * G_2D(ix,iy)
+ ENDDO
+ ENDDO
+
+ END SUBROUTINE convolve_2D_F2R
+
+
+ !! Compute a forward FFT to go back to Fourier space
+ ! - used to save computation after the sum of convolution_2D_F2R in compute_Sapj
+ SUBROUTINE forward_FFT( C_2D )
+
+ USE prec_const
+ USE fourier_grid, ONLY : Nkr, Nkz, Pad
+ USE MKL_DFTI
+
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: C_2D
+ COMPLEX(dp), DIMENSION(Pad*Nkr*Pad*Nkz) :: C_1D
+
+ INTEGER :: ix, iy, Mkr, Mkz
+ INTEGER :: Status, L(2), L_pad(2)
+
+ type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle
+
+ L = (/ Nkr, Nkz /)
+ Mkr = Pad*Nkr; Mkz = Pad*Nkz
+ L_pad = (/ Mkr, Mkz /)
+
+ !WRITE(*,*) 'Reshaping..'
+ DO ix = 1,Mkr
+ DO iy = 1,Mkz
+ IF ( (ix .LE. Nkr) .AND. (iy .LE. Nkz) ) THEN
+ C_1D(Mkr*(iy-1) + ix) = C_2D(ix,iy)
+ ELSE
+ C_1D(Mkr*(iy-1) + ix) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+
+ Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_DOUBLE,&
+ DFTI_COMPLEX, 2, L_pad)
+ Status = DftiCommitDescriptor(My_Desc_Handle)
+ Status = DftiComputeForward (My_Desc_Handle, C_1D)
+ Status = DftiFreeDescriptor (My_Desc_Handle)
+
+ DO ix=1,Nkr
+ DO iy=1,Nkz
+ C_2D(ix,iy) = C_1D(Mkr*(iy-1) + ix)
+ ENDDO
+ ENDDO
+
+ END SUBROUTINE forward_FFT
+
+ !!! Convolution 2D Fourier to Fourier
+ ! - Compute the convolution using the convolution theorem and MKL
+ SUBROUTINE backward_FFT( C_2D )
+
+ USE prec_const
+ USE fourier_grid, ONLY : Nkr, Nkz, Pad
+ USE MKL_DFTI
+
+ IMPLICIT NONE
+
+ COMPLEX(dp), DIMENSION(Nkr,Nkz) :: C_2D
+ COMPLEX(dp), DIMENSION(Pad*Nkr*Pad*Nkz) :: C_1D
+
+ INTEGER :: ix, iy, Mkr, Mkz
+ INTEGER :: Status, L(2), L_pad(2)
+
+ type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle
+
+ L = (/ Nkr, Nkz /)
+ Mkr = Pad*Nkr; Mkz = Pad*Nkz
+ L_pad = (/ Mkr, Mkz /)
+
+ !WRITE(*,*) 'Reshaping..'
+ DO ix = 1,Mkr
+ DO iy = 1,Mkz
+ IF ( (ix .LE. Nkr) .AND. (iy .LE. Nkz) ) THEN
+ C_1D(Mkr*(iy-1) + ix) = C_2D(ix,iy)
+ ELSE
+ C_1D(Mkr*(iy-1) + ix) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+
+ Status = DftiCreateDescriptor(My_Desc_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 2, L_pad)
+ Status = DftiCommitDescriptor(My_Desc_Handle)
+ Status = DftiComputeBackward (My_Desc_Handle, C_1D)
+ Status = DftiFreeDescriptor (My_Desc_Handle)
+
+ DO ix=1,Nkr
+ DO iy=1,Nkz
+ C_2D(ix,iy) = C_1D(Mkr*(iy-1) + ix)/Mkr/Mkz
+ ENDDO
+ ENDDO
+
+ END SUBROUTINE backward_FFT
+
+END MODULE convolution
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 61f8bbd3..a5323291 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -1,8 +1,8 @@
SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
- USE basic
- USE fourier_grid
+ USE basic
+ USE fourier_grid
USE diagnostics_par
USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach
USE model
@@ -10,22 +10,21 @@ SUBROUTINE diagnose(kstep)
USE fields
USE time_integration
- use prec_const
+ USE prec_const
IMPLICIT NONE
-
+
INCLUDE 'srcinfo.h'
INTEGER, INTENT(in) :: kstep
- ! INTEGER, parameter :: BUFSIZE = 20
+ INTEGER, parameter :: BUFSIZE = 2
INTEGER :: rank, dims(1) = (/0/)
CHARACTER(len=256) :: str, fname
- !________________________________________________________________________________
- !
+ !_____________________________________________________________________________
! 1. Initial diagnostics
+
IF (kstep .EQ. 0) THEN
- WRITE(*,'(a)') ' Initial diagnostics'
-
+
! jobnum is either 0 from initialization or some integer values read from chkrst(0)
IF(jobnum .LE. 99) THEN
WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
@@ -43,26 +42,20 @@ SUBROUTINE diagnose(kstep)
WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
call flush(6)
-
+
! Data group
CALL creatg(fidres, "/data", "data")
- !CALL creatg(fidres, "/data/var0d", "0d history arrays")
- !CALL creatg(fidres, "/data/var1d", "1d profiles")
CALL creatg(fidres, "/data/var2d", "2d profiles")
CALL creatg(fidres, "/data/var5d", "5d profiles")
! Initialize counter of number of saves for each category
- !IF (cstep==0) THEN
- ! iframe1d=0
- !END IF
- !CALL attach(fidres,"/data/var1d/" , "frames", iframe1d)
- IF (cstep==0) THEN
+ IF (cstep==0) THEN
iframe2d=0
END IF
CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
- IF (cstep==0) THEN
+ IF (cstep==0) THEN
iframe5d=0
END IF
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
@@ -78,14 +71,14 @@ SUBROUTINE diagnose(kstep)
IF (write_Ni00) THEN
CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
- CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
IF (write_phi) THEN
CALL creatg(fidres, "/data/var2d/phi", "phi")
- CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
- CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
! var5d group (moments)
@@ -98,7 +91,7 @@ SUBROUTINE diagnose(kstep)
CALL putarr(fidres, "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e), "j_e",ionode=0)
CALL putarr(fidres, "/data/var5d/moments_e/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
CALL putarr(fidres, "/data/var5d/moments_e/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
-
+
CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
CALL putarr(fidres, "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i), "p_i",ionode=0)
CALL putarr(fidres, "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i), "j_i",ionode=0)
@@ -106,6 +99,20 @@ SUBROUTINE diagnose(kstep)
CALL putarr(fidres, "/data/var5d/moments_i/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
END IF
+ IF (write_non_lin) THEN
+ CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
+ CALL putarr(fidres, "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e), "j_e",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sepj/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sepj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+
+ CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
+ CALL putarr(fidres, "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i), "p_i",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i), "j_i",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sipj/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+ CALL putarr(fidres, "/data/var5d/Sipj/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+ END IF
+
! Add input namelist variables as attributes of /data/input, defined in srcinfo.h
WRITE(*,*) 'VERSION=', VERSION
WRITE(*,*) 'BRANCH=', BRANCH
@@ -141,7 +148,7 @@ SUBROUTINE diagnose(kstep)
CALL time_integration_outputinputs(fidres, str)
- ! Save STDIN (input file) of this run
+ ! Save STDIN (input file) of this run
IF(jobnum .LE. 99) THEN
WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
ELSE
@@ -150,13 +157,13 @@ SUBROUTINE diagnose(kstep)
INQUIRE(unit=lu_in, name=fname)
CLOSE(lu_in)
-
+
CALL putfile(fidres, TRIM(str), TRIM(fname),ionode=0)
END IF
- !________________________________________________________________________________
+ !_____________________________________________________________________________
! 2. Periodic diagnostics
!
IF (kstep .GE. 0) THEN
@@ -185,13 +192,18 @@ SUBROUTINE diagnose(kstep)
END IF
END IF
- !________________________________________________________________________________
- ! 3. Final diagnostics
+ ! 2.5 Backups
+ !
+ !
+ !
+
+ !_____________________________________________________________________________
+ ! 3. Final diagnostics
ELSEIF (kstep .EQ. -1) THEN
+ CALL cpu_time(finish)
! Close all diagnostic files
CALL closef(fidres)
-
END IF
END SUBROUTINE diagnose
@@ -202,18 +214,15 @@ SUBROUTINE diagnose_0d
USE basic
USE futils, ONLY: append
USE diagnostics_par
- use prec_const
+ USE prec_const
IMPLICIT NONE
WRITE(*,'(a,1x,i7.7,a1,i7.7,20x,a,1pe10.3,10x,a,1pe10.3)') &
'*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time, 'dt =', dt
WRITE(*,*)
-
- ! flush stdout of all ranks. Usually only rank 0 should write, but error messages might be written from other ranks as well
- CALL FLUSH(stdout)
- !CALL append(fidres,"/data/var0d/time" ,time, ionode=0)
- !CALL append(fidres,"/data/var0d/cstep" ,real(cstep,dp), ionode=0)
+ ! flush stdout of all ranks. Usually ONLY rank 0 should write, but error messages might be written from other ranks as well
+ CALL FLUSH(stdout)
END SUBROUTINE diagnose_0d
@@ -225,14 +234,14 @@ SUBROUTINE diagnose_2d
USE fields
USE time_integration
USE diagnostics_par
- use prec_const
+ USE prec_const
IMPLICIT NONE
- CALL append(fidres, "/data/var2d/time", time,ionode=0)
- CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var2d/", "frames",iframe2d)
+ CALL append(fidres, "/data/var2d/time", time,ionode=0)
+ CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var2d/", "frames",iframe2d)
iframe2d=iframe2d+1
- CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
+ CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
IF (write_phi) THEN
CALL write_field2d(phi(:,:), 'phi')
@@ -245,8 +254,8 @@ CONTAINS
SUBROUTINE write_field2d(field, text)
USE futils, ONLY: attach, putarr
- USE fourier_grid, only: ikrs,ikre, ikzs,ikze
- use prec_const
+ USE fourier_grid, ONLY: ikrs,ikre, ikzs,ikze
+ USE prec_const
IMPLICIT NONE
COMPLEX(dp), DIMENSION(ikrs:ikre, ikzs:ikze), INTENT(IN) :: field
@@ -258,69 +267,77 @@ CONTAINS
CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze),ionode=0)
CALL attach(fidres, dset_name, "time", time)
-
+
END SUBROUTINE write_field2d
END SUBROUTINE diagnose_2d
- SUBROUTINE diagnose_5d
+SUBROUTINE diagnose_5d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd
USE fields
- USE fourier_grid, only: ips_e,ipe_e, ips_i, ipe_i, &
+ USE array!, ONLY: Sepj, Sipj
+ USE fourier_grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
ijs_e,ije_e, ijs_i, ije_i
USE time_integration
USE diagnostics_par
- use prec_const
+ USE prec_const
+ USE model, ONLY: NON_LIN
IMPLICIT NONE
-
- CALL append(fidres, "/data/var5d/time", time,ionode=0)
- CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
+
+ CALL append(fidres, "/data/var5d/time", time,ionode=0)
+ CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
+ CALL getatt(fidres, "/data/var5d/", "frames",iframe5d)
iframe5d=iframe5d+1
- CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
-
+ CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+
IF (write_moments) THEN
CALL write_field5d_e(moments_e(:,:,:,:,updatetlevel), 'moments_e')
CALL write_field5d_i(moments_i(:,:,:,:,updatetlevel), 'moments_i')
END IF
-
+
+ IF (write_non_lin .AND. NON_LIN) THEN
+ CALL write_field5d_e(Sepj, 'Sepj')
+ CALL write_field5d_i(Sipj, 'Sipj')
+ END IF
+
CONTAINS
-
+
SUBROUTINE write_field5d_e(field, text)
- USE futils, ONLY: attach, putarr
- USE fourier_grid, only: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
+ USE futils, ONLY: attach, putarr
+ USE fourier_grid, ONLY: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
USE prec_const
IMPLICIT NONE
-
+
COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
-
+
CHARACTER(LEN=50) :: dset_name
-
+
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze),ionode=0)
-
+
CALL attach(fidres, dset_name, "time", time)
-
+
END SUBROUTINE write_field5d_e
SUBROUTINE write_field5d_i(field, text)
- USE futils, ONLY: attach, putarr
- USE fourier_grid, only: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
+ USE futils, ONLY: attach, putarr
+ USE fourier_grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
USE prec_const
IMPLICIT NONE
-
+
COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
-
+
CHARACTER(LEN=50) :: dset_name
-
+
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze),ionode=0)
-
+
CALL attach(fidres, dset_name, "time", time)
-
+
END SUBROUTINE write_field5d_i
- END SUBROUTINE diagnose_5d
\ No newline at end of file
+
+END SUBROUTINE diagnose_5d
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index cce34b04..c0b86ad4 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -5,9 +5,10 @@ MODULE diagnostics_par
IMPLICIT NONE
PRIVATE
- LOGICAL, PUBLIC, PROTECTED :: write_Ni00=.TRUE.
- LOGICAL, PUBLIC, PROTECTED :: write_moments=.TRUE.
+ LOGICAL, PUBLIC, PROTECTED :: write_Ni00=.TRUE.
+ LOGICAL, PUBLIC, PROTECTED :: write_moments=.TRUE.
LOGICAL, PUBLIC, PROTECTED :: write_phi=.TRUE.
+ LOGICAL, PUBLIC, PROTECTED :: write_non_lin=.TRUE.
LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d
@@ -17,9 +18,9 @@ MODULE diagnostics_par
CHARACTER(len=64), PUBLIC :: resfile0 = "results" ! Head of main result file name
CHARACTER(len=64), PUBLIC :: resfile ! Main result file
INTEGER, PUBLIC :: fidres ! FID for resfile
- ! CHARACTER(len=64), PUBLIC :: rstfile0 = "restart" ! Head of restart file name
- ! CHARACTER(len=64), PUBLIC :: rstfile ! Full restart file
- ! INTEGER, PUBLIC :: fidrst ! FID for restart file
+ CHARACTER(len=64), PUBLIC :: rstfile0 = "restart" ! Head of restart file name
+ CHARACTER(len=64), PUBLIC :: rstfile ! Full restart file
+ INTEGER, PUBLIC :: fidrst ! FID for restart file
PUBLIC :: output_par_readinputs, output_par_outputinputs
@@ -34,11 +35,10 @@ CONTAINS
IMPLICIT NONE
NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d
- NAMELIST /OUTPUT_PAR/ write_Ni00, write_moments, write_phi, write_doubleprecision
- NAMELIST /OUTPUT_PAR/ resfile0!, rstfile0
+ NAMELIST /OUTPUT_PAR/ write_Ni00, write_moments, write_phi, write_non_lin, write_doubleprecision
+ NAMELIST /OUTPUT_PAR/ resfile0, rstfile0
- READ(lu_in,output_par)
- WRITE(*,output_par)
+ READ(lu_in,output_par)
END SUBROUTINE output_par_readinputs
@@ -56,12 +56,13 @@ CONTAINS
CALL attach(fidres, TRIM(str), "write_Ni00", write_Ni00)
CALL attach(fidres, TRIM(str), "write_moments", write_moments)
CALL attach(fidres, TRIM(str), "write_phi", write_phi)
+ CALL attach(fidres, TRIM(str), "write_non_lin", write_non_lin)
CALL attach(fidres, TRIM(str), "write_doubleprecision", write_doubleprecision)
CALL attach(fidres, TRIM(str), "nsave_0d", nsave_0d)
CALL attach(fidres, TRIM(str), "nsave_1d", nsave_1d)
CALL attach(fidres, TRIM(str), "nsave_2d", nsave_2d)
CALL attach(fidres, TRIM(str), "nsave_5d", nsave_5d)
- CALL attach(fidres, TRIM(str), "resfile0", resfile0)
+ CALL attach(fidres, TRIM(str), "resfile0", resfile0)
END SUBROUTINE output_par_outputinputs
diff --git a/src/fields_mod.F90 b/src/fields_mod.F90
index 2e97e6b6..58075c2d 100644
--- a/src/fields_mod.F90
+++ b/src/fields_mod.F90
@@ -10,7 +10,6 @@ MODULE fields
!------------------ELECTROSTATIC POTENTIAL------------------
! Normalized electric potential: \hat{\phi} ! (kr,kz)
- COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi
END MODULE fields
-
diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90
index da83b23d..3cbcbc76 100644
--- a/src/fourier_grid_mod.F90
+++ b/src/fourier_grid_mod.F90
@@ -10,12 +10,14 @@ MODULE fourier_grid
INTEGER, PUBLIC, PROTECTED :: jmaxe = 2 ! The maximal electron Laguerre-moment computed
INTEGER, PUBLIC, PROTECTED :: pmaxi = 2 ! The maximal ion Hermite-moment computed
INTEGER, PUBLIC, PROTECTED :: jmaxi = 2 ! The maximal ion Laguerre-moment computed
+ INTEGER, PUBLIC, PROTECTED :: maxj = 2 ! The maximal ion Laguerre-moment computed
INTEGER, PUBLIC, PROTECTED :: nkr = 10 ! Number of total internal grid points in kr
REAL(dp), PUBLIC, PROTECTED :: krmin = 0._dp ! kr coordinate for left boundary
REAL(dp), PUBLIC, PROTECTED :: krmax = 1._dp ! kr coordinate for right boundary
INTEGER, PUBLIC, PROTECTED :: nkz = 10 ! Number of total internal grid points in kz
REAL(dp), PUBLIC, PROTECTED :: kzmin = 0._dp ! kz coordinate for left boundary
REAL(dp), PUBLIC, PROTECTED :: kzmax = 1._dp ! kz coordinate for right boundary
+ INTEGER, PUBLIC, PROTECTED :: Pad = 2 ! Zero padding parameter for convolution
! Indices of s -> start , e-> END
INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e
@@ -27,14 +29,14 @@ MODULE fourier_grid
REAL(dp), PUBLIC, PROTECTED :: deltakz
! Grids containing position in fourier space
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: kzarray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: krarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
! Grid containing the polynomials degrees
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray_e
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: parray_i
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jarray_e
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jarray_i
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i
! Public Functions
PUBLIC :: set_krgrid, set_kzgrid, set_pj
@@ -100,7 +102,7 @@ CONTAINS
ALLOCATE(jarray_i(ijs_i:ije_i))
DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
-
+ maxj = MAX(jmaxi, jmaxe)
END SUBROUTINE set_pj
SUBROUTINE fourier_grid_readinputs
@@ -112,10 +114,11 @@ CONTAINS
IMPLICIT NONE
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- nkr, krmin, krmax, nkz, kzmin, kzmax
+ nkr, krmin, krmax, nkz, kzmin, kzmax, &
+ Pad
READ(lu_in,grid)
- WRITE(*,grid)
+ !WRITE(*,grid)
END SUBROUTINE fourier_grid_readinputs
@@ -140,6 +143,7 @@ CONTAINS
CALL attach(fidres, TRIM(str), "nkz", nkz)
CALL attach(fidres, TRIM(str), "kzmin", kzmin)
CALL attach(fidres, TRIM(str), "kzmax", kzmax)
+ CALL attach(fidres, TRIM(str), "Pad", Pad)
END SUBROUTINE fourier_grid_outputinputs
diff --git a/src/inital.F90 b/src/inital.F90
index 07ff502d..d2e839c0 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,27 +1,40 @@
+!******************************************************************************!
+!!!!!! initialize the moments and load/build coeff tables
+!******************************************************************************!
SUBROUTINE inital
USE basic
- USE model, ONLY : CO
+ USE model, ONLY : CO, NON_LIN
USE prec_const
- implicit none
-
- WRITE(*,'(a/)') '=== Set initial conditions ==='
-
- CALL init_profiles
+ USE coeff
+ USE array, ONLY : Sepj,Sipj
+ implicit none
+ !!!!!! Set the moments arrays Nepj, Nipj !!!!!!
+ IF ( RESTART ) THEN
+ CALL init_moments
+ ELSE
+ CALL load_moments
+ ENDIF
+ !!!!!! Set phi !!!!!!
+ CALL poisson
+ !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
+ IF ( NON_LIN ) THEN;
+ CALL compute_Sapj(Sepj,Sipj)
+ CALL build_dnjs_table
+ ENDIF
+ !!!!!! Load the full coulomb collision operator coefficients !!!!!!
IF (CO .EQ. -1) THEN
- WRITE(*,'(a/)') '=== Load Full Coulomb matrix ==='
-
+ WRITE(*,*) '=== Load Full Coulomb matrix ==='
CALL load_FC_mat
ENDIF
- !
-
END SUBROUTINE inital
+!******************************************************************************!
-
-SUBROUTINE init_profiles
- ! Set initial conditions
-
+!******************************************************************************!
+!!!!!!! Initialize the moments randomly
+!******************************************************************************!
+SUBROUTINE init_moments
USE basic
USE fourier_grid
USE fields
@@ -40,34 +53,101 @@ SUBROUTINE init_profiles
CALL set_updatetlevel(1)
- DO ikr=ikrs,ikre
- DO ikz=ikzs,ikze
+ IF ( only_Na00 ) THEN ! Spike initialization on density only
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- CALL RANDOM_NUMBER(noise)
- !moments_e( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- END DO
+ DO ikr=ikrs,ikre
+ DO ikz=ikzs,ikze
+ moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
END DO
+ END DO
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- CALL RANDOM_NUMBER(noise)
- !moments_i( ip,ij, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ ELSE ! Broad noise initialization
+
+ DO ikr=ikrs,ikre
+ DO ikz=ikzs,ikze
+ DO ip=ips_e,ipe_e
+ DO ij=ijs_e,ije_e
+ CALL RANDOM_NUMBER(noise)
+ moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_e( ip,ij, Nkz-ikz, ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_e( ip,ij, ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_e( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ END DO
+ END DO
+ DO ip=ips_i,ipe_i
+ DO ij=ijs_i,ije_i
+ CALL RANDOM_NUMBER(noise)
+ moments_i( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ moments_i( ip,ij, Nkz-ikz, ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_i( ip,ij, ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ moments_i( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ END DO
END DO
END DO
-
- ! Poke initialization on only Ne00 and Ni00
- moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
-
END DO
- END DO
-
- CALL poisson ! To set phi
-
-END SUBROUTINE init_profiles
+ ENDIF
+END SUBROUTINE init_moments
+!******************************************************************************!
+
+!******************************************************************************!
+!!!!!!! Load moments from a previous save
+!******************************************************************************!
+SUBROUTINE load_moments
+! USE basic
+! USE initial_par
+! USE fields, ONLY : moments_e, moments_i
+! USE futils, ONLY : openf, getarr, closef
+! IMPLICIT NONE
+!
+! INTEGER :: fid
+!
+! CALL openf(backup_file,fid, 'r', 'D');
+! CALL getarr(fid, '/moments_e', moments_e) ! Nepj
+! CALL getarr(fid, '/moments_i', moments_i) ! Nipj
+! CALL getarr(fid, '/time', time) ! time
+! CALL closef(fid)
+END SUBROUTINE
+!******************************************************************************!
+
+!******************************************************************************!
+!!!!!!! Build the Laguerre-Laguerre coupling coefficient table
+!******************************************************************************!
+SUBROUTINE build_dnjs_table
+ USE basic
+ USE array, Only : dnjs
+ USE fourier_grid, Only : jmaxe, jmaxi
+ USE coeff
+ IMPLICIT NONE
+ INTEGER :: in, ij, is, J
+ INTEGER :: n_, j_, s_
+
+ J = max(jmaxe,jmaxi)
+
+ DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
+ n_ = in - 1
+ DO ij = in,J+1
+ j_ = ij - 1
+ DO is = ij,J+1
+ s_ = is - 1
+
+ dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
+ ! By symmetry
+ dnjs(in,is,ij) = dnjs(in,ij,is)
+ dnjs(ij,in,is) = dnjs(in,ij,is)
+ dnjs(ij,is,in) = dnjs(in,ij,is)
+ dnjs(is,ij,in) = dnjs(in,ij,is)
+ dnjs(is,in,ij) = dnjs(in,ij,is)
+ ENDDO
+ ENDDO
+ ENDDO
+END SUBROUTINE
+!******************************************************************************!
+
+!******************************************************************************!
+!!!!!!! Load the Full coulomb coefficient table from COSOlver results
+!******************************************************************************!
SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
USE basic
USE fourier_grid
@@ -80,34 +160,18 @@ SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,
INTEGER :: ip,ij, ip2,ij2
INTEGER :: fid1, fid2, fid3, fid4
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- WRITE(*,*) 'Load self electron collision matrix..'
! get the self electron colision matrix
CALL openf(selfmat_file,fid1, 'r', 'D');
CALL getarr(fid1, '/Caapj/Ceepj', Ceepj) ! get array (moli format)
CALL closef(fid1)
-
! get the Test and Back field electron ion collision matrix
- WRITE(*,*) 'Load test + field electron collision matrix..'
CALL openf(eimat_file,fid2, 'r', 'D');
CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT)
CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF)
CALL closef(fid2)
- write(25,*) eimat_file
- write(25,*) 'Ceepj(3,5)'
- write(25,*) Ceepj(3,5)
- write(25,*) 'CeipjT(12,12)'
- write(25,*) CeipjT(12,12)
- write(25,*) 'CeipjF(12,2)'
- write(25,*) CeipjF(12,2)
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- WRITE(*,*) 'Load self ion collision matrix..'
! get the self electron colision matrix
CALL openf(selfmat_file, fid3, 'r', 'D');
IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
@@ -117,18 +181,9 @@ SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,
ENDIF
CALL closef(fid3)
! get the Test and Back field ion electron collision matrix
- WRITE(*,*) 'Load test + field ion collision matrix..'
CALL openf(iemat_file, fid4, 'r', 'D');
CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT)
CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF)
CALL closef(fid4)
-
- write(25,*) iemat_file
- write(25,*) 'Ciipj(3,5)'
- write(25,*) Ciipj(3,5)
- write(25,*) 'CiepjT(10,10)'
- write(25,*) CiepjT(10,10)
- write(25,*) 'CiepjF(10,2)'
- write(25,*) CiepjF(10,2)
-
END SUBROUTINE load_FC_mat
+!******************************************************************************!
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index 61dbfdb7..50c96ab2 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -5,20 +5,23 @@ MODULE initial_par
IMPLICIT NONE
PRIVATE
+ ! To restart the simulation or look for saved state
+ LOGICAL, PUBLIC, PROTECTED :: RESTART=.true.
+ CHARACTER(len=128), PUBLIC :: backup_file
! Initial background level
- REAL(dp), PUBLIC, PROTECTED :: initback_moments=0._dp
-
+ REAL(dp), PUBLIC, PROTECTED :: initback_moments=0._dp
+ ! Initial background level
+ LOGICAL, PUBLIC, PROTECTED :: only_Na00 = .false.
! Initial background noise amplitude
- REAL(dp), PUBLIC, PROTECTED :: initnoise_moments=1E-6_dp
-
+ REAL(dp), PUBLIC, PROTECTED :: initnoise_moments=1E-6_dp
! Initialization for random number generator
- INTEGER, PUBLIC, PROTECTED :: iseed=42
+ INTEGER, PUBLIC, PROTECTED :: iseed=42
! Parameters of initial smooth sine profiles
REAL(dp), PUBLIC, PROTECTED :: init_nb_oscil_density=2._dp ! Number of oscillations
- REAL(dp), PUBLIC, PROTECTED :: init_nb_oscil_temp=2._dp
+ REAL(dp), PUBLIC, PROTECTED :: init_nb_oscil_temp=2._dp
REAL(dp), PUBLIC, PROTECTED :: init_ampli_density=0.1_dp ! Oscillation amplitude
- REAL(dp), PUBLIC, PROTECTED :: init_ampli_temp=0.1_dp
+ REAL(dp), PUBLIC, PROTECTED :: init_ampli_temp=0.1_dp
CHARACTER(len=128), PUBLIC :: selfmat_file ! COSOlver matrix file names
CHARACTER(len=128), PUBLIC :: iemat_file ! COSOlver matrix file names
@@ -27,15 +30,18 @@ MODULE initial_par
PUBLIC :: initial_outputinputs, initial_readinputs
CONTAINS
-
-
+
+
SUBROUTINE initial_readinputs
! Read the input parameters
- USE basic, ONLY : lu_in
+ USE basic, ONLY : lu_in, RESTART
USE prec_const
IMPLICIT NONE
+ NAMELIST /INITIAL_CON/ RESTART
+ NAMELIST /INITIAL_CON/ backup_file
+ NAMELIST /INITIAL_CON/ only_Na00
NAMELIST /INITIAL_CON/ initback_moments
NAMELIST /INITIAL_CON/ initnoise_moments
NAMELIST /INITIAL_CON/ iseed
@@ -44,7 +50,7 @@ CONTAINS
NAMELIST /INITIAL_CON/ eimat_file
READ(lu_in,initial_con)
- WRITE(*,initial_con)
+ !WRITE(*,initial_con)
END SUBROUTINE initial_readinputs
@@ -58,6 +64,8 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
+ CALL attach(fidres, TRIM(str), "only_Na00", only_Na00)
+
CALL attach(fidres, TRIM(str), "initback_moments", initback_moments)
CALL attach(fidres, TRIM(str), "initnoise_moments", initnoise_moments)
diff --git a/src/memory.F90 b/src/memory.F90
index b695a128..0e1e380b 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -6,7 +6,7 @@ SUBROUTINE memory
USE fields
USE fourier_grid
USE time_integration
- USE model, ONLY: CO
+ USE model, ONLY: CO, NON_LIN
USE prec_const
IMPLICIT NONE
@@ -30,4 +30,12 @@ SUBROUTINE memory
CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1))
CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1))
ENDIF
+
+ ! Non linear terms and dnjs table
+ IF ( NON_LIN ) THEN
+ CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze )
+ CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze )
+ CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
+ ENDIF
+
END SUBROUTINE memory
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index f30447b0..230afb3d 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -5,7 +5,8 @@ MODULE model
IMPLICIT NONE
PRIVATE
- INTEGER(dp), PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
+ INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
+ LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term
REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature
REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp !
@@ -29,11 +30,11 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CO, nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CO, NON_LIN, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, eta_n, eta_T, eta_B, lambdaD
READ(lu_in,model_par)
- WRITE(*,model_par)
+ !WRITE(*,model_par)
! Collision Frequency Normalization ... to match fluid limit
nu = nu*0.532_dp
@@ -51,6 +52,7 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
CALL attach(fidres, TRIM(str), "CO", CO)
+ CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN)
CALL attach(fidres, TRIM(str), "nu", nu)
CALL attach(fidres, TRIM(str), "tau_e", tau_e)
CALL attach(fidres, TRIM(str), "tau_i", tau_i)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index c3c3a9f2..6aa516e4 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -40,9 +40,9 @@ SUBROUTINE moments_eq_rhs
ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxe)
! N_e^{p+2,j} coeff
- xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ xNapp2j = taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
! N_e^{p-2,j} coeff
- xNapm2j = -taue_qe_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
+ xNapm2j = taue_qe_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
factj = 1.0 ! Start of the recursive factorial
@@ -50,11 +50,11 @@ SUBROUTINE moments_eq_rhs
ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxe)
! N_e^{p,j+1} coeff
- xNapjp1 = taue_qe_etaB * (ij_dp + 1._dp)
+ xNapjp1 = -taue_qe_etaB * (ij_dp + 1._dp)
! N_e^{p,j-1} coeff
- xNapjm1 = taue_qe_etaB * ij_dp
+ xNapjm1 = -taue_qe_etaB * ij_dp
! N_e^{pj} coeff
- xNapj = -taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+ xNapj = taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
!! Collision operator pj terms
xCapj = -nu_e*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
@@ -187,11 +187,17 @@ SUBROUTINE moments_eq_rhs
Tphi = 0._dp
ENDIF
- ! Sum of all precomputed terms
+ ! Sum of all linear terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
- -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi)&
+ -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+ TColl
+ ! Adding non linearity
+ IF ( NON_LIN ) THEN
+ moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
+ moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)
+ ENDIF
+
END DO kzloope
END DO krloope
@@ -205,9 +211,9 @@ SUBROUTINE moments_eq_rhs
ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxi)
! x N_i^{p+2,j} coeff
- xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
+ xNapp2j = taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
! x N_i^{p-2,j} coeff
- xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
+ xNapm2j = taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
factj = 1._dp ! Start of the recursive factorial
@@ -215,11 +221,11 @@ SUBROUTINE moments_eq_rhs
ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxi)
! x N_i^{p,j+1} coeff
- xNapjp1 = taui_qi_etaB * (ij_dp + 1._dp)
+ xNapjp1 = -taui_qi_etaB * (ij_dp + 1._dp)
! x N_i^{p,j-1} coeff
- xNapjm1 = taui_qi_etaB * ij_dp
+ xNapjm1 = -taui_qi_etaB * ij_dp
! x N_i^{pj} coeff
- xNapj = -taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
+ xNapj = taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
!! Collision operator pj terms
xCapj = -nu_i*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
@@ -338,8 +344,6 @@ SUBROUTINE moments_eq_rhs
ENDDO jloopie
ENDDO ploopie
- write(25,*) TColl
-
ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein
TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
ENDIF
@@ -354,11 +358,16 @@ SUBROUTINE moments_eq_rhs
Tphi = 0._dp
ENDIF
- ! Sum of all precomputed terms
+ ! Sum of linear terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
- -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi)&
+ -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+ TColl
+ ! Adding non linearity
+ IF ( NON_LIN ) THEN
+ moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
+ moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)
+ ENDIF
END DO kzloopi
END DO krloopi
diff --git a/src/poisson.F90 b/src/poisson.F90
index 1e8eb49d..ee353c0c 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -85,7 +85,12 @@ SUBROUTINE poisson
+ qi2_taui * (1._dp - sum_kernel2_i)
gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
- phi(ikr, ikz) = gammaD_phi/gammaD
+ IF ( (alphaD .EQ. 0._dp) .AND. (gammaD .EQ. 0._dp) ) THEN
+ write(*,*) "Warning : 0/0 occuring"
+ phi(ikr,ikz) = 0._dp
+ ELSE
+ phi(ikr, ikz) = gammaD_phi/gammaD
+ ENDIF
END DO
END DO
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index 645e6641..a7257cc0 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -40,7 +40,7 @@ MODULE prec_const
REAL(dp), PARAMETER :: onesixteenth=0.0625_dp
REAL(dp), PARAMETER :: ninesixteenths=0.5625_dp
REAL(dp), PARAMETER :: thirteentwelfths = 1.083333333333333333333333333333333333333333333_dp
- COMPLEX(dp), PARAMETER :: imagu = (0,-1)
+ COMPLEX(dp), PARAMETER :: imagu = (0,1)
CONTAINS
SUBROUTINE INIT_PREC_CONST
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index cdcb0cc4..879df9c1 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -10,7 +10,7 @@ SUBROUTINE readinputs
USE prec_const
IMPLICIT NONE
- WRITE(*,'(a/)') '=== Define additional data for a new run ==='
+ !WRITE(*,'(a/)') '=== Define additional data for a new run ==='
! The input file must be in the same order as the following read routines commands (eg spacegrid then output then model)
diff --git a/src/stepon.F90 b/src/stepon.F90
index bbd75a2b..678b81e0 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -1,49 +1,46 @@
SUBROUTINE stepon
! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
- USE basic
+ USE basic
USE time_integration
USE fields, ONLY: moments_e, moments_i, phi
- USE array , ONLY: moments_rhs_e, moments_rhs_i
+ USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
USE fourier_grid
USE advance_field_routine, ONLY: advance_time_level, advance_field
USE model
USE utility, ONLY: checkfield
-
use prec_const
IMPLICIT NONE
- INTEGER :: num_step, ip,ij
+ INTEGER :: num_step, ip,ij, ikr, ikz
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
- ! Compute right hand side
- !WRITE (*,*) 'Compute right hand side .. nstep = ', num_step
+ ! Compute right hand side of moments hierarchy equation
CALL moments_eq_rhs
-
- !WRITE (*,*) 'Advance time level .. nstep = ', num_step
- CALL advance_time_level ! Advance from updatetlevel to updatetlevel+1
-
- !WRITE (*,*) 'Advance field .. nstep = ', num_step
-
+ ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
+ CALL advance_time_level
+ ! Update the moments with the hierarchy RHS (step by step)
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
CALL advance_field(moments_e(ip,ij,:,:,:),moments_rhs_e(ip,ij,:,:,:))
ENDDO
ENDDO
-
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
CALL advance_field(moments_i(ip,ij,:,:,:),moments_rhs_i(ip,ij,:,:,:))
ENDDO
ENDDO
-
- ! Solving Poisson equation
- !WRITE (*,*) 'Solving Poisson equation .. nstep = ', num_step
+ ! Update electrostatic potential
CALL poisson
- CALL checkfield_all()
+ ! Update nonlinear term
+ IF ( NON_LIN ) THEN
+ CALL compute_Sapj
+ ENDIF
+ !(The two routines above are called in inital for t=0)
+ CALL checkfield_all()
END DO
CONTAINS
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 8089fe44..4936e6f9 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -32,7 +32,7 @@ CONTAINS
NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme
READ(lu_in,time_integration_par)
- WRITE(*,time_integration_par)
+ !WRITE(*,time_integration_par)
CALL set_numerical_scheme
diff --git a/wk/Test_Non_Lin.m b/wk/Test_Non_Lin.m
index 7f708721..afbec08c 100644
--- a/wk/Test_Non_Lin.m
+++ b/wk/Test_Non_Lin.m
@@ -11,22 +11,24 @@ OUTPUTS.nsave_5d = 1;
OUTPUTS.write_Ni00 = '.true.';
OUTPUTS.write_moments = '.true.';
OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_non_lin = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 2;
-GRID.jmaxe = 2;
-GRID.pmaxi = 2;
-GRID.jmaxi = 2;
-GRID.nkr = 10;
-GRID.krmin = 0.0;
-GRID.krmax = 0.;
-GRID.nkz = 10;
-GRID.kzmin = 0.1;
+GRID.pmaxe = 1;
+GRID.jmaxe = 1;
+GRID.pmaxi = 1;
+GRID.jmaxi = 1;
+GRID.nkr = 2;
+GRID.krmin =-1.0;
+GRID.krmax = 1.0;
+GRID.nkz = 2;
+GRID.kzmin =-1.0;
GRID.kzmax = 1.0;
%% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
+MODEL.NON_LIN = 1; % Non linear term
+MODEL.nu = 0.001; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
@@ -44,9 +46,9 @@ MODEL.eta_B = 0.5; % Magnetic
MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
-BASIC.nrun = 100000;
-BASIC.dt = 0.05;
-BASIC.tmax = 1.0;
+BASIC.nrun = 1;
+BASIC.dt = 0.1;
+BASIC.tmax = 0.1;
INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42;
@@ -66,91 +68,46 @@ INITIAL.iemat_file = ...
%% Write input file and run HeLaZ
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
+MAKE = 'cd ..; make; cd wk';
+system(MAKE)
EXEC = ' ../bin/helaz ';
RUN = ['mpirun -np ' num2str(nproc)];
CMD = [RUN, EXEC, INPUT];
-system(CMD);
+system(CMD); % RUN HeLaZ
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Run MOLI for each grid point and save the moments evolution
-kr_scan = linspace(GRID.krmin,GRID.krmax,GRID.nkr);
-kz_scan = linspace(GRID.kzmin,GRID.kzmax,GRID.nkz);
-Nmtot = (GRID.pmaxe+1) * (GRID.jmaxe+1) + (GRID.pmaxi+1) * (GRID.jmaxi+1);
-Nt = floor(BASIC.tmax/BASIC.dt)+1;
-Napj_2D = zeros(GRID.nkr,GRID.nkz,Nt,Nmtot);
-params.RK4 = 1;
-for ikr = 1:GRID.nkr
- for ikz = 1:GRID.nkz
- kr = kr_scan(ikr); % used by MOLI_time_solver_2D directly
- kz = kz_scan(ikz);
- run ../matlab/MOLI_time_solver_2D
- Napj_2D(ikr,ikz,:,:) = results.Napj;
+%% Load results
+filename = 'results_00.h5';
+[Nepj, pi, ji, kr, kz, time] = load_5D_data(filename,'moments_e');
+[Nipj, pe, je, kr, kz, time] = load_5D_data(filename,'moments_i');
+[Sepj, pe, je, kr, kz, time] = load_5D_data(filename,'Sepj');
+[Sipj, pi, ji, kr, kz, time] = load_5D_data(filename,'Sipj');
+[phi, kr, kz, time] = load_2D_data(filename,'phi');
+
+%% Compute non linear term with a Matlab method
+it = 1; p = 0; j = 0;
+for p = 0:1
+ for j = 0:1
+Sepj_mat = compute_Sapj(p, j, kr, kz, Nepj(:,:,:,:,end), 'e', phi(:,:,end), MODEL, GRID);
end
end
-%% Compute non linear term
-
+Sepj_for = squeeze(Sepj(p+1,j+1,:,:,it+1));
+% Comparison between Matlab and Fortran method
+disp(mean(mean(Sepj_mat - Sepj_for)))
+disp(mean(mean(Sepj_for)))
+disp(mean(mean(Sepj_mat)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Analysis and basic figures
default_plots_options
SAVEFIG = 1;
filename = 'results_00.h5';
default_plots_options
-TITLE = [TITLE,', $k_z=',num2str(GRID.kzmin),'$'];
-if params.RK4; MOLIsolvername = 'RK4';
-else; MOLIsolvername = 'ode15i';
-end
-bare = @(p_,j_) (GRID.jmaxe+1)*p_ + j_ + 1;
-bari = @(p_,j_) bare(GRID.pmaxe, GRID.jmaxe) + (GRID.jmaxi+1)*p_ + j_ + 1;
-
-%% Convert MOLI moments into HeLaZ format
-Nepj_MOLI = zeros(GRID.pmaxe+1, GRID.jmaxe+1,GRID.nkr,GRID.nkz,Nt);
-for ip = 1:GRID.pmaxe+1
- for ij = 1:GRID.jmaxe+1
- Nepj_MOLI(ip,ij,:,:,:) = Napj_2D(:,:,:,bare(ip-1,ij-1));
- end
-end
-Nipj_MOLI = zeros(GRID.pmaxi+1, GRID.jmaxi+1,GRID.nkr,GRID.nkz,Nt);
-for ip = 1:GRID.pmaxi+1
- for ij = 1:GRID.jmaxi+1
- Nipj_MOLI(ip,ij,:,:,:) = Napj_2D(:,:,:,bari(ip-1,ij-1));
- end
-end
-%% Load HeLaZ moments
-
-moment = 'moments_i';
-
-kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
-time = h5read(filename,'/data/var5d/time');
-Nipj_HeLaZ = zeros(GRID.pmaxi+1, GRID.jmaxi+1,numel(kr),numel(kz),numel(time));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
- Nipj_HeLaZ(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
-end
-
-moment = 'moments_e';
-
-kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
-time = h5read(filename,'/data/var5d/time');
-Nepj_HeLaZ = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
- Nepj_HeLaZ(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
-end
-
-%% Check coherence of linear results
-disp('MOLI - HeLaZ differences :');
-disp(sum(sum(sum(sum(sum(Nipj_MOLI-Nipj_HeLaZ)))))...
- / sum(sum(sum(sum(sum(Nipj_MOLI))))));
-
-%% Non linear term computation
-%get phi
-phiHeLaZ = zeros(numel(timephi),numel(kr),numel(kz));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
- phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
-end
-
-compute_Sapj
-
+TITLE = [TITLE,', $Sapj=',num2str(GRID.kzmin),'$'];
+fig = figure;
+subplot(121)
+plot(kz,Sepj_mat)
+hold on
+plot(kz,Sepj_for,'--')
+subplot(122)
+plot(time,squeeze(Sepj(p+1,j+1,1,1,:)))
+FIGNAME = 'out';
+SIMID = BASIC.SIMID; save_figure
\ No newline at end of file
diff --git a/wk/benchmark_kperp_scan.m b/wk/benchmark_kperp_scan.m
index b5432bf4..12658402 100644
--- a/wk/benchmark_kperp_scan.m
+++ b/wk/benchmark_kperp_scan.m
@@ -9,6 +9,7 @@ OUTPUTS.nsave_5d = 0;
OUTPUTS.write_Ni00 = '.true.';
OUTPUTS.write_moments = '.true.';
OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_non_lin = '.false.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
@@ -24,6 +25,7 @@ GRID.kzmin = 0.1;
GRID.kzmax = 1.5;
%% Model parameters
MODEL.CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.NON_LIN = 0; % Non linear term
MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
diff --git a/wk/benchmark_time_solver.m b/wk/benchmark_time_solver.m
index 641f0f0c..34483cef 100644
--- a/wk/benchmark_time_solver.m
+++ b/wk/benchmark_time_solver.m
@@ -1,31 +1,36 @@
+% Run linear/nonlin simulation on a kr,kz grid and compare with linear MOLI
clear all; close all;
-SIMID = 'test_full_coulomb'; % Name of the simulations
+SIMID = 'non_lin_benchmark_time_solver'; % Name of the simulations
addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% outputs options
OUTPUTS.nsave_0d = 0;
OUTPUTS.nsave_1d = 0;
-OUTPUTS.nsave_2d = 1;
-OUTPUTS.nsave_5d = 1;
-OUTPUTS.write_Ni00 = '.true.';
+OUTPUTS.nsave_2d = 20;
+OUTPUTS.nsave_5d = 20;
+OUTPUTS.write_Ni00 = '.false.';
OUTPUTS.write_moments = '.true.';
OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_non_lin = '.false.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
+OUTPUTS.rstfile0 = '''restart''';
%% Grid parameters
-GRID.pmaxe = 15;
-GRID.jmaxe = 6;
-GRID.pmaxi = 15;
-GRID.jmaxi = 6;
-GRID.nkr = 1;
-GRID.krmin = 0.;
-GRID.krmax = 0.;
-GRID.nkz = 1;
-GRID.kzmin = 1.0;
-GRID.kzmax = 1.0;
+GRID.pmaxe = 8;
+GRID.jmaxe = 4;
+GRID.pmaxi = 8;
+GRID.jmaxi = 4;
+GRID.nkr = 8;
+GRID.krmin =-2.0;
+GRID.krmax = 2.0;
+GRID.nkz = 8;
+GRID.kzmin =-2.0;
+GRID.kzmax = 2.0;
+GRID.Pad = 2.0;
%% Model parameters
-MODEL.CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-MODEL.nu = 1.0; % collisionality nu*L_perp/Cs0
+MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.NON_LIN = '.true.'; % Non linear term
+MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0 (~10^2 bigger than Ricci 2006)
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
@@ -36,18 +41,21 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
% gradients L_perp/L_x
-MODEL.eta_n = 0.0; % Density
+MODEL.eta_n = 1.0; % Density
MODEL.eta_T = 0.0; % Temperature
-MODEL.eta_B = 0.0; % Magnetic
+MODEL.eta_B = 0.5; % Magnetic
% Debye length
MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000;
-BASIC.dt = 0.05;
-BASIC.tmax = 10.0;
-INITIAL.initback_moments = 0.01;
-INITIAL.initnoise_moments = 0.;
+BASIC.dt = 0.01;
+BASIC.tmax = 50.0;
+INITIAL.RESTART = '.true.';
+INITIAL.backup_file = '''restart''';
+INITIAL.only_Na00 = '.true.';
+INITIAL.initback_moments = 1.0e-2;
+INITIAL.initnoise_moments = 0.0e-5;
INITIAL.iseed = 42;
INITIAL.selfmat_file = ...
['''../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(GRID.pmaxe),...
@@ -61,20 +69,30 @@ INITIAL.iemat_file = ...
['''../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
'_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Write input file
-INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Run HeLaZ
+INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
+MAKE = 'cd ..; make; cd wk';
+system(MAKE);
EXEC = ' ../bin/helaz ';
RUN = ['mpirun -np ' num2str(nproc)];
CMD = [RUN, EXEC, INPUT];
system(CMD);
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Run MOLI with same input parameters
+% Load results
+filename = [OUTPUTS.resfile0(2:end-1),'_00.h5'];
+[ Nipj, pgrid_i, jgrid_i, kr, kz, time5d ] = load_5D_data( filename, 'moments_i' );
+[ Nepj, pgrid_e, jgrid_e, ~, ~, ~ ] = load_5D_data( filename, 'moments_e' );
+[ phiHeLaZ ,~, ~, time2d ] = load_2D_data( filename, 'phi' );
+if strcmp(OUTPUTS.write_non_lin,'.true.') && strcmp(MODEL.NON_LIN,'.true.')
+[ Sepj, ~, ~, ~, ~, ~ ] = load_5D_data( filename, 'Sepj' );
+[ Sipj, ~, ~, ~, ~, ~ ] = load_5D_data( filename, 'Sipj' );
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%^
+%% Run MOLI with same input parameters and a given position on k grid
+%ikr = ceil(GRID.nkr/2); ikz = ceil(GRID.nkz/2);
+ikr = GRID.nkr; ikz = GRID.nkz;
+kr_MOLI = kr(ikr); kz_MOLI = kz(ikz);
params.RK4 = 1;
run ../matlab/MOLI_time_solver
if params.RK4; MOLIsolvername = 'RK4';
@@ -84,179 +102,163 @@ end
%% Analysis and basic figures
default_plots_options
SAVEFIG = 1;
-filename = 'results_00.h5';
default_plots_options
-TITLE = [TITLE,', $k_z=',num2str(GRID.kzmin),'$'];
-
bare = @(p_,j_) (GRID.jmaxe+1)*p_ + j_ + 1;
bari = @(p_,j_) bare(GRID.pmaxe, GRID.jmaxe) + (GRID.jmaxi+1)*p_ + j_ + 1;
-%% Load moments
-
-moment = 'moments_i';
-
-kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
-time = h5read(filename,'/data/var5d/time');
-Nipj = zeros(GRID.pmaxi+1, GRID.jmaxi+1,numel(kr),numel(kz),numel(time));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
- Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
-end
-
-moment = 'moments_e';
-
-kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
-time = h5read(filename,'/data/var5d/time');
-Nepj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
- Nepj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
-end
-
-
%% Plot moments
fig = figure;
-
-x1 = time;
-x2 = results.time;
-ic = 1;
-
-% Electrons
-subplot(321)
-for ip = 0:1
- for ij = 0:1
- y1 = squeeze(real(Nepj(ip+1,ij+1,1,1,:)));
- plot(x1,y1,'-','DisplayName',...
- ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
- 'color', line_colors(ic,:))
- hold on
- y2 = squeeze(real(results.Napj(:,bare(ip,ij))));
- plot(x2,y2,'--','DisplayName',...
- ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
- 'color', line_colors(ic,:))
- hold on
- ic = ic + 1;
+ tH = time5d; tM = results.time;
+ % Electrons Real%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ subplot(331); ic = 1;
+ for ip = 0:1
+ for ij = 0:1
+ yH = squeeze(real(Nepj(ip+1,ij+1,ikr,ikz,:)));
+ plot(tH,yH,'-','DisplayName',...
+ ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ yM = squeeze(real(results.Napj(:,bare(ip,ij))));
+ plot(tM,yM,'--','DisplayName',...
+ ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', 0.7*line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+ end
+ grid on
+ xlabel('$t$'); ylabel(['Re$(N_e^{pj})$'])
+ % Electrons Imag%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ subplot(334); ic = 1;
+ for ip = 0:1
+ for ij = 0:1
+ yH = squeeze(imag(Nepj(ip+1,ij+1,ikr,ikz,:)));
+ plot(tH,yH,'-','DisplayName',...
+ ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ yM = squeeze(imag(results.Napj(:,bare(ip,ij))));
+ plot(tM,yM,'--','DisplayName',...
+ ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', 0.7*line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+ end
+ grid on
+ xlabel('$t$'); ylabel(['Im$(N_e^{pj})$'])
+ % Error on abs(Ne00)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ yH = squeeze(abs(Nepj(1,1,ikr,ikz,:))); % HeLaZ
+ yM = squeeze(abs(results.Napj(:,bare(0,0))));
+ ERR3 = abs(yM(1:OUTPUTS.nsave_5d:numel(yH)*OUTPUTS.nsave_5d)-yH);
+ subplot(337);
+ yyaxis left
+ semilogy(tH,yH,'-','DisplayName','HeLaZ',...
+ 'color', line_colors(1,:)); hold on;
+ semilogy(tM,yM,'--','DisplayName','MOLI',...
+ 'color', 'k')
+ grid on
+ xlabel('$t$')
+ ylabel('$|N_e^{00}|$')
+ legend('show')
+ yyaxis right
+ semilogy(tH, ERR3,'color', line_colors(2,:));
+ set(gca, 'YScale', 'log')
+ ylabel('$e(N_i^{00})$')
+ % Ions Real%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ subplot(332); ic = 1;
+ for ip = 0:1
+ for ij = 0:1
+ yH = squeeze(real(Nipj(ip+1,ij+1,ikr,ikz,:)));
+ plot(tH,yH,'-','DisplayName',...
+ ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ yM = squeeze(real(results.Napj(:,bari(ip,ij))));
+ plot(tM,yM,'--','DisplayName',...
+ ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', 0.7*line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+ end
+ grid on
+ xlabel('$t$'); ylabel(['Re$(N_i^{pj})$'])
+ % Ions Imag%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ subplot(335); ic = 1;
+ for ip = 0:1
+ for ij = 0:1
+ yH = squeeze(imag(Nipj(ip+1,ij+1,ikr,ikz,:)));
+ plot(tH,yH,'-','DisplayName',...
+ ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', line_colors(ic,:))
+ hold on
+ yM = squeeze(imag(results.Napj(:,bari(ip,ij))));
+ plot(tM,yM,'--','DisplayName',...
+ ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ 'color', 0.7*line_colors(ic,:))
+ hold on
+ ic = ic + 1;
+ end
+ end
+ grid on
+ xlabel('$t$'); ylabel(['Im$(N_i^{pj})$'])
+ % Error on abs(Ni00)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ yH = squeeze(abs(Nipj(1,1,ikr,ikz,:))); % HeLaZ
+ yM = squeeze(abs(results.Napj(:,bari(0,0))));
+ ERR3 = abs(yM(1:OUTPUTS.nsave_5d:numel(yH)*OUTPUTS.nsave_5d)-yH);
+ subplot(338);
+ yyaxis left
+ semilogy(tH,yH,'-','DisplayName','HeLaZ',...
+ 'color', line_colors(1,:)); hold on;
+ semilogy(tM,yM,'--','DisplayName','MOLI',...
+ 'color', 'k')
+ grid on
+ xlabel('$t$')
+ ylabel('$|N_i^{00}|$')
+ legend('show')
+ yyaxis right
+ semilogy(tH, ERR3,'color', line_colors(2,:));
+ set(gca, 'YScale', 'log')
+ ylabel('$e(N_i^{00})$')
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% Phi
+ timephiMOLI = results.time;
+ phiMOLI = zeros(size(timephiMOLI));
+ phiHeLaZ_rz = squeeze(phiHeLaZ(ikr,ikz,:));
+ for it = 1:numel(timephiMOLI)
+ phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
end
-end
-grid on
-xlabel('$t$')
-ylabel(['Re$(N_e^{pj})$'])
-
-% Ions
-ic = 1;
-subplot(322)
-for ip = 0:1
- for ij = 0:1
- y1 = squeeze(real(Nipj(ip+1,ij+1,1,1,:)));
- plot(x1,y1,'-','DisplayName',...
- ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
- 'color', line_colors(ic,:))
+ subplot(333)
+ plot(time2d,real(phiHeLaZ_rz),'-','DisplayName','HeLaZ RK4')
hold on
- y2 = squeeze(real(results.Napj(:,bari(ip,ij))));
- plot(x2,y2,'--','DisplayName',...
- ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
- 'color', line_colors(ic,:))
+ plot(timephiMOLI,real(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
+ grid on
+ xlabel('$t$')
+ ylabel('Re$(\phi)$')
+ %Imag
+ subplot(336)
+ plot(time2d,imag(phiHeLaZ_rz),'-','DisplayName','HeLaZ RK4')
hold on
- ic = ic + 1;
- end
-end
-grid on
-xlabel('$t$')
-ylabel(['Re$(N_i^{pj})$'])
-%suptitle(TITLE);
-%if SAVEFIG; FIGNAME = ['Nipj_kz_',num2str(GRID.kzmin)]; save_figure; end;
-
-% phi
-timephi = h5read(filename,'/data/var2d/time');
-kr = h5read(filename,'/data/var2d/phi/coordkr');
-kz = h5read(filename,'/data/var2d/phi/coordkz');
-phiHeLaZ = zeros(numel(timephi),numel(kr),numel(kz));
-for it = 1:numel(timephi)
- tmp = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
- phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
-end
-
-timephiMOLI = results.time;
-phiMOLI = zeros(size(timephiMOLI));
-for it = 1:numel(timephiMOLI)
- phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
-end
-
-%fig = figure;
-%Real
-subplot(323)
-plot(timephi,real(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
-hold on
-plot(timephiMOLI,real(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
-grid on
-xlabel('$t$')
-ylabel('Re$(\phi)$')
-%Imag
-subplot(324)
-plot(timephi,imag(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
-hold on
-plot(timephiMOLI,imag(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
-grid on
-xlabel('$t$')
-ylabel('Im$(\phi)$')
-%if SAVEFIG; FIGNAME = ['phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
-
-%% phi error
-timephi = h5read(filename,'/data/var2d/time');
-kr = h5read(filename,'/data/var2d/phi/coordkr');
-kz = h5read(filename,'/data/var2d/phi/coordkz');
-phiHeLaZ = zeros(numel(timephi),numel(kr),numel(kz));
-for it = 1:numel(timephi)
- tmp = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
- phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
-end
-
-timephiMOLI = results.time;
-phiMOLI = zeros(size(timephiMOLI));
-for it = 1:numel(timephiMOLI)
- phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
-end
-
-ERR1 = abs(real(phiMOLI(1:numel(timephi)) - phiHeLaZ));
-ERR2 = abs(imag(phiMOLI(1:numel(timephi)) - phiHeLaZ));
-
-%fig = figure;
-subplot(325);
-plot(timephi,ERR1,'-','DisplayName','Real')
-hold on
-plot(timephi,ERR2,'--','DisplayName','Imag')
-%title(TITLE);
-grid on
-xlabel('$t$')
-ylabel('$e(\phi)$')
-legend('show')
-%if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
-
-% Growth rate fit quantity (Ni00)
-subplot(326);
-
-y1 = squeeze(abs(Nipj(1,1,1,1,:))); % HeLaZ
-y2 = squeeze(abs(results.Napj(:,bari(0,0))));
-ERR3 = abs(y2(1:numel(timephi))-y1);
-
-yyaxis left
-semilogy(x1,y1,'-','DisplayName','HeLaZ',...
- 'color', line_colors(1,:)); hold on;
-semilogy(x2,y2,'--','DisplayName','MOLI',...
- 'color', 'k')
-grid on
-xlabel('$t$')
-ylabel('$|N_i^{00}|$')
-legend('show')
-
-yyaxis right
-semilogy(x1, ERR3,'color', line_colors(2,:));
-grid on
-xlabel('$t$')
-ylabel('$e(N_i^{00})$')
-
-% Finish and save
-suptitle(TITLE);
-if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
+ plot(timephiMOLI,imag(phiMOLI),'--','DisplayName',['MOLI ',MOLIsolvername])
+ grid on
+ xlabel('$t$')
+ ylabel('Im$(\phi)$')
+ % phi error
+ ERR1 = abs(real(phiMOLI(1:OUTPUTS.nsave_2d:numel(phiHeLaZ_rz)*OUTPUTS.nsave_2d) - phiHeLaZ_rz));
+ ERR2 = abs(imag(phiMOLI(1:OUTPUTS.nsave_2d:numel(phiHeLaZ_rz)*OUTPUTS.nsave_2d) - phiHeLaZ_rz));
+ subplot(339);
+ semilogy(time2d,ERR1,'-','DisplayName','Real')
+ hold on
+ semilogy(time2d,ERR2,'--','DisplayName','Imag')
+ grid on
+ xlabel('$t$')
+ ylabel('$e(\phi)$')
+ legend('show')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Finish and save
+TITLE = ['$(k_r,k_z)=(',num2str(kr(ikr)),',',num2str(kz(ikz)),')$, ', TITLE];
+suptitle(TITLE);
+if strcmp(MODEL.NON_LIN,'.true.'); LINEARITY = 'nl';
+else; LINEARITY = 'lin'; end;
+if SAVEFIG; FIGNAME = LINEARITY; save_figure; end;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/wk/compute_Sapj.m b/wk/compute_Sapj.m
deleted file mode 100644
index b594fdb1..00000000
--- a/wk/compute_Sapj.m
+++ /dev/null
@@ -1,41 +0,0 @@
-%meshgrid
-[KR, KZ] = meshgrid(kr,kz);
-%BE = sqrt(KR.^2+KZ.^2)*params.mu*sqrt(2);
-BI = sqrt(KR.^2+KZ.^2)*sqrt(2*params.tau);
-
-%non linear term, time evolution
-Sipj = zeros(GRID.nkr, GRID.nkz, numel(time));
-
-%padding
-Mr = 2*GRID.nkr; Mz = 2*GRID.nkz;
-F_pad = zeros(Mr, Mz);
-G_pad = zeros(Mr, Mz);
-
-p = 0; j = 0; %pj moment
-
-for it = 1:numel(time) % time loop
- Sipj_pad = zeros(Mr, Mz);
- for n = 0, GRID.jmaxi % Sum over Laguerre
- %First conv term
- F = (KZ-KR).*phiHeLaZ(it,:,:)*kernel(n,BE);
-
- %Second conv term
- G = 0.*(KZ-KR);
- for s = 0:min(n+j,GRID.jmaxi)
- G = G + dnjs(n,j,s) .* Napj(p,s,:,:,it);
- end
- G = (KZ-KR) .* G;
-
- %Padding
- F_pad(1:GRID.nkr,1:GRID.nkz) = F_;
- G_pad(1:GRID.nkr,1:GRID.nkz) = G_;
-
- %Inverse fourier transform to real space
- f = ifft2(F_pad);
- g = ifft2(G_pad);
-
- %Conv theorem
- Sipj_pad = Sipj_pad + fft2(f.*g);
- end
- Sipj(:,:,it) = Sipj_pad(1:GRID.nkr,1:GRID.nkz);
-end
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index 2a9ffd21..b0220251 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,56 +1,16 @@
&BASIC
nrun=100000
dt=0.01
- tmax=100 ! time normalized to 1/omega_pe
+ tmax=200
/
&GRID
- pmaxe =8 ! Electron Hermite moments
- jmaxe = 4 ! Electron Laguerre moments
- pmaxi = 8 ! Ion Hermite moments
- jmaxi = 4 ! Ion Laguerre moments
- nkr = 16
- krmin = 0.1
- krmax = 1
- nkz = 16
+ pmaxe =5
+ jmaxe = 3
+ pmaxi = 5
+ jmaxi = 3
+ nkr = 1
+ krmin = 0
+ krmax = 0
+ nkz = 1
kzmin = 0.1
- kzmax = 1
-/
-&OUTPUT_PAR
- nsave_0d = 0
- nsave_1d = 0
- nsave_2d = 1
- nsave_5d = 0
- write_Ni00 = .true.
- write_moments = .true.
- write_phi = .true.
- write_doubleprecision = .true.
- resfile0 = 'results'
-/
-&MODEL_PAR
- ! Collisionality
- CO = -2 ! Collision operator (-1:Full Coulomb, 0: Dougherty)
- nu = 0.01 ! Normalized collision frequency normalized to plasma frequency
- tau_e = 1 ! T_e/T_e
- tau_i = 1 ! T_i/T_e temperature ratio
- sigma_e = 0.023338 ! sqrt(m_e/m_i) mass ratio
- sigma_i = 1 ! sqrt(m_i/m_i)
- q_e = -1 ! Electrons charge
- q_i = 1 ! Ions charge
- eta_n = 0 ! L_perp / L_n Density gradient
- eta_T = 0 ! L_perp / L_T Temperature gradient
- eta_B = 1 ! L_perp / L_B Magnetic gradient and curvature
- lambdaD = 0 ! Debye length
-/
-&INITIAL_CON
- ! Background value
- initback_moments=0.01 ! x 1e-3
- ! Noise amplitude
- initnoise_moments=0
- iseed=42
- selfmat_file='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_8_Jmaxe_4_Pmaxi_8_Jmaxi_4_pamaxx_10.h5'
- eimat_file='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_8_Jmaxe_4_Pmaxi_8_Jmaxi_4_pamaxx_10_tau_1.0000_mu_0.0233.h5'
- iemat_file='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_8_Jmaxe_4_Pmaxi_8_Jmaxi_4_pamaxx_10_tau_1.0000_mu_0.0233.h5'
-/
-&TIME_INTEGRATION_PAR
- numerical_scheme='RK4'
-/
\ No newline at end of file
+ kzmax = 0.1
diff --git a/wk/kr_kz_time_evolution.m b/wk/kr_kz_time_evolution.m
index 5be2d573..d7da8ba7 100644
--- a/wk/kr_kz_time_evolution.m
+++ b/wk/kr_kz_time_evolution.m
@@ -1,4 +1,5 @@
-%% Run linear simulation on a kr,kz grid and compute the non linear term afterwards
+% Run linear/nonlin simulation on a kr,kz grid and create gifs on the
+% fourier and real space
clear all; close all;
BASIC.SIMID = 'kr_kz_time_evolution'; % Name of the simulations
addpath(genpath('../matlab')) % ... add
@@ -6,26 +7,30 @@ addpath(genpath('../matlab')) % ... add
%% outputs options
OUTPUTS.nsave_0d = 0;
OUTPUTS.nsave_1d = 0;
-OUTPUTS.nsave_2d = 1;
+OUTPUTS.nsave_2d = 10;
OUTPUTS.nsave_5d = 0;
OUTPUTS.write_Ni00 = '.true.';
OUTPUTS.write_moments = '.true.';
OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_non_lin = '.false.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
+OUTPUTS.rstfile0 = '''restart''';
%% Grid parameters
-GRID.pmaxe = 8;
-GRID.jmaxe = 4;
-GRID.pmaxi = 8;
-GRID.jmaxi = 4;
-GRID.nkr = 16;
-GRID.krmin = 0.1;
-GRID.krmax = 1.;
-GRID.nkz = 16;
-GRID.kzmin = 0.1;
-GRID.kzmax = 1.0;
+GRID.pmaxe = 08; % Electron Hermite moments
+GRID.jmaxe = 04; % Electron Laguerre moments
+GRID.pmaxi = 08; % Ion Hermite moments
+GRID.jmaxi = 04; % Ion Laguerre moments
+GRID.nkr = 128; % kr grid resolution
+GRID.krmin =-2.0; % kr minimum value
+GRID.krmax = 2.0; % kr maximal value
+GRID.nkz = 128; % kz ''
+GRID.kzmin =-2.0; % ''
+GRID.kzmax = 2.0; % ''
+GRID.Pad = 2.0; % Zero padding for dealiasing (Mx = Pad * nkx)
%% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.NON_LIN = '.false.'; % Non linear term
MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
@@ -37,18 +42,21 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
% gradients L_perp/L_x
-MODEL.eta_n = 0.0; % Density
+MODEL.eta_n = 1.0; % Density
MODEL.eta_T = 0.0; % Temperature
-MODEL.eta_B = 1.0; % Magnetic
+MODEL.eta_B = 0.5; % Magnetic
% Debye length
MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000;
BASIC.dt = 0.01;
-BASIC.tmax = 100.0;
-INITIAL.initback_moments = 0.01;
-INITIAL.initnoise_moments = 0.0;
+BASIC.tmax = 200.0; %time normalized to 1/omega_pe
+INITIAL.RESTART = '.true.';
+INITIAL.backup_file = '''restart''';
+INITIAL.only_Na00 = '.false.';
+INITIAL.initback_moments = 1.0e-4;
+INITIAL.initnoise_moments = 5.0e-5;
INITIAL.iseed = 42;
INITIAL.selfmat_file = ...
['''../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(GRID.pmaxe),...
@@ -66,48 +74,46 @@ INITIAL.iemat_file = ...
%% Write input file and run HeLaZ
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
+MAKE = 'cd ..; make; cd wk';
+system(MAKE);
EXEC = ' ../bin/helaz ';
RUN = ['mpirun -np ' num2str(nproc)];
CMD = [RUN, EXEC, INPUT];
-system(CMD);
-
+%system(CMD);
+%% load results
+filename = [OUTPUTS.resfile0(2:end-1),'_00.h5'];
+[Ni00, kr, kz, time] = load_2D_data(filename, 'Ni00');
+[ phi, kr, kz, time] = load_2D_data(filename, 'phi');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Analysis and basic figures
default_plots_options
SAVEFIG = 1;
-filename = 'results_00.h5';
default_plots_options
-TITLE = [TITLE,', $k_z=',num2str(GRID.kzmin),'$'];
-
-%% Load HeLaZ moments
-
-moment = 'Ni00';
-
-kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
-time = h5read(filename,'/data/var2d/time');
-Ni00 = zeros(numel(kr),numel(kz),numel(time));
-for it = 1:numel(time)
- tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
- Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary;
-end
-
%% Gif creation
+if strcmp(MODEL.NON_LIN,'.true.'); LINEARITY = 'nl';
+else; LINEARITY = 'lin'; end;
disp('Creating gif..');
-skip = 50; % To skip some frames
-delay = 0.05; % Speed of the movie (smaller for faster)
+disp('');
+skip = 5; % To skip some frames
+delay = 0.03; % Speed of the movie (smaller for faster)
t = time(1:skip:end);
N = squeeze(Ni00(:,:,1:skip:end));
H = zeros(GRID.nkr, GRID.nkz, numel(t));
h = zeros(GRID.nkr, GRID.nkz, numel(t));
-Pad = 1.0; % Padding
+Pad = 2.0; % Padding
for it = 1:numel(t)
H(:,:,it) = abs(N(:,:,it));
- F = abs(ifft2(H,floor(Pad*GRID.nkr),floor(Pad*GRID.nkz)));
+ F = real(ifft2(N(:,:,it),ceil(Pad*GRID.nkr)+1,ceil(Pad*GRID.nkz)+1));
h(:,:,it) = F(1:GRID.nkr, 1:GRID.nkr);
end
-GIFNAME = 'Ni00';
-create_gif(kr, kz, t, H, BASIC, GRID, MODEL, delay, GIFNAME)
-% GIFNAME = 'Mi00';
-% create_gif(kr, kz, t, h, BASIC, GRID, MODEL, delay, GIFNAME)
+%%
+GIFNAME = [num2str(GRID.nkr),'x',num2str(GRID.nkz),'-',LINEARITY,'-Ni00'];
+create_gif(kr, kz, t, log(H), BASIC, GRID, MODEL, delay, GIFNAME, false)
+
+%%
+Nr = ceil(GRID.nkr/2); Nz = ceil(GRID.nkz/2);
+h_ = h;
+h_(h_==0) = nan;
+GIFNAME = [num2str(GRID.nkr),'x',num2str(GRID.nkz),'-',LINEARITY,'-Mi00'];
+create_gif(kr(Nr:end), kz(Nr:end), t, h_(Nr:end,Nz:end,:), BASIC, GRID, MODEL, delay, GIFNAME, false)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/wk/kz_linear.m b/wk/kz_linear.m
index 0c0036d2..e123f7c0 100644
--- a/wk/kz_linear.m
+++ b/wk/kz_linear.m
@@ -12,6 +12,7 @@ OUTPUTS.write_moments = '.true.';
OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
+OUTPUTS.rstfile0 = ['''../backup/',SIMID,''''];
%% Grid parameters
GRID.pmaxe = 25;
GRID.jmaxe = 9;
diff --git a/wk/launcher.m b/wk/launcher.m
index 9f031c53..c8fa8da5 100644
--- a/wk/launcher.m
+++ b/wk/launcher.m
@@ -13,10 +13,10 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 25;
-GRID.jmaxe = 10;
-GRID.pmaxi = 25;
-GRID.jmaxi = 10;
+GRID.pmaxe = 5;
+GRID.jmaxe = 3;
+GRID.pmaxi = 5;
+GRID.jmaxi = 3;
GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
--
GitLab