main.scmc

;paravec.scmc dm_fv.scmc general_macros.scmc load_gapsio.scmc checkpoint_data_needed.scmc main_iter_macro.scmc

(include< "stdio.h")
(include< "stdlib.h")
(include< "stdint.h")
(include< "string.h")
(include< "math.h")
(include< "time.h")
(include< "assert.h")
(eval-scmc-global (begin (load "pscmc_config_runtime.ss") (load "main_iter_macro.scmc") '()))
(input-scmc "paravec.scmc")
(input-scmc "general_macros.scmc")
(input-scmc "load_gapsio.scmc")
(include- "genrand.h")
(eval-scmc-global (begin (load "checkpoint_data_needed.scmc") '()))
(input-all-pscmc-struct)
(input-all-kernel-and-rt)

(include- "space_filling_curve.h")
(include- "mpifields.h")
(include- "dmshell.h")
(include- "hydroAshell.h")
(include- "init_field3d_mpi.h")
(include< "cgapsio.h")
(include- "mpi_fieldio.h")
(include- "call_curl_kernel.h")
(include- "init_particle.h")
(include- "run_particle.h")
(include- "run_particle_call_fun.h")
(include- "pass_xyzzyx.h")
(include- "sort_particle.h")
(include- "seqfields.h")
(include- "init_adjoint_relation.h")
(include- "init_implicit_particle.h")
(include- "split_shell.h")
(include- "blas_shell.h")
(include- "nonlinear_jfnk.h")
(include- "cfgcst.h")
(typedef Field3D_MPI Field3D_MPI_ALL)

(include- "checkpoint.h")
(define-scmc-global add_the_prefix (lambda (name)
  (concat 'c_ name))
  )
(defmacro sync_ovlp_with-boundary-core (field bound)
  `(begin
     (sync_ovlp_mpi_field ,field)
     (blas_axpy_full_block_Field3D_MPI ,field ,field 1 ,bound)
     )
  )
(define-scmc-global boundary-map (fast-make-single-env-from-var-and-val '((prho_s_vx1 . pboundary_rho_s_vx) (prho_s_vx0 . pboundary_rho_s_vx) (pvA1 . pboundary_vA) (pvA0 . pboundary_vA) (pvA2 . pboundary_vA) (palpha_beta1 . pboundary_alpha_beta1) (palpha_beta0 . pboundary_alpha_beta1))))
(defmacro sync_ovlp_with-boundary (field)
  (define fdvar (fast-find-var-in-single-env field boundary-map))
  (if (eq? fdvar 'VAR-NOT-BOUND)
    (begin
      (write-string (multi-concat "Error: field " field " is unbounded\n"))
      (car 0)
      )
    `(sync_ovlp_with-boundary-core ,field ,fdvar)
    )
  )
(defmacro macro-add-the-prefix (name)
  (add_the_prefix name)
  )
(defmacro gen_cst_specs cstlst
  (begin-map
    (lambda (cst)
      (multi-define name type initfun cst)
      `(begin
	 (declare (,type *) (,name (malloc (* (sizeof ,type ) NUM_SPEC)) ))
	 (for-from-to i 0 NUM_SPEC
	   (vector-set! ,name i (,initfun i))
	   )
	 )
      )
    cstlst
    )
  )
(include< "sys/time.h")
(defun print_max int ((Field3D_MPI* in))
  (define-double tmp0 (blas_sum_Field3D_MPI in in))
  (LOG_RECORD_INFO "ext max=%e\n" tmp0)
  (return 0)
  )
(defmacro dump_field_txt (field name)
  `(block
     (define-FILE* fp (fopen ,name "w"))
     (dump_field ,field fp)
     (fclose fp)
     )
  )
(defmacro dump_field_m filename-field-pairs
  (append 
    (begin-map
      (lambda (x)
	(define filename (car x))
	(define fds (cdr x))
	(append 
	  `(block
	     (decl-var-and-pvar Gaps_IO_DataFile gid)
	     (init_parallel_file_for_mpi_fields ,(car fds) pgid ,filename -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	     )
	  (let loop ((n 0) (fds fds))
	    (cond
	      ((null? fds) '())
	      (else
		(cons 
		  `(mpi_field_write_to_file ,(car fds) pgid ,n) 
		  (loop (+ n 1) (cdr fds)))
		)
	      )
	    )
	  )
	) 
      filename-field-pairs
      )
    `(
       (PS_MPI_Barrier PS_MPI_COMM_WORLD)
       (PS_MPI_Finalize)
       (return 0)
       )
    )
  )
(defmacro gen_ips_info ()
  `(if (and (== rank 0) (eq? M_DISABLE_TS_LOG 0))
     (block
       (define tnow (wclk_now))
       (define tused (- tnow tend))
       (set! tend tnow)
       (LOG_RECORD_INFO "step=%d time used=%fs, ips=%f, allips=%f\n" t tused (/ 1 tused) (/ (+ t 1) (- tend tbeg)))
       ))
  )


(defun wclk_now double ()
  (declare (struct timeval) tv1)
  (declare (struct timezone) tz1)
  (gettimeofday ("&" tv1) ("&" tz1))
  (define-double time_start1 (+ (type-convert double tv1.tv_sec) (* 1e-6 (type-convert double tv1.tv_usec))))
  (return time_start1)
  )
(defun gen_sts_dt double ((int j) (int nt) (double nv))
  (return (runc "1/((-1+nv)*cos((2*j-1.0)/nt*M_PI/2)+1+nv)"))
  )
(eval-scmc-global 
  `(begin 
     ,(begin-map 
	(lambda (x)
	  (gen_loadgapsio x 'double '((long cursp) (double z) (double y) (double x) (long l))))
	'(TEMPERATURE_DIST)
	)
     ,(begin-map
	(lambda (x)
	  (gen_loadgapsio x 'double '((long cursp) (double z) (double y) (double x))))
	'(DENSITY_DIST V0_x V0_y V0_z)
	)
     ))

(input-scmc "dm_fv.scmc")
(defun hydro_bicg_fun int ((Field3D_MPI* out) (Field3D_MPI* in) (hydro_fv* fv))
  (sync_ovlp_mpi_field in)
  (MPI_hydro_push_jac_alpha_beta out in fv->palpha_beta1 fv->prho_s_vx fv->pvA1 fv->pvA1 fv->dt fv->dx fv->dy fv->dz fv->u0 fv->qm0)
  (if 0
    (block
      (class-header-Field3D_Seq out->data)
      (for-from-zero-to i 8
	;(fprintf stderr "r0%d=%e\n" i (GET_FIELD3D_SEQ r0h->data 0 i 0 0 0))
	;(fprintf stderr "r1%d=%e\n" i (GET_FIELD3D_SEQ r0h->data 0 i 0 0 1))
	(fprintf stderr "out%d=%e loc=0x%lx\n" i (GET_FIELD3D_SEQ out->data 0 i 0 0 0) ("&" (GET_FIELD3D_SEQ out->data 0 i 0 0 0)))
	(fprintf stderr "out%d=%e loc=0x%lx\n" i (GET_FIELD3D_SEQ out->data 0 i 0 0 1) ("&" (GET_FIELD3D_SEQ out->data 0 i 0 0 1)))
	)
      ))
  (return 0)
  )
(defun dm_bicg_fun int ((Field3D_MPI* out) (Field3D_MPI* in) (dm_fv* fv))
  (sync_ovlp_mpi_field in)
  (MPI_dm_1st_eqn_core out in fv->pF0 fv->pA0 fv->pA1 fv->pA2 fv->DM_A fv->Q fv->M fv->DT)
  (return 0)
  )
(define-int global_use_dm_core_type 0)
(defmacro xyzmax_ifn1x2 (x)
  `(if (== ,x 1) ,x (* ,x 2))
  )
(defun cayley_dm_new_fun_all int ((Field3D_MPI* out) (Field3D_MPI* in) (dm_fv* fv) (double lftrht))
  ;from \dot{y}=f to x1-dt f(x1/2) = x0 + dt f(x0/2)
  ;(LOG_RECORD_INFO "newfun called, lftrht=%e \n" lftrht)
  ;(LOG_RECORD_INFO "all:")
  ;(print_max in)
  (sync_ovlp_mpi_field in)
  (if (eq? global_use_dm_core_type 0)
    (MPI_dm_cayley_eqn_core out in (structp-ref fv pA0) (structp-ref fv pA1) (structp-ref fv DT) (structp-ref fv DZ) (structp-ref fv DY) (structp-ref fv DX) (structp-ref fv DM_A) (structp-ref fv Q) (structp-ref fv M))
    (MPI_dm_8x8_eqn_psi out in (structp-ref fv pA0) (structp-ref fv pA1) (structp-ref fv DT) (structp-ref fv DZ) (structp-ref fv DY) (structp-ref fv DX) (structp-ref fv DM_A) (structp-ref fv Q) (structp-ref fv M)))
  ;(LOG_RECORD_INFO "Run core:")
  ;(print_max in)
  ;(print_max out)
  (blas_axpby_Field3D_MPI out out 1 in (* lftrht fv->DT))
  ;(LOG_RECORD_INFO "added out:")
  ;(print_max out)
  (return 0)
  )
(defun cayley_dm_new_fun int ((Field3D_MPI* out) (Field3D_MPI* in) (dm_fv* fv))
  (return (cayley_dm_new_fun_all out in fv -0.5))
  )
(defun cayley_dm_new_fun_right int ((Field3D_MPI* out) (Field3D_MPI* in) (dm_fv* fv))
  (cayley_dm_new_fun_all out in fv 0.5)
  ;(print_max out)
  (return 0)
  )
(defun high_order_split_pic void ((void* ppis) (double dt) (int l))
  (cond
    ((< l 0) (FDTD_2_2th_ALL_passes ppis dt))
    ((== l 0) (FDTD_2_4th_ALL_passes ppis dt))
    (else
      (define-double al (runc "1.0/(2-pow(2,l/3.))"))
      (define-double bl (runc "1-2*al"))
      (define-double adt (* al dt))
      (define-double bdt (* bl dt))
      (high_order_split_pic ppis adt (- l 1))
      (high_order_split_pic ppis bdt (- l 1))
      (high_order_split_pic ppis adt (- l 1))
      )
    )
  )
(define-long srand_seed 0)
(defun main int ((int argc) (char** argv))
  (main_proc argc argv)
  (get_constants)
  (PS_MPI_Init ("&" argc) ("&" argv))
  (eval-scmc-global (init-global-gapsio-fun-vars))
  (define-int n NUM_PROCESS)
  (if (eq? NUM_PROCESS 0)
    (begin
      (PS_MPI_Comm_size PS_MPI_COMM_WORLD ("&" n))
      (set! NUM_PROCESS n)
      )
    )
  (define-long num_runtime NUM_RUNTIME)
  (define-long n_hilbert NUM_N_HILBERT)
  (define len_hilbert (shift-l 1 n_hilbert))
  (define numt 1)
  (for-from-to g 0 NUM_N_HILBERT_DIMENSION 
    (set! numt (* numt len_hilbert))
    )
  (define-int rank)
  (PS_MPI_Comm_rank PS_MPI_COMM_WORLD ("&" rank))
  (assert (< rank n))
  (gen-const-vars M_ (USE_TIME_AS_RANDOM_SEED RAND_SEED DELTA_X DELTA_Y DELTA_Z USE_VLO DISABLE_TS_LOG USE_NO_CURRENT USE_REDUCE_DIM USE_NOT_REDUCE_DIM_PARTICLE REDUCE_DIM_X_RAT REDUCE_DIM_Y_RAT REDUCE_DIM_Z_RAT REDUCE_DIM_RANDOM_RATE USE_SAME_RANDOM_SEED USE_HYDRO_DEBUG USE_REL_NCR USE_REL_2ND_SUBSTEP))
  (gen-const-vars G_ (USE_ABC_DIR USE_PEC_DIR USE_DAMP_DIR ))
  (define-int use_rel_ncr M_USE_REL_NCR)
  (dec-array double dims_rat_xyz 6)
  (set! (vrf dims_rat_xyz 0) M_REDUCE_DIM_X_RAT)
  (set! (vrf dims_rat_xyz 1) M_REDUCE_DIM_Y_RAT)
  (set! (vrf dims_rat_xyz 2) M_REDUCE_DIM_Z_RAT)
  (set! (vrf dims_rat_xyz 3) M_REDUCE_DIM_X_RAT)
  (set! (vrf dims_rat_xyz 4) M_REDUCE_DIM_Y_RAT)
  (set! (vrf dims_rat_xyz 5) M_REDUCE_DIM_Z_RAT) 
  (define-int reduce_dim M_USE_REDUCE_DIM)
  (define-double ONE_FORM_REDUCE_DIM_X_RAT M_REDUCE_DIM_X_RAT)
  (define-double ONE_FORM_REDUCE_DIM_Y_RAT M_REDUCE_DIM_Y_RAT)
  (define-double ONE_FORM_REDUCE_DIM_Z_RAT M_REDUCE_DIM_Z_RAT)
  (define-double TWO_FORM_REDUCE_DIM_X_RAT (* M_REDUCE_DIM_Y_RAT M_REDUCE_DIM_Z_RAT))
  (define-double TWO_FORM_REDUCE_DIM_Y_RAT (* M_REDUCE_DIM_X_RAT M_REDUCE_DIM_Z_RAT))
  (define-double TWO_FORM_REDUCE_DIM_Z_RAT (* M_REDUCE_DIM_X_RAT M_REDUCE_DIM_Y_RAT))
  (if (eq? M_DELTA_X 0) (set! M_DELTA_X 1))
  (if (eq? M_DELTA_Y 0) (set! M_DELTA_Y 1))
  (if (eq? M_DELTA_Z 0) (set! M_DELTA_Z 1))
  ;(srand (+ rank 1 (if M_USE_TIME_AS_RANDOM_SEED (time NULL) M_RAND_SEED)))
  (set! srand_seed (+ rank 1))
  (define-long local_seed (+ rank 1 (if M_USE_TIME_AS_RANDOM_SEED (time NULL) M_RAND_SEED)))
  ;(LOG_RECORD_INFO "musrs=%e\n" M_USE_SAME_RANDOM_SEED)
  (if (not M_USE_SAME_RANDOM_SEED) (srand srand_seed) (set! srand_seed local_seed))
  (dec-array long tids (* num_runtime n))
  (define-int64_t* local_tid_array (malloc (* (sizeof int64_t ) numt)))
  (decl-var-and-pvar Field3D_MPI testfield)
  (decl-var-and-pvar Field3D_MPI testfield_2x)
  ;(decl-var-and-pvar Field3D_MPI testfieldB)
  (decl-var-and-pvar Field3D_MPI testfieldSPEC)
  (define-Field3D_Seq fstest)
  (decl-var-and-pvar Field3D_Seq fstest_2x)
  (decl-var-and-pvar Field3D_Seq fstestSPEC)
  ;(decl-var-and-pvar Field3D_Seq fstestB)
  (define-Field3D_Seq* pfstest ("&" fstest))
  ;(init_Field3D_Seq pfstest ppe 4 3 1 cur_tid 2 3 0 len_hilbert len_hilbert 1 rank)
  (dec-array int cd_types NUM_MAX_RUNTIME)
  (dec-array int dev_ids NUM_MAX_RUNTIME)
  (dec-array int cd_performances NUM_MAX_RUNTIME)
  (gen-const-vars G_ (GAPSIO_VERSION GAPSIO_NUM_REDUCEWRITE OVERLAP_LEN USE_G_E USE_PML_ABC_DIR PML_LEVEL PML_SIGMA_MAX USE_DIFFERENT_DEV_PERFORMANCE))
  (if (< G_GAPSIO_VERSION 2) (set! G_GAPSIO_NUM_REDUCEWRITE 0))
  (block
    (define i 0)
    (for i=0 i<NUM_RUNTIME i++
      (vector-set! cd_types i (call_GET_DEV_TYPE i rank))
      (vector-set! dev_ids i (call_GET_DEV_ID i rank))
      (vector-set! cd_performances i (if G_USE_DIFFERENT_DEV_PERFORMANCE (call_CAL_FUN_ONE_PARA "GET_DEV_PERFORMANCE" i) 1))
      )
    )
  (LOG_RECORD_INFO "rank %d init, pid=%d\n" rank (getpid))
  ;(define-long cur_tid (get_cur_num_tid rank numt n tids))
  (memset pfstest 0 (sizeof Field3D_Seq))
  (memset pfstest_2x 0 (sizeof Field3D_Seq))
  (define fieldlen 3)
  (define-int overlap_len 2)
  (if USE_KGM 
    (begin
      (set! fieldlen 10)
      (set! overlap_len 1)
      )
    )
  (if (neq? G_OVERLAP_LEN 0)
    (set! overlap_len G_OVERLAP_LEN)
    )
  (if (eq? rank 0)
    (LOG_RECORD_INFO "overlap=%d\n" overlap_len))
  (dec-array long allxyzmax 3)
  (cond
    ((== NUM_N_HILBERT_DIMENSION 1)
      (dec-array long lengs 3)
      (for-from-to i 0 3 
	(vector-set! lengs i (if (== i HILBERT_DIR) len_hilbert 1))
	)
      (set! (vector-ref allxyzmax 0) (* (vector-ref lengs 0) XMAX))
      (set! (vector-ref allxyzmax 1) (* (vector-ref lengs 1) YMAX))
      (set! (vector-ref allxyzmax 2) (* (vector-ref lengs 2) ZMAX))
      (set_Field3D_Seq pfstest NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len fieldlen 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
      (if_rel_ncr
	(set_Field3D_Seq pfstest_2x NULL M_DELTA_X M_DELTA_Y M_DELTA_Z (xyzmax_ifn1x2 XMAX) (xyzmax_ifn1x2 YMAX) (xyzmax_ifn1x2 ZMAX) 2 overlap_len fieldlen 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
	)
      ;(set_Field3D_Seq pfstestB NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len 3 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
      )
    ((== NUM_N_HILBERT_DIMENSION 2)
      (dec-array long lengs 3)
      (for-from-to i 0 3 
	(vector-set! lengs i (if (== i HILBERT_DIR) 1 len_hilbert))
	)
      (set! (vector-ref allxyzmax 0) (* (vector-ref lengs 0) XMAX))
      (set! (vector-ref allxyzmax 1) (* (vector-ref lengs 1) YMAX))
      (set! (vector-ref allxyzmax 2) (* (vector-ref lengs 2) ZMAX))
      (set_Field3D_Seq pfstest NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len fieldlen 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
      (if_rel_ncr
	(set_Field3D_Seq pfstest_2x NULL M_DELTA_X M_DELTA_Y M_DELTA_Z (xyzmax_ifn1x2 XMAX) (xyzmax_ifn1x2 YMAX) (xyzmax_ifn1x2 ZMAX) 2 overlap_len fieldlen 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
	)
      ;(set_Field3D_Seq pfstestB NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len 3 0 (vector-ref lengs 0) (vector-ref lengs 1) (vector-ref lengs 2) rank)
      )
    ((== NUM_N_HILBERT_DIMENSION 3)
      (set_Field3D_Seq pfstest NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len fieldlen 0 len_hilbert len_hilbert len_hilbert rank)
      (if_rel_ncr
	(set_Field3D_Seq pfstest_2x NULL M_DELTA_X M_DELTA_Y M_DELTA_Z (xyzmax_ifn1x2 XMAX) (xyzmax_ifn1x2 YMAX) (xyzmax_ifn1x2 ZMAX) 2 overlap_len fieldlen 0 len_hilbert len_hilbert len_hilbert rank)
	)
      ;(set_Field3D_Seq pfstestB NULL M_DELTA_X M_DELTA_Y M_DELTA_Z XMAX YMAX ZMAX 2 overlap_len 3 0 len_hilbert len_hilbert len_hilbert rank)
      (set! (vector-ref allxyzmax 0) (* len_hilbert XMAX))
      (set! (vector-ref allxyzmax 1) (* len_hilbert YMAX))
      (set! (vector-ref allxyzmax 2) (* len_hilbert ZMAX))
      )
    )
  (set! fstestSPEC fstest)
  (set! (structp-ref pfstestSPEC num_ele) (* 7 NUM_SPEC))
  (init_Field3D_MPI_ALL ptestfield pfstest n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids cd_performances num_runtime PS_MPI_COMM_WORLD rank n)
  (if_rel_ncr
    (init_Field3D_MPI_ALL ptestfield_2x pfstest_2x n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids cd_performances num_runtime PS_MPI_COMM_WORLD rank n))
  ;(init_Field3D_MPI_ALL ptestfieldB pfstestB n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids num_runtime PS_MPI_COMM_WORLD rank n)
  (init_Field3D_MPI_from_new_num_ele ptestfieldSPEC ptestfield (* 7 NUM_SPEC))
  
  ;(init_Field3D_MPI_ALL ptestfieldSPEC pfstestSPEC n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids num_runtime PS_MPI_COMM_WORLD rank n)
    
  (define tbeg (wclk_now))
  (define tend tbeg)
  (cond
    ((eq? USE_KGM 1) 
      ;(define KGM_M0 (call_GET_VAR "KGM_M0"))
      (gen-const-vars M_ (KGM_M0 KGM_Q0 KGM_DT KGM_PHI0 KGM_AMPX KGM_SGM KGM_EY KGM_FRQ KGM_ASSEMBLE KGM_ASSEMBLE_TIME KGM_MID_LOC KGM_LEN_A0 KGM_SGM_DENS KGM_ASS_E_LOC0 KGM_G_BEG KGM_DX KGM_REFZ0 KGM_EXTG KGM_INIT_SPER KGM_XMID KGM_YMID KGM_ZMID KGM_AVRHO_DUMP_TIMESTEP))
      (decl-var-and-pvar Field3D_MPI F0)
      (decl-var-and-pvar Field3D_MPI F1)
      (decl-var-and-pvar Field3D_MPI F2)
      (decl-var-and-pvar Field3D_MPI extA0)
      (decl-var-and-pvar Field3D_MPI extA1)
      (decl-var-and-pvar Field3D_MPI extAtmp)
      (decl-var-and-pvar Field3D_MPI filter_A)
      ;(decl-var-and-pvar Field3D_MPI exte)
      ;(decl-var-and-pvar Field3D_MPI extb)
      (init_Field3D_MPI_from pF0 ptestfield)
      (init_Field3D_MPI_from pF1 ptestfield)
      ;(init_Field3D_MPI_from pextA0 ptestfieldB)
      ;(init_Field3D_MPI_from pextA1 ptestfieldB)
      (init_Field3D_MPI_from_new_num_ele pextA0 ptestfield 3)
      (init_Field3D_MPI_from_new_num_ele pextA1 ptestfield 3)
      (init_Field3D_MPI_from_new_num_ele pextAtmp ptestfield 3)
      (init_Field3D_MPI_from_new_num_ele pF2 pF0 1)
      (blas_yiszero_Field3D_MPI pF0 pF0)
      (sync_ovlp_mpi_field pF0)
      (blas_yiszero_Field3D_MPI pF1 pF1)
      (sync_ovlp_mpi_field pF1)
      (if M_KGM_ASSEMBLE
	(begin
	  (decl-var-and-pvar Gaps_IO_DataFile ipigs)
	  (GAPS_IO_InitIFile pipigs "tmpKGM_ASSEMBLE")
	  (define-long tread M_KGM_ASSEMBLE_TIME)
	  (read_parallel_file_for_mpi_fields pF0 pipigs tread)
	  (init_kgm_assemble pF0 M_KGM_AMPX (* M_KGM_EY M_KGM_DX) (/ M_KGM_DT M_KGM_DX) M_KGM_SGM M_KGM_FRQ M_KGM_ASS_E_LOC0)
	  )
	(begin
	  (init_kgm_global pF0 M_KGM_PHI0 M_KGM_M0 M_KGM_Q0 M_KGM_AMPX M_KGM_EY M_KGM_DT M_KGM_SGM M_KGM_FRQ M_KGM_MID_LOC M_KGM_LEN_A0 M_KGM_SGM_DENS M_KGM_INIT_SPER M_KGM_XMID M_KGM_YMID M_KGM_ZMID)
	  (init_external_field3d_without_ss_KGM pF0)
	  )
	)
      (decl-var-and-pvar Gaps_IO_DataFile gid)
      (decl-var-and-pvar Gaps_IO_DataFile agid)
      (define-long tsave 0)
      (define-long tsave_rho 0)
      (init_parallel_file_for_mpi_fields ptestfield pgid "tmpKGM" (if tsave tsave -1) G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (if USE_INIT_EXT_EB 
	(begin
	  ;(blas_yiszero_synced_Field3D_MPI pextA0 pextA0)
	  ;(blas_yiszero_synced_Field3D_MPI pextA1 pextA1)
	  (init_external_field3d_E_2d_extend pextA0 M_USE_REDUCE_DIM ONE_FORM_REDUCE_DIM_X_RAT ONE_FORM_REDUCE_DIM_Y_RAT ONE_FORM_REDUCE_DIM_Z_RAT)
	  (init_external_field3d_B_2d_extend pextA1 M_USE_REDUCE_DIM TWO_FORM_REDUCE_DIM_X_RAT TWO_FORM_REDUCE_DIM_Y_RAT TWO_FORM_REDUCE_DIM_Z_RAT)
	  (sync_ovlp_mpi_field pextA0)
	  (sync_ovlp_mpi_field pextA1)
	  )
	(begin
	  (blas_yiszero_synced_Field3D_MPI pextA0 pextA0)
	  (blas_yiszero_synced_Field3D_MPI pextA1 pextA1)
	  )
	)
      (if USE_FILTER
	(begin
	  (init_Field3D_MPI_from_new_num_ele pfilter_A pF0 3)
	  (init_external_field3d_FILTER_B pfilter_A)
	  )
	)
      (if M_KGM_AVRHO_DUMP_TIMESTEP
	(init_parallel_file_for_mpi_fields pF2 pagid "tmpRHO" (if tsave tsave -1) G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE))
      (for-from-to t 0 NUM_TIMESTEP
	(MPI_kgm_eqn_core pF1 pF0 pextA0 pextA1 M_KGM_DT M_KGM_M0 M_KGM_Q0 M_KGM_DX M_KGM_EXTG M_KGM_REFZ0 M_KGM_G_BEG (and t 1))
	;(blas_yisax_Field3D_MPI pextAtmp pextAtmp 1 pextA1)
	;(blas_yisax_Field3D_MPI pextA1 pextA1 2 pextA1)
	;(blas_axpy_Field3D_MPI pextA1 pextA1 -1 pextA0)
	;(blas_yisax_Field3D_MPI pextA0 pextA0 1 pextAtmp)
	(if USE_FILTER
	  (blas_mulxy_numele3_Field3D_MPI pF1 pF1 pfilter_A)
	  )
	(if (== rank 0) (LOG_RECORD_OUT "%d\n" t))
	(if (and 0 (== rank 0)) (begin (sync_main_data_d2h pF1) (incf! (GET_FIELD3D_SEQ pF1->data 0 4 0 0 0) (* M_KGM_AMPX (sin (* t DELTAT)))) (sync_main_data_h2d pF1)))
	(if USE_NP_BOUNDARY
	  (begin
	    ;(gen-const-vars G_ (USE_ABC_DIR USE_PEC_DIR USE_DAMP_DIR))
	    (MPI_Yee_FDTD_MUR_ABC pF1 pF0 M_KGM_DT G_USE_ABC_DIR G_USE_PEC_DIR G_USE_DAMP_DIR 0)
	    )
	  )
	(define-long rdmd (* (eq? M_KGM_ASSEMBLE 1) (type-convert long (* NUM_DUMP_TIMESTEP (/ (* (rand) 1.0) RAND_MAX)))))
	(if (and (type-convert int M_KGM_AVRHO_DUMP_TIMESTEP) (eq? 0 (remainder t (type-convert int M_KGM_AVRHO_DUMP_TIMESTEP))))
	  (begin
	    (MPI_kgm_calc_rho pF2 pF1 M_KGM_DT M_KGM_M0 M_KGM_Q0 M_KGM_DX M_KGM_REFZ0 0.1 0.3 0 0)
	    (mpi_field_write_to_file pF2 pagid tsave_rho)
	    (incf! tsave_rho)
	    (if (== rank 0)
	      (LOG_RECORD_INFO "%d rho done\n" tsave_rho)
	      )
	    )
	  )
	(if (eq? 0 (remainder (+ t rdmd) NUM_DUMP_TIMESTEP))
	  (begin
	    (mpi_field_write_to_file pF1 pgid tsave)
	    (incf! tsave)
	    (if (== rank 0)
	    (LOG_RECORD_INFO "%d done\n" t))
	    )
	  )
	(define-void* tmp pF0)
	(set! pF0 pF1)
	(set! pF1 tmp)
	)
      )
    ((eq? M_USE_HYDRO_DEBUG 1)
      (gen-const-vars M_ (HYDRO_QM0 HYDRO_U0 HYDRO_RHO0 HYDRO_DRHO0 HYDRO_B0 HYDRO_S0 HYDRO_BICGSTAB_ERR HYDRO_SOLVE_ERR))
      (decl-var-and-pvar Field3D_MPI rho_s_vx0)
      (decl-var-and-pvar Field3D_MPI boundary_rho_s_vx)
      (decl-var-and-pvar Field3D_MPI rho_s_vx1)
      (decl-var-and-pvar Field3D_MPI alpha_beta0)
      (decl-var-and-pvar Field3D_MPI alpha_beta1)
      (decl-var-and-pvar Field3D_MPI boundary_alpha_beta1)
      (decl-var-and-pvar Field3D_MPI rest_alpha_beta)
      (decl-var-and-pvar Field3D_MPI v_alpha_beta)
      (decl-var-and-pvar Field3D_MPI vA0)
      (decl-var-and-pvar Field3D_MPI vA1)
      (decl-var-and-pvar Field3D_MPI vA2)
      (decl-var-and-pvar Field3D_MPI boundary_vA)
      (set! (structp-ref pfstest num_ele) 5)
      (init_Field3D_MPI_ALL prho_s_vx0 pfstest n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids cd_performances num_runtime PS_MPI_COMM_WORLD rank n)
      (init_Field3D_MPI_from prho_s_vx1 prho_s_vx0)
      (init_Field3D_MPI_from pboundary_rho_s_vx prho_s_vx0)
      (init_Field3D_MPI_from_new_num_ele palpha_beta0 prho_s_vx0 2)
      (init_Field3D_MPI_from palpha_beta1 palpha_beta0)
      (init_Field3D_MPI_from pboundary_alpha_beta1 palpha_beta0)
      (init_Field3D_MPI_from prest_alpha_beta palpha_beta0)
      (init_Field3D_MPI_from pv_alpha_beta palpha_beta0)
      (init_Field3D_MPI_from_new_num_ele pvA0 prho_s_vx0 3)
      (init_Field3D_MPI_from pvA1 pvA0)
      (init_Field3D_MPI_from pvA2 pvA0)
      (init_Field3D_MPI_from pboundary_vA pvA0)
      
      ;init alpha beta
      ;(init_hydro_alpha_beta palpha_beta0 1e-2 1e-2)
      (init_hydro_alpha_beta palpha_beta0 0 0)
      (init_hydro_init_alpha_beta pboundary_alpha_beta1 (- M_HYDRO_B0) (* M_HYDRO_RHO0 M_HYDRO_QM0))
      (blas_axpy_Field3D_MPI palpha_beta0 palpha_beta0 1 pboundary_alpha_beta1)
      (blas_yisax_Field3D_MPI palpha_beta1 palpha_beta1 1 palpha_beta0)
      (init_hydro_boundary_alpha_beta pboundary_alpha_beta1 (- M_HYDRO_B0) (* M_HYDRO_RHO0 M_HYDRO_QM0))

      ;init s0
      (init_hydro_init_s0 pboundary_rho_s_vx 1 1)
      (init_hydro_rho_s_vx prho_s_vx0 M_HYDRO_RHO0 M_HYDRO_DRHO0 M_HYDRO_S0 0 0 0)
      (blas_axpy_full_block_Field3D_MPI prho_s_vx0 prho_s_vx0 1 pboundary_rho_s_vx)
      (init_hydro_boundary_s0 pboundary_rho_s_vx 1 1)
      ;init vA
      (init_hydro_init_A0y pvA0 M_HYDRO_B0)
      (init_hydro_init_A0y pvA1 M_HYDRO_B0)
      (init_hydro_init_A0y pvA2 M_HYDRO_B0)
      (init_hydro_boundary_A0y pboundary_vA M_HYDRO_B0)
      (sync_ovlp_mpi_field pvA0)
      (sync_ovlp_mpi_field pvA1)
      (sync_ovlp_mpi_field pvA2)
      (blas_axpy_full_block_Field3D_MPI pvA0 pvA0 1 pboundary_vA)
      (blas_axpy_full_block_Field3D_MPI pvA1 pvA1 1 pboundary_vA)
      (blas_axpy_full_block_Field3D_MPI pvA2 pvA2 1 pboundary_vA)
      ;(blas_yiszero_synced_Field3D_MPI pvA0 pvA0)
      ;(blas_yiszero_synced_Field3D_MPI pvA1 pvA1)
      ;(blas_yiszero_synced_Field3D_MPI pvA2 pvA2)

      (decl-var-and-pvar Gaps_IO_DataFile gid)
      (decl-var-and-pvar Gaps_IO_DataFile gaid)
      (decl-var-and-pvar Gaps_IO_DataFile grid)
      (init_parallel_file_for_mpi_fields prho_s_vx1 pgrid "tmpRHOSV" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (init_parallel_file_for_mpi_fields pvA0 pgaid "tmpA" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (init_parallel_file_for_mpi_fields palpha_beta0 pgid "tmpALPHABETA" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      ;calculate v0
      (block
	(blas_yisax_Field3D_MPI prho_s_vx1 prho_s_vx1 1 prho_s_vx0)
	(sync_ovlp_with-boundary prho_s_vx0)
	(sync_ovlp_with-boundary prho_s_vx1)
	(sync_ovlp_with-boundary palpha_beta1)
	(if 0
	  (begin
	    (dump_field_txt pvA1 "Aall_txt")
	    (dump_field_txt prho_s_vx1 "rhoall_bef_txt")
	    ;(blas_yiszero_synced_Field3D_MPI pvA1 pvA1)
	    (MPI_hydro_push_vx prho_s_vx0 palpha_beta1 palpha_beta1 prho_s_vx1 pvA1 pvA1 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	    (dump_field_txt prho_s_vx0 "rhoall_txt")
	    (dump_field_txt pboundary_rho_s_vx "rhobound_txt")
	    (dump_field_m ("tmpRSV_dbg" prho_s_vx0 prho_s_vx1) ("tmpA_dbg" pvA1) ("tmpAB_dbg" palpha_beta1))
	    )
	  (MPI_hydro_push_vx prho_s_vx0 palpha_beta1 palpha_beta1 prho_s_vx1 pvA1 pvA1 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	  )
	)
      ;newton_method

      (block
	(decl-var-and-pvar bicg_space bs)
	(decl-var-and-pvar hydro_fv fv) ;hydro_fv is defined in dm_fv.scmc
	(eval-scmc-global 
	  (begin-map
	    (lambda (xy)
	      (multi-define x y xy)
	      `(set! (structp-ref pfv ,x) ,y)
	      )
	    `((palpha_beta1 palpha_beta1)
	      (prho_s_vx prho_s_vx1)
	      (pvA1 pvA1)
	      (dt DELTAT)
	      (dx M_DELTA_X)
	      (dy M_DELTA_Y)
	      (dz M_DELTA_Z)
	      (u0 M_HYDRO_U0)
	      (qm0 M_HYDRO_QM0)
	      )
	    ))
	(mpi_init_bicg pbs hydro_bicg_fun palpha_beta1 1000 M_HYDRO_BICGSTAB_ERR pfv)
	;(init_hydro_ext_s0 pext_rho_s_vx 1 1)
	;(init_hydro_ext_A0y pvAext M_HYDRO_B0)


	;(init_external_field3d_without_ss_EXT_FLUID_RHO_S_VX0 prho_s_vx0)

	(for-from-zero-to t NUM_TIMESTEP
	  (if (eq? 0 (remainder t NUM_DUMP_TIMESTEP))
	    (begin
	      (define-long ts (/ t NUM_DUMP_TIMESTEP))
	      (mpi_field_write_to_file prho_s_vx1 pgrid ts)
	      (mpi_field_write_to_file palpha_beta1 pgid ts)
	      (mpi_field_write_to_file pvA1 pgaid (* ts 2))
	      (mpi_field_write_to_file pvA2 pgaid (+ (* ts 2) 1))
	      (if (== rank 0) (LOG_RECORD_INFO "output fields done\n"))
	      )
	    )
	  (block
	    (define-double rem_beta)
	    (if (eq? rank 0)
	      (begin
		(sync_field_d2h palpha_beta0)
		(set! rem_beta (GET_FIELD3D_SEQ palpha_beta0->data 0 0 0 0 1))
		)
	      )
	    (PS_MPI_Bcast ("&" rem_beta) 1 PS_MPI_DOUBLE 0 PS_MPI_COMM_WORLD)
	    (LOG_RECORD_INFO "rem_beta=%e\n" rem_beta)
	    (init_hydro_alpha_beta palpha_beta1 0 rem_beta)
	    (blas_axpy_Field3D_MPI palpha_beta0 palpha_beta0 -1 palpha_beta1)
	    ;(blas_yisconst_Field3D_MPI palpha_beta1->data)
	    (set_hydro_s_0 prho_s_vx0)
	    (set_hydro_s_0 prho_s_vx1)
	    )
	  (sync_ovlp_with-boundary prho_s_vx0)
	  ;(dump_field_txt pboundary_rho_s_vx "rhoall_bef_txt")
	  (sync_ovlp_with-boundary pvA1)
	  ;(dump_field_txt pvA1 "Aall_txt")
	  (sync_ovlp_with-boundary palpha_beta0)
	  ;Doc: Eulerian symplectic
	  ;push rho and s
	  ;(dump_field_txt prho_s_vx0 "rhoall_bef_txt")
	  (MPI_hydro_push_rho_s prho_s_vx1 palpha_beta0 palpha_beta0 prho_s_vx0 pvA1 pvA1 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	  ;(dump_field_txt prho_s_vx1 "rhoall_txt")
	  (blas_yisax_Field3D_MPI palpha_beta1 palpha_beta1 1 palpha_beta0)
	  (sync_ovlp_with-boundary palpha_beta1)
	  ;(sync_ovlp_with-boundary palpha_beta0)

	  ;(sync_ovlp_mpi_field palpha_beta1)
	  ;(blas_axpby_Field3D_MPI palpha_beta1 palpha_beta1 1 pboundary_alpha_beta1)
	  (sync_ovlp_with-boundary prho_s_vx1)
	  ;(sync_ovlp_mpi_field prho_s_vx1)
	  ;(dump_field_txt prho_s_vx0 "rhoall_bef_txt")
	  ;(blas_axpy_full_block_Field3D_MPI prho_s_vx1 prho_s_vx1 1 pboundary_rho_s_vx)
	  ;push alpha and beta
	  (for-from-zero-to z 100
	    (sync_ovlp_with-boundary palpha_beta1)
	    (MPI_hydro_push_alpha_beta prest_alpha_beta palpha_beta0 palpha_beta1 prho_s_vx1 pvA1 pvA1 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	    (if 0
	      (block 
		(decl-var-and-pvar Gaps_IO_DataFile gid)
		(decl-var-and-pvar Gaps_IO_DataFile grid)
		;(LOG_RECORD_INFO "OK here\n")
		(init_parallel_file_for_mpi_fields prho_s_vx1 pgrid "tmpRHO_dbg" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
		;(dump_field_txt prho_s_vx1 "rhoall_txt")
		;(LOG_RECORD_INFO "OK aft init rho\n")
		(init_parallel_file_for_mpi_fields prest_alpha_beta pgid "tmpField_dbg" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
		(mpi_field_write_to_file prho_s_vx1 pgrid 0)
		(mpi_field_write_to_file prho_s_vx0 pgrid 1)
		(mpi_field_write_to_file prest_alpha_beta pgid 0)
		(mpi_field_write_to_file palpha_beta1 pgid 1)
		(mpi_field_write_to_file palpha_beta0 pgid 2)
		(mpi_field_write_to_file pboundary_alpha_beta1 pgid 3)
		(PS_MPI_Barrier PS_MPI_COMM_WORLD)
		(PS_MPI_Finalize)
		(return 0)
		))
	    ;(LOG_RECORD_INFO "bef findmax\n")
	    (define-double newton_err (blas_findmax_Field3D_MPI prest_alpha_beta prest_alpha_beta))
	    (LOG_RECORD_INFO "newton_err z=%d rest_ab=%e\n" z newton_err)
	    (if (< newton_err M_HYDRO_SOLVE_ERR) break)
	    (mpi_simple_bicgstab pbs pv_alpha_beta prest_alpha_beta)
	    (blas_axpy_Field3D_MPI palpha_beta1 palpha_beta1 -1 pv_alpha_beta)
	    )
	  (sync_ovlp_with-boundary palpha_beta1)
	  (sync_ovlp_with-boundary prho_s_vx1)
	  ;push vx
	  (if 0 
	    (begin
	      (dump_field_txt prho_s_vx1 "rhoall_txt")
	      (dump_field_txt palpha_beta1 "alpbeta_txt")
	      (dump_field_txt pvA1 "Aall_txt")
	      (MPI_Abort MPI_COMM_WORLD 0)
	      ))
	  (MPI_hydro_push_vx prho_s_vx1 palpha_beta1 palpha_beta1 prho_s_vx1 pvA1 pvA1 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	  ;push vA
	  (sync_ovlp_with-boundary prho_s_vx1)
	  (MPI_hydro_push_vA pvA2 palpha_beta1 palpha_beta1 prho_s_vx1 pvA1 pvA0 DELTAT M_DELTA_X M_DELTA_Y M_DELTA_Z M_HYDRO_U0 M_HYDRO_QM0)
	  ;(LOG_RECORD_INFO "t=%d done\n" t)
	  (gen_ips_info)
	  (if 0 
	    (block 
	      (decl-var-and-pvar Gaps_IO_DataFile gid)
	      (decl-var-and-pvar Gaps_IO_DataFile grid)
	      (LOG_RECORD_INFO "OK here\n")
	      (init_parallel_file_for_mpi_fields prho_s_vx1 pgrid "tmpRHO_dbg" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	      (LOG_RECORD_INFO "OK aft init rho\n")
	      (init_parallel_file_for_mpi_fields prest_alpha_beta pgid "tmpField_dbg" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	      (mpi_field_write_to_file prho_s_vx1 pgrid 0)
	      (mpi_field_write_to_file prho_s_vx0 pgrid 1)
	      (mpi_field_write_to_file prest_alpha_beta pgid 0)
	      (mpi_field_write_to_file palpha_beta1 pgid 1)
	      (mpi_field_write_to_file palpha_beta0 pgid 2)
	      (PS_MPI_Barrier PS_MPI_COMM_WORLD)
	      (PS_MPI_Finalize)
	      (return 0)
	      ))
	  (block
	    (blas_yisax_Field3D_MPI pvA0 pvA0 1 pvA1)
	    (sync_ovlp_with-boundary pvA0)
	    (blas_yisax_Field3D_MPI pvA1 pvA1 1 pvA2)
	    (sync_ovlp_with-boundary pvA1)
	    (blas_yisax_Field3D_MPI palpha_beta0 palpha_beta0 1 palpha_beta1)
	    (sync_ovlp_with-boundary palpha_beta0)
	    (blas_yisax_Field3D_MPI prho_s_vx0 prho_s_vx0 1 prho_s_vx1)
	    (sync_ovlp_with-boundary prho_s_vx0)
	    )
	  )
	)
      )
    ((eq? USE_DM 1)
      (gen-const-vars M_ (DM_M0 DM_Q0 DM_LAMBDA DM_DT DM_PHI0 DM_PZ_I DM_PZ_R DM_Z_OFFSET DM_AMPX DM_SGM DM_EY DM_YZ DM_FRQ DM_SOLVE_ERR DM_A DM_PHI34V DM_PZ_R_GAUSS_RAND DM_USE_DISP_INIT_CONDITION DM_DISP_RAND_AMP DM_DUAL_FIELD DM_AMP_PSI DM_NUMP DM_USE_8x8 DM_USE_SINGLE_PSI))
      (if (eq? M_DM_LAMBDA 0) (set! M_DM_LAMBDA 1))
      (if M_DM_USE_8x8 (set! global_use_dm_core_type 1))
      (decl-var-and-pvar Field3D_MPI F0)
      (decl-var-and-pvar Field3D_MPI F1)
      (decl-var-and-pvar Field3D_MPI Ff0)
      (decl-var-and-pvar Field3D_MPI Ff1)
      (decl-var-and-pvar Field3D_MPI b_AF0)
      (decl-var-and-pvar Field3D_MPI b_AFf0)
      (decl-var-and-pvar Field3D_MPI A0)
      (decl-var-and-pvar Field3D_MPI A1)
      (decl-var-and-pvar Field3D_MPI A2)
      ;(decl-var-and-pvar Field3D_MPI Atmp)
      (set! (structp-ref pfstest num_ele) (if M_DM_USE_8x8 16 8))
      (init_Field3D_MPI_ALL pb_AF0 pfstest n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids cd_performances num_runtime PS_MPI_COMM_WORLD rank n)
      (init_Field3D_MPI_from pF0 pb_AF0)
      (init_Field3D_MPI_from pF1 pb_AF0)
      (if M_DM_DUAL_FIELD
	(begin
	  (init_Field3D_MPI_from pb_AFf0 pb_AF0)
	  (init_Field3D_MPI_from pFf0 pb_AF0)
	  (init_Field3D_MPI_from pFf1 pb_AF0)))
      (set! (structp-ref pfstest num_ele) 3)
      (init_Field3D_MPI_ALL ptestfield pfstest n_hilbert NUM_N_HILBERT_DIMENSION 0 tids local_tid_array cd_types dev_ids cd_performances num_runtime PS_MPI_COMM_WORLD rank n)
      (init_Field3D_MPI_from pA0 ptestfield)
      (init_Field3D_MPI_from pA1 ptestfield)
      (init_Field3D_MPI_from pA2 ptestfield)
      ;(init_Field3D_MPI_from pAtmp ptestfield)
      (decl-var-and-pvar bicg_space bs)
      (decl-var-and-pvar bicg_space bsf)
      (decl-var-and-pvar dm_fv fv)
      (set! (structp-ref pfv Q) M_DM_Q0)
      (set! (structp-ref pfv M) M_DM_M0)
      (set! (structp-ref pfv DT) M_DM_DT)
      (set! (structp-ref pfv DM_A) M_DM_A)
      (set! (structp-ref pfv DX) M_DELTA_X)
      (set! (structp-ref pfv DY) M_DELTA_Y)
      (set! (structp-ref pfv DZ) M_DELTA_Z)


      (decl-var-and-pvar Gaps_IO_DataFile gid)
      (decl-var-and-pvar Gaps_IO_DataFile ghd)
      (decl-var-and-pvar Gaps_IO_DataFile gad)
      (define-long tsave 0)
      (init_parallel_file_for_mpi_fields pF0 pgid "tmpDM" (if tsave tsave -1) G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (init_parallel_file_for_mpi_fields pA0 pghd "tmpH" (if tsave tsave -1) G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (init_parallel_file_for_mpi_fields pA0 pgad "tmpA" (if tsave tsave -1) G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
      (define-double gauss_rand_tmp (maxwell_dist 0 1))
      (if M_DM_USE_DISP_INIT_CONDITION
	(init_dm_phi_global_rand pF0 M_DM_PHI0 M_DM_DISP_RAND_AMP M_DM_USE_SINGLE_PSI)
	(init_dm_phi_global pF0 M_DM_M0 (* M_DM_PZ_R (if M_DM_PZ_R_GAUSS_RAND gauss_rand_tmp 1)) M_DM_PZ_I M_DM_Z_OFFSET M_DM_LAMBDA M_DM_A (sqrt (- 1 (* M_DM_A M_DM_A))) 1 M_DELTA_Z 0 M_DM_PHI34V 0 M_DM_AMP_PSI))
      (if M_DM_DUAL_FIELD
	(begin
	  ;(init_dm_phi_global pFf0 M_DM_M0 (* M_DM_PZ_R (if M_DM_PZ_R_GAUSS_RAND gauss_rand_tmp 1)) M_DM_PZ_I M_DM_Z_OFFSET M_DM_LAMBDA M_DM_A (sqrt (- 1 (* M_DM_A M_DM_A))) 1 M_DELTA_Z 0 M_DM_PHI34V M_PI M_DM_AMP_PSI)
	  (if M_DM_USE_DISP_INIT_CONDITION
	    (init_dm_phi_global_rand pFf0 M_DM_PHI0 M_DM_DISP_RAND_AMP M_DM_USE_SINGLE_PSI)
	    (init_dm_dual_phi_global pF0 pFf0 M_DM_AMP_PSI M_DM_SGM M_DM_NUMP M_DM_M0 M_DELTA_Z (vector-ref allxyzmax 2)))
	  (init_external_field3d_without_ss_DMf pFf0)
	  )
      )
      (init_external_field3d_without_ss_DM pF0)
      ;(write_parallel_file_for_mpi_fields pF0 pgid 0)
      ;(exit 0)
      (define-int use_old_unstable_alg 0)
      (init_dm_A0_global pA0 0 0)
      (init_dm_A1_global pA1 M_DM_AMPX (if use_old_unstable_alg M_DM_EY M_DM_YZ) (if use_old_unstable_alg M_DM_DT 1))
      ;if not use_old_unstable_alg, then pA0 is A, pA1 is Y, pA2 is J

      (if use_old_unstable_alg
	(mpi_init_bicg pbs dm_bicg_fun pF0 1000 M_DM_SOLVE_ERR pfv)
	(begin
	  (mpi_init_bicg pbs cayley_dm_new_fun pF0 1000 M_DM_SOLVE_ERR pfv)
	  (if M_DM_DUAL_FIELD
	    (mpi_init_bicg pbsf cayley_dm_new_fun pFf0 1000 M_DM_SOLVE_ERR pfv)
	    )
	  )
	)

      (define-double beg_time_dm (wclk_now))
      (define-double beg_time_dm_lst beg_time_dm)
      (for-from-to t 0 NUM_TIMESTEP 
	(eval-scmc-global 
	  (begin-map 
	    (lambda (x)
	      `(set! (structp-ref pfv ,x) ,x)
	      )
	    '(pA0 pA1 pA2 pF0)
	    )
	  )
	;(init_dm_phi_global pF0 M_DM_PHI0 0 0 0 0 0 0 0 M_DM_SGM M_DM_FRQ)
	(sync_ovlp_mpi_field pF0)
	(if M_DM_DUAL_FIELD
	  (sync_ovlp_mpi_field pFf0)
	  )
	(sync_ovlp_mpi_field pA0)
	(sync_ovlp_mpi_field pA1)
	(if (eq? 0 (remainder t NUM_DUMP_TIMESTEP))
	  (begin
	    (if M_DM_DUAL_FIELD
	      (begin
		(mpi_field_write_to_file pF0 pgid (* 2 tsave))
		(mpi_field_write_to_file pFf0 pgid (+ 1 (* 2 tsave)))
		(MPI_dm_calc_hamt_dual pFf0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
		)
	      (begin
		(mpi_field_write_to_file pF0 pgid tsave)
		(if M_DM_USE_8x8
		  (MPI_dm_8x8_cal_H pb_AF0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
		  (MPI_dm_calc_hamt pb_AF0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0))
		)
	      )
	    (mpi_field_write_to_file pA2 pghd tsave)
	    (mpi_field_write_to_file pA0 pgad (* tsave 2))
	    (mpi_field_write_to_file pA1 pgad (+ (* tsave 2) 1))
	    (incf! tsave)
	    (define-double tnow (wclk_now))
	    (define-double tall (- tnow beg_time_dm))
	    (define-double tonestep (- tnow beg_time_dm_lst))
	    (LOG_RECORD_INFO "%d done, ips=%f, ipsall=%f\n" t (/ NUM_DUMP_TIMESTEP tonestep) (/ t tall))
	    (set! beg_time_dm_lst tnow)
	    )
	  )

	(cond
	  (use_old_unstable_alg 
	    (MPI_dm_1st_eqn_right pb_AF0 pF1 pF0 pA0 pA1 pA2 M_DM_A M_DM_Q0 M_DM_M0 M_DM_DT)
	    (sync_ovlp_mpi_field pb_AF0)
	    (mpi_simple_bicgstab pbs pF1 pb_AF0)
	    (MPI_dm_1st_eqn_fdtd pb_AF0 pF1 pF0 pA0 pA1 pA2 M_DM_A M_DM_Q0 M_DM_M0 M_DM_DT)
	    )
	  (else
	    ;(sync_ovlp_mpi_field pA0)
	    (cayley_dm_new_fun_right pb_AF0 pF0 pfv)
	    (sync_ovlp_mpi_field pb_AF0)
	    (mpi_simple_bicgstab pbs pF1 pb_AF0)
	    (blas_axpby_Field3D_MPI pF0 pF0 0.5 pF1 0.5)
	    (blas_yiszero_synced_Field3D_MPI pA2 pA2)
	    ;from here pF0 is (pF1+pF0)/2
	    (if M_DM_DUAL_FIELD
	      (begin
		(set! (structp-ref pfv pF0) pFf0)
		(cayley_dm_new_fun_right pb_AFf0 pFf0 pfv)
		;(define-double tmp0 (blas_findmax_Field3D_MPI pb_AFf0 pb_AFf0))
		;(LOG_RECORD_INFO "pbAFf0 ext max=%e\n" tmp0)
		(sync_ovlp_mpi_field pb_AFf0)
		(mpi_simple_bicgstab pbsf pFf1 pb_AFf0)
		(blas_axpby_Field3D_MPI pFf0 pFf0 0.5 pFf1 0.5)
		(MPI_dm_bihamt_dual_psi_eqn_J pFf0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)

		)
	      (if M_DM_USE_8x8
		(MPI_dm_8x8_eqn_J pF0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
		(MPI_dm_cayley_eqn_J pF0 pF0 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0))
	      )
	    ;pA2 is J, pA1 is Y
	    (blas_axpy_Field3D_MPI pA1 pA1 M_DM_DT pA2)
	    ;above is \dot{F}={F,H^1_d}

	    (if M_DM_USE_8x8
	      (MPI_dm_8x8_eqn_psi_m pF0 pF1 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
	      (MPI_dm_exact_eqn_m pF0 pF1 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
	      )
	    (blas_yisax_Field3D_MPI pF1 pF1 1 pF0)
	    (if M_DM_DUAL_FIELD
	      (begin 
		(MPI_dm_exact_eqn_m pFf0 pFf1 pA0 pA2 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
		(blas_yisax_Field3D_MPI pFf1 pFf1 1 pFf0)
		)
	      )
	    ;\dot{F}={F,H^2_d}

	    (blas_axpy_Field3D_MPI pA0 pA0 M_DM_DT pA1)
	    ;\dot{F}={F,H^3_{dY}}

	    (sync_ovlp_mpi_field pA0)
	    (MPI_dm_bihamt_eqn_dydt pF0 pF1 pA0 pA1 M_DM_DT M_DELTA_Z M_DELTA_Y M_DELTA_X M_DM_A M_DM_Q0 M_DM_M0)
	    ;\dot{F}={F,H^3_{dA}}
	    )
	  )

	(define-void* ftmp)
	(if use_old_unstable_alg
	  (begin
	    (set! ftmp pF0)
	    (set! pF0 pF1)
	    (set! pF1 ftmp)
	    (set! ftmp pA0)
	    (set! pA0 pA1)
	    (set! pA1 pA2)
	    (set! pA2 ftmp)))

	)
      ;(MPI_dm_1st_eqn_core pF2 pF0 pF1 pA0 pA1 pA2 M_DM_Q0 M_DM_M0 M_DM_DT)
      ;(init_dm_global )
      )
    (else 
      (gen_cst_specs 
	(pnpm double call_GET_NPM)
	(pchg double call_GET_CHARGE)
	(pmass double call_GET_MASS)
	(pgcache long call_GET_GRID_CACHE_LEN)
	(pcucache long call_GET_CU_CACHE_LEN)
	)
      (decl-var-and-pvar Particle_in_Cell_MPI pis)
      (gen-const-vars M_ (USE_TYPE3_KERNEL USE_SMALL_NUM_GRIDS USE_USER_DEFINED_PARTICLE_DISTRIBUTION USE_INIT_PARTICLE_FILE USE_DUMP_PARTICLE_FILE USE_ITG_MODE ITG_ENE_CONS ITG_CONST_NE0 USE_SMALL_TIMESTEP_MODE USE_MULTI_STEP USE_MULTI_SUBSTEP USE_SPLIT_MULTI_SUBSTEP USE_EOUT_SCHEME ITG_GC USE_DEBUG NUM_MULTI_STEP_SORT USE_MIDP_IMPLICIT USE_IMPLICIT_2ND USE_FULLY_IMPLICIT_ELECTRON PARTICLE_SOLVE_ERR JFNK_SOLVE_ERR JFNK_BICGSTAB_SOLVE_ERR JFNK_EPSILON))
      (set! M_JFNK_EPSILON (if (eq? M_JFNK_EPSILON 0) 1e-7 M_JFNK_EPSILON))
      (set! M_JFNK_BICGSTAB_SOLVE_ERR (if (eq? M_JFNK_BICGSTAB_SOLVE_ERR 0) 1e-10 M_JFNK_BICGSTAB_SOLVE_ERR))
      (set! M_JFNK_SOLVE_ERR (if (eq? M_JFNK_SOLVE_ERR 0) 1e-7 M_JFNK_SOLVE_ERR))
      (define-long num_substep 4)
      (define-long num_sort_e_substep 4)
      (if (or M_USE_MULTI_SUBSTEP M_USE_SPLIT_MULTI_SUBSTEP)
	(begin
	  (gen-const-vars M_ (NUM_SUBSTEP SORTE_NUM_SUBSTEP))
	  (if M_NUM_SUBSTEP (set! num_substep M_NUM_SUBSTEP))
	  (if M_SORTE_NUM_SUBSTEP (set! num_sort_e_substep M_SORTE_NUM_SUBSTEP))
	  )
	)
      (define-int num_multi_step_sort (if (eq? M_NUM_MULTI_STEP_SORT 0) 1 M_NUM_MULTI_STEP_SORT))
      (assert (> num_multi_step_sort 0))

      (gen-const-vars M_ (ITG_B0 ITG_MIN_R0 ITG_MAJ_R0 ITG_Q0 ITG_ZMID) 1)
      (gen-const-vars G_ (USE_REL USE_NON_UNI_IO_FOR_EACH_SPEC))

      (if (and (neq? M_USE_ITG_MODE 0) (neq? USE_TORI 0)) 
	(gen-const-vars M_ (ITG_B0 ITG_MIN_R0 ITG_MAJ_R0 ITG_Q0 ITG_ZMID) 2))
      (gen-const-vars T_ (TORI_X0 TORI_SOLVE_ERR) 1)
      (if USE_TORI (gen-const-vars T_ (TORI_X0 TORI_SOLVE_ERR) 2))
      (gen-const-vars G_ (REL_SOLVE_ERR NR_SOLVE_ERR) 1)
      (if (or G_USE_REL use_rel_ncr) (gen-const-vars G_ (REL_SOLVE_ERR) 2))
      (if M_USE_IMPLICIT_2ND (gen-const-vars G_ (NR_SOLVE_ERR) 2))
      ;(LOG_RECORD_INFO "tx=%e\n" T_TORI_X0) (exit 0)
      (define-int* p_particle_type (TYPE_MALLOC int NUM_SPEC))
      (for-from-zero-to i NUM_SPEC
	(if (or M_USE_SMALL_TIMESTEP_MODE M_USE_MULTI_SUBSTEP M_USE_FULLY_IMPLICIT_ELECTRON)
	  (vector-set! p_particle_type i (call_CAL_FUN_ONE_PARA "GET_PARTICLE_TYPE" i))
	  (vector-set! p_particle_type i 0)
	  )
	)

      (init_global_particles ppis ptestfield ptestfield_2x ptestfieldSPEC M_USE_SMALL_NUM_GRIDS G_USE_PML_ABC_DIR G_PML_LEVEL G_PML_SIGMA_MAX DELTAT NUM_SPEC allxyzmax pmass pchg pnpm pgcache pcucache p_particle_type M_USE_REDUCE_DIM M_REDUCE_DIM_X_RAT M_REDUCE_DIM_Y_RAT M_REDUCE_DIM_Z_RAT M_REDUCE_DIM_RANDOM_RATE M_USE_VLO use_rel_ncr)
      (if (and M_USE_ITG_MODE USE_TORI)
	(sync_ovlp_mpi_field ("&" (structp-ref ppis MPI_fieldB_ext)))
	)
      (if M_USE_DEBUG
	(begin
	  (LOG_RECORD_INFO "dumping debug field\n")
	  (decl-var-and-pvar Gaps_IO_DataFile gid)
	  (init_parallel_file_for_mpi_fields ptestfield pgid "tmpField_dbg" -1 G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	  (mpi_field_write_to_file ("&" (structp-ref ppis MPI_fieldE)) pgid 0)
	  (mpi_field_write_to_file ("&" (structp-ref ppis MPI_fieldB)) pgid 1)
	  (mpi_field_write_to_file ("&" (structp-ref ppis MPI_fieldE_ext)) pgid 2)
	  (mpi_field_write_to_file ("&" (structp-ref ppis MPI_fieldB_ext)) pgid 3)
	  )
	)
      (define-jfnk_newton_space jns)
      (if M_USE_SMALL_TIMESTEP_MODE
	(begin
	  (init_implicit_particle_mpi ppis)
	  (if M_USE_MIDP_IMPLICIT
	    (init_jfnk_newton_space ("&" jns) ("&" ppis->MPI_fieldE1) nonlin_fun_midp_vlasov_maxwell M_JFNK_SOLVE_ERR M_JFNK_BICGSTAB_SOLVE_ERR 18 8 M_JFNK_EPSILON ppis)
	    (init_jfnk_newton_space ("&" jns) ("&" ppis->MPI_fieldE) nonlin_fun_cur_min_curlB M_JFNK_SOLVE_ERR M_JFNK_BICGSTAB_SOLVE_ERR 18 8 M_JFNK_EPSILON ppis)
	    ) ;newton_error bicg_error dx_epsl
	  ;(init_jfnk_newton_space ("&" jns) ("&" ppis->MPI_fieldE) one_step_calc_current 1e-12 1e-5 4 4 1e-6 ppis)
	  (jfnk_newton_init_parameters ("&" jns) 1 1 pmass pchg)
	  ;(if USE_INIT_EB0 (begin (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) 1 ("&" ppis->MPI_fieldE)) (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldB))))
	  )
	)
      (if M_USE_FULLY_IMPLICIT_ELECTRON
	(begin
	  (init_implicit_particle_mpi ppis)
	  (init_jfnk_newton_space ("&" jns) ("&" ppis->MPI_fieldE) nonlin_fun_shell_for_mask M_JFNK_SOLVE_ERR M_JFNK_BICGSTAB_SOLVE_ERR 18 8 M_JFNK_EPSILON ppis)
	  (jfnk_newton_init_parameters ("&" jns) 1 1 pmass pchg)
	  )
	)

      (if M_USE_ITG_MODE 
	;(blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldE_ext) ("&" ppis->MPI_fieldE_ext) -1 ("&" ppis->MPI_fieldE_ext))
	(blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) -1 ("&" ppis->MPI_fieldE)) ; B1 is -Te
	)
      (cond
	(0
	  (init_single_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)))
	  )
	(M_USE_USER_DEFINED_PARTICLE_DISTRIBUTION
	  (init_user_defined_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)))
	  (call_particle_sort_mpi ppis 0 0)
	  (call_particle_sort_mpi ppis 1 0)
	  (call_particle_sort_mpi ppis 2 0)
	  )
	(M_USE_INIT_PARTICLE_FILE
	  (decl-var-and-pvar Gaps_IO_DataFile gid_grid)
	  (decl-var-and-pvar Gaps_IO_DataFile gid_cu)
	  (init_parallel_file_particle_for_mpi_fields_V0 ("&" (structp-ref ppis MPI_fieldE)) pgid_grid pgid_cu pgcache pcucache "GRID_PARTICLE_file" "CU_PARTICLE_file" 0 1)
	  (LOG_RECORD_INFO "rank=%d, Loading particles ...\n" rank)
	  (if (and M_USE_REDUCE_DIM (not M_USE_NOT_REDUCE_DIM_PARTICLE))
	    (read_particle_parallel_file_for_mpi_fields_V0 ("&" (structp-ref ppis MPI_fieldE)) pgid_grid pgid_cu pgcache pcucache 0 M_USE_REDUCE_DIM M_REDUCE_DIM_X_RAT M_REDUCE_DIM_Y_RAT M_REDUCE_DIM_Z_RAT G_USE_NON_UNI_IO_FOR_EACH_SPEC)
	    (read_particle_parallel_file_for_mpi_fields_V0 ("&" (structp-ref ppis MPI_fieldE)) pgid_grid pgid_cu pgcache pcucache 0 0 1 1 1 G_USE_NON_UNI_IO_FOR_EACH_SPEC)
	    )
	  (GAPS_IO_DeleteDataInfo pgid_cu)
	  (GAPS_IO_DeleteDataInfo pgid_grid)
	  )
	(USE_NON_UNI_DENSITY 
	  (gen-const-vars M_ (INIT_VMAX))
	  (if (eq? M_INIT_VMAX 0) (set! M_INIT_VMAX 1))
	  (init_non_uni_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)) M_INIT_VMAX))
	;(USE_TORI (init_single_particle_fmpi ("&" (structp-ref ppis MPI_fieldE))))
	;(M_USE_SMALL_TIMESTEP_MODE (init_dist_particle_nogrid_fmpi ("&" (structp-ref ppis MPI_fieldE)) 1))
	(else
	  ;(init_single_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)))
	  ;(LOG_RECORD_INFO "OK here\n")
	  (init_maxwell_dist_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)))
	  (init_dist_particle_nogrid_fmpi ("&" (structp-ref ppis MPI_fieldE)) 1)
	  )
	)
      ;(init_single_particle_fmpi ("&" (structp-ref ppis MPI_fieldE)))
      ;(generate_random_particle_mpi ptestfield)
      (dec-array char fname1 100)
      (sprintf fname1 "dbg1%d" rank)
      ;(dump_particles ptestfield (fopen fname1 "w"))
      ;(call_particle_sort_mpi ptestfield 0)
      ;(call_particle_sort_mpi ptestfield 1)
      ;(call_particle_sort_mpi ptestfield 2)
      ;(sprintf fname1 "dbg2%d" rank)
      ;(dump_particles ptestfield (fopen fname1 "w"))
      ;(delete_particle_mpi ptestfield)
      (decl-var-and-pvar Gaps_IO_DataFile gid)
      (decl-var-and-pvar Gaps_IO_DataFile gid2x)
      (decl-var-and-pvar Gaps_IO_DataFile gidj)
      (decl-var-and-pvar Gaps_IO_DataFile gide)

      (decl-var-and-pvar Field3D_MPI aver_E)
      (decl-var-and-pvar Field3D_MPI aver_B)
      (decl-var-and-pvar Field3D_MPI aver_J)
      (decl-var-and-pvar Field3D_MPI aver_EN)
      (decl-var-and-pvar Field3D_MPI ext_krook)
      (decl-var-and-pvar Field3D_MPI ext_rho_dist)

      ;(decl-var-and-pvar Field3D_MPI ext_add_par)
      (init_Field3D_MPI_from paver_E ptestfield)
      (init_Field3D_MPI_from paver_B ptestfield)
      (init_Field3D_MPI_from paver_J ptestfield)
      (init_Field3D_MPI_from pext_rho_dist ptestfield)
      (init_Field3D_MPI_from_new_num_ele pext_krook ptestfield 6)
      ;(init_Field3D_MPI_from pext_add_par ptestfield)
      (init_Field3D_MPI_from paver_EN ptestfieldSPEC)
      (blas_yiszero_Field3D_MPI paver_E paver_E)
      (blas_yiszero_Field3D_MPI paver_B paver_B)
      (blas_yiszero_Field3D_MPI paver_J paver_J)
      (blas_yiszero_Field3D_MPI paver_EN paver_EN)

      (gen-const-vars G_ (CHECKPOINT_TIME CHECKPOINT_NUM_SAVE CHECKPOINT_NUM_REDUCEWRITE USE_KROOK_FIELD KROOK_MU_FREQ USE_BORIS_YEE KROOK_SPEED_LIMIT STS_N))
      (if (eq? G_CHECKPOINT_NUM_REDUCEWRITE 0) (set! G_CHECKPOINT_NUM_REDUCEWRITE 1))
      (define-int cur_proc rank)
      (define-int num_max_proc)
      (PS_MPI_Comm_size PS_MPI_COMM_WORLD ("&" num_max_proc))
      (define-int num_one_group G_CHECKPOINT_NUM_REDUCEWRITE)


      (gen-const-vars M_ (USE_AVER_OUTPUT USE_1ST USE_OUTPUT_PERFORMANCE AVER_NUM_REDUCE))
      (if (eq? M_AVER_NUM_REDUCE 0)
	(set! M_AVER_NUM_REDUCE 1)
	)
      (if (or G_KROOK_MU_FREQ G_USE_KROOK_FIELD)
	;(init_external_field3d_FILTER_KROOK pext_krook)
	(init_external_field3d_FILTER_KROOK_2d_extend pext_krook reduce_dim 1 1 1)
	)
      (define-double sts_nv 1e-5)
      (if G_STS_N
	(begin
	  (gen-const-vars G_ (STS_NV))
	  (set! sts_nv G_STS_NV)
	  )
	)
      ;(block (define-FILE* fpout (fopen "dbgEB" "w")) (test_set_mainland ptestfield 1.) ;(test_set_mainland ptestfieldB 0.) ;(set! (GET_FIELD3D_SEQ (structp-ref ptestfield data) 0 0 0 0 1) 1.0) (sync_main_data_h2d ptestfield) ;(sync_ovlp_mpi_field ptestfield) (class-header-Field3D_MPI ptestfield))

      (block
	(test_set_mainland ("&" ppis->MPI_FoutJ) 0.)
	(sync_main_data_h2d ("&" ppis->MPI_FoutJ))
	(define t 0)
	(define tsave 0)

	(define-int current_check 0)
	(if USE_CHECKPOINT 
	  (begin
	    (LOG_RECORD_INFO "rank=%d, Loading...\n" rank)
	    (define-FILE* fp (load_identity_state_and_alloc_file_pointer rank num_one_group))
	    (define-int ifexistfp fp)
	    (PS_MPI_Bcast ("&" ifexistfp) 1 PS_MPI_INT 0 PS_MPI_COMM_WORLD)
	    (if ifexistfp
	      (begin 
		(eval-scmc-global (gen-save-load-data-needed 'load))
		(if fp (fclose fp))
		(set! current_check 1)
		)
	      (begin
		(if (== rank 0)
		  (LOG_RECORD_INFO "Warning: unable to locate the saved checkpoints, will creating new one\n")
		  )
		)
	      )
	    (LOG_RECORD_INFO "rank=%d, Done\n" rank)
	    (block
	      (define-long tn_new)
	      (if (== rank 0)
		(PS_MPI_Bcast ("&" t) 1 PS_MPI_LONG_LONG 0 PS_MPI_COMM_WORLD)
		(begin
		  (PS_MPI_Bcast ("&" tn_new) 1 PS_MPI_LONG_LONG 0 PS_MPI_COMM_WORLD)
		  (if (neq? t tn_new)
		    (begin
		      (LOG_RECORD_INFO "Error: the checkpoints are not saved in the same time step! rank=%d t=%ld, t should be %ld" rank t tn_new)
		      (PS_MPI_Abort PS_MPI_COMM_WORLD -1)
		      )
		    )
		  )
		)
	      (PS_MPI_Barrier PS_MPI_COMM_WORLD)
	      )
	    ;(LOG_RECORD_INFO "tcs=%e %e\n" G_CHECKPOINT_TIME G_CHECKPOINT_NUM_SAVE)
	    (if G_CHECKPOINT_TIME
	      (set! t G_CHECKPOINT_TIME)
	      )
	    (if G_CHECKPOINT_NUM_SAVE
	      (set! tsave G_CHECKPOINT_NUM_SAVE)
	      )
	    ))
	(block
	  (define-long newtsave (if tsave tsave -1))
	  (LOG_RECORD_INFO "tsave=%d\n" tsave)
	  (init_parallel_file_for_mpi_fields ptestfield pgid "tmpEB" newtsave G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	  (if use_rel_ncr
	    (init_parallel_file_for_mpi_fields ptestfield_2x pgid2x "tmpEB2X" newtsave G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE))
	  (init_parallel_file_for_mpi_fields ptestfield pgidj "tmpJ" newtsave G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE)
	  (init_parallel_file_for_mpi_fields ptestfieldSPEC pgide "tmpEN" newtsave G_GAPSIO_VERSION G_GAPSIO_NUM_REDUCEWRITE) 
	  )

	(define tbeg (wclk_now))
	(define tend tbeg)
	(gen-const-vars G_ (USE_ECRH_INPUT NUM_HIGH_ORDER_L USE_MULTIPLE_PARTICLE_DUMP USE_FIXB_DIR USE_CALC_EQ_B USE_INNER_SPLIT NUM_J_SOURCE))
	(define-double external_J_FREQ (if G_NUM_J_SOURCE (call_CAL_FUN_ONE_PARA "GET_EXTERNAL_J_FREQUENCY" 0) 0))
	(define-double external_J_AMPLITUTE (if G_NUM_J_SOURCE (call_CAL_FUN_ONE_PARA "GET_EXTERNAL_J_AMPLITUTE" 0) 0))
	(decl-var-and-pvar Field3D_MPI external_J_SOURCE)
	(if G_NUM_J_SOURCE 
	  (begin 
	    (init_Field3D_MPI_from ("&" external_J_SOURCE) ptestfield)
	    (init_external_field3d_without_ss_J_SOURCE ("&" external_J_SOURCE))))
	(if G_USE_REL
	  (gen-const-vars T_ (TORI_X0 TORI_SOLVE_ERR) 2))
	(define-double curt_profile_only (wclk_now))
	(if M_USE_ITG_MODE
	  (init_external_field3d_DENSITY_DIST pext_rho_dist)
	  )
	(define-int ext_field_set 0)
	(define t0 t)
	(for 0 t<NUM_TIMESTEP t++
	  (if (and USE_CHECKPOINT (== 0 (remainder t NUM_CHECKPOINT_TIMESTEP)))
	    (begin
	      (LOG_RECORD_INFO "rank=%d, Saving...\n" rank)
	      (if (not current_check)
		(begin
		  (define-FILE* fp (save_identity_state_and_alloc_file_pointer rank num_one_group))
		  (eval-scmc-global (gen-save-load-data-needed 'save))
		  (if fp (fclose fp)))
		(set! current_check 0)
		)
	      (if M_USE_DUMP_PARTICLE_FILE
		(begin
		  (define-long numdump (/ t NUM_CHECKPOINT_TIMESTEP))
		  (if (not (type-convert int G_USE_MULTIPLE_PARTICLE_DUMP))
		    (set! numdump 0)
		    )
		  (decl-var-and-pvar Gaps_IO_DataFile gid_grid)
		  (decl-var-and-pvar Gaps_IO_DataFile gid_cu)
		  (init_parallel_file_particle_for_mpi_fields_V0 ("&" (structp-ref ppis MPI_fieldE)) pgid_grid pgid_cu pgcache pcucache "tmpGRID_PARTICLE_file" "tmpCU_PARTICLE_file" numdump (if numdump 0 0))
		  (block
		    (define-int msg 0)
		    (define-MPI_Status st)
		    (if (neq? rank 0)
		      (PS_MPI_Recv ("&" msg) 1 PS_MPI_INT (- rank 1) (- rank 1) PS_MPI_COMM_WORLD ("&" st))
		      )
		    (dump_particle_parallel_file_for_mpi_fields_V0 ("&" (structp-ref ppis MPI_fieldE)) pgid_grid pgid_cu pgcache pcucache numdump G_USE_NON_UNI_IO_FOR_EACH_SPEC)
		    (if (neq? rank (- NUM_PROCESS 1))
		      (PS_MPI_Send ("&" msg) 1 PS_MPI_INT (+ rank 1) rank PS_MPI_COMM_WORLD)
		      )
		    )
		  (LOG_RECORD_INFO "rank=%d, Saving particles ...\n" rank)
		  (GAPS_IO_DeleteDataInfo pgid_cu)
		  (GAPS_IO_DeleteDataInfo pgid_grid)
		  (if (eq? rank 0) 
		    (begin
		      (define-FILE* fp (fopen "CHECKPOINT_TIME.ss" "wb"))
		      (fprintf fp "(define CHECKPOINT_TIME %d)\n" t)
		      (fprintf fp "(define CHECKPOINT_NUM_SAVE %d)\n" tsave)
		      (fclose fp)
		      )
		    )
		  ;(exit 0)
		  )
		)
	      (LOG_RECORD_INFO "rank=%d, Done\n" rank)
	      )
	    )
	  (if (and 0 (== rank 0))
	    (begin
	      (LOG_RECORD_OUT "%d\n" t)
	      (sync_main_data_d2h ("&" ppis->MPI_fieldE))
	      (incf! (GET_FIELD3D_SEQ ppis->MPI_fieldE.data 0 7 7 0 0) (sin (* t DELTAT)))
	      (sync_main_data_h2d ("&" ppis->MPI_fieldE))
	      )
	    )
	  (define-long num_av_red M_AVER_NUM_REDUCE)
	  (if (and (== 0 (remainder t num_av_red)) M_USE_AVER_OUTPUT)
	    (begin
	      (blas_axpby_Field3D_MPI paver_E paver_E (/ M_AVER_NUM_REDUCE NUM_DUMP_TIMESTEP) ("&" ppis->MPI_fieldE) 1.)
	      (blas_axpby_Field3D_MPI paver_B paver_B (/ M_AVER_NUM_REDUCE NUM_DUMP_TIMESTEP) ("&" ppis->MPI_fieldB) 1.)
	      (blas_axpby_Field3D_MPI paver_J paver_J (/ M_AVER_NUM_REDUCE NUM_DUMP_TIMESTEP) ppis->pMPI_FoutJ 1.)
	      (blas_yiszero_Field3D_MPI ptestfieldSPEC ptestfieldSPEC)
	      (dump_ene_num_mpi ppis 0)
	      (sync_main_data_d2h ptestfieldSPEC)
	      (blas_axpby_Field3D_MPI paver_EN paver_EN (/ M_AVER_NUM_REDUCE NUM_DUMP_TIMESTEP) ptestfieldSPEC 1.)
	      )
	    )
	  ;(sync_main_data_d2h ptestfield)
	  (if (== 0 (remainder t NUM_DUMP_TIMESTEP))
	    (begin
	      (if (== rank 0) (LOG_RECORD_INFO "outputing..."))
	      (if M_USE_OUTPUT_PERFORMANCE
		(PS_MPI_Barrier PS_MPI_COMM_WORLD)
		)
	      (define-double tbeg (wclk_now))
	      (if M_USE_AVER_OUTPUT
		(begin
		  (sync_main_data_d2h paver_E)
		  (sync_main_data_d2h paver_B)
		  (sync_main_data_d2h paver_J)
		  (sync_main_data_d2h paver_EN)
		  (mpi_field_write_to_file paver_E pgid (* tsave 2))
		  (mpi_field_write_to_file paver_B pgid (+ (* tsave 2) 1))
		  (mpi_field_write_to_file paver_J pgidj tsave)
		  (mpi_field_write_to_file paver_EN pgide tsave)
		  (blas_yiszero_Field3D_MPI paver_E paver_E)
		  (blas_yiszero_Field3D_MPI paver_B paver_B)
		  (blas_yiszero_Field3D_MPI paver_J paver_J)
		  (blas_yiszero_Field3D_MPI paver_EN paver_EN)
		  )
		(begin
		  (sync_main_data_d2h ("&" ppis->MPI_fieldB))
		  (sync_main_data_d2h ("&" ppis->MPI_fieldE))
		  (sync_main_data_d2h ppis->pMPI_FoutJ)
		  (mpi_field_write_to_file ("&" ppis->MPI_fieldE) pgid (* tsave 2))
		  (mpi_field_write_to_file ("&" ppis->MPI_fieldB) pgid (+ (* tsave 2) 1))
		  (if use_rel_ncr
		    (begin
		      (mpi_field_write_to_file ("&" ppis->MPI_fieldE01_2x) pgid2x (* tsave 2))
		      (mpi_field_write_to_file ("&" ppis->MPI_fieldB0_2x) pgid2x (+ (* tsave 2) 1))))
		  (mpi_field_write_to_file ppis->pMPI_FoutJ pgidj tsave)
		  ;(define-double themax (blas_findmax_Field3D_MPI ppis->pMPI_FoutJ ppis->pMPI_FoutJ))
		  ;(define-double thesum (blas_sum_Field3D_MPI ppis->pMPI_FoutJ ppis->pMPI_FoutJ))
		  ;(define-double thedot (blas_dot_Field3D_MPI ppis->pMPI_FoutJ ppis->pMPI_FoutJ ("&" ppis->MPI_fieldE)))
		  (blas_yiszero_Field3D_MPI ptestfieldSPEC ptestfieldSPEC)
		  ;(LOG_RECORD_INFO "bef dump ok\n")
		  ;(PS_MPI_Barrier PS_MPI_COMM_WORLD)
		  (dump_ene_num_mpi ppis 0)
		  ;(PS_MPI_Barrier PS_MPI_COMM_WORLD)
		  ;(LOG_RECORD_INFO "%d m=%.15e d=%.15e d=%.15e\n" rank themax thedot thesum)
		  (sync_main_data_d2h ptestfieldSPEC)
		  (mpi_field_write_to_file ptestfieldSPEC pgide tsave)
		  )
		)
	      (GAPS_IO_FileFlush pgid)
	      (GAPS_IO_FileFlush pgidj)
	      (GAPS_IO_FileFlush pgide)
	      (if M_USE_OUTPUT_PERFORMANCE
		(PS_MPI_Barrier PS_MPI_COMM_WORLD)
		)
	      (set! tbeg (- (wclk_now) tbeg))
	      (if (== rank 0) (LOG_RECORD_INFO "done, time used=%fs\n" tbeg))
	      (incf! tsave))
	    )

	  (if USE_LHCD_INPUT
	    (begin
	      (gen-const-vars L_ (LHCD_FREQ LHCD_K LHCD_SOURCE_LOCATION LHCD_AMPX LHCD_AMPY))
	      (set_lhcd_fields ("&" ppis->MPI_fieldE) (* t L_LHCD_FREQ DELTAT) L_LHCD_K L_LHCD_SOURCE_LOCATION  L_LHCD_AMPX L_LHCD_AMPY)
	      )
	    )
	  (if (type-convert int G_USE_ECRH_INPUT)
	    (begin
	      (gen-const-vars E_ 
		(ECRH_FREQ ECRH_AMPX ECRH_AMPY ECRH_AMPZ 
		  ECRH_KY ECRH_KZ ECRH_LOC_X0 ECRH_LOC_Y0 ECRH_LOC_Y1 ECRH_LOC_Z0 ECRH_LOC_Z1 ECRH_IKT ECRH_SQUARE_NUM ECRH_PHASE_0))
	      (set_ecrh_fields (if M_USE_ITG_MODE ("&" ppis->MPI_fieldE_ext) ("&" ppis->MPI_fieldE)) (* t E_ECRH_FREQ DELTAT )
		E_ECRH_AMPX E_ECRH_AMPY E_ECRH_AMPZ E_ECRH_KY E_ECRH_KZ 
		E_ECRH_IKT E_ECRH_LOC_X0 E_ECRH_LOC_Y0 E_ECRH_LOC_Y1 
		E_ECRH_LOC_Z0 E_ECRH_LOC_Z1 E_ECRH_SQUARE_NUM E_ECRH_PHASE_0) 
	      )
	    )

	  ;;;;here comes the main iteration step
	  (block
	    (define-double outter_dt DELTAT)
	    (define-int sts_n G_STS_N)
	    (if (== sts_n 0)
	      (set! sts_n 1)
	      )
	    (define-double total_dt 0)
	    (for-from-to i 0 sts_n
	      (incf! total_dt (gen_sts_dt (+ i 1) sts_n sts_nv))
	      )
	    (define-double sts_basic_dt (/ outter_dt total_dt))
	    ;(define-long num_substep 4)
	    (define-long e_numstep 0)
	    (define-double dt_substep (/ DELTAT num_substep))
	    (for-from-to i 0 sts_n
	      (define-double DELTAT (* sts_basic_dt (gen_sts_dt (+ i 1) sts_n sts_nv)))
	      
	      (if G_NUM_J_SOURCE
		(begin
		  (define-double ampoff (sin (* external_J_FREQ DELTAT t)))
		  (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) (* -1 external_J_AMPLITUTE ampoff DELTAT) ("&" external_J_SOURCE))

		  )
		)

	      ;(LOG_RECORD_INFO "STS_DT=%e\n" DELTAT)
	      (eval-scmc-global 
		(gen-all-milp
		  `(
		     ((and (or USE_TORI G_USE_REL) (not M_USE_ITG_MODE))
		       (cond
			 (M_USE_NO_CURRENT
			   (if (and USE_INIT_EXT_EB (not ext_field_set))
			     (begin
			       (incf! ext_field_set)
			       (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) 1 ("&" ppis->MPI_fieldE_ext))
			       (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldB_ext))
			       (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) 1 ("&" ppis->MPI_fieldB_ext))
			       (sync_ovlp_mpi_field ("&" ppis->MPI_fieldE))
			       (sync_ovlp_mpi_field ("&" ppis->MPI_fieldB))
			       (sync_ovlp_mpi_field ("&" ppis->MPI_fieldB1))
			       )
			     )
			   )
			 (else
			   (comment "set B1=B")
			   (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) 1.0 ("&" ppis->MPI_fieldB))
			   (comment "PML or not")
			   (if G_USE_PML_ABC_DIR
			     (begin
			       (class-header-Particle_in_Cell_MPI ppis)
			       (class-header-Field3D_Seq (struct-ref MPI_fieldE data))
			       (sync_ovlp_mpi_field ("&" ppis->MPI_fieldPMLE))
			       (comment "Push PML version of FDTD dB/dt solver")
			       (MPI_PML_FDTD_CURL_FWD ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldPMLB) ("&" ppis->MPI_fieldPMLE) DELTAT 0 0 delta_x delta_y delta_z G_USE_PML_ABC_DIR G_PML_LEVEL 3 G_PML_SIGMA_MAX allxmax allymax allzmax) 
			       )
			     ;(MPI_Yee_FDTD_Curl_E ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldE) DELTAT))
			     (begin
			       (comment "Push normal version of FDTD dB/dt solver")
			     (MPI_YEE_CURL_R ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldE) DELTAT)))
			   (if G_USE_FIXB_DIR
			     (MPI_Yee_FDTD_MUR_ABC ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB) DELTAT 0 0 0 G_USE_FIXB_DIR)
			     )
			   (PERFORMANCE_TIME "CURLE")
			   (if (eq? USE_FILTER 2)
			     (begin
			       (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldEtmp1) ("&" ppis->MPI_fieldEtmp1) 1 ("&" ppis->MPI_fieldE))
			       (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldBtmp1) ("&" ppis->MPI_fieldBtmp1) 1 ("&" ppis->MPI_fieldB))
			       )
			     )
			   (if USE_FILTER
			     (begin
			       (blas_mulxy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE_filter))
			       ;(sync_ovlp_mpi_field ("&" ppis->MPI_fieldE))
			       (blas_mulxy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB_filter))
			       ;(sync_ovlp_mpi_field ("&" ppis->MPI_fieldB))
			       (if (eq? USE_FILTER 2)
				 (begin
				   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldEtmp1) ("&" ppis->MPI_fieldEtmp1) -1 ("&" ppis->MPI_fieldE))
				   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldBtmp1) ("&" ppis->MPI_fieldBtmp1) -1 ("&" ppis->MPI_fieldB))
				   )
				 )
			       )
			     )

			   (if USE_INIT_EXT_EB
			     (begin
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) 1 ("&" ppis->MPI_fieldE_ext))
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldB_ext))
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) 1 ("&" ppis->MPI_fieldB_ext)
				 )) 
			     )
			   (define-double tnow)
			   ;(if USE_PROFILE (begin (MPI_Barrier MPI_COMM_WORLD) (define-double now1 (wclk_now)) (fprintf "%fs used") (set! tnow now1)))
			   (sync_ovlp_mpi_field ("&" ppis->MPI_fieldE))
			   (sync_ovlp_mpi_field ("&" ppis->MPI_fieldB))
			   (sync_ovlp_mpi_field ("&" ppis->MPI_fieldB1))
			   (PERFORMANCE_TIME "sync E, B, B1")

			   ;(LOG_RECORD_INFO "OK here\n")
			   (if (not G_USE_CALC_EQ_B)
			     (blas_yiszero_synced_Field3D_MPI ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_FoutJ)))
			   ;(sync_ovlp_mpi_field ppis->pMPI_FoutJ)
			   (PERFORMANCE_TIME "sync FoutJ")

			   (define-FILE* fp)))

		       (if USE_TORI
			 (if (not (and G_USE_CALC_EQ_B (> t t0)))
			   (cond
			     ((not G_USE_REL)
			       (class-header-Particle_in_Cell_MPI ppis)
			       (define-Field3D_MPI* pMPI_fieldE ("&" MPI_fieldE))
			       (define-Field3D_MPI* pMPI_fieldB ("&" MPI_fieldB))
			       (define-double curt (wclk_now))
			       ;(if 1 (block (define-FILE* fp (fopen "dbg00" "wb")) (dump_particles ("&" ppis->MPI_fieldE) fp) (fclose fp)))
			       (comment "Push particle, use vlo version if USE_VLO is set to 1")
			       (if M_USE_VLO
				 (MPI_geo_nr_Bfield_pushJ_vlo ppis pMPI_fieldE pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) T_TORI_X0 0 0 0 0 0)
				 (MPI_geo_nr_Bfield_pushJ ppis pMPI_fieldE pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) T_TORI_X0 0 0 0 0 0)
				 )
			       )
			     ((eq? M_USE_SMALL_NUM_GRIDS 1)
			       (define-double curt (wclk_now))
			       (blas_yiszero_Field3D_MPI ("&" ppis->MPI_LFoutJ) ("&" ppis->MPI_LFoutJ))
			       ;(sync_ovlp_mpi_field ("&" ppis->MPI_LFoutJ))
			       (PERFORMANCE_TIME "sync ovlp")
			       (MPI_geo_rel_1st_fwd_small_grids ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_LFoutJ) pmass pchg DELTAT T_TORI_X0 T_TORI_SOLVE_ERR)
			       (PERFORMANCE_TIME "geo_rel_sg")
			       ;(define-double* foj (structp-ref (type-convert c_pscmc_mem* (structp-ref (struct-ref ppis->MPI_LFoutJ data) main_data)) d_data))
			       ;(LOG_RECORD_DEBUG "foj=0x%lx foj[27951]=%e\n" foj (vector-ref foj 27951))
			       ;(LOG_RECORD_DEBUG "foj[32343]=%e\n" (vector-ref foj 32343))
			       (MPI_merge_current ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_LFoutJ))
			       )
			     ((eq? M_USE_SMALL_NUM_GRIDS 2)
			       (define-double curt (wclk_now))
			       (blas_yiszero_synced_Field3D_MPI ("&" ppis->MPI_LFoutJ) ("&" ppis->MPI_LFoutJ))
			       ;(sync_ovlp_mpi_field ("&" ppis->MPI_LFoutJ))
			       (PERFORMANCE_TIME "sync ovlp")
			       (MPI_geo_rel_1st_fwd_sg2_small_grids ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_LFoutJ) pmass pchg DELTAT T_TORI_X0 T_TORI_SOLVE_ERR)
			       (PERFORMANCE_TIME "geo_rel")
			       (MPI_merge_current_2 ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_LFoutJ))
			       )
			     (else
			       (MPI_geo_rel_1st_fwd ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT T_TORI_X0 T_TORI_SOLVE_ERR))
			     ))
			 (cond
			   ((eq? M_USE_SMALL_NUM_GRIDS 1)
			     (blas_yiszero_Field3D_MPI ("&" ppis->MPI_LFoutJ) ("&" ppis->MPI_LFoutJ))
			     (PERFORMANCE_TIME "sync ovlp")
			     (MPI_relng_1st_small_grids ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_LFoutJ) pmass pchg DELTAT T_TORI_X0 T_TORI_SOLVE_ERR)
			     (PERFORMANCE_TIME "geo_rel")
			     (MPI_merge_current ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_LFoutJ))
			     )
			   ((eq? M_USE_SMALL_NUM_GRIDS 2)
			     (blas_yiszero_synced_Field3D_MPI ("&" ppis->MPI_LFoutJ) ("&" ppis->MPI_LFoutJ))
			     (PERFORMANCE_TIME "sync ovlp")
			     (MPI_relng_1st_sg2_small_grids ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_LFoutJ) pmass pchg DELTAT T_TORI_X0 T_TORI_SOLVE_ERR)
			     (PERFORMANCE_TIME "geo_rel")
			     (MPI_merge_current_2 ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_LFoutJ))
			     )
			   (M_USE_TYPE3_KERNEL 
			     ;(LOG_RECORD_INFO "TYPE3_KERNEL\n")
			     (MPI_relng_1st_goto ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT (* t T_TORI_X0) T_TORI_SOLVE_ERR)
			     )
			   (else
			     (MPI_relng_1st ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT (* t T_TORI_X0) T_TORI_SOLVE_ERR)
			     ))
			 )

		       ;(if 0 (begin (set! fp (fopen "dbg01" "wb")) (dump_particles ("&" ppis->MPI_fieldE) fp) (fclose fp) (exit 0)))
		       (PERFORMANCE_TIME "merge_current/relng/geo_rel")
		       ;(dump_particles ("&" ppis->MPI_fieldE) stderr)
		       (comment "sort particle if it is the sort time-step")
		       (if (eq? (remainder t num_multi_step_sort) 0)
			 (begin
			   (comment "sort x direction")
			   (call_particle_sort_mpi ppis 0 0)
			   (PERFORMANCE_TIME "sort x particle")
			   ;(LOG_RECORD_INFO "bef 0\n")
			   ;(exit 0)
			   ;(dump_particles ("&" ppis->MPI_fieldE) stderr)
			   (comment "sort y direction")
			   (call_particle_sort_mpi ppis 1 0)
			   (PERFORMANCE_TIME "sort y particle")
			   ;(LOG_RECORD_INFO "bef 2\n")
			   ;(dump_particles ("&" ppis->MPI_fieldE) stderr)
			   (comment "sort z direction")
			   (call_particle_sort_mpi ppis 2 0)
			   ;(LOG_RECORD_INFO "aft 2\n")
			   ;(dump_particles ("&" ppis->MPI_fieldE) stderr)
			   (PERFORMANCE_TIME "sort z particle")))
		       (if (not M_USE_NO_CURRENT)
			 (begin
			   (if (not (and G_USE_CALC_EQ_B (> t t0)))
			     (merge_ovlp_mpi_field ("&" ppis->MPI_FoutJ)))
			   (PERFORMANCE_TIME "merge ovlp")
			   ;(LOG_RECORD_INFO "end cal 1st rel, time used=%fs\n" (- (wclk_now) tnow))

			   (if USE_INIT_EXT_EB
			     (begin
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) -1 ("&" ppis->MPI_fieldE_ext))
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) -1 ("&" ppis->MPI_fieldB_ext))
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) -1 ("&" ppis->MPI_fieldB_ext)
				 ))
			     )
			   (PERFORMANCE_TIME "yisax E B B1")
			   (if (eq? USE_FILTER 2)
			     (begin
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) 1 ("&" ppis->MPI_fieldEtmp1))
			       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldBtmp1))
			       )
			     )
			   (comment "add current to Electric field E")
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) -1 ("&" ppis->MPI_FoutJ))
			   (PERFORMANCE_TIME "CAL DEDT=-J")
			   (if USE_NP_BOUNDARY
			     (begin
			       (define-Field3D_MPI* pMPI_fieldE ("&" ppis->MPI_fieldE))
			       (define-Field3D_MPI* pMPI_fieldB ("&" ppis->MPI_fieldB))
			       (define-Field3D_MPI* pMPI_fieldEtmp ("&" ppis->MPI_fieldEtmp))
			       (if 0
				 (begin
				   (sync_main_data_d2h pMPI_fieldE)
				   (incf! (GET_FIELD3D_SEQ pMPI_fieldE->data 0 63 0 63 0) (sin (* t DELTAT)))
				   (incf! (GET_FIELD3D_SEQ pMPI_fieldE->data 0 63 0 63 2) (sin (* t DELTAT)))
				   (sync_main_data_h2d pMPI_fieldE)))
			       (blas_yisax_Field3D_MPI pMPI_fieldEtmp pMPI_fieldEtmp 1. pMPI_fieldE)
			       (if USE_TORI
				 ;(MPI_GEO_YEE_CURL_L ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) T_TORI_X0 DELTAT)
				 (begin
				   (comment "Use tori version of FDTD dE/dt")
				   (MPI_GEO_YEE_CURL_L ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) T_TORI_X0 M_DELTA_X M_DELTA_Y M_DELTA_Z DELTAT))
				 ;(MPI_Yee_FDTD_Curl_B ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) DELTAT)
				 (if G_USE_PML_ABC_DIR
				   (begin
				     (class-header-Particle_in_Cell_MPI ppis)
				     ;(class-header-Field3D_Seq pMPI_fieldE->data)
				     (class-header-Field3D_Seq (struct-ref MPI_fieldE data))
				     (sync_ovlp_mpi_field ("&" ppis->MPI_fieldPMLB))
				     (comment "Use PML version of FDTD dE/dt")
				     (MPI_PML_FDTD_CURL_BWD ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldPMLE) ("&" ppis->MPI_fieldPMLB) DELTAT 0 0 delta_x delta_y delta_z G_USE_PML_ABC_DIR G_PML_LEVEL 3 G_PML_SIGMA_MAX allxmax allymax allzmax) 
				     )
				   (begin
				     (comment "Otherwise use normal version of FDTD dE/dt")
				   (MPI_RECT_YEE_CURL_L ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) M_DELTA_Z M_DELTA_Y M_DELTA_X DELTAT)))
				 )
			       (PERFORMANCE_TIME "CURLB")
			       (comment "Apply Mur absorbing bondary condition")
			       (MPI_Yee_FDTD_MUR_ABC pMPI_fieldE pMPI_fieldEtmp DELTAT G_USE_ABC_DIR G_USE_PEC_DIR G_USE_DAMP_DIR 0)
			       (PERFORMANCE_TIME "MURABC")
			       )
			     (begin 
			       (comment "Following are corresponding non-NP bondary condition version of FDTD dE/dt")
			       (if USE_TORI
				 (MPI_GEO_YEE_CURL_L ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) T_TORI_X0 M_DELTA_X M_DELTA_Y M_DELTA_Z DELTAT)
				 ;(MPI_Yee_FDTD_Curl_B ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) DELTAT)
				 (if G_USE_PML_ABC_DIR
				   (begin
				     (class-header-Particle_in_Cell_MPI ppis)
				     ;(class-header-Field3D_Seq pMPI_fieldE->data)
				     (class-header-Field3D_Seq (struct-ref MPI_fieldE data))
				     (sync_ovlp_mpi_field ("&" ppis->MPI_fieldPMLB))
				     (MPI_PML_FDTD_CURL_BWD ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldPMLE) ("&" ppis->MPI_fieldPMLB) DELTAT 0 0 delta_x delta_y delta_z G_USE_PML_ABC_DIR G_PML_LEVEL 3 G_PML_SIGMA_MAX allxmax allymax allzmax) 
				     )
				   (MPI_RECT_YEE_CURL_L ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB1) M_DELTA_Z M_DELTA_Y M_DELTA_X DELTAT))
				 )
			       (PERFORMANCE_TIME "CURLB"))
			     )
			   (comment "set B=B1")
			   (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldB1))))
		       )
		     (G_USE_BORIS_YEE
		       (sync_ovlp_mpi_field  ("&" ppis->MPI_fieldE))
		       (MPI_Yee_FDTD_Curl_E ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldE) (/ DELTAT 2))

		       (blas_yiszero_Field3D_MPI ("&" ppis->MPI_FoutJ) ("&" ppis->MPI_FoutJ))
		       (sync_ovlp_mpi_field ("&" ppis->MPI_FoutJ))
		       (if USE_INIT_EXT_EB
			 (begin
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) 1 ("&" ppis->MPI_fieldE_ext))
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) 1 ("&" ppis->MPI_fieldB_ext))
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) 1 ("&" ppis->MPI_fieldB_ext)
			     )) 
			 )
		       (sync_ovlp_mpi_field  ("&" ppis->MPI_fieldE))
		       (sync_ovlp_mpi_field  ("&" ppis->MPI_fieldB))
		       (MPI_boris_yee ppis ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldB) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT)
		       (if USE_INIT_EXT_EB
			 (begin
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) -1 ("&" ppis->MPI_fieldE_ext))
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldB) -1 ("&" ppis->MPI_fieldB_ext))
			   (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldB1) ("&" ppis->MPI_fieldB1) -1 ("&" ppis->MPI_fieldB_ext)
			     )) 
			 )
		       (merge_ovlp_mpi_field ("&" ppis->MPI_FoutJ))

		       (MPI_Yee_FDTD_Curl_E ("&" ppis->MPI_fieldB) ("&" ppis->MPI_fieldE) (/ DELTAT 2))
		       (call_particle_sort_mpi ppis 0 0)
		       (call_particle_sort_mpi ppis 1 0)
		       (call_particle_sort_mpi ppis 2 0)
		       (sync_ovlp_mpi_field  ("&" ppis->MPI_fieldB))
		       (pass_GeqB ppis DELTAT)
		       (blas_axpy_Field3D_MPI ("&" ppis->MPI_fieldE) ("&" ppis->MPI_fieldE) (- DELTAT) ("&" ppis->MPI_FoutJ))
		       )
		     (G_NUM_HIGH_ORDER_L
		       (high_order_split_pic ppis DELTAT G_NUM_HIGH_ORDER_L)
		       )
		     ;(G_USE_MULTI_STEP_LOW_FREQ_MODE)
		     (M_USE_1ST
		       (split_1st_all_passes ppis DELTAT)) 
		     (M_USE_ITG_MODE
		       ;(LOG_RECORD_INFO "tr=%d OK here\n" USE_TORI)
		       (cond
			 (USE_TORI
			   (class-header-Particle_in_Cell_MPI ppis)
			   (define-Field3D_MPI* pMPI_fieldE ("&" MPI_fieldE))
			   (define-Field3D_MPI* pMPI_fieldB ("&" MPI_fieldB))
			   ;(blas_yisconst_Field3D_MPI pMPI_fieldB pMPI_fieldB 0)
			   (if 1 
			     (begin
			       (blas_yiszero_synced_Field3D_MPI pMPI_FoutJ pMPI_FoutJ)
			       (calculate_rho_mpi ppis 0)
			       (merge_ovlp_mpi_field pMPI_FoutJ)
			       (if USE_NON_UNI_DENSITY
				 (blas_yisax_Field3D_MPI pMPI_fieldE pMPI_fieldE M_ITG_CONST_NE0 pext_rho_dist)
				 (blas_yisconst_Field3D_MPI pMPI_fieldE pMPI_fieldE M_ITG_CONST_NE0))
			       (blas_yisax_Field3D_MPI pMPI_fieldB pMPI_fieldB 1 pMPI_FoutJ)
			       (blas_get_ITG_Potential_Field3D_MPI pMPI_FoutJ pMPI_FoutJ ("&" MPI_fieldB1) pMPI_fieldE (/ 1.0 (vector-ref pchg 0))) ;B1 is -Te, E is n e0
			       (sync_ovlp_mpi_field pMPI_FoutJ)
			       (MPI_Yee_FDTD_Grad_FWD pMPI_fieldE pMPI_FoutJ 0)
			       (if USE_INIT_EXT_EB
				 (begin
				   ;(LOG_RECORD_INFO "OK Here\n")
				   ;(blas_yiszero_synced_Field3D_MPI pMPI_fieldE pMPI_fieldE)
				   (blas_axpy_Field3D_MPI ("&" MPI_fieldE ) ("&" MPI_fieldE ) 1. ("&" MPI_fieldE_ext))

				   ;(decl-var-and-pvar Gaps_IO_DataFile gad)
				   ;(init_parallel_file_for_mpi_fields ("&" MPI_fieldE) pgad "tmpEout" -1 0)
				   ;(mpi_field_write_to_file ("&" MPI_fieldE) pgad 0)
				   ;(exit 0)
				   )
				 )
			       (sync_ovlp_mpi_field pMPI_fieldE)
			       ()
			       )
			     (blas_yiszero_synced_Field3D_MPI pMPI_fieldE pMPI_fieldE)
			     )
			   ;(LOG_RECORD_INFO "ok bef nr tx=%e\n" T_TORI_X0)
			   (if 0
			     (if M_USE_VLO
			       (MPI_geo_nr_vlo ppis pMPI_fieldE pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg DELTAT T_TORI_X0 M_ITG_MAJ_R0 M_ITG_MIN_R0 M_ITG_Q0 M_ITG_B0 M_ITG_ZMID)
			       (MPI_geo_nr ppis pMPI_fieldE pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg DELTAT T_TORI_X0 M_ITG_MAJ_R0 M_ITG_MIN_R0 M_ITG_Q0 M_ITG_B0 M_ITG_ZMID)
			       ))
			   (if M_USE_VLO 
			     (MPI_geo_nr_Bfield_pushJ_vlo ppis pMPI_fieldE ("&" MPI_fieldB_ext) ("&" MPI_fieldB_ext) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) T_TORI_X0 M_ITG_MAJ_R0 M_ITG_MIN_R0 M_ITG_Q0 M_ITG_B0 M_ITG_ZMID)
			     (MPI_geo_nr_Bfield_pushJ ppis pMPI_fieldE ("&" MPI_fieldB_ext) ("&" MPI_fieldB_ext) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) T_TORI_X0 M_ITG_MAJ_R0 M_ITG_MIN_R0 M_ITG_Q0 M_ITG_B0 M_ITG_ZMID)
			     )
			   ;(LOG_RECORD_INFO "OK here\n")
			   (call_particle_sort_mpi ppis 0 0)
			   (call_particle_sort_mpi ppis 1 0)
			   (call_particle_sort_mpi ppis 2 0)
			   ;(LOG_RECORD_INFO "OK here 2\n")
			   (if USE_INIT_EXT_EB
			     (begin
			       (blas_axpy_Field3D_MPI ("&" MPI_fieldE ) ("&" MPI_fieldE ) -1. ("&" MPI_fieldE_ext))
			       )
			     )
			   )
			 (M_ITG_ENE_CONS
			   (ITG_split_2nd_all_passes_enecons ppis pext_rho_dist DELTAT (/ 1. (vector-ref pchg 0)) M_ITG_CONST_NE0 M_ITG_GC)
			   )
			 (else
			   (ITG_split_2nd_all_passes ppis pext_rho_dist DELTAT (/ 1. (vector-ref pchg 0)) M_ITG_CONST_NE0 M_ITG_GC)
			   )))
		     (M_USE_SMALL_TIMESTEP_MODE
		       ;(split_small_timestep ppis DELTAT pmass pchg 4 1 1)
		       ;(LOG_RECORD_INFO "OK here\n")
		       (if M_USE_MIDP_IMPLICIT
			 (begin
			   (blas_yisax_Field3D_MPI ("&" ppis->MPI_fieldE1) ("&" ppis->MPI_fieldE1) 1 ("&" ppis->MPI_fieldE))
			   (jfnk_newton_simple_mpi ("&" jns) ("&" ppis->MPI_fieldE1))
			   )
			 (begin
			   (blas_yiszero_Field3D_MPI ("&" ppis->MPI_fieldE_filter) ("&" ppis->MPI_fieldE_filter))
			   (jfnk_newton_simple_mpi ("&" jns) ("&" ppis->MPI_fieldE_filter))
			   )
			 )
		       (block
			 (class-header-Particle_in_Cell_MPI ppis)
			 (define-Field3D_MPI* pMPI_fieldE ("&" MPI_fieldE))
			 (define-Field3D_MPI* pMPI_fieldB ("&" MPI_fieldB))
			 (define-Field3D_MPI* pMPI_fieldB1 ("&" MPI_fieldB1))
			 (if M_USE_MIDP_IMPLICIT 
			   (begin
			     (define-Field3D_MPI* pMPI_fieldE1 ("&" MPI_fieldE1))
			     ;(MPI_Yee_FDTD_Curl_E pMPI_fieldB1 pMPI_fieldE (* DELTAT 0.5))
			     ;(MPI_Yee_FDTD_Curl_E pMPI_fieldB1 pMPI_fieldE1 (* DELTAT 0.5))
			     ;(blas_yiszero_synced_Field3D_MPI pMPI_FoutJ pMPI_FoutJ)
			     (sync_ovlp_mpi_field pMPI_fieldE1)
			     (sync_ovlp_mpi_field pMPI_fieldB1)
			     (if M_USE_VLO
			       (MPI_slab_midp_nr_Bfield_pushJ_vlo ppis pMPI_fieldE pMPI_fieldE1 pMPI_fieldB pMPI_fieldB1 pMPI_FoutJ pmass pchg DELTAT 0 0 0 0 0 0)
			       (MPI_slab_midp_nr_Bfield_pushJ ppis pMPI_fieldE pMPI_fieldE1 pMPI_fieldB pMPI_fieldB1 pMPI_FoutJ pmass pchg DELTAT 0 0 0 0 0 0)
			       )
			     ;(merge_ovlp_mpi_field pMPI_FoutJ)
			     ;(sync_ovlp_mpi_field pMPI_FoutJ)
			     (blas_yisax_Field3D_MPI pMPI_fieldE pMPI_fieldE 1 pMPI_fieldE1)
			     )
			   (begin
			     (blas_yisax_Field3D_MPI ("&" MPI_fieldB1) ("&" MPI_fieldB1) 1 pMPI_fieldB)
			     (MPI_Yee_FDTD_Curl_E ("&" MPI_fieldB1) ("&" ppis->MPI_fieldE_filter) DELTAT)
			     (MPI_slab_nr_Bfield_pushJ ppis ("&" ppis->MPI_fieldE_filter) pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) 0 0 0 0 0 0)
			     )
			   )
			 (blas_yisax_Field3D_MPI pMPI_fieldB pMPI_fieldB 1 pMPI_fieldB1) 
			 (call_particle_sort_mpi ppis 0 0)
			 (call_particle_sort_mpi ppis 1 0)
			 (call_particle_sort_mpi ppis 2 0)
			 )
		       ;(MPI_geo_nr_Bfield_pushJ_vlo ppis pMPI_fieldE pMPI_fieldB ("&" MPI_fieldB1) pMPI_FoutJ pmass pchg (* 0.5 DELTAT) T_TORI_X0 0 0 0 0 0)
		       ;(MPI_geo_nr_Bfield_pushJ_)
		       )
		     (M_USE_MULTI_STEP
		       ;(LOG_RECORD_INFO "multi step\n")
		       ;(split_multi_step_passes_test ppis DELTAT 1 3 pmass pchg)
		       (split_2nd_all_passes_multi_step ppis DELTAT G_USE_G_E 2 1 pmass pchg)
		       )
		     (M_USE_EOUT_SCHEME
		       (split_2nd_all_passes_Eout ppis DELTAT G_USE_G_E) 
		       )
		     (else
		       (split_2nd_all_passes ppis DELTAT G_USE_G_E (if M_USE_SPLIT_MULTI_SUBSTEP num_substep 0)))
		     )))))
	  ;;;; main iteration step ends
	  (if (eq? (remainder t num_multi_step_sort) 0)
	    (begin
	      (if G_KROOK_MU_FREQ
		(MPI_krook_collision_test ppis pext_krook ("&" ppis->MPI_fieldB) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT G_KROOK_MU_FREQ)
		)
	      (if G_USE_KROOK_FIELD
		(MPI_krook_collision_test ppis pext_krook ("&" ppis->MPI_fieldB) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT 0)
		)
	      (if G_KROOK_SPEED_LIMIT
		(MPI_krook_collision_remove_small_speed ppis pext_krook ("&" ppis->MPI_fieldB) ("&" ppis->MPI_FoutJ) pmass pchg DELTAT G_KROOK_SPEED_LIMIT)
		)
	      )
	    )
	  (gen_ips_info)
	  ;(split_1st_all_passes ppis 0.5)
	  ;(MPI_Yee_FDTD_Curl_E ptestfield ptestfieldB 0.5)
	  ;(MPI_Yee_FDTD_Curl_B ptestfieldB ptestfield 0.5)
	  )
	(GAPS_IO_DeleteDataInfo pgid)
	(if (eq? rank 0)
	  (begin
	    (LOG_RECORD_INFO "task complete: ")
	    (for-from-to i 0 argc 
	      (LOG_RECORD_INFO "%s " (vector-ref argv i))
	      )
	    (LOG_RECORD_INFO "\n")
	    )
	  ;(LOG_RECORD_INFO "task complete: %s %s\n" (vector-ref argv 0))
	  )

	;(if 0 (block (dec-array char fname 100) (sprintf fname "dbgF%03d" t) (sync_main_data_d2h ptestfield) (mpi_field_2_outfile ptestfield fname)))
	)))
  ;(fffor_xyz_from_zero_to 6 6 2 (LOG_RECORD_INFO "%e \n" (GET_FIELD3D_SEQ (structp-ref ptestfield data) 2 (- xyzx 2) (- xyzy 2) xyzz 0)))
  ;(LOG_RECORD_INFO "sync...\n")
  ;(sync_ovlp_mpi_field ptestfield)
  ;(mpi_field_2_outfile ptestfield "dbgfile")
  ;(eval-scmc-global `(decl-var-and-pvar ,(add_the_prefix 'pscmc_env) pe))
  ;((macro-add-the-prefix pscmc_env_init) ppe 0 NULL)
  (define i 0)
  ;(sync_main_data_d2h ptestfield)
  ;(choose_the_runtime (pscmc_mem_sync_h2d (vector-ref (+ (structp-ref ptestfield data) i))))
  ;(for i=0 i<n i++ (LOG_RECORD_OUT "rank %d, i=%d tid[i]=%d\n" rank i (vector-ref tids i)))
  ;(LOG_RECORD_OUT "rank %d, cur_tid=%d\n" rank cur_tid)
  ;(init_adjoint_relations pfstest n_hilbert 0 n tids)
  
  ;(LOG_RECORD_OUT "rank %d, adjid=%d\n" rank (vector-ref (structp-ref pfstest adj_processes) 13))
  ;(LOG_RECORD_OUT "%d\n" (vector-ref tids 4))
  ;(for i=0 (< i NUM_SYNC_LAYER) i++ 0 (LOG_RECORD_OUT "rank %d: i=%ld pid=%ld tid=%ld\n" rank i (vector-ref (structp-ref pfstest adj_processes) i) (vector-ref (structp-ref pfstest adj_ids) i)))
  ;(fffor_xyz_from_zero_to 1 6 6 (LOG_RECORD_INFO "0x%x 0x%x %e xyz=[%d %d %d]\n" ("&" (GET_FIELD3D_SEQ (structp-ref ptestfield data) 0 -2 -2 -2 0)) ("&" (GET_FIELD3D_SEQ (structp-ref ptestfield data) 3 0 (- xyzy 2) (- xyzz 2) 0)) (GET_FIELD3D_SEQ (structp-ref ptestfield data) 3 0 (- xyzy 2) (- xyzz 2) 0) xyzx xyzy xyzz))
  ;(LOG_RECORD_INFO "0x%x 0x%x\n" ("&" (GET_FIELD3D_SEQ (structp-ref ptestfield data) 0 -2 -2 -2 0)) ("&" (GET_FIELD3D_SEQ (structp-ref ptestfield data) 3 0 0 0 0)))
  ;(delete_Field3D_MPI ptestfield)
  (if (eq? rank 0)
    (LOG_RECORD_INFO "Calling Finalize\n"))
  (PS_MPI_Finalize)
  (return 0)
  )