Time-step dependent force interpolation scheme

epoch1d/Makefile

@@ -206,6 +206,9 @@ DEFINES := $(DEFINE)
 # Use fifth order particle weighting (default is third order).
 #DEFINES += $(D)PARTICLE_SHAPE_BSPLINE3
 
+# Use WT scheme for field interpolation
+#DEFINES += $(D)WT_INTERPOLATION
+
 # Include a unique global particle ID. The first flag defines the ID using
 # an 8-byte integer, the second uses 4-bytes.
 #DEFINES += $(D)PARTICLE_ID

epoch1d/src/constants.F90

@@ -151,6 +151,7 @@ MODULE constants
   INTEGER, PARAMETER :: c_maxwell_solver_lehe_z = 4
   INTEGER, PARAMETER :: c_maxwell_solver_cowan = 5
   INTEGER, PARAMETER :: c_maxwell_solver_pukhov = 6
+  INTEGER, PARAMETER :: c_maxwell_solver_m4 = 7
 
   ! Particle distribution type codes
   INTEGER, PARAMETER :: c_ic_df_thermal = 1
@@ -266,6 +267,7 @@ MODULE constants
   INTEGER(i8), PARAMETER :: c_def_use_isatty = 2**24
   INTEGER(i8), PARAMETER :: c_def_use_mpi3 = 2**25
   INTEGER(i8), PARAMETER :: c_def_bremsstrahlung = 2**26
+  INTEGER(i8), PARAMETER :: c_def_wt_interpolation = 2**27
 
   ! Stagger types
   INTEGER, PARAMETER :: c_stagger_ex = c_stagger_face_x
@@ -448,11 +450,12 @@ MODULE constants
   INTEGER, PARAMETER :: c_const_maxwell_solver_lehe_z = 104
   INTEGER, PARAMETER :: c_const_maxwell_solver_cowan = 105
   INTEGER, PARAMETER :: c_const_maxwell_solver_pukhov = 106
-  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_m4 = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 108
 
-  INTEGER, PARAMETER :: c_const_px = 108
-  INTEGER, PARAMETER :: c_const_py = 109
-  INTEGER, PARAMETER :: c_const_pz = 110
+  INTEGER, PARAMETER :: c_const_px = 109
+  INTEGER, PARAMETER :: c_const_py = 110
+  INTEGER, PARAMETER :: c_const_pz = 111
 
   ! Custom constants
   INTEGER, PARAMETER :: c_const_deck_lowbound = 4096

epoch1d/src/deck/deck_control_block.F90

@@ -367,6 +367,7 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
           .AND. maxwell_solver /= c_maxwell_solver_lehe_x &
           .AND. maxwell_solver /= c_maxwell_solver_cowan &
           .AND. maxwell_solver /= c_maxwell_solver_pukhov &
+          .AND. maxwell_solver /= c_maxwell_solver_m4 &
           .AND. maxwell_solver /= c_maxwell_solver_custom) THEN
         errcode = c_err_bad_value
       END IF

epoch1d/src/fields.f90

@@ -59,6 +59,11 @@ SUBROUTINE set_maxwell_solver
       dx_cdt = dx / (c * dt)
       deltax = 0.25_num * (1.0_num - dx_cdt**2 * SIN(0.5_num * pi / dx_cdt)**2)
       alphax = 1.0_num - 3.0_num * deltax
+
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+      deltax = ((c * dt / dx)**2 - 1.0) / 12.0_num
+      alphax = 1.0_num - 3.0_num * deltax
     END IF
 
     IF (rank == 0 .AND. maxwell_solver /= c_maxwell_solver_yee) THEN

epoch1d/src/housekeeping/setup.F90

@@ -582,6 +582,8 @@ SUBROUTINE set_dt        ! sets CFL limited step
       ! R. Lehe, PhD Thesis (2014)
       dt = dx / c
 
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      dt = dx / c
     END IF
 
     IF (any_open) THEN
@@ -621,6 +623,16 @@ SUBROUTINE set_dt        ! sets CFL limited step
 
     dt = dt_multiplier * dt
 
+#ifdef WT_INTERPOLATION
+    IF (c * dt / dx > 0.5_num) THEN
+      IF (rank == 0) THEN
+        PRINT*, '*** ERROR ***'
+        PRINT*, 'Cannot use WT beacause c*dt/dx>0.5'
+      END IF
+      CALL abort_code(c_err_bad_setup)
+    END IF
+#endif
+
     IF (.NOT. any_average) RETURN
 
     DO io = 1, n_io_blocks

epoch1d/src/housekeeping/welcome.F90

@@ -99,6 +99,9 @@ SUBROUTINE compiler_directives
 #ifdef PARTICLE_SHAPE_TOPHAT
     found = .TRUE.
 #endif
+#ifdef WT_INTERPOLATION
+    found = .TRUE.
+#endif
 #ifdef PER_SPECIES_WEIGHT
     found = .TRUE.
 #endif
@@ -186,6 +189,10 @@ SUBROUTINE compiler_directives
     defines = IOR(defines, c_def_particle_shape_tophat)
     WRITE(*,*) 'Top-hat particle shape -DPARTICLE_SHAPE_TOPHAT'
 #endif
