diff --git a/Makefile b/Makefile index d104f73..6a4c8ee 100644 --- a/Makefile +++ b/Makefile @@ -1,144 +1,144 @@ 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) run: all (cd wk; $(EXEC);) dirs: mkdir -p $(BINDIR) mkdir -p $(OBJDIR) mkdir -p $(MODDIR) - mkdir -p $(BACKUPDIR) + mkdir -p $(CHCKPTDIR) src/srcinfo.h: ( cd src/srcinfo; $(MAKE);) clean: cleanobj cleanmod @rm -f src/srcinfo.h @rm -f src/srcinfo/srcinfo.h cleanobj: @rm -f $(OBJDIR)/*o cleanmod: @rm -f $(MODDIR)/*mod @rm -f *.mod cleanbin: @rm -f $(EXEC) $(OBJDIR)/diagnose.o : src/srcinfo.h FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \ -$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/convolution_mod.o \ +$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_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)/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)/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 $(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)/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 $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@ $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(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 + $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/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)/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)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/mkl_dfti.o $(OBJDIR)/grid_mod.o + $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_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 + $(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)/grid_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 $(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 $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@ $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@ - $(OBJDIR)/fourier_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)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o + $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@ - $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o + $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/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 $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@ $(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@ - $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/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)/grid_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -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)/model_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)/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)/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 $(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)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@ $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90 $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@ - $(OBJDIR)/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)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.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)/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 $(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 $(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 + $(OBJDIR)/utility_mod.o : src/utility_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/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 cafb307..6ff72f9 100644 --- a/local/dirs.inc +++ b/local/dirs.inc @@ -1,10 +1,12 @@ #Local code, binaries, pputils library -PREFIX = $(HOME)/Documents/HeLaZ#write the path to the root of your MONOLIZ distro here +PREFIX = $(HOME)/Documents/HeLaZ#write the path to the root of your HeLaZ distrib here SRCDIR = $(PREFIX)/src 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 +FFTW3DIR = /usr/local/fftw3 +CHCKPTDIR = $(PREFIX)/checkpoint + +# Naming ideas : HeLaZ, MoNoLiT (Moment Non Linear Torroidal) diff --git a/local/make.inc b/local/make.inc index 5a62151..b1943c2 100644 --- a/local/make.inc +++ b/local/make.inc @@ -1,201 +1,207 @@ ################################################################ # # Section I: Preprocessor options # ################################################################ #0/1 = turn pardiso off/on PARDISO=0 #0/1 = turn MUMPS off/on ifdef MUMPS USEMUMPS=1 else USEMUMPS=0 endif #0/1 = turn parallel MG solver off/on USEPARMG=1 #0/1 = turn OpenMP compilation off/on USEOPENMP=0 #0/1 = turn verification off/on USEVERIFICATION=0 #Please set precompiler flags here FPPFLAGS=-DPARDISO=$(PARDISO) -DMUMPS=$(USEMUMPS) -DPARMG=$(USEPARMG) -DVERIFICATION=$(USEVERIFICATION) ################################################################ # # Section II: Compiler options # ################################################################ #Default optimization level (O=optimized, g=debug) OPTLEVEL = O F90FLAGS = CFLAGS = ifeq ($(OPTLEVEL), O) #optimized ifeq ($(COMPTYPE), i) #intel F90FLAGS += -O3 -xHOST endif ifeq ($(COMPTYPE), g) #gnu F90FLAGS += -ffree-line-length-0 -O3 endif ifeq ($(COMPTYPE), c) #cray F90FLAGS += endif endif ifeq ($(OPTLEVEL), g) #debug ifeq ($(COMPTYPE), i) #intel F90FLAGS += -g -traceback -CB endif ifeq ($(COMPTYPE), g) #gnu F90FLAGS += -ffree-line-length-0 -g -fbacktrace -fcheck=all -pedantic -Wall endif ifeq ($(COMPTYPE), c) #cray F90FLAGS += -g -O0 endif endif ifeq ($(USEOPENMP), 1) ifeq ($(COMPTYPE), i) #intel F90FLAGS += -qopenmp CFLAGS += -qopenmp endif ifeq ($(COMPTYPE), g) #gnu F90FLAGS += -fopenmp CFLAGS += -fopenmp endif endif ################################################################ # # Section III: Libraries and where to find them # ################################################################ 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 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 # ifdef MKL_LIB # common libraries LIBS += -lmkl_scalapack_lp64 -lmkl_core -lmkl_blacs_intelmpi_lp64 -lm # compiler dependent libraries ifeq ($(COMPTYPE), i) #intel LIBS += -lmkl_intel_lp64 endif ifeq ($(COMPTYPE), g) #gnu LIBS += -lmkl_gf_lp64 endif # different library depending on threaded or sequential ifeq ($(USEOPENMP), 1) LIBS += -lmkl_intel_thread else LIBS += -lmkl_sequential endif # MKL needs intel omp library when compiled with gnu ifeq ($(USEOPENMP), 1) ifeq ($(COMPTYPE), g) #gnu LIBS += -liomp5 endif endif LDIRS += -L$(MKL_LIB) IDIRS += -I$(MKL_LIB)/../../include endif +ifdef FFTW3DIR + LIBS += -lfftw3 + LDIRS += -L/usr/local/fftw3/lib + IDIRS += -I/usr/local/fftw3/include +endif + # # Add mumps libraries if preprocessor set # ifeq ($(USEMUMPS),1) LIBS += -lzmumps -ldmumps -lmumps_common -lpord IDIRS += -I$(MUMPS)/include LDIRS += -L$(MUMPS)/lib ifdef PARMETIS LIBS += -lparmetis -lmetis LDIRS += -L$(PARMETIS)/lib -L$(METIS)/lib endif ifdef SCOTCH LIBS += -lesmumps -lscotch -lscotcherr LDIRS += -L$(SCOTCH)/lib endif endif # # Add parallel multigrid libraries if preprocessor set # #ifeq ($(USEPARMG),1) # LIBS += -lparmgrid #endif ################################################################ # # Section IV: Set up compiler and compiler flags # ################################################################ CC = $(SPC_MPICC) FC = $(SPC_MPIF90) F90 = $(SPC_MPIF90) GBS_EXEC = $(BINDIR)/gbs_$(PLAT)_$(COMP) LDFLAGS = $(F90FLAGS) # # Intel + Dora + OpenMP needs dynamic linking # ifeq ($(PLAT), dora) ifeq ($(COMPTYPE), i) #intel ifeq ($(USEOPENMP),1) LDFLAGS += -dynamic endif endif endif ################################################################ # # Section V: Set up inclusion of modules and libraries during # compiling / linking phase # ################################################################ #Flag for finding external modules in MODDIR ifeq ($(COMPTYPE), i) #intel EXTMOD = -module $(MODDIR) endif ifeq ($(COMPTYPE), g) #gnu EXTMOD = -J $(MODDIR) endif ifeq ($(COMPTYPE), c) #cray EXTMOD = -em -J $(MODDIR) endif #Flag for finding external libraries in LDIR EXTLIBS = $(LDIRS) -Wl,--start-group $(LIBS) -Wl,--end-group #Flag for finding external include files EXTINC = $(IDIRS) diff --git a/src/advance_field.F90 b/src/advance_field.F90 index 709686d..c615d87 100644 --- a/src/advance_field.F90 +++ b/src/advance_field.F90 @@ -1,56 +1,56 @@ MODULE advance_field_routine USE prec_const implicit none CONTAINS SUBROUTINE advance_time_level USE basic USE time_integration use prec_const IMPLICIT NONE call set_updatetlevel(mod(updatetlevel,ntimelevel)+1) END SUBROUTINE advance_time_level SUBROUTINE advance_field( f, f_rhs ) USE basic USE time_integration USE array - USE fourier_grid + USE grid use prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f_rhs - INTEGER :: ikr, ikz, istage + INTEGER :: istage - SELECT CASE (updatetlevel) + SELECT CASE (updatetlevel) CASE(1) - DO ikz=ikzs,ikze + DO ikz=ikzs,ikze DO ikr=ikrs,ikre DO istage=1,ntimelevel f(ikr,ikz,1) = f(ikr,ikz,1) + dt*b_E(istage)*f_rhs(ikr,ikz,istage) END DO END DO END DO CASE DEFAULT - DO ikz=ikzs,ikze + DO ikz=ikzs,ikze DO ikr=ikrs,ikre f(ikr,ikz,updatetlevel) = f(ikr,ikz,1); DO istage=1,updatetlevel-1 f(ikr,ikz,updatetlevel) = f(ikr,ikz,updatetlevel) + & dt*A_E(updatetlevel,istage)*f_rhs(ikr,ikz,istage) END DO END DO END DO END SELECT END SUBROUTINE advance_field END MODULE advance_field_routine diff --git a/src/auxval.F90 b/src/auxval.F90 index 412809b..85ef9e9 100644 --- a/src/auxval.F90 +++ b/src/auxval.F90 @@ -1,19 +1,19 @@ subroutine auxval ! Set auxiliary values, at beginning of simulation USE basic - USE fourier_grid + USE grid USE array ! use mumps_bsplines, only: putrow_mumps => putrow, updtmat_mumps => updtmat, factor_mumps => factor, bsolve_mumps => bsolve ! from BSPLINES use prec_const IMPLICIT NONE - + INTEGER :: irows,irowe, irow, icol WRITE(*,*) '=== Set auxiliary values ===' call set_krgrid call set_kzgrid call set_pj CALL memory ! Allocate memory for global arrays END SUBROUTINE auxval diff --git a/src/basic_mod.F90 b/src/basic_mod.F90 index 859e81a..735a3b4 100644 --- a/src/basic_mod.F90 +++ b/src/basic_mod.F90 @@ -1,239 +1,238 @@ MODULE basic ! Basic module for time dependent problems use prec_const IMPLICIT none - 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) - INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose) - INTEGER :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose) + 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 :: RESTART = .FALSE. ! Signal end of run + 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) + INTEGER :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose) + INTEGER :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose) ! List of logical file units 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 MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4 MODULE PROCEDURE allocate_array_l1,allocate_array_l2,allocate_array_l3,allocate_array_l4 END INTERFACE allocate_array CONTAINS !================================================================================ SUBROUTINE basic_data ! Read basic data for input file use prec_const IMPLICIT NONE - NAMELIST /BASIC/ nrun, dt, tmax + NAMELIST /BASIC/ nrun, dt, tmax, RESTART READ(lu_in,basic) - !WRITE(*,basic) END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! Print date and time use prec_const IMPLICIT NONE CHARACTER(len=*) , INTENT(in) :: str CHARACTER(len=16) :: d, t, dat, time !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) time=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), time(1:12) ! END SUBROUTINE daytim !================================================================================ ! To allocate arrays of doubles, integers, etc. at run time SUBROUTINE allocate_array_dp1(a,is1,ie1) IMPLICIT NONE real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0.0_dp END SUBROUTINE allocate_array_dp1 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 ALLOCATE(a(is1:ie1,is2:ie2)) a=0.0_dp END SUBROUTINE allocate_array_dp2 SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0.0_dp END SUBROUTINE allocate_array_dp3 SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0.0_dp END SUBROUTINE allocate_array_dp4 SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a 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.0_dp END SUBROUTINE allocate_array_dp5 !======================================== SUBROUTINE allocate_array_dc1(a,is1,ie1) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc1 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 ALLOCATE(a(is1:ie1,is2:ie2)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc2 SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc3 SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc4 SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE DOUBLE COMPLEX, DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a 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=CMPLX(0.0_dp,0.0_dp) END SUBROUTINE allocate_array_dc5 !======================================== SUBROUTINE allocate_array_i1(a,is1,ie1) IMPLICIT NONE INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=0 END SUBROUTINE allocate_array_i1 SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2) IMPLICIT NONE INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=0 END SUBROUTINE allocate_array_i2 SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=0 END SUBROUTINE allocate_array_i3 SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=0 END SUBROUTINE allocate_array_i4 SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a 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 !======================================== SUBROUTINE allocate_array_l1(a,is1,ie1) IMPLICIT NONE LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1 ALLOCATE(a(is1:ie1)) a=.false. END SUBROUTINE allocate_array_l1 SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2) IMPLICIT NONE LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2 ALLOCATE(a(is1:ie1,is2:ie2)) a=.false. END SUBROUTINE allocate_array_l2 SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3)) a=.false. END SUBROUTINE allocate_array_l3 SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4) IMPLICIT NONE LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4 ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4)) a=.false. END SUBROUTINE allocate_array_l4 SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5) IMPLICIT NONE real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a 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=.false. END SUBROUTINE allocate_array_l5 END MODULE basic diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90 index ca2787e..f6c2fd2 100644 --- a/src/compute_Sapj.F90 +++ b/src/compute_Sapj.F90 @@ -1,109 +1,109 @@ SUBROUTINE compute_Sapj(Sepj, Sipj) USE array, ONLY : dnjs USE basic - USE convolution + USE fourier USE fields!, ONLY : phi, moments_e, moments_i - USE fourier_grid + USE 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 + INTEGER :: in, is 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/convolution_mod.F90 b/src/convolution_mod.F90 index b50e917..7ca4cb8 100644 --- a/src/convolution_mod.F90 +++ b/src/convolution_mod.F90 @@ -1,221 +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 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 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 + !! 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 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 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 a532329..47e9dac 100644 --- a/src/diagnose.F90 +++ b/src/diagnose.F90 @@ -1,343 +1,382 @@ SUBROUTINE diagnose(kstep) ! Diagnostics, writing simulation state to disk USE basic - USE fourier_grid + USE grid USE diagnostics_par - USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach + USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf USE model USE initial_par USE fields USE time_integration USE prec_const IMPLICIT NONE INCLUDE 'srcinfo.h' INTEGER, INTENT(in) :: kstep INTEGER, parameter :: BUFSIZE = 2 INTEGER :: rank, dims(1) = (/0/) CHARACTER(len=256) :: str, fname !_____________________________________________________________________________ ! 1. Initial diagnostics - IF (kstep .EQ. 0) THEN - + IF ((kstep .EQ. 0)) THEN ! 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' ELSE WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' END IF ! 1.1 Initial run IF (write_doubleprecision) THEN CALL creatf(resfile, fidres, real_prec='d') ELSE CALL creatf(resfile, fidres) END IF WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' call flush(6) ! Data group CALL creatg(fidres, "/data", "data") 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 - iframe2d=0 - END IF + iframe2d=0 + ENDIF CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) IF (cstep==0) THEN - iframe5d=0 + iframe5d=0 END IF CALL attach(fidres,"/data/var5d/" , "frames", iframe5d) ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) ! var2d group (electro. pot., Ni00 moment) rank = 0 CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R") CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number") 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 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) 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 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) END IF ! var5d group (moments) rank = 0 CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R") CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number") IF (write_moments) THEN - CALL creatg(fidres, "/data/var5d/moments_e", "moments_e") - CALL putarr(fidres, "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0) - 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) - CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0) - CALL putarr(fidres, "/data/var5d/moments_i/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0) + CALL creatg(fidres, "/data/var5d/moments_e", "moments_e") + CALL putarr(fidres, "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e), "p_e",ionode=0) + 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) + CALL putarr(fidres, "/data/var5d/moments_i/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0) + 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) + 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 WRITE(*,*) 'AUTHOR=', AUTHOR WRITE(*,*) 'HOST=', HOST - IF(jobnum .LE. 99) THEN - WRITE(str,'(a,i2.2)') "/data/input.",jobnum - ELSE - WRITE(str,'(a,i3.2)') "/data/input.",jobnum - END IF + WRITE(str,'(a,i2.2)') "/data/input" + rank=0 CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters') CALL attach(fidres, TRIM(str), "version", VERSION) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "branch", BRANCH) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "author", AUTHOR) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "host", HOST) !defined in srcinfo.h CALL attach(fidres, TRIM(str), "start_time", time) - CALL attach(fidres, TRIM(str), "start_cstep", cstep) + CALL attach(fidres, TRIM(str), "start_cstep", cstep-1) CALL attach(fidres, TRIM(str), "dt", dt) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) - CALL fourier_grid_outputinputs(fidres, str) + CALL grid_outputinputs(fidres, str) CALL output_par_outputinputs(fidres, str) CALL model_outputinputs(fidres, str) CALL initial_outputinputs(fidres, str) CALL time_integration_outputinputs(fidres, str) ! Save STDIN (input file) of this run IF(jobnum .LE. 99) THEN WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum ELSE WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum END IF INQUIRE(unit=lu_in, name=fname) CLOSE(lu_in) CALL putfile(fidres, TRIM(str), TRIM(fname),ionode=0) - END IF + ELSEIF((kstep .EQ. 0)) THEN + + IF(jobnum .LE. 99) THEN + WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' + ELSE + WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5' + END IF + CALL openf(resfile,fidres, 'D'); + + ENDIF !_____________________________________________________________________________ ! 2. Periodic diagnostics ! IF (kstep .GE. 0) THEN ! 2.1 0d history arrays IF (nsave_0d .NE. 0) THEN IF ( MOD(cstep, nsave_0d) == 0 ) THEN CALL diagnose_0d END IF END IF ! 2.2 1d profiles ! empty in our case ! 2.3 2d profiles IF (nsave_2d .NE. 0) THEN IF (MOD(cstep, nsave_2d) == 0) THEN CALL diagnose_2d + IF (MOD(cstep, nsave_2d*10) == 0) THEN + WRITE(*,"(F4.0,A,F4.0)") time,"/",tmax + ENDIF END IF END IF ! 2.4 3d profiles IF (nsave_5d .NE. 0) THEN IF (MOD(cstep, nsave_5d) == 0) THEN CALL diagnose_5d END IF END IF ! 2.5 Backups - ! - ! - ! + IF (MOD(cstep, nsave_cp) == 0) THEN + CALL checkpoint_save + ENDIF !_____________________________________________________________________________ ! 3. Final diagnostics ELSEIF (kstep .EQ. -1) THEN CALL cpu_time(finish) ! Close all diagnostic files CALL closef(fidres) END IF END SUBROUTINE diagnose SUBROUTINE diagnose_0d USE basic USE futils, ONLY: append USE diagnostics_par 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) END SUBROUTINE diagnose_0d SUBROUTINE diagnose_2d USE basic USE futils, ONLY: append, getatt, attach, putarrnd USE fields USE time_integration USE diagnostics_par 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) iframe2d=iframe2d+1 CALL attach(fidres,"/data/var2d/" , "frames", iframe2d) IF (write_phi) THEN CALL write_field2d(phi(:,:), 'phi') END IF IF (write_Ni00) THEN CALL write_field2d(moments_i(1,1,:,:,updatetlevel), 'Ni00') END IF CONTAINS SUBROUTINE write_field2d(field, text) USE futils, ONLY: attach, putarr - USE fourier_grid, ONLY: ikrs,ikre, ikzs,ikze + USE grid, ONLY: ikrs,ikre, ikzs,ikze USE prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ikrs:ikre, ikzs:ikze), INTENT(IN) :: field CHARACTER(*), INTENT(IN) :: text CHARACTER(LEN=50) :: dset_name WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d 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 USE basic USE futils, ONLY: append, getatt, attach, putarrnd USE fields 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 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 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) iframe5d=iframe5d+1 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 + IF (write_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 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 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 + +SUBROUTINE checkpoint_save + USE basic + USE grid + USE diagnostics_par + USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf + USE model + USE initial_par + USE fields + USE time_integration + IMPLICIT NONE + + WRITE(rstfile,'(a,a3)') TRIM(rstfile0),'.h5' + + CALL creatf(rstfile, fidrst, real_prec='d') + CALL creatg(fidrst, '/Basic', 'Basic data') + CALL attach(fidrst, '/Basic', 'cstep', cstep) + CALL attach(fidrst, '/Basic', 'time', time) + CALL attach(fidrst, '/Basic', 'jobnum', jobnum) + CALL attach(fidrst, '/Basic', 'dt', dt) + CALL attach(fidrst, '/Basic', 'iframe2d', iframe2d) + CALL attach(fidrst, '/Basic', 'iframe5d', iframe5d) + + ! Write state of system to restart file + CALL putarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,& + ikrs:ikre,ikzs:ikze,1),ionode=0) + CALL putarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,& + ikrs:ikre,ikzs:ikze,1),ionode=0) + CALL closef(fidrst) + WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved!" + +END SUBROUTINE checkpoint_save diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90 index c0b86ad..a9005ea 100644 --- a/src/diagnostics_par_mod.F90 +++ b/src/diagnostics_par_mod.F90 @@ -1,70 +1,70 @@ MODULE diagnostics_par ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE 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 - + INTEGER, PUBLIC, PROTECTED :: nsave_cp = 1e3 ! HDF5 file 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 + 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 + INTEGER, PUBLIC :: fidrst ! FID for restart file PUBLIC :: output_par_readinputs, output_par_outputinputs CONTAINS SUBROUTINE output_par_readinputs ! Read the input parameters USE basic, ONLY : lu_in USE prec_const IMPLICIT NONE NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d NAMELIST /OUTPUT_PAR/ write_Ni00, write_moments, write_phi, write_non_lin, write_doubleprecision NAMELIST /OUTPUT_PAR/ resfile0, rstfile0 READ(lu_in,output_par) END SUBROUTINE output_par_readinputs SUBROUTINE output_par_outputinputs(fidres, str) ! ! Write the input parameters to the results_xx.h5 file ! USE prec_const USE futils, ONLY: attach IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str 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) END SUBROUTINE output_par_outputinputs END MODULE diagnostics_par diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90 deleted file mode 100644 index 3cbcbc7..0000000 --- a/src/fourier_grid_mod.F90 +++ /dev/null @@ -1,163 +0,0 @@ -MODULE fourier_grid - ! Grid module for spatial discretization - - USE prec_const - IMPLICIT NONE - PRIVATE - - ! GRID Namelist - INTEGER, PUBLIC, PROTECTED :: pmaxe = 2 ! The maximal electron Hermite-moment computed - 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 - INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i - INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze - - ! Toroidal direction - REAL(dp), PUBLIC, PROTECTED :: deltakr - REAL(dp), PUBLIC, PROTECTED :: deltakz - - ! Grids containing position in fourier space - REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray - REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray - - ! Grid containing the polynomials degrees - 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 - PUBLIC :: fourier_grid_readinputs, fourier_grid_outputinputs - PUBLIC :: bare, bari - -CONTAINS - - SUBROUTINE set_krgrid - USE prec_const - IMPLICIT NONE - INTEGER :: ikr - ! Start and END indices of grid - ikrs = 1 - ikre = nkr - ! Grid spacings, precompute some inverses - IF ( nkr .GT. 1 ) THEN ! To avoid case with 0 intervals - deltakr = (krmax - krmin) / REAL(nkr-1,dp) - ELSE - deltakr = 0; - ENDIF - ! Discretized kr positions - ALLOCATE(krarray(ikrs:ikre)) - DO ikr = ikrs,ikre - krarray(ikr) = krmin + REAL(ikr-1,dp) * deltakr - END DO - END SUBROUTINE set_krgrid - - SUBROUTINE set_kzgrid - USE prec_const - IMPLICIT NONE - INTEGER :: ikz - ! Start and END indices of grid - ikzs = 1 - ikze = nkz - ! Grid spacings, precompute some inverses - IF ( nkz .GT. 1 ) THEN - deltakz = (kzmax - kzmin) / REAL(nkz-1,dp) - ELSE - deltakz = 0; - ENDIF - ! Discretized kz positions - ALLOCATE(kzarray(ikzs:ikze)) - DO ikz = ikzs,ikze - kzarray(ikz) = kzmin + REAL(ikz-1,dp) * deltakz - END DO - END SUBROUTINE set_kzgrid - - SUBROUTINE set_pj - USE prec_const - IMPLICIT NONE - INTEGER :: ip, ij - ips_e = 1; ipe_e = pmaxe + 1 - ips_i = 1; ipe_i = pmaxi + 1 - ALLOCATE(parray_e(ips_e:ipe_e)) - ALLOCATE(parray_i(ips_i:ipe_i)) - DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO - DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO - - ijs_e = 1; ije_e = jmaxe + 1 - ijs_i = 1; ije_i = jmaxi + 1 - ALLOCATE(jarray_e(ijs_e:ije_e)) - 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 - ! Read the input parameters - - USE basic, ONLY : lu_in - - USE prec_const - IMPLICIT NONE - - NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, & - nkr, krmin, krmax, nkz, kzmin, kzmax, & - Pad - - READ(lu_in,grid) - !WRITE(*,grid) - - END SUBROUTINE fourier_grid_readinputs - - - SUBROUTINE fourier_grid_outputinputs(fidres, str) - ! Write the input parameters to the results_xx.h5 file - - USE futils, ONLY: attach - - USE prec_const - IMPLICIT NONE - - INTEGER, INTENT(in) :: fidres - CHARACTER(len=256), INTENT(in) :: str - CALL attach(fidres, TRIM(str), "pmaxe", pmaxe) - CALL attach(fidres, TRIM(str), "jmaxe", jmaxe) - CALL attach(fidres, TRIM(str), "pmaxi", pmaxi) - CALL attach(fidres, TRIM(str), "jmaxi", jmaxi) - CALL attach(fidres, TRIM(str), "nkr", nkr) - CALL attach(fidres, TRIM(str), "krmin", krmin) - CALL attach(fidres, TRIM(str), "krmax", krmax) - 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 - - - FUNCTION bare(p_,j_) - IMPLICIT NONE - INTEGER :: bare, p_, j_ - bare = (jmaxe+1)*p_ + j_ + 1 - END FUNCTION - - FUNCTION bari(p_,j_) - IMPLICIT NONE - INTEGER :: bari, p_, j_ - bari = (jmaxi+1)*p_ + j_ + 1 - END FUNCTION - - -END MODULE fourier_grid diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90 new file mode 100644 index 0000000..7bb0bc2 --- /dev/null +++ b/src/fourier_mod.F90 @@ -0,0 +1,129 @@ +MODULE fourier + + USE prec_const + USE grid + use, intrinsic :: iso_c_binding + implicit none + + INCLUDE 'fftw3.f03' + + PRIVATE + PUBLIC :: fft_r2cc + PUBLIC :: ifft_cc2r + PUBLIC :: fft2_r2cc + PUBLIC :: ifft2_cc2r + PUBLIC :: convolve_2D_F2F + PUBLIC :: set_descriptors, free_descriptors ! Temporary to switch easily with MKL DFTI + + CONTAINS + + SUBROUTINE fft_r2cc(fx_in, Fk_out) + + IMPLICIT NONE + REAL(dp), DIMENSION(Nr), INTENT(IN) :: fx_in + COMPLEX(dp), DIMENSION(Nkr), INTENT(OUT):: Fk_out + integer*8 plan + + call dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE) + call dfftw_execute_dft_r2c(plan, fx_in, Fk_out) + call dfftw_destroy_plan(plan) + + END SUBROUTINE fft_r2cc + + + + SUBROUTINE ifft_cc2r(Fk_in, fx_out) + + IMPLICIT NONE + COMPLEX(dp), DIMENSION(Nkr), INTENT(IN) :: Fk_in + REAL(dp), DIMENSION(Nr), INTENT(OUT):: fx_out + integer*8 plan + + call dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE) + call dfftw_execute_dft_c2r(plan, Fk_in, fx_out) + call dfftw_destroy_plan(plan) + fx_out = fx_out/Nr + + END SUBROUTINE ifft_cc2r + + + + SUBROUTINE fft2_r2cc( ffx_in, FFk_out ) + + IMPLICIT NONE + REAL(dp), DIMENSION(Nr,Nz), INTENT(IN) :: ffx_in + COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out + integer*8 plan + + !!! 2D Forward FFT ________________________! + call dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE) + call dfftw_execute_dft_r2c(plan,ffx_in,FFk_out) + call dfftw_destroy_plan(plan) + + END SUBROUTINE fft2_r2cc + + + + SUBROUTINE ifft2_cc2r( FFk_in, ffx_out ) + + IMPLICIT NONE + COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: FFk_in + REAL(dp), DIMENSION(Nr,Nz), INTENT(OUT) :: ffx_out + COMPLEX(dp), DIMENSION(Nkr,Nkz) :: tmp_c + integer*8 plan + + tmp_c = FFk_in + !!! 2D Backward FFT ________________________! + call dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE) + call dfftw_execute_dft_c2r(plan,tmp_c,ffx_out) + call dfftw_destroy_plan(plan) + ffx_out = ffx_out/Nr/Nz + + END SUBROUTINE ifft2_cc2r + + + + !!! 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 grid, ONLY : AA_r, AA_z + IMPLICIT NONE + COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: F_2D, G_2D ! input fields + COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT) :: C_2D ! output convolutioned field + REAL(dp), DIMENSION(Nr,Nz) :: ff, gg ! iFFT of input fields + REAL(dp), DIMENSION(Nr,Nz) :: ffgg ! will receive the product of f*g in Real + INTEGER :: ikr, ikz + REAL :: a_r + + ! 2D inverse Fourier transform + CALL ifft2_cc2r(F_2D,ff); + CALL ifft2_cc2r(G_2D,gg); + + ! Product in physical space + ffgg = ff * gg; + + ! 2D Fourier tranform + CALL fft2_r2cc(ffgg,C_2D); + + ! Anti aliasing (2/3 rule) + DO ikr = 1,Nkr + a_r = AA_r(ikr) + DO ikz = 1,Nkz + C_2D(ikr,ikz) = C_2D(ikr,ikz) * a_r * AA_z(ikz) + ENDDO + ENDDO + +END SUBROUTINE convolve_2D_F2F + + +! Empty set/free routines to switch easily with MKL DFTI +SUBROUTINE set_descriptors + +END SUBROUTINE set_descriptors + +SUBROUTINE free_descriptors + +END SUBROUTINE free_descriptors + +END MODULE fourier diff --git a/src/grid_mod.F90 b/src/grid_mod.F90 new file mode 100644 index 0000000..c239c8d --- /dev/null +++ b/src/grid_mod.F90 @@ -0,0 +1,240 @@ +MODULE grid + ! Grid module for spatial discretization + USE prec_const + IMPLICIT NONE + PRIVATE + + ! GRID Namelist + INTEGER, PUBLIC, PROTECTED :: pmaxe = 1 ! The maximal electron Hermite-moment computed + INTEGER, PUBLIC, PROTECTED :: jmaxe = 1 ! The maximal electron Laguerre-moment computed + INTEGER, PUBLIC, PROTECTED :: pmaxi = 1 ! The maximal ion Hermite-moment computed + INTEGER, PUBLIC, PROTECTED :: jmaxi = 1 ! The maximal ion Laguerre-moment computed + INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal ion Laguerre-moment computed + INTEGER, PUBLIC, PROTECTED :: Nr = 16 ! Number of total internal grid points in r + REAL(dp), PUBLIC, PROTECTED :: Lr = 1._dp ! horizontal length of the spatial box + INTEGER, PUBLIC, PROTECTED :: Nz = 16 ! Number of total internal grid points in z + REAL(dp), PUBLIC, PROTECTED :: Lz = 1._dp ! vertical length of the spatial box + INTEGER, PUBLIC, PROTECTED :: Nkr = 8 ! Number of total internal grid points in kr + REAL(dp), PUBLIC, PROTECTED :: Lkr = 1._dp ! horizontal length of the fourier box + INTEGER, PUBLIC, PROTECTED :: Nkz = 16 ! Number of total internal grid points in kz + REAL(dp), PUBLIC, PROTECTED :: Lkz = 1._dp ! vertical length of the fourier box + + ! For Orszag filter + REAL(dp), PUBLIC, PROTECTED :: two_third_krmax + REAL(dp), PUBLIC, PROTECTED :: two_third_kzmax + + ! 1D Antialiasing arrays (2/3 rule) + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_r + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_z + + ! Grids containing position in physical space + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: rarray + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray + REAL(dp), PUBLIC, PROTECTED :: deltar, deltaz + INTEGER, PUBLIC, PROTECTED :: irs, ire, izs, ize + INTEGER, PUBLIC :: ir,iz ! counters + + ! Grids containing position in fourier space + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray + REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray + REAL(dp), PUBLIC, PROTECTED :: deltakr, deltakz + INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze, a + INTEGER, PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin + INTEGER, PUBLIC :: ikr, ikz, ip, ij ! counters + + ! Grid containing the polynomials degrees + INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e + INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i + INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e + INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i + INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e + INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i + + ! Public Functions + PUBLIC :: set_pj + PUBLIC :: set_rgrid, set_zgrid + PUBLIC :: set_krgrid, set_kzgrid + PUBLIC :: grid_readinputs, grid_outputinputs + PUBLIC :: bare, bari +CONTAINS + + SUBROUTINE set_pj + USE prec_const + IMPLICIT NONE + INTEGER :: ip, ij + ips_e = 1; ipe_e = pmaxe + 1 + ips_i = 1; ipe_i = pmaxi + 1 + ALLOCATE(parray_e(ips_e:ipe_e)) + ALLOCATE(parray_i(ips_i:ipe_i)) + DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO + DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO + + ijs_e = 1; ije_e = jmaxe + 1 + ijs_i = 1; ije_i = jmaxi + 1 + ALLOCATE(jarray_e(ijs_e:ije_e)) + 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 set_rgrid + USE prec_const + IMPLICIT NONE + INTEGER :: ir + ! Start and END indices of grid + irs = 1 + ire = Nr + ! Grid spacing + IF ( Nr .GT. 1 ) THEN ! To avoid case with 0 intervals + deltar = Lr / REAL(Nr-1,dp) + ELSE + deltar = 0; + ENDIF + ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1) + ALLOCATE(rarray(irs:ire)) + DO ir = irs,ire + rarray(ir) = deltar*(MODULO(ir-1,Nr/2)-Nr/2*FLOOR(2.*real(ir-1)/real(Nr))) + END DO + rarray(Nr/2+1) = -rarray(Nr/2+1) + + END SUBROUTINE set_rgrid + + SUBROUTINE set_zgrid + USE prec_const + IMPLICIT NONE + INTEGER :: iz + ! Start and END indices of grid + izs = 1 + ize = Nr + ! Grid spacing + IF ( Nz .GT. 1 ) THEN ! To avoid case with 0 intervals + deltaz = Lz / REAL(Nz-1,dp) + ELSE + deltaz = 0; + ENDIF + ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1) + ALLOCATE(zarray(irs:ire)) + DO iz = izs,ize + zarray(iz) = deltaz*(MODULO(iz-1,Nz/2)-Nz/2*FLOOR(2.*real(iz-1)/real(Nz))) + END DO + zarray(Nz/2+1) = -zarray(Nz/2+1) + + END SUBROUTINE set_zgrid + + SUBROUTINE set_krgrid + USE prec_const + IMPLICIT NONE + INTEGER :: ikr + + Nkr = Nr/2+1 ! Defined only on positive kr since fields are real + ! Start and END indices of grid + ikrs = 1 + ikre = Nkr + ! Grid spacings + deltakr = 2._dp*PI/Nr/deltar + + ! Discretized kr positions ordered as dk*(0 1 2) + ALLOCATE(krarray(ikrs:ikre)) + DO ikr = ikrs,ikre + krarray(ikr) = REAL(ikr-1,dp) * deltakr + if (krarray(ikr) .EQ. 0) THEN + ikr_0 = ikr + ENDIF + END DO + + ! Orszag 2/3 filter + two_third_krmax = 2._dp/3._dp*deltakr*Nkr + ALLOCATE(AA_r(ikrs:ikre)) + DO ikr = ikrs,ikre + IF ( (krarray(ikr) .GT. -two_third_krmax) .AND. (krarray(ikr) .LT. two_third_krmax) ) THEN + AA_r(ikr) = 1._dp; + ELSE + AA_r(ikr) = 0._dp; + ENDIF + END DO + END SUBROUTINE set_krgrid + + SUBROUTINE set_kzgrid + USE prec_const + IMPLICIT NONE + + Nkz = Nz; + ! Start and END indices of grid + ikzs = 1 + ikze = Nkz + ! Grid spacings + deltakz = 2._dp*PI/Nkz/deltaz + + ! Discretized kz positions ordered as dk*(0 1 2 -3 -2 -1) + ALLOCATE(kzarray(ikzs:ikze)) + DO ikz = ikzs,ikze + kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz))) + if (kzarray(ikz) .EQ. 0) THEN + ikz_0 = ikz + ENDIF + END DO + kzarray(Nz/2+1) = -kzarray(Nz/2+1) + + ! Orszag 2/3 filter + two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2); + ALLOCATE(AA_z(ikzs:ikze)) + DO ikz = ikzs,ikze + IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN + AA_z(ikz) = 1._dp; + ELSE + AA_z(ikz) = 0._dp; + ENDIF + END DO + END SUBROUTINE set_kzgrid + + SUBROUTINE grid_readinputs + ! Read the input parameters + USE prec_const + IMPLICIT NONE + INTEGER :: lu_in = 90 ! File duplicated from STDIN + + NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, & + Nr, Lr, Nz, Lz + READ(lu_in,grid) + + END SUBROUTINE grid_readinputs + + + SUBROUTINE grid_outputinputs(fidres, str) + ! Write the input parameters to the results_xx.h5 file + + USE futils, ONLY: attach + + USE prec_const + IMPLICIT NONE + + INTEGER, INTENT(in) :: fidres + CHARACTER(len=256), INTENT(in) :: str + CALL attach(fidres, TRIM(str), "pmaxe", pmaxe) + CALL attach(fidres, TRIM(str), "jmaxe", jmaxe) + CALL attach(fidres, TRIM(str), "pmaxi", pmaxi) + CALL attach(fidres, TRIM(str), "jmaxi", jmaxi) + CALL attach(fidres, TRIM(str), "nr", nr) + CALL attach(fidres, TRIM(str), "Lr", Lr) + CALL attach(fidres, TRIM(str), "nz", nz) + CALL attach(fidres, TRIM(str), "Lz", Lz) + CALL attach(fidres, TRIM(str), "nkr", nkr) + CALL attach(fidres, TRIM(str), "Lkr", Lkr) + CALL attach(fidres, TRIM(str), "nkz", nkz) + CALL attach(fidres, TRIM(str), "Lkz", Lkz) + END SUBROUTINE grid_outputinputs + + FUNCTION bare(p_,j_) + IMPLICIT NONE + INTEGER :: bare, p_, j_ + bare = (jmaxe+1)*p_ + j_ + 1 + END FUNCTION + + FUNCTION bari(p_,j_) + IMPLICIT NONE + INTEGER :: bari, p_, j_ + bari = (jmaxi+1)*p_ + j_ + 1 + END FUNCTION + +END MODULE grid diff --git a/src/inital.F90 b/src/inital.F90 index d2e839c..69312a5 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,189 +1,188 @@ !******************************************************************************! !!!!!! initialize the moments and load/build coeff tables !******************************************************************************! SUBROUTINE inital USE basic USE model, ONLY : CO, NON_LIN USE prec_const 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(*,*) '=== Load Full Coulomb matrix ===' CALL load_FC_mat ENDIF END SUBROUTINE inital !******************************************************************************! !******************************************************************************! !!!!!!! Initialize the moments randomly !******************************************************************************! SUBROUTINE init_moments USE basic - USE fourier_grid + USE grid USE fields USE initial_par USE time_integration USE prec_const IMPLICIT NONE - INTEGER :: ip,ij, ikr,ikz INTEGER, DIMENSION(12) :: iseedarr REAL(dp) :: noise ! Seed random number generator iseedarr(:)=iseed CALL RANDOM_SEED(PUT=iseedarr) CALL set_updatetlevel(1) IF ( only_Na00 ) THEN ! Spike initialization on density only 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 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 END DO 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 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 + USE grid USE initial_par USE array use model USE futils, ONLY : openf, getarr, closef IMPLICIT NONE - INTEGER :: ip,ij, ip2,ij2 + INTEGER :: ip2,ij2 INTEGER :: fid1, fid2, fid3, fid4 !!!!!!!!!!!! Electron matrices !!!!!!!!!!!! ! 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 CALL openf(eimat_file,fid2, 'r', 'D'); CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT) CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF) CALL closef(fid2) !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!! ! 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 CALL getarr(fid3, '/Caapj/Ceepj', Ciipj) ! get array (moli format) ELSE CALL getarr(fid3, '/Caapj/Ciipj', Ciipj) ! get array (moli format) ENDIF CALL closef(fid3) ! get the Test and Back field ion electron collision matrix CALL openf(iemat_file, fid4, 'r', 'D'); CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT) CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF) CALL closef(fid4) END SUBROUTINE load_FC_mat !******************************************************************************! diff --git a/src/memory.F90 b/src/memory.F90 index 0e1e380..58b8b39 100644 --- a/src/memory.F90 +++ b/src/memory.F90 @@ -1,41 +1,41 @@ SUBROUTINE memory ! Allocate arrays (done dynamically otherwise size is unknown) USE array USE basic USE fields - USE fourier_grid + USE grid USE time_integration USE model, ONLY: CO, NON_LIN USE prec_const IMPLICIT NONE ! Moments and moments rhs CALL allocate_array( moments_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel ) ! Electrostatic potential CALL allocate_array(phi, ikrs,ikre, ikzs,ikze) ! Collision matrix IF (CO .EQ. -1) THEN CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1)) CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1)) CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1)) CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1)) 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/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90 index 6aa516e..0b85811 100644 --- a/src/moments_eq_rhs.F90 +++ b/src/moments_eq_rhs.F90 @@ -1,377 +1,377 @@ SUBROUTINE moments_eq_rhs USE basic USE time_integration USE array USE fields - USE fourier_grid + USE grid USE model USE prec_const IMPLICIT NONE - INTEGER :: ip,ij, ikr,ikz, ip2, ij2 ! loops indices + INTEGER :: ip2, ij2 ! loops indices REAL(dp) :: ip_dp, ij_dp REAL(dp) :: kr, kz, kperp2 REAL(dp) :: taue_qe_etaB, taui_qi_etaB REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop REAL(dp) :: xphij, xphijp1, xphijm1 ! ESpot. factors depending on the pj loop REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency !write(*,*) '----------------------------------------' !Precompute species dependant factors taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling taui_qi_etaB = tau_i/q_i * eta_B sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp nu_e = nu ! electron-ion collision frequency (where already multiplied by 0.532) nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ. nu_ee = nu_e/SQRT2 ! e-e coll. frequ. nu_ie = nu*sigma_e**2 ! i-e coll. frequ. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Electrons moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1 ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxe) ! N_e^{p+2,j} coeff xNapp2j = taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp)) ! N_e^{p-2,j} coeff xNapm2j = taue_qe_etaB * SQRT(ip_dp * (ip_dp - 1._dp)) factj = 1.0 ! Start of the recursive factorial jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1 ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxe) ! N_e^{p,j+1} coeff xNapjp1 = -taue_qe_etaB * (ij_dp + 1._dp) ! N_e^{p,j-1} coeff xNapjm1 = -taue_qe_etaB * ij_dp ! N_e^{pj} coeff xNapj = taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp) !! Collision operator pj terms xCapj = -nu_e*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis ! Dougherty part IF ( CO .EQ. -2) THEN IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20 xCa20 = nu_e * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01 xCa20 = -nu_e * SQRT2 * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN ! kronecker pj10 xCa20 = 0._dp xCa01 = 0._dp xCa10 = nu_e ELSE xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp ENDIF ENDIF !! Electrostatic potential pj terms IF (ip .EQ. 1) THEN ! kronecker p0 xphij = (eta_n + 2.*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp) ) xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp) xphijm1 = -(eta_T - eta_B)* ij_dp ELSE IF (ip .EQ. 3) THEN ! kronecker p2 xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphijp1 = 0._dp; xphijm1 = 0._dp ELSE xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp ENDIF ! Recursive factorial IF (ij_dp .GT. 0) THEN factj = factj * ij_dp ELSE factj = 1._dp ENDIF ! Loop on kspace krloope : DO ikr = ikrs,ikre kzloope : DO ikz = ikzs,ikze kr = krarray(ikr) ! Poloidal wavevector kz = kzarray(ikz) ! Toroidal wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector b_e2 = kperp2 * sigmae2_taue_o2 ! Bessel argument !! Compute moments and mixing terms ! term propto N_e^{p,j} TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel) ! term propto N_e^{p+2,j} IF (ip+2 .LE. pmaxe+1) THEN ! OoB check TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel) ELSE TNapp2j = 0._dp ENDIF ! term propto N_e^{p-2,j} IF (ip-2 .GE. 1) THEN ! OoB check TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel) ELSE TNapm2j = 0._dp ENDIF ! xterm propto N_e^{p,j+1} IF (ij+1 .LE. jmaxe+1) THEN ! OoB check TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel) ELSE TNapjp1 = 0._dp ENDIF ! term propto N_e^{p,j-1} IF (ij-1 .GE. 1) THEN ! OoB check TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel) ELSE TNapjm1 = 0._dp ENDIF !! Collision IF (CO .EQ. -2) THEN ! Dougherty Collision terms IF ( (pmaxe .GE. 2) ) THEN ! OoB check TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel) ELSE TColl20 = 0._dp ENDIF IF ( (jmaxe .GE. 1) ) THEN ! OoB check TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel) ELSE TColl01 = 0._dp ENDIF IF ( (pmaxe .GE. 1) ) THEN ! OoB check TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel) ELSE TColl10 = 0._dp ENDIF ! Total collisional term TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)& + TColl20 + TColl01 + TColl10 ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix) TColl = 0._dp ! Initialization ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms jloopee: DO ij2 = 1,jmaxe+1 TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) & *( nu_e * CeipjT(bare(ip-1,ij-1), bare(ip2-1,ij2-1)) & +nu_ee * Ceepj (bare(ip-1,ij-1), bare(ip2-1,ij2-1))) ENDDO jloopee ENDDO ploopee ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms jloopei: DO ij2 = 1,jmaxi+1 TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) & *(nu_e * CeipjF(bare(ip-1,ij-1), bari(ip2-1,ij2-1))) END DO jloopei ENDDO ploopei ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel) ENDIF !! Electrical potential term IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2 kernelj = b_e2**(ij-1) * exp(-b_e2)/factj kerneljp1 = kernelj * b_e2 /(ij_dp + 1._dp) kerneljm1 = kernelj * ij_dp / b_e2 Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz) ELSE Tphi = 0._dp ENDIF ! Sum of all linear terms moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = & -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 END DO jloope END DO ploope !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ploopi : DO ip = ips_i, ipe_i ! This loop is from 1 to pmaxi+1 ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxi) ! x N_i^{p+2,j} coeff xNapp2j = taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp)) ! x N_i^{p-2,j} coeff xNapm2j = taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1._dp)) factj = 1._dp ! Start of the recursive factorial jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxi) ! x N_i^{p,j+1} coeff xNapjp1 = -taui_qi_etaB * (ij_dp + 1._dp) ! x N_i^{p,j-1} coeff xNapjm1 = -taui_qi_etaB * ij_dp ! x N_i^{pj} coeff xNapj = taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp) !! Collision operator pj terms xCapj = -nu_i*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis ! Dougherty part IF ( CO .EQ. -2) THEN IF ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20 xCa20 = nu_i * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01 xCa20 = -nu_i * SQRT2 * 2._dp/3._dp xCa01 = -SQRT2 * xCa20 xCa10 = 0._dp ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN xCa20 = 0._dp xCa01 = 0._dp xCa10 = nu_i ELSE xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp ENDIF ENDIF !! Electrostatic potential pj terms IF (ip .EQ. 1) THEN ! krokecker p0 xphij = (eta_n + 2._dp*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp)) xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp) xphijm1 = -(eta_T - eta_B)* ij_dp ELSE IF (ip .EQ. 3) THEN !krokecker p2 xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphijp1 = 0._dp; xphijm1 = 0._dp ELSE xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp ENDIF ! Recursive factorial IF (ij_dp .GT. 0) THEN factj = factj * ij_dp ELSE factj = 1._dp ENDIF ! Loop on kspace krloopi : DO ikr = ikrs,ikre kzloopi : DO ikz = ikzs,ikze kr = krarray(ikr) ! Poloidal wavevector kz = kzarray(ikz) ! Toroidal wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector b_i2 = kperp2 * sigmai2_taui_o2 ! Bessel argument !! Compute moments and mixing terms ! term propto N_i^{p,j} TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel) ! term propto N_i^{p+2,j} IF (ip+2 .LE. pmaxi+1) THEN ! OoB check TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel) ELSE TNapp2j = 0._dp ENDIF ! term propto N_i^{p-2,j} IF (ip-2 .GE. 1) THEN ! OoB check TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel) ELSE TNapm2j = 0._dp ENDIF ! xterm propto N_i^{p,j+1} IF (ij+1 .LE. jmaxi+1) THEN ! OoB check TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel) ELSE TNapjp1 = 0._dp ENDIF ! term propto N_i^{p,j-1} IF (ij-1 .GE. 1) THEN ! OoB check TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel) ELSE TNapjm1 = 0._dp ENDIF !! Collision IF (CO .EQ. -2) THEN ! Dougherty Collision terms IF ( (pmaxi .GE. 2) ) THEN ! OoB check TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel) ELSE TColl20 = 0._dp ENDIF IF ( (jmaxi .GE. 1) ) THEN ! OoB check TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel) ELSE TColl01 = 0._dp ENDIF IF ( (pmaxi .GE. 1) ) THEN ! OoB check TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel) ELSE TColl10 = 0._dp ENDIF ! Total collisional term TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)& + TColl20 + TColl01 + TColl10 ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb for ions (COSOlver matrix) !!! TColl = 0._dp ! Initialization ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms jloopii: DO ij2 = 1,jmaxi+1 TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) & *( nu_ie * CiepjT(bari(ip-1,ij-1), bari(ip2-1,ij2-1)) & +nu_i * Ciipj (bari(ip-1,ij-1), bari(ip2-1,ij2-1))) ENDDO jloopii ENDDO ploopii ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms jloopie: DO ij2 = 1,jmaxe+1 TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) & *(nu_ie * CiepjF(bari(ip-1,ij-1), bare(ip2-1,ij2-1))) ENDDO jloopie ENDDO ploopie ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel) ENDIF !! Electrical potential term IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2 kernelj = b_i2**(ij-1) * exp(-b_i2)/factj kerneljp1 = kernelj * b_i2 /(ij_dp + 1._dp) kerneljm1 = kernelj * ij_dp / b_i2 Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz) ELSE Tphi = 0._dp ENDIF ! Sum of linear terms moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = & -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 END DO jloopi END DO ploopi END SUBROUTINE moments_eq_rhs diff --git a/src/poisson.F90 b/src/poisson.F90 index ee353c0..abfa726 100644 --- a/src/poisson.F90 +++ b/src/poisson.F90 @@ -1,100 +1,100 @@ SUBROUTINE poisson ! Solve poisson equation to get phi USE time_integration, ONLY: updatetlevel USE array USE fields - USE fourier_grid + USE grid use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD USE prec_const IMPLICIT NONE - INTEGER :: ikr,ikz, ini,ine, i0j + INTEGER :: ini,ine, i0j REAL(dp) :: ini_dp, ine_dp REAL(dp) :: kr, kz, kperp2 REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build REAL(dp) :: b_e2, b_i2, alphaD REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation REAL(dp) :: sum_kernel2_e, sum_kernel2_i ! Store sum Kn^2 COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store sum Kn*Napn REAL(dp) :: gammaD COMPLEX(dp) :: gammaD_phi !Precompute species dependant factors sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2 ! = kperp^2 tau_a sigma_a^2/2 qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum qi2_taui = (q_i**2)/tau_i DO ikr=ikrs,ikre DO ikz=ikzs,ikze kr = krarray(ikr) kz = kzarray(ikz) kperp2 = kr**2 + kz**2 !! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2) b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)) ! Initialization for n = 0 (ine = 1) Kne = exp(-b_e2) sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel) sum_kernel2_e = Kne**2 !write(*,*) 'K0e = ', Kne ! loop over n only if the max polynomial degree is 1 or more if (jmaxe .GT. 0) then DO ine=2,jmaxe+1 ! ine = n+1 ine_dp = REAL(ine-1,dp) Kne = Kne * b_e2/ine_dp ! We update iteratively the kernel functions (to spare factorial computations) sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel) sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ... !write(*,*) 'K',ine-1,'e = ', Kne END DO endif !! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2) b_i2 = kperp2 * sigmai2_taui_o2 ! Initialization for n = 0 (ini = 1) Kni = exp(-b_i2) sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel) sum_kernel2_i = Kni**2 !write(*,*) 'K0i = ', Kni ! loop over n only if the max polynomial degree is 1 or more if (jmaxi .GT. 0) then DO ini=2,jmaxi+1 ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax) Kni = Kni * b_i2/ini_dp ! We update iteratively the kernel functions this time for ions ... sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel) sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ... !write(*,*) 'K',ini-1,'i = ', Kni END DO endif !!! Assembling the poisson equation alphaD = kperp2 * lambdaD**2 gammaD = alphaD + qe2_taue * (1._dp - sum_kernel2_e) & ! Called Poisson_ in MOLI + qi2_taui * (1._dp - sum_kernel2_i) gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i 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 !write(*,*) 'gammaD =',gammaD END SUBROUTINE poisson diff --git a/src/readinputs.F90 b/src/readinputs.F90 index 879df9c..44bfcf8 100644 --- a/src/readinputs.F90 +++ b/src/readinputs.F90 @@ -1,33 +1,32 @@ SUBROUTINE readinputs ! Additional data specific for a new run - USE fourier_grid, ONLY: fourier_grid_readinputs + USE grid, ONLY: grid_readinputs USE diagnostics_par, ONLY: output_par_readinputs USE model, ONLY: model_readinputs USE initial_par, ONLY: initial_readinputs USE time_integration, ONLY: time_integration_readinputs USE prec_const IMPLICIT NONE !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) ! Load grid data from input file - CALL fourier_grid_readinputs + CALL grid_readinputs ! Load diagnostic options from input file CALL output_par_readinputs ! Load model parameters from input file CALL model_readinputs ! Load initial condition parameters from input file CALL initial_readinputs ! Load parameters for time integration from input file CALL time_integration_readinputs END SUBROUTINE readinputs - diff --git a/src/stepon.F90 b/src/stepon.F90 index 678b81e..47d0014 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,65 +1,65 @@ SUBROUTINE stepon ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme) USE basic USE time_integration USE fields, ONLY: moments_e, moments_i, phi USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj - USE fourier_grid + USE 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, ikr, ikz + INTEGER :: num_step DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4 ! Compute right hand side of moments hierarchy equation CALL moments_eq_rhs ! 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 ! Update electrostatic potential CALL poisson ! 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 SUBROUTINE checkfield_all ! Check all the fields for inf or nan IF(.NOT.nlend) THEN nlend=nlend .or. checkfield(phi,' phi') DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e nlend=nlend .or. checkfield(moments_e(ip,ij,:,:,updatetlevel),' moments_e') ENDDO ENDDO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i nlend=nlend .or. checkfield(moments_i(ip,ij,:,:,updatetlevel),' moments_i') ENDDO ENDDO ENDIF END SUBROUTINE checkfield_all END SUBROUTINE stepon diff --git a/src/utility_mod.F90 b/src/utility_mod.F90 index 24445b8..16dd278 100644 --- a/src/utility_mod.F90 +++ b/src/utility_mod.F90 @@ -1,74 +1,74 @@ MODULE utility USE basic use prec_const IMPLICIT NONE CONTAINS FUNCTION is_nan(x,str) RESULT(isn) USE time_integration use prec_const IMPLICIT NONE real(dp), INTENT(IN) :: x CHARACTER(LEN=*), INTENT(IN) :: str LOGICAL :: isn isn=.FALSE. IF(x .NE. x) THEN isn=.TRUE. END IF IF((isn).AND.(str.NE.'')) THEN WRITE(*,'(a20,a25,i6.6,a20,i1)') str,' = NaN at timestep',cstep, ' and substep',updatetlevel CALL FLUSH(stdout) END IF END FUNCTION is_nan FUNCTION is_inf(x,str) RESULT(isi) USE time_integration use prec_const IMPLICIT NONE real(dp), INTENT(IN) :: x CHARACTER(LEN=*), INTENT(IN) :: str LOGICAL :: isi isi=.FALSE. IF(x+1.0== x) THEN isi=.TRUE. END IF IF((isi).AND.(str.NE.'')) THEN WRITE(*,'(a20,a25,i6.6,a20,i1)') str,' = Inf at timestep',cstep, ' and substep',updatetlevel CALL FLUSH(stdout) END IF END FUNCTION is_inf FUNCTION checkfield(field,str) RESULT(mlend) - USE fourier_grid + USE grid use prec_const IMPLICIT NONE COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(IN) :: field CHARACTER(LEN=*), INTENT(IN) :: str LOGICAL :: mlend COMPLEX(dp) :: sumfield sumfield=SUM(field) mlend= is_nan( REAL(sumfield),str).OR.is_inf( REAL(sumfield),str) & .OR. is_nan(AIMAG(sumfield),str).OR.is_inf(AIMAG(sumfield),str) END FUNCTION checkfield END MODULE utility