diff --git a/Makefile b/Makefile index 1da9fec..cc3344b 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 $(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)/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)/model_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \ $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \ $(OBJDIR)/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)/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)/grid_mod.o + $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(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)/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)/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)/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)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@ $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@ $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@ $(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)/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)/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)/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)/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)/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)/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)/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 2685432..e5471f0 100644 --- a/local/dirs.inc +++ b/local/dirs.inc @@ -1,14 +1,13 @@ #Local code, binaries, pputils library 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 -#FFTW3DIR = /usr/local/fftw3 FFTW3DIR = /home/ahoffman/lib/fftw-3.3.8 CHCKPTDIR = $(PREFIX)/checkpoint FUTILS = /home/ahoffman/lib/futils_1.4/src # Naming ideas : HeLaZ, MoNoLiT (Moment Non Linear Torroidal) diff --git a/local/make.inc b/local/make.inc index 4a88e77..d293ec8 100644 --- a/local/make.inc +++ b/local/make.inc @@ -1,209 +1,95 @@ ################################################################ # -# 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=0 - -#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 +# Section I: Compiler options # ################################################################ #Default optimization level (O=optimized, g=debug) -OPTLEVEL = O +OPTLEVEL = g 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 +# Section II: Libraries and where to find them # ################################################################ IDIRS := -I$(SPC_LOCAL)/include/$(OPTLEVEL) -#IDIRS := -I$(FUTILS)/include/$(OPTLEVEL) -I$(HDF5)/include -I$(ZLIB)/include -#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 -#LDIRS := -L$(FUTILS)/lib -L$(HDF5)/lib -L$(ZLIB)/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 + LIBS += -lfftw3 -lfftw3_mpi LDIRS += -L$(FFTW3DIR)/lib64 IDIRS += -I$(FFTW3DIR)/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/auxval.F90 b/src/auxval.F90 index 9520579..50efa00 100644 --- a/src/auxval.F90 +++ b/src/auxval.F90 @@ -1,26 +1,24 @@ subroutine auxval ! Set auxiliary values, at beginning of simulation USE basic USE grid USE array - ! use mumps_bsplines, only: putrow_mumps => putrow, updtmat_mumps => updtmat, factor_mumps => factor, bsolve_mumps => bsolve ! from BSPLINES + USE fourier, ONLY: initialize_FFT use prec_const IMPLICIT NONE INTEGER :: irows,irowe, irow, icol IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ===' CALL set_krgrid - CALL mpi_barrier(MPI_COMM_WORLD, ierr) - CALL set_kzgrid - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + CALL set_kzgrid ! Distributed dimension CALL set_pj - CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL memory ! Allocate memory for global arrays - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + + CALL initialize_FFT ! intialization of FFTW plans END SUBROUTINE auxval diff --git a/src/basic_mod.F90 b/src/basic_mod.F90 index 7499fd8..e71b4d0 100644 --- a/src/basic_mod.F90 +++ b/src/basic_mod.F90 @@ -1,280 +1,281 @@ MODULE basic ! Basic module for time dependent problems - + use, intrinsic :: iso_c_binding use prec_const IMPLICIT none ! INCLUDE 'fftw3-mpi.f03' 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 :: my_id ! identification number of current process INTEGER :: num_procs ! number of MPI processes - INTEGER :: grp_world ! world communicator - INTEGER :: comm_world + integer(C_INTPTR_T) :: alloc_local ! total local data allocation size (mpi process) + integer(C_INTPTR_T) :: local_nkr ! local size of parallelized kz direction + integer(C_INTPTR_T) :: local_kr_start ! local start of parallelized kz direction 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 ! To measure computation time real :: start, finish real(dp) :: t0_rhs, t0_adv_field, t0_poisson, t0_Sapj, t0_diag, t0_checkfield, t0_step real(dp) :: t1_rhs, t1_adv_field, t1_poisson, t1_Sapj, t1_diag, t1_checkfield, t1_step real(dp) :: tc_rhs, tc_adv_field, tc_poisson, tc_Sapj, tc_diag, tc_checkfield, tc_step 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, RESTART READ(lu_in,basic) !Init cumulative timers tc_rhs = 0.;tc_adv_field = 0.; tc_poisson = 0. tc_Sapj = 0.; tc_diag = 0.; tc_checkfield = 0. 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 !================================================================================ SUBROUTINE display_h_min_s(time) real :: time integer :: days, hours, mins, secs days = FLOOR(time/24./3600.); hours= FLOOR(time/3600.); mins = FLOOR(time/60.); secs = FLOOR(time); IF ( days .GT. 0 ) THEN !display day h min s hours = (time/3600./24. - days) * 24 mins = (time/3600. - days*24. - hours) * 60 secs = (time/60. - days*24.*60 - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]' ELSEIF ( hours .GT. 0 ) THEN !display h min s mins = (time/3600. - hours) * 60 secs = (time/60. - hours*60 - mins) * 60 WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]' ELSEIF ( mins .GT. 0 ) THEN !display min s secs = (time/60. - mins) * 60 WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]' ELSE ! display s WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]' ENDIF END SUBROUTINE display_h_min_s !================================================================================ ! 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 43693ae..b4d0251 100644 --- a/src/compute_Sapj.F90 +++ b/src/compute_Sapj.F90 @@ -1,190 +1,207 @@ SUBROUTINE compute_Sapj ! This routine is meant to compute the non linear term for each specie and degree !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes !! Sapj = Sum_n (ikr Kn phi)#(ikz Sum_s d_njs Naps) - (ikz Kn phi)#(ikr Sum_s d_njs Naps) !! where # denotes the convolution. USE array, ONLY : dnjs, Sepj, Sipj USE basic USE fourier USE fields!, ONLY : phi, moments_e, moments_i 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 + COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: F_, G_, CONV INTEGER :: in, is REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2 ! Execution time start CALL cpu_time(t0_Sapj) !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!! sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln 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 krloope1: DO ikr = ikrs,ikre ! Loop over kr kzloope1: DO ikz = ikzs,ikze ! Loop over kz kr = krarray(ikr) kz = kzarray(ikz) be_2 = sigmae2_taue_o2 * (kr**2 + kz**2) kerneln = be_2**(in-1)/factn * EXP(-be_2) ! First convolution term IF ( NON_LIN ) THEN F_(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln ELSE F_(ikr,ikz) = 0._dp ENDIF ! Background sinusoidal electrostatic potential phi_0 for KH inst. IF ( q_e .NE. 0._dp ) THEN ! If electron are not removed IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0 IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0 F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz) ENDIF ENDIF ENDIF ! 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) = imagu*kz*G_(ikr,ikz) ENDDO kzloope1 ENDDO krloope1 CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier - Sepj(ip,ij,:,:) = Sepj(ip,ij,:,:) + CONV ! Add it to Sepj + + DO ikr = ikrs,ikre ! Loop over kr + DO ikz = ikzs,ikze ! Loop over kz + Sepj(ip,ij,ikr,ikz) = Sepj(ip,ij,ikr,ikz) + CONV(ikr,ikz) ! Add it to Sepj + ENDDO + ENDDO IF ( NON_LIN ) THEN ! Fully non linear term ikz phi * ikr Napj krloope2: DO ikr = ikrs,ikre ! Loop over kr kzloope2: DO ikz = ikzs,ikze ! Loop over kz kr = krarray(ikr) kz = kzarray(ikz) be_2 = sigmae2_taue_o2 * (kr**2 + kz**2) kerneln = be_2**(in-1)/factn * EXP(-be_2) ! First convolution term F_(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln ! 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) = imagu*kr*G_(ikr,ikz) ENDDO kzloope2 ENDDO krloope2 CALL convolve_2D_F2F( F_, G_, CONV ) - Sepj(ip,ij,:,:) = Sepj(ip,ij,:,:) - CONV ! substract it to Sepj + + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + Sepj(ip,ij,ikr,ikz) = Sepj(ip,ij,ikr,ikz) - CONV(ikr,ikz) ! substract it to Sepj + ENDDO + ENDDO + ENDIF IF ( in + 1 .LE. jmaxe+1 ) THEN factn = real(in,dp) * factn ! compute (n+1)! ENDIF ENDDO nloope 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 krloopi1: DO ikr = ikrs,ikre ! Loop over kr kzloopi1: DO ikz = ikzs,ikze ! Loop over kz kr = krarray(ikr) kz = kzarray(ikz) bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2) kerneln = bi_2**(in-1)/factn * EXP(-bi_2) ! First convolution term IF ( NON_LIN ) THEN F_(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln ELSE F_(ikr,ikz) = 0._dp ENDIF ! Background sinusoidal electrostatic potential phi_0 for KH inst. IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0 IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0 F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz) ENDIF ENDIF ! Second convolution term 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) = imagu*kz*G_(ikr,ikz) ENDDO kzloopi1 ENDDO krloopi1 CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier - Sipj(ip,ij,:,:) = Sipj(ip,ij,:,:) + CONV ! Add it to Sipj + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + Sipj(ip,ij,ikr,ikz) = Sipj(ip,ij,ikr,ikz) + CONV(ikr,ikz) ! add it to Sipj + ENDDO + ENDDO krloopi2: DO ikr = ikrs,ikre ! Loop over kr kzloopi2: DO ikz = ikzs,ikze ! Loop over kz kr = krarray(ikr) kz = kzarray(ikz) bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2) kerneln = bi_2**(in-1)/factn * EXP(-bi_2) ! First convolution term F_(ikr,ikz) = imagu*kz*phi(ikr,ikz) * kerneln ! Second convolution term G_(ikr,ikz) = 0._dp ! initialization of the sum IF ( NON_LIN ) THEN 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) = imagu*kr*G_(ikr,ikz) ENDIF ! Background sinusoidal ion denstiy n_i0 for KH inst. IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0 IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=+-kr0 G_(ikr,ikz) = -A0KH*kr0KH**2/2._dp + G_(ikr,ikz) ENDIF ENDIF ENDDO kzloopi2 ENDDO krloopi2 CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve and back to Fourier - Sipj(ip,ij,:,:) = Sipj(ip,ij,:,:) - CONV! substract it to Sipj + DO ikr = ikrs,ikre + DO ikz = ikzs,ikze + Sipj(ip,ij,ikr,ikz) = Sipj(ip,ij,ikr,ikz) - CONV(ikr,ikz) ! substract it to Sepj + ENDDO + ENDDO IF ( in + 1 .LE. jmaxi+1 ) THEN factn = real(in,dp) * factn ! factorial(n+1) ENDIF ENDDO nloopi ENDDO jloopi ENDDO ploopi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Execution time end CALL cpu_time(t1_Sapj) tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj) END SUBROUTINE compute_Sapj diff --git a/src/control.F90 b/src/control.F90 index 3a78b91..b2c710a 100644 --- a/src/control.F90 +++ b/src/control.F90 @@ -1,94 +1,94 @@ SUBROUTINE control ! Control the run use basic use prec_const IMPLICIT NONE CALL cpu_time(start) !________________________________________________________________________________ ! 1. Prologue ! 1.1 Initialize the parallel environment - IF (my_id .EQ. 0) WRITE(*,*) 'Initialize MPI...' + IF (my_id .EQ. 1) WRITE(*,*) 'Initialize MPI...' CALL ppinit - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...MPI initialized' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...MPI initialized' CALL daytim('Start at ') ! 1.2 Define data specific to run - IF (my_id .EQ. 0) WRITE(*,*) 'Load basic data...' + IF (my_id .EQ. 1) WRITE(*,*) 'Load basic data...' CALL basic_data - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...basic data loaded.' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...basic data loaded.' - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! 1.3 Read input parameters from input file - IF (my_id .EQ. 0) WRITE(*,*) 'Read input parameters...' + IF (my_id .EQ. 1) WRITE(*,*) 'Read input parameters...' CALL readinputs - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...input parameters read' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...input parameters read' - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! 1.4 Set auxiliary values (allocate arrays, set grid, ...) - IF (my_id .EQ. 0) WRITE(*,*) 'Calculate auxval...' + IF (my_id .EQ. 1) WRITE(*,*) 'Calculate auxval...' CALL auxval - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...auxval calculated' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...auxval calculated' - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! 1.5 Initial conditions - IF (my_id .EQ. 0) WRITE(*,*) 'Create initial state...' + IF (my_id .EQ. 1) WRITE(*,*) 'Create initial state...' CALL inital - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial state created' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...initial state created' - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! 1.6 Initial diagnostics - IF (my_id .EQ. 0) WRITE(*,*) 'Initial diagnostics...' + IF (my_id .EQ. 1) WRITE(*,*) 'Initial diagnostics...' CALL diagnose(0) - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...initial diagnostics done' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...initial diagnostics done' - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL FLUSH(stdout) !________________________________________________________________________________ - IF (my_id .EQ. 0) WRITE(*,*) 'Time integration loop..' + IF (my_id .EQ. 1) WRITE(*,*) 'Time integration loop..' !________________________________________________________________________________ ! 2. Main loop DO CALL cpu_time(t0_step) ! Measuring time step = step + 1 cstep = cstep + 1 - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL stepon CALL mpi_barrier(MPI_COMM_WORLD, ierr) time = time + dt CALL tesend IF( nlend ) EXIT ! exit do loop CALL cpu_time(t0_diag) ! Measuring time CALL diagnose(step) CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag) - CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL cpu_time(t1_step); tc_step = tc_step + (t1_step - t0_step) END DO - IF (my_id .EQ. 0) WRITE(*,'(a/)') '...time integration done' + IF (my_id .EQ. 1) WRITE(*,'(a/)') '...time integration done' !________________________________________________________________________________ ! 9. Epilogue CALL diagnose(-1) CALL endrun - IF (my_id .EQ. 0) CALL daytim('Done at ') + IF (my_id .EQ. 1) CALL daytim('Done at ') CALL ppexit END SUBROUTINE control diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90 index a3dcca3..03978ff 100644 --- a/src/fourier_mod.F90 +++ b/src/fourier_mod.F90 @@ -1,192 +1,118 @@ MODULE fourier USE prec_const USE grid + USE basic use, intrinsic :: iso_c_binding implicit none ! INCLUDE 'fftw3.f03' INCLUDE 'fftw3-mpi.f03' PRIVATE - PUBLIC :: fft_r2cc - PUBLIC :: ifft_cc2r - PUBLIC :: fft2_r2cc - PUBLIC :: ifft2_cc2r + type(C_PTR) :: planf, planb, real_ptr, cmpx_ptr + real(C_DOUBLE), pointer :: real_data_1(:,:), real_data_2(:,:) + complex(C_DOUBLE_complex), pointer :: cmpx_data(:,:) + integer (C_INTPTR_T) :: nx, ny, nh + + PUBLIC :: initialize_FFT, finalize_FFT PUBLIC :: convolve_2D_F2F 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 ) - + SUBROUTINE initialize_FFT IMPLICIT NONE - REAL(dp), DIMENSION(Nr,Nz), INTENT(IN) :: ffx_in - COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out - integer*8 plan + nx = Nr; ny = Nz; nh = ( nx / 2 ) + 1 - !!! 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) + IF ( my_id .EQ. 1)write(*,*) 'Initialize FFT..' + CALL fftw_mpi_init - END SUBROUTINE fft2_r2cc + IF ( my_id .EQ. 1)write(*,*) 'distribution of the array along kr..' + alloc_local = fftw_mpi_local_size_2d(nh, ny, MPI_COMM_WORLD, local_nkr, local_kr_start) + write(*,*) 'local_nkr', local_nkr, 'local_kr_start', local_kr_start + real_ptr = fftw_alloc_real(2 * alloc_local) + cmpx_ptr = fftw_alloc_complex(alloc_local) + call c_f_pointer(real_ptr, real_data_1, [2*local_nkr, nx]) + call c_f_pointer(real_ptr, real_data_2, [2*local_nkr, nx]) + call c_f_pointer(cmpx_ptr, cmpx_data, [ local_nkr, nx]) - SUBROUTINE ifft2_cc2r( FFk_in, ffx_out ) + IF ( my_id .EQ. 1)write(*,*) 'setting FFT iFFT 2D plans..' + ! Backward FFT plan + planb = fftw_mpi_plan_dft_c2r_2d(ny, nx, cmpx_data, real_data_1, MPI_COMM_WORLD, & + FFTW_PATIENT) + ! Forward FFT plan + planf = fftw_mpi_plan_dft_r2c_2d(ny, nx, real_data_1, cmpx_data, MPI_COMM_WORLD, & + FFTW_PATIENT) + END SUBROUTINE initialize_FFT + SUBROUTINE finalize_FFT 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 + IF ( my_id .EQ. 0)write(*,*) 'finalize FFTW' + CALL dfftw_destroy_plan(planf) + CALL dfftw_destroy_plan(planb) + call fftw_mpi_cleanup + call fftw_free(real_ptr) + call fftw_free(cmpx_ptr) - END SUBROUTINE ifft2_cc2r + END SUBROUTINE finalize_FFT !!! Convolution 2D Fourier to Fourier - ! - Compute the convolution using the convolution theorem and MKL + ! - Compute the convolution using the convolution theorem SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D ) - USE basic - USE grid, ONLY : AA_r, AA_z, Lr, Lz 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 - + COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(IN) :: F_2D, G_2D ! input fields + COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze), INTENT(OUT) :: C_2D ! output convolutioned field + + ! initialize data to some function my_function(i,j) + do ikr = ikrs,ikre + do ikz = ikzs,ikze + cmpx_data(ikr-local_kr_start, ikz) = F_2D(ikr, ikz) + end do + end do + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + + ! 2D backward Fourier transform + ! write(*,*) my_id, 'iFFT(F) ..' + call fftw_mpi_execute_dft_c2r(planb, cmpx_data, real_data_1) + + ! initialize data to some function my_function(i,j) + do ikr = ikrs,ikre + do ikz = ikzs,ikze + cmpx_data(ikr-local_kr_start, ikz) = G_2D(ikr, ikz) + end do + end do ! 2D inverse Fourier transform - IF ( num_procs .EQ. 1 ) THEN - CALL ifft2_cc2r(F_2D,ff) - CALL ifft2_cc2r(G_2D,gg) - ELSE - CALL ifft2_cc2r_mpi(F_2D,ff) - CALL ifft2_cc2r_mpi(G_2D,gg) - ENDIF + ! write(*,*) my_id, 'iFFT(G) ..' + call fftw_mpi_execute_dft_c2r(planb, cmpx_data, real_data_2) ! Product in physical space - ffgg = ff * gg; + ! write(*,*) 'multply..' + real_data_1 = real_data_1 * real_data_2 ! 2D Fourier tranform - IF ( num_procs .EQ. 1 ) THEN - CALL fft2_r2cc(ffgg,C_2D) - ELSE - CALL fft2_r2cc_mpi(ffgg,C_2D) - ENDIF + ! write(*,*) my_id, 'FFT(f*g) ..' + call fftw_mpi_execute_dft_r2c(planf, real_data_1, cmpx_data) + + ! COPY TO OUTPUT ARRAY + ! write(*,*) my_id, 'out ..' + do ikr = ikrs,ikre + do ikz = ikzs,ikze + cmpx_data(ikr-local_kr_start,ikz) = C_2D(ikr,ikz) + end do + end do ! Anti aliasing (2/3 rule) DO ikr = ikrs,ikre - a_r = AA_r(ikr) DO ikz = ikzs,ikze - C_2D(ikr,ikz) = C_2D(ikr,ikz) * a_r * AA_z(ikz) + C_2D(ikr,ikz) = C_2D(ikr,ikz) * AA_r(ikr) * AA_z(ikz) ENDDO ENDDO END SUBROUTINE convolve_2D_F2F -!! MPI routines -SUBROUTINE fft2_r2cc_mpi( ffx_in, FFk_out ) - - IMPLICIT NONE - REAL(dp), DIMENSION(Nr,Nz), INTENT(IN) :: ffx_in - COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out - REAL(dp), DIMENSION(Nkr,Nkz) :: tmp_r - type(C_PTR) :: plan - integer(C_INTPTR_T) :: L - integer(C_INTPTR_T) :: M - - ! L = Nr; M = Nz - ! - ! tmp_r = ffx_in - ! !!! 2D Forward FFT ________________________! - ! plan = fftw_mpi_plan_dft_r2c_2d(L, M, tmp_r, FFk_out, MPI_COMM_WORLD, FFTW_ESTIMATE) - ! - ! if ((.not. c_associated(plan))) then - ! write(*,*) "plan creation error!!" - ! stop - ! end if - ! - ! CALL fftw_mpi_execute_dft_r2c(plan,tmp_r,FFk_out) - ! CALL fftw_destroy_plan(plan) - -END SUBROUTINE fft2_r2cc_mpi - - - -SUBROUTINE ifft2_cc2r_mpi( 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 - type(C_PTR) :: plan - integer(C_INTPTR_T) :: L - integer(C_INTPTR_T) :: M - - ! L = Nr; M = Nz - ! - ! tmp_c = FFk_in - ! !!! 2D Backward FFT ________________________! - ! plan = fftw_mpi_plan_dft_c2r_2d(L, M, tmp_c, ffx_out, MPI_COMM_WORLD, FFTW_ESTIMATE) - ! - ! if ((.not. c_associated(plan))) then - ! write(*,*) "plan creation error!!" - ! stop - ! end if - ! - ! CALL fftw_mpi_execute_dft_c2r(plan,tmp_c,ffx_out) - ! CALL fftw_destroy_plan(plan) - ! - ! ffx_out = ffx_out/Nr/Nz - -END SUBROUTINE ifft2_cc2r_mpi - -! 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 index 12f9adf..a3c0286 100644 --- a/src/grid_mod.F90 +++ b/src/grid_mod.F90 @@ -1,229 +1,221 @@ MODULE grid ! Grid module for spatial discretization USE prec_const USE basic 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 REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component LOGICAL, PUBLIC, PROTECTED :: CANCEL_ODD_P = .false. ! To cancel odd Hermite polynomials INTEGER, PUBLIC, PROTECTED :: pskip = 0 ! variable to skip p degrees or not ! 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_krgrid, set_kzgrid PUBLIC :: grid_readinputs, grid_outputinputs PUBLIC :: bare, bari CONTAINS SUBROUTINE set_pj USE prec_const IMPLICIT NONE INTEGER :: ip, ij IF (CANCEL_ODD_P) THEN pskip = 1 ELSE pskip = 0 ENDIF ips_e = 1; ipe_e = pmaxe/(1+pskip) + 1 ips_i = 1; ipe_i = pmaxi/(1+pskip) + 1 ALLOCATE(parray_e(ips_e:ipe_e)) ALLOCATE(parray_i(ips_i:ipe_i)) DO ip = ips_e,ipe_e; parray_e(ip) = (1+pskip)*(ip-1); END DO DO ip = ips_i,ipe_i; parray_i(ip) = (1+pskip)*(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_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 - IF ( Nkr .GT. 1 ) THEN - ikrs = my_id * (Nkr-1)/num_procs + 1 - ikre = (my_id+1) * (Nkr-1)/num_procs - ELSE - ikrs = 1; ikre = Nkr - ENDIF + ikrs = my_id * Nkr/num_procs + 1 + ikre = (my_id+1) * Nkr/num_procs + + IF(my_id .EQ. num_procs) ikre = Nkr - IF (my_id .EQ. num_procs-1) THEN - ikre = Nkr - ENDIF WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre ! Grid spacings IF (Lr .GT. 0) THEN deltakr = 2._dp*PI/Lr ELSE deltakr = 1.0 ENDIF ! 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 USE model, ONLY : kr0KH, ikz0KH IMPLICIT NONE Nkz = Nz; ! Start and END indices of grid ikzs = 1 ikze = Nkz - ! ikzs = my_id * Nz/num_procs + 1 - ! ikze = (my_id+1) * Nz/num_procs WRITE(*,*) 'ID = ',my_id,' ikzs = ', ikzs, ' ikze = ', ikze ! Grid spacings deltakz = 2._dp*PI/Lz ! 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 + if (kzarray(ikz) .EQ. 0) ikz_0 = ikz + if (ikz .EQ. Nz/2+1) kzarray(ikz) = -kzarray(ikz) END DO - kzarray(Nz/2+1) = -kzarray(Nz/2+1) + IF (my_id .EQ. 1) WRITE(*,*) 'ID = ',my_id,' kz = ',kzarray ! 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 ! Put kz0RT to the nearest grid point on kz ikz0KH = NINT(kr0KH/deltakr)+1 - kr0KH = kzarray(ikz0KH) + IF ( (ikz0KH .GE. ikzs) .AND. (ikz0KH .LE. ikze) ) kr0KH = kzarray(ikz0KH) 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, kpar, CANCEL_ODD_P 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 79d0068..7f50dbb 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -1,282 +1,281 @@ !******************************************************************************! !!!!!! initialize the moments and load/build coeff tables !******************************************************************************! SUBROUTINE inital USE basic USE model, ONLY : CO, NON_LIN USE prec_const USE coeff USE time_integration USE array, ONLY : Sepj,Sipj implicit none real :: start_init, end_init, time_estimation CALL set_updatetlevel(1) !!!!!! Set the moments arrays Nepj, Nipj !!!!!! ! WRITE(*,*) 'Init moments' IF ( RESTART ) THEN CALL load_cp ELSE CALL init_moments ENDIF !!!!!! Set phi !!!!!! - ! WRITE(*,*) 'Init phi' + IF (my_id .EQ. 1) WRITE(*,*) 'Init phi' CALL poisson !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!! IF ( NON_LIN .OR. (A0KH .NE. 0)) THEN; - ! WRITE(*,*) 'Init Sapj' + IF (my_id .EQ. 1) WRITE(*,*) 'Init Sapj' CALL compute_Sapj - ! WRITE(*,*) 'Building Dnjs table' CALL build_dnjs_table ENDIF !!!!!! Load the full coulomb collision operator coefficients !!!!!! IF (CO .EQ. -1) THEN IF (my_id .EQ. 0) WRITE(*,*) '=== Load Full Coulomb matrix ===' CALL load_FC_mat IF (my_id .EQ. 0) WRITE(*,*) '..done' ENDIF END SUBROUTINE inital !******************************************************************************! !******************************************************************************! !!!!!!! Initialize the moments randomly !******************************************************************************! SUBROUTINE init_moments USE basic USE grid USE fields USE prec_const USE utility, ONLY: checkfield USE initial_par IMPLICIT NONE REAL(dp) :: noise REAL(dp) :: kr, kz, sigma, gain, kz_shift INTEGER, DIMENSION(12) :: iseedarr ! Seed random number generator iseedarr(:)=iseed CALL RANDOM_SEED(PUT=iseedarr) 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 DO ikz=2,Nkz/2 !symmetry at kr = 0 CALL RANDOM_NUMBER(noise) moments_e( 1,1,1,ikz, :) = moments_e( 1,1,1,Nkz+2-ikz, :) moments_i( 1,1,1,ikz, :) = moments_i( 1,1,1,Nkz+2-ikz, :) END DO ELSE !**** Gaussian initialization (Hakim 2017) ********************************* ! sigma = 5._dp ! Gaussian sigma ! gain = 0.5_dp ! Gaussian mean ! kz_shift = 0.5_dp ! Gaussian z shift ! moments_i = 0; moments_e = 0; ! DO ikr=ikrs,ikre ! kr = krarray(ikr) ! DO ikz=ikzs,ikze ! kz = kzarray(ikz) ! moments_i( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp) ! moments_e( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp) ! END DO ! END DO !**** Broad noise initialization ******************************************* DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e DO ikr=ikrs,ikre DO ikz=ikzs,ikze CALL RANDOM_NUMBER(noise) moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) END DO END DO - ! IF ( ikrs .EQ. 1 ) THEN + IF ( ikrs .EQ. 1 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 CALL RANDOM_NUMBER(noise) moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :) END DO - ! ENDIF + ENDIF END DO END DO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i DO ikr=ikrs,ikre DO ikz=ikzs,ikze CALL RANDOM_NUMBER(noise) moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) END DO END DO - ! IF ( ikrs .EQ. 1 ) THEN + IF ( ikrs .EQ. 1 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 CALL RANDOM_NUMBER(noise) moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :) END DO - ! ENDIF + ENDIF END DO END DO ENDIF END SUBROUTINE init_moments !******************************************************************************! !******************************************************************************! !!!!!!! Load moments from a previous save !******************************************************************************! SUBROUTINE load_cp USE basic USE futils, ONLY: openf, closef, getarr, getatt, isgroup, isdataset USE grid USE fields USE diagnostics_par USE time_integration IMPLICIT NONE INTEGER :: rank, sz_, n_ INTEGER :: dims(1) = (/0/) CHARACTER(LEN=50) :: dset_name WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5' WRITE(*,'(3x,a)') "Resume from previous run" CALL openf(rstfile, fidrst) IF (isgroup(fidrst,'/Basic/moments_e')) THEN n_ = 0 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ DO WHILE (isdataset(fidrst, dset_name)) n_ = n_ + 1 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ ENDDO n_ = n_ - 1 WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_e", n_ CALL getatt(fidrst, dset_name, 'cstep', cstep) CALL getatt(fidrst, dset_name, 'time', time) CALL getatt(fidrst, dset_name, 'jobnum', jobnum) jobnum = jobnum+1 CALL getatt(fidrst, dset_name, 'iframe2d',iframe2d) CALL getatt(fidrst, dset_name, 'iframe5d',iframe5d) iframe2d = iframe2d-1; iframe5d = iframe5d-1 ! Read state of system from restart file CALL getarr(fidrst, dset_name, moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0) WRITE(dset_name, "(A, '/', i6.6)") "/Basic/moments_i", n_ CALL getarr(fidrst, dset_name, moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0) ELSE CALL getatt(fidrst, '/Basic', 'cstep', cstep) CALL getatt(fidrst, '/Basic', 'time', time) CALL getatt(fidrst, '/Basic', 'jobnum', jobnum) jobnum = jobnum+1 CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d) CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d) iframe2d = iframe2d-1; iframe5d = iframe5d-1 ! Read state of system from restart file CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0) CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0) ENDIF CALL closef(fidrst) WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!" END SUBROUTINE load_cp !******************************************************************************! !******************************************************************************! !!!!!!! Build the Laguerre-Laguerre coupling coefficient table !******************************************************************************! SUBROUTINE build_dnjs_table USE basic USE array, Only : dnjs 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 grid USE initial_par USE array use model USE futils, ONLY : openf, getarr, closef IMPLICIT NONE 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/poisson.F90 b/src/poisson.F90 index 0c8db4f..fbc1a3e 100644 --- a/src/poisson.F90 +++ b/src/poisson.F90 @@ -1,102 +1,102 @@ SUBROUTINE poisson ! Solve poisson equation to get phi USE basic, ONLY: t0_poisson, t1_poisson, tc_poisson USE time_integration, ONLY: updatetlevel USE array USE fields USE grid use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, DK USE prec_const IMPLICIT NONE 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) :: be_2, bi_2, 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 ! Execution time start CALL cpu_time(t0_poisson) !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 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)) be_2 = kperp2 * sigmae2_taue_o2 bi_2 = kperp2 * sigmai2_taui_o2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2) !! Sum(kernel * Moments_0n) ! Initialization for n = 0 (ine = 1) Kne = exp(-be_2) sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel) sum_kernel2_e = Kne**2 ! 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) ! We update iteratively the kernel functions (to spare factorial computations) Kne = Kne * be_2/ine_dp 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 ... END DO ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2) !! Sum(kernel * Moments_0n) ! Initialization for n = 0 (ini = 1) Kni = exp(-bi_2) sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel) sum_kernel2_i = Kni**2 ! 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) ! We update iteratively to spare factorial computations Kni = Kni * bi_2/ini_dp 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 ... 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 ( (gammaD .EQ. 0) .AND. (abs(kr)+abs(kz) .NE. 0._dp) ) THEN WRITE(*,*) 'Warning gammaD = 0', sum_kernel2_i ENDIF phi(ikr, ikz) = gammaD_phi/gammaD END DO END DO ! Cancel origin singularity -phi(ikr_0,ikz_0) = 0 +IF ((ikr_0 .GE. ikrs) .AND. (ikr_0 .LE. ikre)) phi(ikr_0,ikz_0) = 0 ! Execution time end CALL cpu_time(t1_poisson) tc_poisson = tc_poisson + (t1_poisson - t0_poisson) END SUBROUTINE poisson diff --git a/src/ppexit.F90 b/src/ppexit.F90 index 0c0ddb7..9995c91 100644 --- a/src/ppexit.F90 +++ b/src/ppexit.F90 @@ -1,13 +1,15 @@ SUBROUTINE ppexit ! Exit parallel environment USE basic + USE fourier, ONLY : finalize_FFT use prec_const IMPLICIT NONE - + + CALL finalize_FFT CALL mpi_barrier(MPI_COMM_WORLD, ierr) CALL mpi_finalize(ierr) - + END SUBROUTINE ppexit diff --git a/src/ppinit.F90 b/src/ppinit.F90 index c80bf87..57d7dc3 100644 --- a/src/ppinit.F90 +++ b/src/ppinit.F90 @@ -1,24 +1,18 @@ SUBROUTINE ppinit ! Parallel environment USE basic USE model, ONLY : NON_LIN use prec_const IMPLICIT NONE INTEGER :: version_prov=-1 CALL MPI_INIT(ierr) ! CALL MPI_INIT_THREAD(MPI_THREAD_SINGLE,version_prov,ierr) ! CALL MPI_INIT_THREAD(MPI_THREAD_FUNNELED,version_prov,ierr) CALL MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr) CALL MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr) - CALL MPI_COMM_GROUP(MPI_COMM_WORLD,grp_world,ierr) - ! CALL MPI_COMM_CREATE_GROUP(MPI_COMM_WORLD,grp_world,0,comm_self,ierr) - - ! IF (NON_LIN .AND. (num_procs .GT. 1)) THEN - ! CALL fftw_mpi_init - ! ENDIF - + END SUBROUTINE ppinit diff --git a/src/srcinfo.h b/src/srcinfo.h index e6f1a84..d9f601e 100644 --- a/src/srcinfo.h +++ b/src/srcinfo.h @@ -1,10 +1,10 @@ character(len=40) VERSION character(len=40) BRANCH character(len=20) AUTHOR character(len=40) EXECDATE character(len=40) HOST -parameter (VERSION='8775464-dirty') -parameter (BRANCH='master') +parameter (VERSION='003c315-dirty') +parameter (BRANCH='MPI') parameter (AUTHOR='ahoffman') -parameter (EXECDATE='Mon Nov 9 13:42:43 CET 2020') +parameter (EXECDATE='Fri Nov 20 14:14:05 CET 2020') parameter (HOST ='spcpc606') diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h index e6f1a84..d9f601e 100644 --- a/src/srcinfo/srcinfo.h +++ b/src/srcinfo/srcinfo.h @@ -1,10 +1,10 @@ character(len=40) VERSION character(len=40) BRANCH character(len=20) AUTHOR character(len=40) EXECDATE character(len=40) HOST -parameter (VERSION='8775464-dirty') -parameter (BRANCH='master') +parameter (VERSION='003c315-dirty') +parameter (BRANCH='MPI') parameter (AUTHOR='ahoffman') -parameter (EXECDATE='Mon Nov 9 13:42:43 CET 2020') +parameter (EXECDATE='Fri Nov 20 14:14:05 CET 2020') parameter (HOST ='spcpc606') diff --git a/src/stepon.F90 b/src/stepon.F90 index 4826a63..7f16f1a 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,98 +1,110 @@ 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 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 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 + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme) CALL advance_time_level + CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Update the moments with the hierarchy RHS (step by step) ! Execution time start CALL cpu_time(t0_adv_field) 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 + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! Execution time end CALL cpu_time(t1_adv_field) tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field) + CALL mpi_barrier(MPI_COMM_WORLD, ierr) ! Update electrostatic potential CALL poisson + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + ! Update nonlinear term IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN CALL compute_Sapj ENDIF + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + !(The two routines above are called in inital for t=0) CALL checkfield_all() + CALL mpi_barrier(MPI_COMM_WORLD, ierr) + END DO CONTAINS SUBROUTINE checkfield_all ! Check all the fields for inf or nan ! Execution time start CALL cpu_time(t0_checkfield) - CALL enforce_symetry() ! Enforcing symmetry on kr = 0 + IF ( ikrs .EQ. 1 ) CALL enforce_symetry() ! Enforcing symmetry on kr = 0 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 ! Execution time end CALL cpu_time(t1_checkfield) tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield) END SUBROUTINE checkfield_all SUBROUTINE enforce_symetry DO ip=ips_e,ipe_e DO ij=ijs_e,ije_e DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :)) END DO END DO END DO DO ip=ips_i,ipe_i DO ij=ijs_i,ije_i DO ikz=2,Nkz/2 !symmetry at kr = 0 moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :)) END DO END DO END DO DO ikz=2,Nkz/2 !symmetry at kr = 0 phi(1,ikz) = phi(1,Nkz+2-ikz) END DO END SUBROUTINE enforce_symetry END SUBROUTINE stepon diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m index 71fab61..20d1a91 100644 --- a/wk/analysis_2D.m +++ b/wk/analysis_2D.m @@ -1,379 +1,379 @@ %% Load results % JOBNUM = 0; load_results; % JOBNUM = 1; load_results; compile_results %% Retrieving max polynomial degree and sampling info Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe); Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi); Ns5D = numel(Ts5D); Ns2D = numel(Ts2D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts2D = Ts2D'; if strcmp(OUTPUTS.write_non_lin,'.true.') Si00 = squeeze(Sipj(1,1,:,:,:)); Se00 = squeeze(Sepj(1,1,:,:,:)); end %% Build grids Nkr = numel(kr); Nkz = numel(kz); [KZ,KR] = meshgrid(kz,kr); Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz); dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1); Lk = max(Lkr,Lkz); dr = 2*pi/Lk; dz = 2*pi/Lk; Nr = max(Nkr,Nkz); Nz = Nr; r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r); z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z); [ZZ,RR] = meshgrid(z,r); %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nr,Nz,Ns2D); ni00 = zeros(Nr,Nz,Ns2D); si00 = zeros(Nr,Nz,Ns5D); phi = zeros(Nr,Nz,Ns2D); drphi = zeros(Nr,Nz,Ns2D); dzphi = zeros(Nr,Nz,Ns2D); for it = 1:numel(Ts2D) NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz))); ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz))); phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz))); drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz))); dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz))); end if strcmp(OUTPUTS.write_non_lin,'.true.') for it = 1:numel(Ts5D) si00(:,:,it) = real(fftshift(ifft2(squeeze(Si00(:,:,it)),Nr,Nz))); end end % Post processing disp('- post processing') E_pot = zeros(1,Ns2D); % Potential energy n^2 E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2 ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi| Flux_ri = zeros(1,Ns2D); % Particle flux Gamma = Flux_zi = zeros(1,Ns2D); % Particle flux Gamma = Flux_re = zeros(1,Ns2D); % Particle flux Gamma = Flux_ze = zeros(1,Ns2D); % Particle flux Gamma = Ne_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj Ni_norm = zeros(Npi,Nji,Ns5D);% . if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj Si_norm = zeros(Npi,Nji,Ns5D);% . Sne00_norm = zeros(1,Ns2D); % Time evol. of the amp of e nonlin term Sni00_norm = zeros(1,Ns2D); % end Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2; for it = 1:numel(Ts2D) % Loop over 2D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); E_pot(it) = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id E_kin(it) = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr; ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz)))); Flux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; Flux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz; Flux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz; end E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2); dEdt = diff(E_pot+E_kin)./dt2D; for it = 1:numel(Ts5D) % Loop over 5D arrays NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it); Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz; if strcmp(OUTPUTS.write_non_lin,'.true.') Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz; Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz; Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz; Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz; end end %% Compute growth rate disp('- growth rate') tend = Ts2D(end); tstart = 0.6*tend; g_ = zeros(Nkr,Nkz); [~,ikr0KH] = min(abs(kr-KR0KH)); for ikr = 1:Nkr for ikz = 1:Nkz g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend); end end % gkr0kz_Ni00 = max(real(g_(:,:)),[],1); gkr0kz_Ni00 = real(g_(ikr0KH,:)); %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 0 %% Time evolutions fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)]; set(gcf, 'Position', [100, 100, 900, 800]) subplot(221); for ip = 1:Npe for ij = 1:Nje plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$']; semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); subplot(222) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$']; semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); subplot(223) plot(Ts2D,Flux_ri,'-','DisplayName','$\Gamma_{ri}$'); hold on; plot(Ts2D,Flux_zi,'-','DisplayName','$\Gamma_{zi}$'); hold on; plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$') plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$') grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma$'); %legend('show'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224) for ip = 1:Npi for ij = 1:Nji plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$S_i^{',num2str(ip-1),num2str(ij-1),'}$']; semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on; end end grid on; xlabel('$t c_s/R$'); ylabel('$\sum_{k_r,k_z}|S_i^{pj}|$'); %legend('show'); else %% Growth rate subplot(224) [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz(1:Nz/2),gkr0kz_Ni00(1:Nz/2),... 'DisplayName',['J = ',num2str(JMAXI)]); xlabel('$k_z$'); ylabel('$\gamma_{Ni}$'); xlim([0. 1.0]); %ylim([0. 0.04]); end suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]); save_figure end %% if 0 %% Photomaton : real space tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['photo_real',sprintf('_t=%.0f',Ts2D(it))]; set(gcf, 'Position', [100, 100, 1500, 500]) subplot(131); plt = @(x) (((x))); pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$'); legend('$n_i$'); subplot(132); plt = @(x) ((x)); DATA = plt(ni00(:,:,it))-plt(ne00(:,:,it)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); legend('$n_i-n_e$'); set(gca,'ytick',[]); subplot(133); plt = @(x) ((x)); DATA = plt(phi(:,:,it)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); set(gca,'ytick',[]); legend('$\phi$'); % if strcmp(OUTPUTS.write_non_lin,'.true.') % subplot(133); plt = @(x) fftshift((abs(x)),2); % pclr = pcolor((RR),(ZZ),plt(si00(:,:,it5D))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) % xlabel('$r/\rho_s$'); legend('$|S_i^{00}|$'); set(gca,'ytick',[]) % end suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB), sprintf(', $t c_s/R=%.0f$',Ts2D(it))]); save_figure end if 1 %% Photomaton : real space % FIELD = ni00; FNAME = 'ni'; FIELD = phi; FNAME = 'phi'; tf = 60; [~,it1] = min(abs(Ts2D-tf)); tf = 65; [~,it2] = min(abs(Ts2D-tf)); tf = 200; [~,it3] = min(abs(Ts2D-tf)); tf = 400; [~,it4] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position', [100, 100, 1500, 400]) subplot(141); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it1)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it1))); subplot(142); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it2)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); subplot(143); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it3)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it3))); subplot(144); plt = @(x) ((x)); DATA = plt(FIELD(:,:,it4)); pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1]) xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); title(sprintf('$t c_s/R=%.0f$',Ts2D(it4))); % suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),... % ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']); save_figure end %% -t0 = 50; +t0 = 0; skip_ = 1; DELAY = 0.01*skip_; FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D); %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% Density ion GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_i$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Density electron GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Density ion - electron GIFNAME = ['ni-ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1; FIELD = real(ni00+ne00); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$n_i-n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$'; create_gif end if 0 %% Phi GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1; FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% Density ion frequency GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% Density ion frequency @ kr = 0 GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.'; X = (kz); T = Ts2D; YMIN = -.1; YMAX = 1.1; XMIN = min(kz); XMAX = max(kz); FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r\rho_s$'; create_gif_1D end %% if 1 %% Space time diagramm (fig 11 Ivanov 2020) t0 = 1; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D)); fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position', [100, 100, 1200, 400]) subplot(211) plot(Ts2D(it0:it1),Flux_ri(it0:it1)); ylabel('$\Gamma_r$'); grid on % title(['$\eta=',num2str(ETAB),'\quad',... % '\nu_{',CONAME,'}=',num2str(NU),'$']) % legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) ylim([-0.01,1.1*max(Flux_ri(it0:it1))]); subplot(212) [TY,TX] = meshgrid(r,Ts2D(it0:it1)); pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,it0:it1),2))'); set(pclr, 'edgecolor','none'); %colorbar; xlabel('$t c_s/R$'), ylabel('$r/\rho_s$') legend('$\langle\partial_r \phi\rangle_z$') save_figure end if 0 %% Flow skip = 10; plt = @(x) x(1:skip:end,1:skip:end); tf = 120; [~,it] = min(abs(Ts2D-tf)); fig = figure; FIGNAME = 'space_time_drphi'; quiver(plt(RR),plt(ZZ),plt(dzphi(:,:,it)),plt(-drphi(:,:,it)),0); xlabel('$r$'), ylabel('$z$') title(sprintf('$v_E$, $t=%.0f c_s/R$',Ts2D(it))); end if 0 %% Growth rate fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)]; [~,ikr0KH] = min(abs(kr-KR0KH)); plot(kz/kr(ikr0KH),gkr0kz_Ni00/(kr(ikr0KH)*A0KH),'DisplayName',['J = ',num2str(JMAXI)]) xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$'); xlim([0. 1.0]); ylim([0. 0.6]); save_figure end %% if 0 %% Mode time evolution [~,ik00] = min(abs(kz)); [~,idk] = min(abs(kz-dkz)); [~,ik50] = min(abs(kz-0.5*KR0KH)); [~,ik75] = min(abs(kz-0.7*KR0KH)); [~,ik10] = min(abs(kz-1.00*KR0KH)); plt = @(x) abs(squeeze(x)); fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)]; semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik00)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(idk)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik50)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik75)/KR0KH)]); hold on semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ... ['$k_z/k_{r0} = $',num2str(kz(ik10)/KR0KH)]); hold on xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show'); semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts2D-Ts2D(end))),'k--') semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts2D-Ts2D(end))),'k--') plot(tstart*[1 1],ylim,'k-','LineWidth',0.5); plot(tend*[1 1],ylim,'k-','LineWidth',0.5); save_figure end %% if 0 %% Show frame in kspace tf = 20; [~,it] = min(abs(Ts5D-tf)); fig = figure; FIGNAME = ['krkz_frame',sprintf('_%.2d',JOBNUM)]; subplot(221); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat\phi|$'); subplot(222); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_i^{00}|$'); subplot(223); plt = @(x) fftshift(abs(x),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_e^{00}|$'); if strcmp(OUTPUTS.write_non_lin,'.true.') subplot(224); plt = @(x) fftshift((abs(x)),2); pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar; xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$'); end save_figure end diff --git a/wk/flux_results.m b/wk/flux_results.m index 1cefbb8..d1c37d7 100644 --- a/wk/flux_results.m +++ b/wk/flux_results.m @@ -1,92 +1,92 @@ default_plots_options if 0 %% Compute time average and std of the mean flow t0 = 180; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D)); range = it0:it1; avg = mean(Flux_ri(range)) stdev = std(Flux_ri(range))^(.5) figure hist(Flux_ri(range),20); xlabel('$\Gamma$') end %% Handwritten results for nu = 0.01 % High definition results (256x128) Results_256x128.Gamma = [0.02, 0.03, 0.20, 0.037, 2.7, 2.25, 4, 5e-4, 2e-3, 0.03]; Results_256x128.L = [ 66, 66, 66, 50, 66, 66, 66, 66, 66, 66]; Results_256x128.P = [ 2, 3, 4, 5, 2, 3, 4, 2, 3, 4]; Results_256x128.J = [ 1, 2, 2, 3, 1, 2, 2, 1, 2, 2]; Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6]; Results_256x128.mrkr = [ 'v', '>', '^', 'o', 'v', '>', '^', 'v', '>', '^']; Results_256x128.clr = [ 'k', 'k', 'k', 'r', 'r', 'r', 'r', 'k', 'k', 'k']; % Low definition results (128x64) % Results_128x64.Gamma = [0.29, 0.05, 7e-4, 0.31, 3.7, 2e-3]; % Results_128x64.L = [ 25, 25, 25, 33 33, 33]; % Results_128x64.P = [ 2, 2, 2, 2, 2, 2]; % Results_128x64.J = [ 1, 1, 1, 1, 1, 1]; % Results_128x64.NU = [0.01, 0.1, 0.01, 0.01, 0.01, 0.01]; % Results_128x64.etaB = [ 0.5, 0.5, 0.67, 0.5 0.4, 0.6]; % Results_128x64.mrkr = [ 's', 's', 's', 's', 's', 's']; % Results_128x64.clr = [ 'b', 'b', 'b', 'r', 'r', 'r']; % Ricci_Rogers.Gamma = [2.5 1 1e-2]; % Ricci_Rogers.etaB = [0.4 0.5 1.0]; Ricci_Rogers.Gamma = [10 1e-1]; Ricci_Rogers.etaB = [0.5 1.0]; if 1 % Fig 3 of Ricci Rogers 2006 fig = figure; semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6); hold on; res = Results_256x128; for i = 1:numel(res.Gamma) semilogy(res.etaB(i),res.Gamma(i),... res.mrkr(i),'DisplayName','256x128', 'color', res.clr(i)); hold on; end % res = Results_128x64; % for i = 1:numel(res.Gamma) % if res.NU(i) == 0.01 % semilogy(res.etaB(i),res.Gamma(i),... % res.mrkr(i),'DisplayName','128x64', 'color', res.clr(i)); % end % hold on; % end xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$') end grid on; title('$\nu = 0.01$') legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$','$P=5$, $J=3$') -FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png']; +FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png']; ylim([1e-4, 10]) saveas(fig,FIGNAME); disp(['Figure saved @ : ',FIGNAME]) %% Handwritten results for nu = 0.1 Results_256x128.Gamma = [0.026,0.026, 1e-2, 1, 1, 1, 2e-2, 1, 0.15, 3e-3]; Results_256x128.P = [ 2, 3, 4, 2, 3, 4, 2, 3, 4, 4]; Results_256x128.J = [ 1, 2, 2, 1, 2, 2, 1, 2, 2, 2]; Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.7]; Results_256x128.mrkr = [ 'v', '>', '^', 'v', '>', '^', 'v', '>', '^', '^']; Results_256x128.clr = [ 'k', 'k', 'k', 'b', 'b', 'b', 'r', 'b', 'r', 'r']; % Ricci_Rogers.Gamma = [2 1e-1]; % Ricci_Rogers.etaB = [0.5 1.0]; Ricci_Rogers.Gamma = [10 1e-1]; Ricci_Rogers.etaB = [0.5 1.0]; if 1 % Fig 3 of Ricci Rogers 2006 fig = figure; semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6); hold on; res = Results_256x128; for i = 1:numel(res.Gamma) semilogy(res.etaB(i),res.Gamma(i),... res.mrkr(i),'DisplayName','256x128', 'color', res.clr(i)); hold on; end xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$') end grid on; title('$\nu = 0.1$') legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$') FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png']; saveas(fig,FIGNAME); disp(['Figure saved @ : ',FIGNAME]) \ No newline at end of file diff --git a/wk/fort.90 b/wk/fort.90 index d823444..412c03c 100644 --- a/wk/fort.90 +++ b/wk/fort.90 @@ -1,64 +1,64 @@ &BASIC nrun = 100000000 - dt = 0.004 - tmax = 400 - RESTART = .true. + dt = 0.01 + tmax = 1 + RESTART = .false. / &GRID - pmaxe =5 - jmaxe = 3 - pmaxi = 5 - jmaxi = 3 - Nr = 256 - Lr = 66 - Nz = 256 - Lz = 66 + pmaxe =2 + jmaxe = 1 + pmaxi = 2 + jmaxi = 1 + Nr = 64 + Lr = 20 + Nz = 64 + Lz = 20 kpar = 0 CANCEL_ODD_P = .false. / &OUTPUT_PAR - nsave_0d = 250 + nsave_0d = 1 nsave_1d = 0 - nsave_2d = 250 - nsave_5d = 2500 + nsave_2d = 1 + nsave_5d = 1 write_Na00 = .true. write_moments = .true. write_phi = .true. write_non_lin = .true. write_doubleprecision = .true. - resfile0 = '../results/ZP/256x128_L_66_Pe_5_Je_3_Pi_5_Ji_3_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/outputs' - rstfile0 = '../results/ZP/256x128_L_66_Pe_5_Je_3_Pi_5_Ji_3_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/checkpoint' + resfile0 = '../results/test_parallel/64x32_L_20_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/outputs' + rstfile0 = '../results/test_parallel/64x32_L_20_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-06/checkpoint' job2load = 0 / &MODEL_PAR ! Collisionality - CO = -1 + CO = -2 DK = .false. NON_LIN = .true. - mu = 0.0005 - nu = 0.1 + mu = 5e-06 + nu = 0.01 tau_e = 1 tau_i = 1 sigma_e = 0.023338 sigma_i = 1 q_e = -1 q_i = 1 eta_n = 1 eta_T = 0 eta_B = 0.5 lambdaD = 0 kr0KH = 0 A0KH = 0 / &INITIAL_CON only_Na00 =.false. initback_moments =0 initnoise_moments =1e-05 iseed =42 - selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10.h5' - eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10_tau_1.0000_mu_0.0233.h5' - iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_5_Jmaxe_3_Pmaxi_5_Jmaxi_3_pamaxx_10_tau_1.0000_mu_0.0233.h5' + selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5' + eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' + iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5' / &TIME_INTEGRATION_PAR numerical_scheme='RK4' / \ No newline at end of file