+#ifdef WT_INTERPOLATION
+    defines = IOR(defines, c_def_wt_interpolation)
+    WRITE(*,*) 'WT interpolation scheme -DWT_INTERPOLATION'
+#endif
 #ifdef PER_SPECIES_WEIGHT
     defines = IOR(defines, c_def_per_particle_weight)
     WRITE(*,*) 'Per species weighting -DPER_SPECIES_WEIGHT'

epoch1d/src/include/bspline3/gx_wt.inc

@@ -0,0 +1,15 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/24
+        wt_var1 = wt_facx * (wt_dtx + cell_frac_x)**4
+        wt_var2 = wt_facx * (wt_dtx - cell_frac_x)**4
+        wt_var3 = SIGN(wt_var1, wt_dtx + cell_frac_x) &
+            + SIGN(wt_var2, wt_dtx - cell_frac_x)
+        wt_var4 = ((3.0_num - cell_frac_x) * cell_frac_x + 3.0_num - wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        wt_var5 = ((3.0_num + cell_frac_x) * cell_frac_x - 3.0_num + wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        gx(-2) = wt_var3 + wt_var1 - wt_var2
+        gx(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gx( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gx( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gx( 2) = wt_var3 + wt_var2 - wt_var1

epoch1d/src/include/tophat/gx_wt.inc

@@ -0,0 +1,4 @@
+        wt_var1 = wt_facx * (ABS(wt_dtx + cell_frac_x) &
+            - ABS(wt_dtx - cell_frac_x))
+        gx( 0) = 0.5_num + wt_var1
+        gx( 1) = 0.5_num - wt_var1

epoch1d/src/include/triangle/gx_wt.inc

@@ -0,0 +1,8 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/2
+        wt_var1 = wt_facx * ( &
+            ABS(wt_dtx + cell_frac_x) * (wt_dtx + cell_frac_x) &
+            + ABS(wt_dtx - cell_frac_x) * (wt_dtx - cell_frac_x))
+        gx(-1) = wt_var1 + cell_frac_x
+        gx( 0) = 2.0_num - 2.0_num * wt_var1
+        gx( 1) = wt_var1 - cell_frac_x

epoch1d/src/parser/evaluator_blocks.F90

@@ -634,6 +634,11 @@ SUBROUTINE do_constant(opcode, simplify, parameters, err)
       RETURN
     END IF
 
+    IF (opcode == c_const_maxwell_solver_m4) THEN
+      CALL push_on_eval(REAL(c_maxwell_solver_m4, num))
+      RETURN
+    END IF
+
     IF (opcode == c_const_maxwell_solver_custom) THEN
       CALL push_on_eval(REAL(c_maxwell_solver_custom, num))
       RETURN

epoch1d/src/parser/tokenizer_blocks.f90

@@ -302,6 +302,9 @@ FUNCTION as_constant(name)
     ELSE IF (str_cmp(name, 'pukhov')) THEN
       as_constant = c_const_maxwell_solver_pukhov
 
+    ELSE IF (str_cmp(name, 'm4')) THEN
+      as_constant = c_const_maxwell_solver_m4
+
     ELSE IF (str_cmp(name, 'lehe_x')) THEN
       as_constant = c_const_maxwell_solver_lehe_x
 

epoch1d/src/particles.F90

@@ -82,6 +82,14 @@ SUBROUTINE push_particles
     ! This is to deal with the grid stagger
     REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx
 
+    ! For WT scheme, the weight factors used to weight particle properties onto
+    ! grid to calculate J is same as bspline3, triangle or tophat.  But the
+    ! particle weight factors to weight fields onto particles are not bspline3,
+    ! triangle or tophat. We use hx0 to store the former weight factors at Xi0*.
+#ifdef WT_INTERPOLATION
+    REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx0
+#endif
+
     ! Fields at particle location
     REAL(num) :: ex_part, ey_part, ez_part, bx_part, by_part, bz_part
 
@@ -125,6 +133,15 @@ SUBROUTINE push_particles
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 #ifdef DELTAF_METHOD
     REAL(num) :: weight_back
@@ -141,6 +158,9 @@ SUBROUTINE push_particles
     jz = 0.0_num
 
     gx = 0.0_num
+#ifdef WT_INTERPOLATION
+    hx0 = 0.0_num
+#endif
 
     ! Unvarying multiplication factors
 
@@ -153,6 +173,14 @@ SUBROUTINE push_particles
     idtf = idt * fac
     idxf = idx * fac
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     DO ispecies = 1, n_species
       current => species_list(ispecies)%attached_list%head
       IF (species_list(ispecies)%immobile) CYCLE
@@ -267,6 +295,21 @@ SUBROUTINE push_particles
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+
+        ! For WT scheme, use hx0 to store the weight factors used to weigh
+        ! particle properties onto grid at Xi0*. And calculate gx for WT.
+        ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+        ! NOTE: These weights require an additional multiplication factor!
+#ifdef WT_INTERPOLATION
+        hx0 = gx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -431,6 +474,9 @@ SUBROUTINE push_particles
 
           ! Now change Xi1* to be Xi1*-Xi0*. This makes the representation of
           ! the current update much simpler
+#ifdef WT_INTERPOLATION
+          gx = hx0
+#endif
           hx = hx - gx
 
           ! Remember that due to CFL condition particle can never cross more