Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_time_propagator_pic_vm_1d2v_helper.F90
Go to the documentation of this file.
1 module sll_m_time_propagator_pic_vm_1d2v_helper
2  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 #include "sll_assert.h"
4 #include "sll_errors.h"
5 #include "sll_memory.h"
6 #include "sll_working_precision.h"
7 
8  use sll_m_collective, only: &
9  sll_f_get_collective_rank, &
10  sll_o_collective_allreduce, &
11  sll_v_world_collective
12 
13  use sll_m_control_variate, only: &
14  sll_t_control_variates
15 
16  use sll_m_filter_base_1d, only: &
17  sll_c_filter_base_1d
18 
19  use sll_m_linear_operator_particle_mass_1d, only : &
20  sll_t_linear_operator_particle_mass_1d
21 
22  use sll_m_linear_operator_particle_mass_smooth_1d, only : &
23  sll_t_linear_operator_particle_mass_smooth_1d
24 
25  use sll_m_linear_operator_particle_mass_cl_1d, only : &
26  sll_t_linear_operator_particle_mass_cl_1d
27 
28  use sll_m_linear_operator_schur_ev_1d, only : &
29  sll_t_linear_operator_schur_ev_1d
30 
31  use sll_m_linear_operator_schur_phiv_1d, only : &
32  sll_t_linear_operator_schur_phiv_1d
33 
34  use sll_m_linear_solver_cg, only : &
35  sll_t_linear_solver_cg
36 
37  use sll_m_maxwell_1d_base, only: &
38  sll_c_maxwell_1d_base
39 
40  use mpi, only: &
41  mpi_sum, &
42  mpi_max
43 
44  use sll_m_particle_mass_1d_base, only: &
45  sll_c_particle_mass_1d_base
46 
47  use sll_m_particle_group_base, only: &
48  sll_t_particle_array
49 
50  use sll_m_particle_mesh_coupling_base_1d, only: &
51  sll_c_particle_mesh_coupling_1d
52 
53  use sll_m_particle_mesh_coupling_spline_1d, only: &
54  sll_t_particle_mesh_coupling_spline_1d
55 
56  use sll_m_particle_mesh_coupling_spline_smooth_1d, only: &
57  sll_t_particle_mesh_coupling_spline_smooth_1d
58 
59  use sll_m_particle_mesh_coupling_spline_strong_1d, only: &
60  sll_t_particle_mesh_coupling_spline_strong_1d
61 
62  use sll_m_time_propagator_base, only: &
63  sll_c_time_propagator_base
64 
65  use sll_m_time_propagator_pic_vm_3d3v_cl_helper, only: &
66  sll_p_boundary_particles_periodic, &
67  sll_p_boundary_particles_singular, &
68  sll_p_boundary_particles_reflection, &
69  sll_p_boundary_particles_absorption
70 
71 
72  implicit none
73 
74  public :: &
75  sll_t_time_propagator_pic_vm_1d2v_helper
76 
77  private
78  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79 
81  type :: sll_t_time_propagator_pic_vm_1d2v_helper
82  class(sll_c_maxwell_1d_base), pointer :: maxwell_solver
83  class(sll_c_particle_mesh_coupling_1d), pointer :: kernel_smoother_0
84  class(sll_c_particle_mesh_coupling_1d), pointer :: kernel_smoother_1
85  class(sll_t_particle_array), pointer :: particle_group
86  class( sll_c_particle_mass_1d_base ), allocatable :: particle_mass_1
87  class( sll_c_particle_mass_1d_base ), allocatable :: particle_mass_0
88 
89  type( sll_t_linear_operator_schur_ev_1d ) :: linear_operator_1
90  type(sll_t_linear_solver_cg) :: linear_solver_1
91  type( sll_t_linear_operator_schur_ev_1d ) :: linear_operator_0
92  type(sll_t_linear_solver_cg) :: linear_solver_0
93 
94  type( sll_t_linear_operator_schur_phiv_1d ) :: linear_op_schur_phiv
95  type( sll_t_linear_solver_cg ) :: linear_solver_schur_phiv
96 
97  sll_int32 :: spline_degree
98  sll_real64 :: lx
99  sll_real64 :: x_min
100  sll_real64 :: x_max
101  sll_real64 :: delta_x
102  sll_int32 :: n_cells
103  sll_int32 :: n_dofs0
104  sll_int32 :: n_dofs1
105 
106  logical :: smooth = .false.
107  sll_int32 :: size_particle_mass_0
108  sll_int32 :: size_particle_mass_1
109  logical :: build_particle_mass_loc = .true.
110 
111  sll_int32 :: boundary_particles = 100
112  sll_int32 :: counter_left = 0
113  sll_int32 :: counter_right = 0
114  logical :: boundary = .false.
115 
116  sll_real64, pointer :: phi_dofs(:)
117  sll_real64, pointer :: efield_dofs(:,:)
118  sll_real64, pointer :: bfield_dofs(:)
119  sll_real64, allocatable :: j_dofs(:,:)
120  sll_real64, allocatable :: j_dofs_local(:,:)
121 
122  sll_real64, allocatable :: particle_mass_0_local(:,:)
123  sll_real64, allocatable :: particle_mass_1_local(:,:)
124  sll_real64, allocatable :: residual(:)
125 
126  sll_real64, allocatable :: xnew(:,:)
127  sll_real64, allocatable :: vnew(:,:,:)
128  sll_real64, allocatable :: efield_dofs_new(:,:)
129  sll_real64, allocatable :: efield_dofs_work(:,:)
130 
131 
132  sll_int32 :: niter(100000)
133  sll_int32 :: nevaliter(100000)
134  sll_real64, allocatable :: nsubiter(:,:)
135  sll_int32 :: iter_counter = 0
136  sll_int32, allocatable :: sub_iter_counter(:)
137  sll_int32 :: eval_counter = 0
138  sll_int32 :: max_iter = 10
139  sll_int32 :: max_iter_sub = 100
140  sll_real64 :: iter_tolerance = 1d-12
141  sll_real64 :: iter_tolerance_sub = 1d-10
142  sll_int32 :: n_failed = 0
143 
144  sll_real64 :: betar(2) = 1.0_f64
145  sll_real64 :: force_sign = 1._f64
146  logical :: jmean = .false.
147  logical :: adiabatic_electrons = .false.
148  sll_real64 :: solver_tolerance = 1d-12
149 
150  sll_int32 :: n_sub_iter(2) = [8,2]
151 
152  character(len=256) :: output_prefix = "disgrad"
153 
154  sll_int32 :: n_particles_max ! Helper variable for size of temporary particle arrays
155 
156 
157  sll_real64, allocatable :: efield_filter(:,:)
158  sll_real64, allocatable :: bfield_filter(:)
159 
160  sll_real64, allocatable :: efield_to_val(:,:)
161  sll_real64, allocatable :: bfield_to_val(:)
162 
163  class(sll_c_filter_base_1d), pointer :: filter
164 
165  ! For version with control variate
166  class(sll_t_control_variates), pointer :: control_variate => null()
167  sll_int32 :: i_weight
168 
169  contains
170  procedure :: advect_x => advect_x_pic_vm_1d2v_helper
171  procedure :: advect_vb => advect_vb_pic_vm_1d2v_helper
172  procedure :: advect_eb => advect_eb_pic_vm_1d2v_helper
173  procedure :: advect_e => advect_e_pic_vm_1d2v_helper
174  procedure :: advect_ex => advect_ex_pic_vm_1d2v_helper
175  procedure :: advect_e_sub => advect_e_sub_pic_vm_1d2v_helper
176 
177  procedure :: init => initialize_pic_vm_1d2v_helper
178  procedure :: init_from_file => initialize_file_pic_vm_1d2v_helper
179  procedure :: free => delete_pic_vm_1d2v_helper
180 
181  procedure :: reinit_fields
182 
183  end type sll_t_time_propagator_pic_vm_1d2v_helper
184 
185 contains
186 
187  !---------------------------------------------------------------------------!
189  subroutine advect_x_pic_vm_1d2v_helper ( self, dt )
190  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
191  sll_real64, intent(in) :: dt
192  !local variables
193  sll_int32 :: i_part, i_sp
194  sll_real64 :: xi(3), xnew(3), vi(3), wp(3)
195 
196  do i_sp = 1, self%particle_group%n_species
197  do i_part = 1, self%particle_group%group(i_sp)%n_particles
198  ! Read out particle position and velocity
199  xi = self%particle_group%group(i_sp)%get_x(i_part)
200  vi = self%particle_group%group(i_sp)%get_v(i_part)
201 
202  xnew(1) = xi(1) + dt * vi(1)
203 
204  call compute_particle_boundary( self, xi(1), xnew(1), vi(1) )
205  call self%particle_group%group(i_sp)%set_x ( i_part, xnew )
206  call self%particle_group%group(i_sp)%set_v ( i_part, vi )
207  if (self%particle_group%group(i_sp)%n_weights == 3) then
208  ! Update weights if control variate
209  wp = self%particle_group%group(i_sp)%get_weights(i_part)
210  wp(3) = self%control_variate%cv(i_sp)%update_df_weight( xnew(1:1), vi(1:2), 0.0_f64, wp(1), wp(2))
211  call self%particle_group%group(i_sp)%set_weights(i_part, wp)
212  end if
213  end do
214  end do
215 
216  end subroutine advect_x_pic_vm_1d2v_helper
217 
218 
220  subroutine compute_particle_boundary( self, xold, xnew, vi )
221  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent( inout ) :: self
222  sll_real64, intent( inout ) :: xold
223  sll_real64, intent( inout ) :: xnew
224  sll_real64, intent( inout ) :: vi
225  !local variables
226  sll_real64 :: xmid, xbar, dx
227 
228  if(xnew < -self%x_max .or. xnew > 2._f64*self%x_max ) then
229  print*, xnew
230  sll_error("particle boundary", "particle out of bound")
231  else if(xnew < self%x_min .or. xnew > self%x_max )then
232  if(xnew < self%x_min )then
233  xbar = self%x_min
234  self%counter_left = self%counter_left+1
235  else if(xnew > self%x_max)then
236  xbar = self%x_max
237  self%counter_right = self%counter_right+1
238  end if
239  dx = (xbar- xold)/(xnew-xold)
240  xmid = xold + dx * (xnew-xold)
241  xmid = xbar
242 
243  select case(self%boundary_particles)
244  case(sll_p_boundary_particles_reflection)
245  vi = -vi
246  xnew = 2._f64*xbar-xnew
247  case(sll_p_boundary_particles_absorption)
248  case( sll_p_boundary_particles_periodic)
249  xnew = self%x_min + modulo(xnew-self%x_min, self%Lx)
250  case default
251  xnew = self%x_min + modulo(xnew-self%x_min, self%Lx)
252  end select
253  end if
254 
255  end subroutine compute_particle_boundary
256 
257 
258  !---------------------------------------------------------------------------!
261  subroutine advect_vb_pic_vm_1d2v_helper ( self, dt )
262  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
263  sll_real64, intent(in) :: dt
264  !local variables
265  sll_int32 :: i_part, i_sp
266  sll_real64 :: qmdt
267  sll_real64 :: vi(3), v_new(3), xi(3), wp(3)
268  sll_real64 :: bfield, cs, sn
269 
270 
271  call self%maxwell_solver%transform_dofs( self%bfield_filter, self%bfield_to_val, 0 )
272 
273  do i_sp = 1, self%particle_group%n_species
274  qmdt = self%particle_group%group(i_sp)%species%q_over_m()*dt;
275  do i_part = 1, self%particle_group%group(i_sp)%n_particles
276  vi = self%particle_group%group(i_sp)%get_v(i_part)
277  xi = self%particle_group%group(i_sp)%get_x(i_part)
278  call self%kernel_smoother_1%evaluate &
279  (xi(1), self%bfield_to_val, bfield)
280 
281  bfield = qmdt*bfield
282 
283  cs = cos(bfield)
284  sn = sin(bfield)
285 
286  v_new(1) = cs * vi(1) + sn * vi(2)
287  v_new(2) = -sn* vi(1) + cs * vi(2)
288 
289  call self%particle_group%group(i_sp)%set_v( i_part, v_new )
290  if (self%particle_group%group(i_sp)%n_weights == 3) then
291  ! Update weights if control variate
292  wp = self%particle_group%group(i_sp)%get_weights(i_part)
293  wp(3) = self%control_variate%cv(i_sp)%update_df_weight( xi(1:1), vi(1:2), 0.0_f64, wp(1), wp(2))
294  call self%particle_group%group(i_sp)%set_weights(i_part, wp)
295  end if
296  end do
297  end do
298 
299 
300  end subroutine advect_vb_pic_vm_1d2v_helper
301 
302  !---------------------------------------------------------------------------!
307  subroutine advect_eb_pic_vm_1d2v_helper ( self, dt )
308  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
309  sll_real64, intent(in) :: dt
310 
311 
312  call self%maxwell_solver%compute_curl_part( dt, self%efield_dofs(:,2), self%bfield_dofs, self%betar(1) )
313 
314 
315  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
316  call self%filter%apply( self%bfield_dofs, self%bfield_filter )
317 
318  end subroutine advect_eb_pic_vm_1d2v_helper
319 
320 
321  !---------------------------------------------------------------------------!
325  subroutine advect_e_pic_vm_1d2v_helper ( self, dt )
326  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
327  sll_real64, intent(in) :: dt
328  ! local variables
329  sll_int32 :: i_part, i_sp
330  sll_real64 :: vi(3), xi(3), wi(1), wp(3)
331  sll_real64 :: qoverm, factor
332  sll_real64 :: efield(2)
333  sll_real64 :: rhs0(self%n_dofs0)
334  sll_real64 :: rhs1(self%n_dofs1)
335 
336 
337  ! Set to zero
338  self%j_dofs_local = 0.0_f64
339  self%particle_mass_1_local = 0.0_f64
340  self%particle_mass_0_local = 0.0_f64
341 
342  ! First particle loop
343  do i_sp=1,self%particle_group%n_species
344  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
345  factor = dt**2*0.25_f64*self%betar(2)*qoverm
346  do i_part = 1,self%particle_group%group(i_sp)%n_particles
347  vi = self%particle_group%group(i_sp)%get_v(i_part)
348  xi = self%particle_group%group(i_sp)%get_x(i_part)
349 
350  ! Get charge for accumulation of j
351  wi = self%particle_group%group(i_sp)%get_charge(i_part)
352 
353  ! Accumulate the particle mass matrix diagonals
354  call self%kernel_smoother_1%add_particle_mass_full( xi(1), &
355  wi(1) * factor, &
356  self%particle_mass_1_local )
357  call self%kernel_smoother_0%add_particle_mass_full( xi(1), &
358  wi(1) * factor, &
359  self%particle_mass_0_local )
360 
361  ! Accumulate jx
362  call self%kernel_smoother_1%add_charge( xi(1), wi(1)*vi(1), &
363  self%j_dofs_local(1:self%n_dofs1,1) )
364  ! Accumulate jy
365  call self%kernel_smoother_0%add_charge( xi(1), wi(1)*vi(2), &
366  self%j_dofs_local(:,2) )
367 
368  ! Evaulate E_x at particle position and propagate v a half step
369  call self%kernel_smoother_1%evaluate &
370  (xi(1), self%efield_filter(1:self%n_dofs1,1), efield(1) )
371  ! Evaulate E_y at particle position and propagate v a half step
372  call self%kernel_smoother_0%evaluate &
373  (xi(1), self%efield_filter(:,2), efield(2))
374 
375  vi(1:2) = vi(1:2) + dt* 0.5_f64* qoverm * efield
376  call self%particle_group%group(i_sp)%set_v( i_part, vi )
377 
378  if (self%particle_group%group(i_sp)%n_weights == 3) then
379  ! Update weights if control variate
380  wp = self%particle_group%group(i_sp)%get_weights(i_part)
381  wp(3) = self%control_variate%cv(i_sp)%update_df_weight( xi(1:1), vi(1:2), 0.0_f64, wp(1), wp(2))
382  call self%particle_group%group(i_sp)%set_weights(i_part, wp)
383  end if
384  end do
385  end do
386 
387  self%j_dofs = 0.0_f64
388  self%particle_mass_0%particle_mass = 0.0_f64
389  self%particle_mass_1%particle_mass = 0.0_f64
390  ! MPI to sum up contributions from each processor
391  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(1:self%n_dofs1,1), &
392  self%n_dofs1, mpi_sum, self%j_dofs(1:self%n_dofs1,1) )
393  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
394  self%n_dofs0, mpi_sum, self%j_dofs(:,2) )
395 
396  self%particle_mass_1%particle_mass = 0.0_f64
397  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_1_local, &
398  self%size_particle_mass_1, mpi_sum, self%particle_mass_1%particle_mass)
399  self%particle_mass_0%particle_mass = 0.0_f64
400  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_0_local, &
401  self%size_particle_mass_0, mpi_sum, self%particle_mass_0%particle_mass)
402 
403  call self%filter%apply_inplace( self%j_dofs(:,1) )
404  call self%filter%apply_inplace( self%j_dofs(:,2) )
405 
406  if ( self%jmean ) then
407  self%j_dofs(:,1) = self%j_dofs(:,1) - sum(self%j_dofs(1:self%n_dofs1,1))/real(self%n_dofs1, f64)
408  end if
409 
410  if( self%adiabatic_electrons ) then
411  ! Compute rhs
412  self%linear_op_schur_phiv%sign = -1._f64
413  call self%linear_op_schur_phiv%dot(self%phi_dofs, rhs0)
414  call self%maxwell_solver%multiply_gt( self%j_dofs(1:self%n_dofs1,1), self%j_dofs(:,2) )
415  rhs0 = rhs0 + dt * self%j_dofs(:,2)
416 
417  ! Solve the Schur complement
418  self%linear_op_schur_phiv%sign = 1._f64
419  call self%linear_solver_schur_phiv%set_guess(self%phi_dofs)
420  call self%linear_solver_schur_phiv%solve( rhs0, self%phi_dofs )
421  call self%maxwell_solver%multiply_g( self%phi_dofs, self%efield_dofs(1:self%n_dofs1,1) )
422  self%efield_dofs(:,1) = -self%efield_dofs(:,1)
423  else
424 
425  ! Compute rhs
426  self%linear_operator_1%sign = -self%force_sign
427  call self%linear_operator_1%dot(self%efield_dofs(1:self%n_dofs1,1) , rhs1 )
428  rhs1 = rhs1 - dt * self%force_sign * self%j_dofs(1:self%n_dofs1,1)
429  ! Solve the Schur complement
430  self%linear_operator_1%sign = self%force_sign
431  call self%linear_solver_1%set_guess(self%efield_dofs(1:self%n_dofs1,1))
432  call self%linear_solver_1%solve( rhs1 , self%efield_dofs(1:self%n_dofs1,1) )
433 
434  ! Compute rhs
435  self%linear_operator_0%sign = -self%force_sign
436  call self%linear_operator_0%dot(self%efield_dofs(:,2) , rhs0 )
437  rhs0 = rhs0 - self%force_sign * dt * self%j_dofs(:,2)
438 
439  ! Solve the Schur complement
440  self%linear_operator_0%sign = self%force_sign
441  call self%linear_solver_0%set_guess(self%efield_dofs(:,2))
442  call self%linear_solver_0%solve( rhs0 , self%efield_dofs(:,2) )
443  end if
444 
445  if( self%boundary ) then !perfect conductor boundary
446  self%efield_dofs(1,2) = 0._f64
447  self%efield_dofs(self%n_dofs0,2) = 0._f64
448  end if
449 
450  call self%filter%apply( self%efield_dofs(:,1), self%efield_filter(:,1) )
451  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
452  ! Second particle loop (second half step of particle propagation)
453  do i_sp=1,self%particle_group%n_species
454  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
455  do i_part = 1, self%particle_group%group(i_sp)%n_particles
456 
457  vi = self%particle_group%group(i_sp)%get_v(i_part)
458  xi = self%particle_group%group(i_sp)%get_x(i_part)
459  ! Evaulate E_x at particle position and propagate v a half step
460  call self%kernel_smoother_1%evaluate &
461  (xi(1), self%efield_filter(1:self%n_dofs1,1), efield(1))
462  ! Evaulate E_y at particle position and propagate v a half step
463  call self%kernel_smoother_0%evaluate &
464  (xi(1), self%efield_filter(:,2), efield(2))
465 
466  vi(1:2) = vi(1:2) + dt* 0.5_f64* qoverm * efield
467  call self%particle_group%group(i_sp)%set_v( i_part, vi )
468 
469  if (self%particle_group%group(i_sp)%n_weights == 3) then
470  ! Update weights if control variate
471  wp = self%particle_group%group(i_sp)%get_weights(i_part)
472  wp(3) = self%control_variate%cv(i_sp)%update_df_weight( xi(1:1), vi(1:2), 0.0_f64, wp(1), wp(2))
473  call self%particle_group%group(i_sp)%set_weights(i_part, wp)
474  end if
475  end do
476  end do
477 
478  end subroutine advect_e_pic_vm_1d2v_helper
479 
480 
481  !---------------------------------------------------------------------------!
483  subroutine advect_ex_pic_vm_1d2v_helper ( self, dt )
484  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
485  sll_real64, intent(in) :: dt
486  ! local variables
487  sll_int32 :: i_part, i_sp
488  sll_int32 :: niter
489  sll_real64 :: vi(3), xi(3), wi(1), xnew(3), vbar(2)
490  sll_real64 :: qoverm
491  sll_real64 :: efield(2)
492  sll_real64 :: residual, residuals(4), residuals_all(4)
493  sll_real64 :: xmid(1), xbar
494 
495  residuals = 0.0_f64
496 
497  ! Instead of a first particle loop as in the standard Picard iteration, we now
498  ! perform the advect_x and advect_e part of the DISGRADE scheme
499  self%efield_dofs_new = self%efield_dofs
500  !call disgradE_for_start ( self, dt, self%efield_dofs_new )
501  do i_sp = 1, self%particle_group%n_species
502  do i_part = 1, self%particle_group%group(i_sp)%n_particles
503  vi = self%particle_group%group(i_sp)%get_v(i_part)
504  xi = self%particle_group%group(i_sp)%get_x(i_part)
505  self%xnew(i_sp,i_part) = xi(1)
506  self%vnew(i_sp,:,i_part) = vi(1:2)
507  end do
508  end do
509 
510  residual = 1.0_f64
511  niter = 0
512  !print*, 'Iteration', self%iter_counter+1
513  do while ( (residual > self%iter_tolerance) .and. niter < self%max_iter )
514  niter = niter + 1
515  residuals = 0.0_f64
516 
517  call self%filter%apply( self%efield_dofs_new(:,1), self%efield_dofs_work(:,1) )
518  call self%filter%apply( self%efield_dofs_new(:,2), self%efield_dofs_work(:,2) )
519 
520  self%efield_dofs_work = 0.5_f64*( self%efield_filter + self%efield_dofs_work )
521 
522  call self%maxwell_solver%transform_dofs( self%efield_dofs_work(:,1), self%efield_to_val(:,1), 0 )
523  call self%maxwell_solver%transform_dofs( self%efield_dofs_work(:,2), self%efield_to_val(:,2), 1 )
524 
525  ! Set to zero
526  self%j_dofs_local = 0.0_f64
527  ! First particle loop
528  do i_sp = 1, self%particle_group%n_species
529  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
530  do i_part = 1,self%particle_group%group(i_sp)%n_particles
531  vi = self%particle_group%group(i_sp)%get_v(i_part)
532  xi = self%particle_group%group(i_sp)%get_x(i_part)
533 
534  ! Get charge for accumulation of j
535  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
536 
537  vbar = 0.5_f64 * (self%vnew(i_sp,:, i_part)+vi(1:2))
538  xnew(1) = xi(1) + dt * vbar(1)
539 
540  if(xnew(1) < -self%x_max .or. xnew(1) > 2._f64*self%x_max ) then
541  print*, xnew
542  sll_error("particle boundary", "particle out of bound")
543  else if(xnew(1) < self%x_min .or. xnew(1) > self%x_max )then
544  if(xnew(1) < self%x_min )then
545  xbar = self%x_min
546  else if(xnew(1) > self%x_max)then
547  xbar = self%x_max
548  end if
549 
550  xmid = xbar
551  if ( abs(vbar(1)) > 1.0d-16 ) then
552  call self%kernel_smoother_1%add_current_evaluate &
553  ( xi(1), xmid(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
554  efield(1) )
555 
556  call self%kernel_smoother_0%add_current_evaluate &
557  ( xi(1), xmid(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
558  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
559  efield(2) )
560  else
561  call self%kernel_smoother_1%evaluate &
562  (xi(1), self%efield_to_val(1:self%n_dofs1,1), efield(1) )
563  efield(1) = efield(1)*dt
564  call self%kernel_smoother_0%add_charge( xi(1), &
565  wi(1)* vbar(2)*dt, self%j_dofs_local(:,2) )
566  call self%kernel_smoother_0%evaluate &
567  (xi(1), self%efield_to_val(:,2), efield(2) )
568  efield(2) = efield(2)*dt
569  end if
570  vi(1:2) = vi(1:2) + qoverm * efield(1:2)
571 
572  select case(self%boundary_particles)
573  case(sll_p_boundary_particles_reflection)
574  xnew(1) = 2._f64*xbar-xnew(1)
575  vi(1) = -vi(1)
576  vbar(1) = -vbar(1)
577  case(sll_p_boundary_particles_absorption)
578  case( sll_p_boundary_particles_periodic)
579  xnew(1) = self%x_min + modulo(xnew(1)-self%x_min, self%Lx)
580  xmid = self%x_max+self%x_min-xbar
581  case default
582  xnew(1) = self%x_min + modulo(xnew(1)-self%x_min, self%Lx)
583  xmid = self%x_max+self%x_min-xbar
584  end select
585 
586  if ( abs(vbar(1)) > 1.0d-16 ) then
587  call self%kernel_smoother_1%add_current_evaluate &
588  ( xmid(1), xnew(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
589  efield(1) )
590  call self%kernel_smoother_0%add_current_evaluate &
591  ( xmid(1), xnew(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
592  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
593  efield(2) )
594  vi(1:2) = vi(1:2) + qoverm * efield(1:2)
595  end if
596  if(self%boundary_particles == sll_p_boundary_particles_reflection) then
597  vi(1) = - vi(1)
598  end if
599  else
600  if ( abs(vbar(1)) > 1.0d-16 ) then
601  call self%kernel_smoother_1%add_current_evaluate &
602  ( xi(1), xnew(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
603  efield(1) )
604 
605  call self%kernel_smoother_0%add_current_evaluate &
606  ( xi(1), xnew(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
607  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
608  efield(2) )
609  else
610  call self%kernel_smoother_1%evaluate &
611  (xi(1), self%efield_to_val(1:self%n_dofs1,1), efield(1) )
612  efield(1) = efield(1)*dt
613  call self%kernel_smoother_0%add_charge( xi(1), &
614  wi(1)* vbar(2)*dt, self%j_dofs_local(:,2) )
615  call self%kernel_smoother_0%evaluate &
616  (xi(1), self%efield_to_val(:,2), efield(2) )
617  efield(2) = efield(2)*dt
618  end if
619  vi(1:2) = vi(1:2) + qoverm * efield(1:2)
620  end if
621  ! Here yo can change the residual criterion
622  residuals(1) = residuals(1) + (xnew(1)-self%xnew(i_sp,i_part))**2*abs(wi(1))!max( residuals(1), abs(xnew(1)-self%xnew(i_sp,i_part)) )
623  residuals(2) = residuals(2) + (vi(1)-self%vnew(i_sp,1,i_part))**2*abs(wi(1))!max( residuals(2), abs(self%vnew(i_sp,1,i_part)-vi(1)) )
624  residuals(3) = residuals(3) + (vi(2)-self%vnew(i_sp,2,i_part))**2*abs(wi(1))!max( residuals(3), abs(self%vnew(i_sp,2,i_part)-vi(2)) )
625  self%xnew(i_sp,i_part) = xnew(1)
626  self%vnew(i_sp,:,i_part) = vi(1:2)
627  end do
628  end do
629 
630  ! Update d_n
631  self%j_dofs = 0.0_f64
632  ! MPI to sum up contributions from each processor
633  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(1:self%n_dofs1,1), &
634  self%n_dofs1, mpi_sum, self%j_dofs(1:self%n_dofs1,1) )
635  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
636  self%n_dofs0, mpi_sum, self%j_dofs(:,2) )
637 
638  call self%filter%apply_inplace( self%j_dofs(:,1) )
639  call self%filter%apply_inplace( self%j_dofs(:,2) )
640 
641  ! Update the electric field.
642  self%efield_dofs_work = self%efield_dofs
643  call self%maxwell_solver%compute_E_from_j(self%j_dofs(1:self%n_dofs1,1), 1, self%efield_dofs_work(1:self%n_dofs1,1) )
644  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, self%efield_dofs_work(:,2) )
645 
646 !!$ if ( self%maxwell_solver%strong_ampere .eqv. .true. ) then
647 !!$ residuals(4) = self%maxwell_solver%l2norm_squared( self%efield_dofs_work(1:self%n_dofs1,1)-self%efield_dofs_new(1:self%n_dofs1,1), self%maxwell_solver%s_deg_0 ) + &
648 !!$ self%maxwell_solver%l2norm_squared( efield_dofs_old(:,2)-&
649 !!$ efield_dofs(:,2), self%maxwell_solver%s_deg_0 -1 )
650 !!$ else
651 !!$ residuals(4) = self%maxwell_solver%l2norm_squared( efield_dofs_old(:,1)-&
652 !!$ efield_dofs(:,1), self%maxwell_solver%s_deg_0 - 1 ) + &
653 !!$ self%maxwell_solver%l2norm_squared( self%efield_dofs_work(:,2)-self%efield_dofs_new(:,2), self%maxwell_solver%s_deg_0 )
654 !!$ end if
655  ! Compute residual based on e
656  ! Here you can change the residual norm
657  residuals(4) = (sum((self%efield_dofs_work-self%efield_dofs_new)**2))*self%delta_x!maxval( self%efield_dofs_work-self%efield_dofs_new )
658 
659 
660  call sll_o_collective_allreduce( sll_v_world_collective, residuals, 4, mpi_max, residuals_all )
661 
662  residuals_all = sqrt(residuals_all)
663  residual = residuals_all(4)!max(residuals_all(4)), maxval(residuals_all(1:3))*0.01_f64)!residuals_all(4)!maxval(residuals_all)!residuals_all(4)! maxval(residuals_all)
664  !call filter( efield_dofs, self%j_dofs_local, self%kernel_smoother_1%n_dofs )
665  !efield_dofs = self%j_dofs_local
666 
667  self%efield_dofs_new = self%efield_dofs_work
668  end do
669 
670  ! OLD VERSION TAKING FAILED ITERATION
671  if ( residual > self%iter_tolerance ) then
672  print*, 'Warning: Iteration no.', self%iter_counter+1 ,'did not converge.', residuals_all(4), niter
673  self%n_failed = self%n_failed+1
674  end if
675 
676  call self%filter%apply( self%efield_dofs_new(:,1), self%efield_dofs_work(:,1) )
677  call self%filter%apply( self%efield_dofs_new(:,2), self%efield_dofs_work(:,2) )
678 
679  self%efield_dofs_work = 0.5_f64*( self%efield_filter + self%efield_dofs_work )
680 
681  call self%maxwell_solver%transform_dofs( self%efield_dofs_work(:,1), self%efield_to_val(:,1), 0 )
682  call self%maxwell_solver%transform_dofs( self%efield_dofs_work(:,2), self%efield_to_val(:,2), 1 )
683  self%j_dofs_local = 0.0_f64
684 
685  ! Particle loop
686  do i_sp=1,self%particle_group%n_species
687  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
688  do i_part = 1,self%particle_group%group(i_sp)%n_particles
689  vi = self%particle_group%group(i_sp)%get_v(i_part)
690  xi = self%particle_group%group(i_sp)%get_x(i_part)
691 
692  ! Get charge for accumulation of j
693  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
694 
695  vbar = 0.5_f64 * (self%vnew(i_sp,:, i_part)+vi(1:2))
696 
697  xnew(1) = xi(1) + dt * vbar(1)
698  call compute_particle_boundary_current_evaluate( self, xi, xnew, vi, vbar, wi, qoverm, dt )
699 
700  call self%particle_group%group(i_sp)%set_x( i_part, xnew )
701  call self%particle_group%group(i_sp)%set_v( i_part, vi )
702  end do
703  end do
704 
705  self%j_dofs = 0.0_f64
706  ! MPI to sum up contributions from each processor
707  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(1:self%n_dofs1,1), &
708  self%n_dofs1, mpi_sum, self%j_dofs(1:self%n_dofs1,1) )
709  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
710  self%n_dofs0, mpi_sum, self%j_dofs(:,2) )
711 
712  call self%filter%apply_inplace( self%j_dofs(:,1) )
713  call self%filter%apply_inplace( self%j_dofs(:,2) )
714 
715  ! Update the electric field.
716  call self%maxwell_solver%compute_E_from_j(self%j_dofs(1:self%n_dofs1,1), 1, self%efield_dofs(1:self%n_dofs1,1) )
717  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, self%efield_dofs(:,2) )
718 
719  call self%filter%apply( self%efield_dofs(:,1), self%efield_filter(:,1) )
720  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
721 
722  self%iter_counter = self%iter_counter + 1
723  self%niter(self%iter_counter) = niter
724 
725  end subroutine advect_ex_pic_vm_1d2v_helper
726 
727 
729  subroutine compute_particle_boundary_current_evaluate( self, xi, xnew, vi, vbar, wi, qoverm, dt )
730  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent( inout ) :: self
731  sll_real64, intent( in ) :: xi(1)
732  sll_real64, intent( inout ) :: xnew(1)
733  sll_real64, intent( inout ) :: vi(2)
734  sll_real64, intent( inout ) :: vbar(2)
735  sll_real64, intent( in ) :: wi(1)
736  sll_real64, intent( in ) :: qoverm
737  sll_real64, intent( in ) :: dt
738  !local variables
739  sll_real64 :: xmid(1), xbar
740  sll_real64 :: efield(2)
741 
742  if(xnew(1) < -self%x_max .or. xnew(1) > 2._f64*self%x_max ) then
743  print*, xnew
744  sll_error("particle boundary", "particle out of bound")
745  else if(xnew(1) < self%x_min .or. xnew(1) > self%x_max )then
746  if(xnew(1) < self%x_min )then
747  xbar = self%x_min
748  self%counter_left = self%counter_left+1
749  else if(xnew(1) > self%x_max)then
750  xbar = self%x_max
751  self%counter_right = self%counter_right+1
752  end if
753 
754  xmid = xbar
755  if ( abs(vbar(1)) > 1.0d-16 ) then
756  call self%kernel_smoother_1%add_current_evaluate &
757  ( xi(1), xmid(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
758  efield(1) )
759 
760  call self%kernel_smoother_0%add_current_evaluate &
761  ( xi(1), xmid(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
762  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
763  efield(2) )
764  else
765  call self%kernel_smoother_1%evaluate &
766  (xi(1), self%efield_to_val(1:self%n_dofs1,1), efield(1) )
767  efield(1) = efield(1)*dt
768  call self%kernel_smoother_0%add_charge( xi(1), &
769  wi(1)* vbar(2)*dt, self%j_dofs_local(:,2) )
770  call self%kernel_smoother_0%evaluate &
771  (xi(1), self%efield_to_val(:,2), efield(2) )
772  efield(2) = efield(2)*dt
773  end if
774  vi = vi + qoverm * efield
775 
776  select case(self%boundary_particles)
777  case(sll_p_boundary_particles_reflection)
778  xnew(1) = 2._f64*xbar-xnew(1)
779  vi(1) = -vi(1)
780  vbar(1) = -vbar(1)
781  case(sll_p_boundary_particles_absorption)
782  case( sll_p_boundary_particles_periodic)
783  xnew(1) = self%x_min + modulo(xnew(1)-self%x_min, self%Lx)
784  xmid = self%x_max+self%x_min-xbar
785  case default
786  xnew(1) = self%x_min + modulo(xnew(1)-self%x_min, self%Lx)
787  xmid = self%x_max+self%x_min-xbar
788  end select
789 
790  if ( abs(vbar(1)) > 1.0d-16 ) then
791  call self%kernel_smoother_1%add_current_evaluate &
792  ( xmid(1), xnew(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
793  efield(1) )
794  call self%kernel_smoother_0%add_current_evaluate &
795  ( xmid(1), xnew(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
796  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
797  efield(2) )
798  vi = vi + qoverm * efield
799  end if
800  else
801  if ( abs(vbar(1)) > 1.0d-16 ) then
802  call self%kernel_smoother_1%add_current_evaluate &
803  ( xi(1), xnew(1), wi(1), vbar(1), self%efield_to_val(1:self%n_dofs1,1), self%j_dofs_local(1:self%n_dofs1,1), &
804  efield(1) )
805 
806  call self%kernel_smoother_0%add_current_evaluate &
807  ( xi(1), xnew(1), wi(1)*vbar(2)/vbar(1), vbar(1), &
808  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
809  efield(2) )
810  else
811  call self%kernel_smoother_1%evaluate &
812  (xi(1), self%efield_to_val(1:self%n_dofs1,1), efield(1) )
813  efield(1) = efield(1)*dt
814  call self%kernel_smoother_0%add_charge( xi(1), &
815  wi(1)* vbar(2)*dt, self%j_dofs_local(:,2) )
816  call self%kernel_smoother_0%evaluate &
817  (xi(1), self%efield_to_val(:,2), efield(2) )
818  efield(2) = efield(2)*dt
819  end if
820  vi = vi + qoverm * efield
821  end if
822 
823  end subroutine compute_particle_boundary_current_evaluate
824 
825 
826 
828  subroutine advect_e_sub_pic_vm_1d2v_helper ( self, dt )
829  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
830  sll_real64, intent(in) :: dt
831  ! local variables
832  sll_real64 :: residuals(4)
833  sll_real64 :: residuals_all(4)
834  sll_real64 :: residual
835  sll_real64 :: qoverm
836  sll_real64 :: dt_sub(self%particle_group%n_species)
837  sll_real64 :: xi(3), vi(3), wi(1)
838  sll_real64 :: efield_dofs(self%kernel_smoother_1%n_dofs,2)
839  sll_real64 :: efield_dofs_old(self%kernel_smoother_1%n_dofs,2)
840  sll_int32 :: i_part, i_sub, niter, i_sp
841 
842 
843 
844  do i_sp = 1, self%particle_group%n_species
845  dt_sub(i_sp) = dt/real( self%n_sub_iter(i_sp), f64 )
846  end do
847  ! Set to zero
848  self%j_dofs_local = 0.0_f64
849 
850  niter = 0
851 
852  ! Set initial value for electric field
853  ! Two options to initialize the efield: either by old time step value or by value from disgradE propagator
854  !call disgradE_for_start ( self, dt, efield_dofs )
855  efield_dofs = self%efield_dofs
856  do i_sp = 1, self%particle_group%n_species
857  do i_part = 1, self%particle_group%group(i_sp)%n_particles
858  vi = self%particle_group%group(i_sp)%get_v(i_part)
859  xi = self%particle_group%group(i_sp)%get_x(i_part)
860  self%xnew(i_sp,i_part) = xi(1)
861  self%vnew(i_sp,:,i_part) = vi(1:2)
862  end do
863  end do
864  efield_dofs_old = efield_dofs
865  residual = 1.0_f64
866  do while ( (residual > self%iter_tolerance) .and. niter < self%max_iter )
867  niter = niter + 1
868  residuals = 0.0_f64
869  efield_dofs = 0.5_f64*( self%efield_dofs + efield_dofs )
870  ! Set to zero
871  self%j_dofs_local = 0.0_f64
872  do i_sp = 1, self%particle_group%n_species
873  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
874  do i_part = 1,self%particle_group%group(i_sp)%n_particles
875  vi = self%particle_group%group(i_sp)%get_v(i_part)
876  xi = self%particle_group%group(i_sp)%get_x(i_part)
877 
878  ! Get charge for accumulation of j
879  wi = self%particle_group%group(i_sp)%get_charge(i_part)
880 
881  do i_sub = 1, self%n_sub_iter(i_sp)
882  call subcycle_xv( self, dt_sub(i_sp), qoverm, efield_dofs, wi, xi(1:1), &
883  vi(1:2), self%sub_iter_counter(i_sp) )
884  end do
885 
886  residuals(1) = residuals(1) + (xi(1)-self%xnew(i_sp,i_part))**2*abs(wi(1))
887  residuals(2) = residuals(2) + (vi(1)-self%vnew(i_sp,1,i_part))**2*abs(wi(1))
888  residuals(3) = residuals(3) + (vi(2)-self%vnew(i_sp,2,i_part))**2*abs(wi(1))
889 
890  self%xnew(i_sp,i_part) = xi(1)
891  self%vnew(i_sp,:,i_part) = vi(1:2)
892  end do
893  end do
894 
895  ! Update d_n
896  self%j_dofs = 0.0_f64
897  ! MPI to sum up contributions from each processor
898  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,1), &
899  self%kernel_smoother_1%n_dofs, mpi_sum, self%j_dofs(:,1) )
900  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
901  self%kernel_smoother_1%n_dofs, mpi_sum, self%j_dofs(:,2) )
902 
903  ! Update the electric field.
904  efield_dofs = self%efield_dofs
905  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,1), 1, efield_dofs(:,1) )
906  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, efield_dofs(:,2) )
907  residuals(4) = (sum((efield_dofs_old-efield_dofs)**2))*self%delta_x!maxval( efield_dofs_old - efield_dofs ) ! Here the error norm can be changed
908 
909  call sll_o_collective_allreduce( sll_v_world_collective, residuals, 4, mpi_max, residuals_all )
910  residuals_all = sqrt(residuals_all)
911  residual = residuals_all(4)!maxval(residuals_all) ! Here the error criterion for the iteration can be changed
912  efield_dofs_old = efield_dofs
913  end do
914  self%efield_dofs = efield_dofs
915  do i_sp =1,self%particle_group%n_species
916  do i_part = 1, self%particle_group%group(i_sp)%n_particles
917  vi(1:2) = self%vnew(i_sp,:,i_part)
918  xi(1) = modulo(self%xnew(i_sp, i_part), self%Lx )
919  call self%particle_group%group(i_sp)%set_v( i_part, vi )
920  call self%particle_group%group(i_sp)%set_x( i_part, xi )
921  end do
922  end do
923 
924  self%iter_counter = self%iter_counter + 1
925  self%niter(self%iter_counter) = niter
926  do i_sp = 1,self%particle_group%n_species
927  self%nsubiter(i_sp,self%iter_counter) = real(self%sub_iter_counter(i_sp), f64)/ &
928  real( self%particle_group%group(i_sp)%n_particles * niter * self%n_sub_iter(i_sp),f64)
929  end do
930  self%sub_iter_counter = 0
931 
932  end subroutine advect_e_sub_pic_vm_1d2v_helper
933 
934 
935 
937  subroutine subcycle_xv( self, dt, qoverm, efield_dofs, wi, position, velocity, sub_iter_counter )
938  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
939  sll_real64, intent(in) :: dt
940  sll_real64, intent(in) :: qoverm
941  sll_real64, intent(in) :: efield_dofs(:,:)
942  sll_real64, intent(in ) :: wi(1)
943  sll_real64, intent(inout) :: position(1)
944  sll_real64, intent(inout) :: velocity(2)
945  sll_int32, intent(inout) :: sub_iter_counter
946 
947  ! local variables
948  sll_real64 :: xnew(1), vnew(2), vbar, xold(1), vold(2)
949  sll_real64 :: efield(2)
950  sll_int32 :: niter
951  sll_real64 :: residuals(3), residuals_all(3), residual
952 
953  !xnew = position + dt * velocity(1)
954  vnew = velocity
955  xnew = position
956 
957  niter = 0
958  residual = 1.0_f64
959  do while ( (residual > self%iter_tolerance_sub) .and. niter < self%max_iter_sub )
960  niter = niter + 1
961  residuals = 0.0_f64
962  vbar = 0.5_f64*(velocity(1)+vnew(1))
963  vold = vnew
964  xold = xnew
965  xnew = position + dt*vbar
966 
967  if ( abs(vbar) > 1.0d-16 ) then
968 
969  select type ( q=> self%kernel_smoother_1 )
970  type is (sll_t_particle_mesh_coupling_spline_1d )
971  call q%evaluate_int &
972  ( position, xnew, vbar, efield_dofs(:,1), &
973  efield(1) )
974  type is (sll_t_particle_mesh_coupling_spline_smooth_1d )
975  call q%evaluate_int &
976  ( position, xnew, vbar, efield_dofs(:,1), &
977  efield(1) )
978  end select
979  select type ( q0=> self%kernel_smoother_0 )
980  type is (sll_t_particle_mesh_coupling_spline_1d )
981  call q0%evaluate_int &
982  ( position, xnew, vbar, &
983  efield_dofs(:,2), &
984  efield(2) )
985  type is (sll_t_particle_mesh_coupling_spline_smooth_1d )
986  call q0%evaluate_int &
987  ( position, xnew, vbar, &
988  efield_dofs(:,2), &
989  efield(2) )
990  end select
991  else
992  call self%kernel_smoother_1%evaluate &
993  (position, efield_dofs(:,1), efield(1) )
994  efield(1) = efield(1)*dt
995 
996  call self%kernel_smoother_0%evaluate &
997  (position, efield_dofs(:,2), efield(2) )
998  efield(2) = efield(2)*dt
999  end if
1000 
1001  vnew = velocity + qoverm * efield(1:2)
1002 
1003  residuals(1) = abs(xnew(1)-xold(1))
1004  residuals(2) = abs(vnew(1)-vold(1))
1005  residuals(3) = abs(vnew(2)-vold(2))
1006  call sll_o_collective_allreduce( sll_v_world_collective, residuals, 3, mpi_max, residuals_all )
1007  residual = maxval(residuals_all)
1008 
1009  end do
1010 
1011  if ( residual > self%iter_tolerance_sub ) then
1012  print*, 'Warning: Inner iteration of iteration no.', self%iter_counter+1 ,'did not converge.', residuals_all
1013  !self%n_failed = self%n_failed+1
1014  end if
1015 
1016  vbar = 0.5_f64*(velocity(1)+vnew(1))
1017  xnew = position + dt*vbar
1018  if ( abs(vbar) > 1.0d-16 ) then
1019 
1020  select type ( q=> self%kernel_smoother_1 )
1021  type is (sll_t_particle_mesh_coupling_spline_1d )
1022  call q%add_current_evaluate &
1023  ( position, xnew, wi(1), vbar, efield_dofs(:,1), self%j_dofs_local(:,1), &
1024  efield(1) )
1025  type is (sll_t_particle_mesh_coupling_spline_smooth_1d )
1026  call q%add_current_evaluate &
1027  ( position, xnew, wi(1), vbar, efield_dofs(:,1), self%j_dofs_local(:,1), &
1028  efield(1) )
1029  end select
1030  select type ( q0=> self%kernel_smoother_0 )
1031  type is (sll_t_particle_mesh_coupling_spline_1d )
1032  call q0%add_current_evaluate &
1033  ( position, xnew, wi(1)*0.5_f64*(velocity(2)+vnew(2))/vbar, vbar, &
1034  efield_dofs(:,2), self%j_dofs_local(:,2), &
1035  efield(2) )
1036  type is (sll_t_particle_mesh_coupling_spline_smooth_1d )
1037  call q0%add_current_evaluate &
1038  ( position, xnew, wi(1)*0.5_f64*(velocity(2)+vnew(2))/vbar, vbar, &
1039  efield_dofs(:,2), self%j_dofs_local(:,2), &
1040  efield(2) )
1041  end select
1042  else
1043  call self%kernel_smoother_1%evaluate &
1044  (position, efield_dofs(:,1), efield(1) )
1045  efield(1) = efield(1)*dt
1046 
1047  call self%kernel_smoother_0%add_charge( position, &
1048  wi(1)*0.5_f64*(velocity(2)+vnew(2))*dt, self%j_dofs_local(:,2) )
1049  call self%kernel_smoother_0%evaluate &
1050  (position, efield_dofs(:,2), efield(2) )
1051  efield(2) = efield(2)*dt
1052  end if
1053 
1054  vnew = velocity + qoverm*efield(1:2)
1055 
1056  velocity = vnew
1057  position = modulo(xnew, self%Lx)
1058 
1059  sub_iter_counter = sub_iter_counter + niter + 1
1060 
1061  end subroutine subcycle_xv
1062 
1063 
1064 !!$
1065 !!$ subroutine filter( field_in, field_out, n_dofs )
1066 !!$ sll_real64, intent(in) :: field_in(:,:) !< array for the coefficients of the efields
1067 !!$ sll_real64, intent(out) :: field_out(:,:) !< array for the coefficients of the efields
1068 !!$ sll_int32, intent(in) :: n_dofs
1069 !!$
1070 !!$ sll_int32 :: j
1071 !!$
1072 !!$ ! Filter
1073 !!$ field_out(1,:) = 0.25_f64*( field_in(1,:)*2.0_f64 + &
1074 !!$ field_in(n_dofs,:)+field_in(2,:))
1075 !!$ do j=2,n_dofs-1
1076 !!$ field_out(j,:) = 0.25_f64*( field_in(j,:)*2.0_f64 + &
1077 !!$ field_in(j-1,:)+field_in(j+1,:))
1078 !!$ end do
1079 !!$ field_out(n_dofs,:) = 0.25_f64*( field_in(n_dofs,:)*2.0_f64 + &
1080 !!$ field_in(n_dofs-1,:)+field_in(1,:))
1081 !!$
1082 !!$ end subroutine filter
1083 
1084 
1085 
1086  !---------------------------------------------------------------------------!
1088  subroutine initialize_pic_vm_1d2v_helper(&
1089  self, &
1090  maxwell_solver, &
1091  kernel_smoother_0, &
1092  kernel_smoother_1, &
1093  particle_group, &
1094  phi_dofs, &
1095  efield_dofs, &
1096  bfield_dofs, &
1097  x_min, &
1098  Lx, &
1099  filter, &
1100  build_particle_mass, &
1101  boundary_particles, &
1102  solver_tolerance, &
1103  iter_tolerance, max_iter, &
1104  force_sign, &
1105  control_variate, &
1106  i_weight, &
1107  betar, &
1108  jmean)
1109  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
1110  class(sll_c_maxwell_1d_base), target, intent(in) :: maxwell_solver
1111  class(sll_c_particle_mesh_coupling_1d), target, intent(in) :: kernel_smoother_0
1112  class(sll_c_particle_mesh_coupling_1d), target, intent(in) :: kernel_smoother_1
1113  class(sll_t_particle_array), target, intent(in) :: particle_group
1114  sll_real64, target, intent(in) :: phi_dofs(:)
1115  sll_real64, target, intent(in) :: efield_dofs(:,:)
1116  sll_real64, target, intent(in) :: bfield_dofs(:)
1117  sll_real64, intent(in) :: x_min
1118  sll_real64, intent(in) :: lx
1119  class( sll_c_filter_base_1d ), intent( in ), target :: filter
1120  logical, optional, intent(in) :: build_particle_mass
1121  sll_int32, optional, intent( in ) :: boundary_particles
1122  sll_real64, intent(in), optional :: solver_tolerance
1123  sll_real64, optional, intent( in ) :: iter_tolerance
1124  sll_int32, optional, intent( in ) :: max_iter
1125  sll_real64, optional, intent(in) :: force_sign
1126  class(sll_t_control_variates), optional, target, intent(in) :: control_variate
1127  sll_int32, optional, intent(in) :: i_weight
1128  sll_real64, optional, intent( in ) :: betar(2)
1129  logical, optional, intent(in) :: jmean
1130  !local variables
1131  sll_int32 :: ierr, i_sp
1132  sll_int32 :: n_span_particle_mass_0, n_span_particle_mass_1
1133 
1134  if (present(boundary_particles) ) then
1135  self%boundary_particles = boundary_particles
1136  end if
1137 
1138  if (present(solver_tolerance) ) then
1139  self%solver_tolerance = solver_tolerance
1140  end if
1141 
1142  if (present(iter_tolerance) ) then
1143  self%iter_tolerance = iter_tolerance
1144  end if
1145 
1146  if (present(max_iter) ) then
1147  self%max_iter = max_iter
1148  end if
1149 
1150  if( present(force_sign) )then
1151  self%force_sign = force_sign
1152  end if
1153 
1154  if (self%force_sign == 1._f64) then
1155  if( particle_group%group(1)%species%q > 0._f64) self%adiabatic_electrons = .true.
1156  end if
1157 
1158  if (present(betar)) then
1159  self%betar = betar!32.89_f64
1160  end if
1161 
1162  if (present(jmean)) then
1163  self%jmean = jmean
1164  end if
1165 
1166  if ( present( build_particle_mass )) self%build_particle_mass_loc = build_particle_mass
1167 
1168  if( self%boundary_particles < 5 ) self%boundary = .true.
1169 
1170  self%maxwell_solver => maxwell_solver
1171  self%kernel_smoother_0 => kernel_smoother_0
1172  self%kernel_smoother_1 => kernel_smoother_1
1173  self%particle_group => particle_group
1174  self%efield_dofs => efield_dofs
1175  self%bfield_dofs => bfield_dofs
1176 
1177  ! Check that n_dofs is the same for both kernel smoothers.
1178  sll_assert( self%kernel_smoother_0%n_cells == self%kernel_smoother_1%n_cells )
1179  sll_assert( self%kernel_smoother_0%n_cells == self%maxwell_solver%n_cells )
1180 
1181  self%n_cells = self%maxwell_solver%n_cells
1182  self%n_dofs0 = self%maxwell_solver%n_dofs0
1183  self%n_dofs1 = self%maxwell_solver%n_dofs1
1184 
1185  self%spline_degree = self%kernel_smoother_0%spline_degree
1186  self%x_min = x_min
1187  self%Lx = lx
1188  self%x_max = x_min +lx
1189  self%delta_x = self%maxwell_solver%delta_x
1190 
1191 
1192  sll_allocate( self%j_dofs(self%n_dofs0,2), ierr )
1193  sll_allocate( self%j_dofs_local(self%n_dofs0,2), ierr )
1194  sll_allocate(self%residual(self%n_dofs0*2), ierr)
1195 
1196 
1197 
1198  if ( self%maxwell_solver%s_deg_0 > -1 .and. self%build_particle_mass_loc .eqv. .true. ) then
1199  n_span_particle_mass_0 = 2*self%spline_degree+1
1200  n_span_particle_mass_1 = 2*self%spline_degree-1
1201  select type( pm=>kernel_smoother_0)
1202  type is ( sll_t_particle_mesh_coupling_spline_smooth_1d )
1203  self%smooth = .true.
1204  n_span_particle_mass_0 = 2*self%spline_degree+5
1205  n_span_particle_mass_1 = 2*self%spline_degree+3
1206  end select
1207  self%size_particle_mass_0 = n_span_particle_mass_0* self%n_dofs0
1208  self%size_particle_mass_1 = n_span_particle_mass_1* self%n_dofs1
1209 
1210  sll_allocate( self%particle_mass_1_local( n_span_particle_mass_1, self%n_dofs1), ierr )
1211  sll_allocate( self%particle_mass_0_local( n_span_particle_mass_0, self%n_dofs0), ierr )
1212  self%particle_mass_1_local = 0.0_f64
1213  self%particle_mass_0_local = 0.0_f64
1214  if( self%n_cells+self%spline_degree == self%maxwell_solver%n_dofs0 ) then
1215  allocate( sll_t_linear_operator_particle_mass_cl_1d :: self%particle_mass_1 )
1216  select type ( q => self%particle_mass_1 )
1217  type is ( sll_t_linear_operator_particle_mass_cl_1d )
1218  call q%create( self%spline_degree-1, self%n_dofs1 )
1219  end select
1220 
1221  allocate( sll_t_linear_operator_particle_mass_cl_1d :: self%particle_mass_0 )
1222  select type ( q => self%particle_mass_0 )
1223  type is ( sll_t_linear_operator_particle_mass_cl_1d )
1224  call q%create( self%spline_degree, self%n_dofs0 )
1225  end select
1226  else if ( self%n_cells == self%maxwell_solver%n_dofs0 ) then
1227  if( self%smooth ) then
1228  allocate( sll_t_linear_operator_particle_mass_smooth_1d :: self%particle_mass_1 )
1229  select type ( q => self%particle_mass_1 )
1230  type is ( sll_t_linear_operator_particle_mass_smooth_1d )
1231  call q%create( self%spline_degree-1, self%n_dofs1 )
1232  end select
1233 
1234  allocate( sll_t_linear_operator_particle_mass_smooth_1d :: self%particle_mass_0 )
1235  select type ( q => self%particle_mass_0 )
1236  type is ( sll_t_linear_operator_particle_mass_smooth_1d )
1237  call q%create( self%spline_degree, self%n_dofs0 )
1238  end select
1239  else
1240  allocate( sll_t_linear_operator_particle_mass_1d :: self%particle_mass_1 )
1241  select type ( q => self%particle_mass_1 )
1242  type is ( sll_t_linear_operator_particle_mass_1d )
1243  call q%create( self%spline_degree-1, self%n_dofs1 )
1244  end select
1245 
1246  allocate( sll_t_linear_operator_particle_mass_1d :: self%particle_mass_0 )
1247  select type ( q => self%particle_mass_0 )
1248  type is ( sll_t_linear_operator_particle_mass_1d )
1249  call q%create( self%spline_degree, self%n_dofs0 )
1250  end select
1251  end if
1252  end if
1253 
1254  if( self%adiabatic_electrons ) then
1255  call self%linear_op_schur_phiv%create( self%maxwell_solver, self%particle_mass_1 )
1256  call self%linear_solver_schur_phiv%create( self%linear_op_schur_phiv )
1257  self%linear_solver_schur_phiv%atol = self%solver_tolerance
1258  else
1259  call self%linear_operator_1%create( self%maxwell_solver, self%particle_mass_1, self%n_dofs1, self%maxwell_solver%s_deg_1 )
1260  call self%linear_solver_1%create( self%linear_operator_1 )
1261  self%linear_solver_1%atol = self%solver_tolerance
1262  !self%linear_solver_1%verbose = .true.
1263 
1264  call self%linear_operator_0%create( self%maxwell_solver, self%particle_mass_0, self%n_dofs0, self%maxwell_solver%s_deg_0 )
1265  call self%linear_solver_0%create( self%linear_operator_0 )
1266  self%linear_solver_0%atol = self%solver_tolerance
1267  !self%linear_solver_0%verbose = .true.
1268  end if
1269  end if
1270 
1271  self%n_particles_max = self%particle_group%group(1)%n_particles
1272  do i_sp = 2,self%particle_group%n_species
1273  self%n_particles_max = max(self%n_particles_max, self%particle_group%group(i_sp)%n_particles )
1274  end do
1275 
1276  sll_allocate( self%xnew(self%particle_group%n_species,self%n_particles_max), ierr )
1277  sll_allocate( self%vnew(self%particle_group%n_species,2,self%n_particles_max), ierr )
1278  sll_allocate( self%efield_dofs_new(self%n_dofs0,2), ierr )
1279  sll_allocate( self%efield_dofs_work(self%n_dofs0,2), ierr )
1280 
1281  allocate(self%sub_iter_counter( self%particle_group%n_species ) )
1282  self%sub_iter_counter = 0
1283  allocate(self%nsubiter(self%particle_group%n_species,100000 ) )
1284 
1285  ! For binomial filter
1286  self%filter => filter
1287  sll_allocate(self%efield_filter(self%n_dofs0,2), ierr)
1288  sll_allocate(self%bfield_filter(self%n_dofs1), ierr)
1289  sll_allocate(self%efield_to_val(self%n_dofs0,2), ierr)
1290  sll_allocate(self%bfield_to_val(self%n_dofs1), ierr)
1291 
1292  call self%filter%apply( self%efield_dofs(:,1), self%efield_filter(:,1) )
1293  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
1294  call self%filter%apply( self%bfield_dofs, self%bfield_filter )
1295 
1296  self%i_weight = 1
1297  if (present(i_weight)) self%i_weight = i_weight
1298  if(present(control_variate)) then
1299  allocate(self%control_variate )
1300  allocate(self%control_variate%cv(self%particle_group%n_species) )
1301  self%control_variate => control_variate
1302  end if
1303 
1304  end subroutine initialize_pic_vm_1d2v_helper
1305 
1307  subroutine initialize_file_pic_vm_1d2v_helper(&
1308  self, &
1309  maxwell_solver, &
1310  kernel_smoother_0, &
1311  kernel_smoother_1, &
1312  particle_group, &
1313  phi_dofs, &
1314  efield_dofs, &
1315  bfield_dofs, &
1316  x_min, &
1317  Lx, &
1318  filter, &
1319  filename, &
1320  build_particle_mass, &
1321  boundary_particles, &
1322  force_sign, &
1323  control_variate, &
1324  i_weight, &
1325  betar, &
1326  jmean)
1327  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(out) :: self
1328  class(sll_c_maxwell_1d_base), target, intent(in) :: maxwell_solver
1329  class(sll_c_particle_mesh_coupling_1d), target, intent(in) :: kernel_smoother_0
1330  class(sll_c_particle_mesh_coupling_1d), target, intent(in) :: kernel_smoother_1
1331  class(sll_t_particle_array), target, intent(in) :: particle_group
1332  sll_real64, target, intent(in) :: phi_dofs(:)
1333  sll_real64, target, intent(in) :: efield_dofs(:,:)
1334  sll_real64, target, intent(in) :: bfield_dofs(:)
1335  sll_real64, intent(in) :: x_min
1336  sll_real64, intent(in) :: lx
1337  class( sll_c_filter_base_1d ), intent( in ), target :: filter
1338  character(len=*), intent(in) :: filename
1339  logical, optional, intent(in) :: build_particle_mass
1340  sll_int32, optional, intent( in ) :: boundary_particles
1341  sll_real64, optional, intent(in) :: force_sign
1342  class(sll_t_control_variates), optional, target, intent(in) :: control_variate
1343  sll_int32, optional, intent(in) :: i_weight
1344  sll_real64, optional, intent( in ) :: betar(2)
1345  logical, optional, intent(in) :: jmean
1346  !local variables
1347  character(len=256) :: file_prefix
1348  sll_int32 :: input_file, rank
1349  sll_int32 :: io_stat, io_stat0, io_stat1, file_id, boundary_particles_set
1350  sll_real64 :: force_sign_set, betar_set(2)
1351  sll_real64 :: iter_tolerance, sub_iter_tolerance, maxwell_tolerance
1352  sll_int32 :: max_iter, max_iter_sub, n_sub_iter(2)
1353  logical :: jmean_set
1354 
1355 
1356  namelist /output/ file_prefix
1357  namelist /time_solver/ maxwell_tolerance
1358  namelist /time_iterate/ iter_tolerance, sub_iter_tolerance, max_iter, max_iter_sub, n_sub_iter
1359 
1360  rank = sll_f_get_collective_rank(sll_v_world_collective)
1361 
1362 
1363  if (present(boundary_particles)) then
1364  boundary_particles_set = boundary_particles
1365  else
1366  boundary_particles_set = 100
1367  end if
1368 
1369  if( present(force_sign) )then
1370  force_sign_set = force_sign
1371  else
1372  force_sign_set = 1._f64
1373  end if
1374 
1375  if( present(betar) )then
1376  betar_set = betar
1377  else
1378  betar_set = 1._f64
1379  end if
1380 
1381  if (present(jmean)) then
1382  jmean_set = jmean
1383  else
1384  jmean_set = .false.
1385  end if
1386 
1387  ! Read tolerance and max iter from nml file
1388 
1389  open(newunit = input_file, file=filename, status='old',iostat=io_stat)
1390  if (io_stat /= 0) then
1391  print*, 'initialize_pic_vm_1d2v_helper() failed to open nml-file.'
1392  else
1393 
1394  close(input_file)
1395  end if
1396 
1397  if (present(control_variate) ) then
1398  ! Read in solver tolerance
1399  open(newunit = input_file, file=filename, status='old',iostat=io_stat)
1400  if (io_stat /= 0) then
1401  if (rank == 0 ) then
1402  print*, 'sll_m_time_propagator_pic_vm_1d2v_helper: Input file does not exist. Set default tolerance.'
1403  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1404  write(file_id, *) 'solver tolerance:', 1d-12
1405  write(file_id, *) 'force_sign:', force_sign_set
1406  write(file_id, *) 'betar:', betar_set
1407  close(file_id)
1408  end if
1409  call self%init( maxwell_solver, &
1410  kernel_smoother_0, &
1411  kernel_smoother_1, &
1412  particle_group, &
1413  phi_dofs, &
1414  efield_dofs, &
1415  bfield_dofs, &
1416  x_min, &
1417  lx,&
1418  filter, &
1419  build_particle_mass, &
1420  boundary_particles = boundary_particles_set, &
1421  force_sign=force_sign_set, &
1422  control_variate = control_variate, &
1423  i_weight=i_weight, &
1424  betar=betar_set, &
1425  jmean=jmean_set)
1426  else
1427  read(input_file, output, iostat=io_stat0)
1428  read(input_file, time_solver,iostat=io_stat)
1429  read(input_file, time_iterate,iostat=io_stat1)
1430 
1431  self%output_prefix = trim(file_prefix)
1432  self%solver_tolerance = maxwell_tolerance
1433  self%iter_tolerance = iter_tolerance
1434  self%iter_tolerance_sub = sub_iter_tolerance
1435  self%max_iter = max_iter
1436  self%max_iter_sub = max_iter_sub
1437  self%n_sub_iter = n_sub_iter
1438  if (io_stat /= 0) then
1439  if (rank == 0 ) then
1440  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1441  write(file_id, *) 'solver tolerance:', 1d-12
1442  write(file_id, *) 'force_sign:', force_sign_set
1443  write(file_id, *) 'betar:', betar_set
1444  close(file_id)
1445  end if
1446  call self%init( maxwell_solver, &
1447  kernel_smoother_0, &
1448  kernel_smoother_1, &
1449  particle_group, &
1450  phi_dofs, &
1451  efield_dofs, &
1452  bfield_dofs, &
1453  x_min, &
1454  lx,&
1455  filter, &
1456  build_particle_mass, &
1457  boundary_particles = boundary_particles_set, &
1458  force_sign=force_sign_set, &
1459  control_variate = control_variate, &
1460  i_weight=i_weight, &
1461  betar=betar_set, &
1462  jmean=jmean_set)
1463  else
1464  if (rank == 0 ) then
1465  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1466  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1467  write(file_id, *) 'force_sign:', force_sign_set
1468  write(file_id, *) 'betar:', betar_set
1469  close(file_id)
1470  end if
1471  call self%init( maxwell_solver, &
1472  kernel_smoother_0, &
1473  kernel_smoother_1, &
1474  particle_group, &
1475  phi_dofs, &
1476  efield_dofs, &
1477  bfield_dofs, &
1478  x_min, &
1479  lx, &
1480  filter, &
1481  build_particle_mass, &
1482  boundary_particles = boundary_particles_set, &
1483  solver_tolerance=maxwell_tolerance, &
1484  force_sign=force_sign_set, &
1485  control_variate = control_variate, &
1486  i_weight=i_weight, &
1487  betar=betar_set, &
1488  jmean=jmean_set)
1489 
1490  end if
1491  close(input_file)
1492  end if
1493  else
1494  ! Read in solver tolerance
1495  open(newunit = input_file, file=filename, status='old',iostat=io_stat)
1496  if (io_stat /= 0) then
1497  if (rank == 0 ) then
1498  print*, 'sll_m_time_propagator_pic_vm_1d2v_helper: Input file does not exist. Set default tolerance.'
1499  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1500  write(file_id, *) 'solver tolerance:', 1d-12
1501  write(file_id, *) 'force_sign:', force_sign_set
1502  write(file_id, *) 'betar:', betar_set
1503  close(file_id)
1504  end if
1505  call self%init( maxwell_solver, &
1506  kernel_smoother_0, &
1507  kernel_smoother_1, &
1508  particle_group, &
1509  phi_dofs, &
1510  efield_dofs, &
1511  bfield_dofs, &
1512  x_min, &
1513  lx,&
1514  filter, &
1515  build_particle_mass, &
1516  boundary_particles = boundary_particles_set, &
1517  force_sign=force_sign_set, &
1518  betar=betar_set, &
1519  jmean=jmean_set)
1520  else
1521  read(input_file, output, iostat=io_stat0)
1522  read(input_file, time_solver,iostat=io_stat)
1523  read(input_file, time_iterate,iostat=io_stat1)
1524 
1525  self%output_prefix = trim(file_prefix)
1526  self%solver_tolerance = maxwell_tolerance
1527  self%iter_tolerance = iter_tolerance
1528  self%iter_tolerance_sub = sub_iter_tolerance
1529  self%max_iter = max_iter
1530  self%max_iter_sub = max_iter_sub
1531  self%n_sub_iter = n_sub_iter
1532 
1533  if (io_stat /= 0) then
1534  if (rank == 0 ) then
1535  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1536  write(file_id, *) 'solver tolerance:', 1d-12
1537  write(file_id, *) 'force_sign:', force_sign_set
1538  write(file_id, *) 'betar:', betar_set
1539  close(file_id)
1540  end if
1541  call self%init( maxwell_solver, &
1542  kernel_smoother_0, &
1543  kernel_smoother_1, &
1544  particle_group, &
1545  phi_dofs, &
1546  efield_dofs, &
1547  bfield_dofs, &
1548  x_min, &
1549  lx,&
1550  filter, &
1551  build_particle_mass, &
1552  boundary_particles = boundary_particles_set, &
1553  force_sign=force_sign_set, &
1554  betar=betar_set, &
1555  jmean=jmean_set)
1556  else
1557  if (rank == 0 ) then
1558  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1559  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1560  write(file_id, *) 'force_sign:', force_sign_set
1561  write(file_id, *) 'betar:', betar_set
1562  close(file_id)
1563  end if
1564  call self%init( maxwell_solver, &
1565  kernel_smoother_0, &
1566  kernel_smoother_1, &
1567  particle_group, &
1568  phi_dofs, &
1569  efield_dofs, &
1570  bfield_dofs, &
1571  x_min, &
1572  lx, &
1573  filter, &
1574  build_particle_mass, &
1575  boundary_particles = boundary_particles_set, &
1576  solver_tolerance=maxwell_tolerance, &
1577  force_sign=force_sign_set, &
1578  betar=betar_set, &
1579  jmean=jmean_set)
1580 
1581  end if
1582  close(input_file)
1583  end if
1584  end if
1585 
1586 
1587  end subroutine initialize_file_pic_vm_1d2v_helper
1588 
1589 
1590 
1591  !---------------------------------------------------------------------------!
1593  subroutine delete_pic_vm_1d2v_helper(self)
1594  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent( inout ) :: self
1595  !local variables
1596  sll_int32 :: file_id, j
1597 
1598  if( self%boundary ) then
1599  print*, 'left boundary', self%counter_left
1600  print*, 'right boundary', self%counter_right
1601  end if
1602 
1603  open(newunit=file_id, file=trim(self%output_prefix)//"-outer-iters.dat", action='write')
1604  do j=1,self%iter_counter
1605  write(file_id,*) self%niter(j)
1606  end do
1607  close(file_id)
1608  open(newunit=file_id, file=trim(self%output_prefix)//"-inner-iters.dat", action='write')
1609  do j=1, self%iter_counter
1610  write(file_id,*) self%nsubiter(:,j)
1611  end do
1612  close(file_id)
1613  print*, 'No. of failed iterations:', self%n_failed
1614 
1615  if ( self%maxwell_solver%s_deg_0 > -1 .and. self%build_particle_mass_loc .eqv. .true. ) then
1616  call self%linear_solver_1%free()
1617  call self%linear_operator_1%free()
1618  call self%particle_mass_1%free()
1619 
1620  call self%linear_solver_0%free()
1621  call self%linear_operator_0%free()
1622  call self%particle_mass_0%free()
1623 
1624  deallocate( self%particle_mass_1_local )
1625  deallocate( self%particle_mass_0_local )
1626 
1627  if(self%adiabatic_electrons) then
1628  call self%linear_solver_schur_phiv%free()
1629  call self%linear_op_schur_phiv%free()
1630  end if
1631  end if
1632 
1633  deallocate(self%j_dofs)
1634  deallocate(self%j_dofs_local)
1635  self%maxwell_solver => null()
1636  self%kernel_smoother_0 => null()
1637  self%kernel_smoother_1 => null()
1638  self%particle_group => null()
1639  self%efield_dofs => null()
1640  self%bfield_dofs => null()
1641 
1642  end subroutine delete_pic_vm_1d2v_helper
1643 
1644 
1645 
1646 
1648  subroutine advect_e_start_disgrade_pic_vm_1d2v_helper ( self, dt )
1649  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
1650  sll_real64, intent(in) :: dt
1651 
1652  ! local variables
1653  sll_int32 :: i_part, i_sp
1654  sll_int32 :: niter
1655  sll_real64 :: vi(3), xi(3), wi(1), xnew(3), vnew(3), vbar
1656  sll_real64 :: qoverm
1657  sll_real64 :: efield(2)
1658  sll_real64 :: residual, residuals(4), residuals_all(4)
1659  sll_real64 :: efield_dofs(self%kernel_smoother_1%n_dofs,2)
1660  sll_real64 :: efield_dofs_old(self%kernel_smoother_1%n_dofs,2)
1661 
1662  residuals = 0.0_f64
1663 
1664  ! Instead of a first particle loop as in the standard Picard iteration, we now
1665  ! perform the advect_x and advect_e part of the DISGRADE scheme
1666  efield_dofs = self%efield_dofs
1667  !call disgradE_for_start ( self, dt, efield_dofs )
1668  do i_sp = 1, self%particle_group%n_species
1669  do i_part = 1, self%particle_group%group(i_sp)%n_particles
1670  vnew = self%particle_group%group(i_sp)%get_v(i_part)
1671  xnew = self%particle_group%group(i_sp)%get_x(i_part)
1672  self%xnew(i_sp,i_part) = xnew(1)
1673  self%vnew(i_sp,:,i_part) = vnew(1:2)
1674  end do
1675  end do
1676 
1677  efield_dofs_old = efield_dofs
1678  residual = 1.0_f64
1679  niter = 1
1680 
1681  !print*, 'Iteration', self%iter_counter+1
1682  do while ( (residual > self%iter_tolerance) .and. niter < self%max_iter )
1683  niter = niter + 1
1684  residuals = 0.0_f64
1685  call self%filter%apply( efield_dofs_old(:,1), efield_dofs(:,1) )
1686  call self%filter%apply( efield_dofs_old(:,2), efield_dofs(:,2) )
1687  efield_dofs = 0.5_f64*(self%efield_filter + efield_dofs )
1688 
1689  call self%maxwell_solver%transform_dofs( efield_dofs(:,1), self%efield_to_val(:,1), 0 )
1690  call self%maxwell_solver%transform_dofs( efield_dofs(:,2), self%efield_to_val(:,2), 1 )
1691 
1692  ! Set to zero
1693  self%j_dofs_local = 0.0_f64
1694  ! First particle loop
1695  do i_sp = 1, self%particle_group%n_species
1696  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
1697  do i_part = 1,self%particle_group%group(i_sp)%n_particles
1698  vi = self%particle_group%group(i_sp)%get_v(i_part)
1699  xi = self%particle_group%group(i_sp)%get_x(i_part)
1700 
1701  ! Get charge for accumulation of j
1702  wi = self%particle_group%group(i_sp)%get_charge(i_part)
1703 
1704  vnew(1:2) = self%vnew(i_sp,1:2, i_part)
1705  xnew(1) = xi(1) + dt * 0.5_f64 * ( vi(1) + vnew(1) )
1706 
1707 
1708  vbar = 0.5_f64*(vi(1)+vnew(1))
1709  if ( abs(vbar) > 1.0d-16 ) then
1710  call self%kernel_smoother_1%add_current_evaluate &
1711  ( xi(1), xnew(1), wi(1), vbar, self%efield_to_val(:,1), self%j_dofs_local(:,1), &
1712  efield(1) )
1713 
1714  call self%kernel_smoother_0%add_current_evaluate &
1715  ( xi(1), xnew(1), wi(1)*(vi(2)+vnew(2))/(vi(1)+vnew(1)), vbar, &
1716  self%efield_to_val(:,2), self%j_dofs_local(:,2), &
1717  efield(2) )
1718  else
1719  call self%kernel_smoother_1%evaluate &
1720  (xi(1), self%efield_to_val(:,1), efield(1) )
1721  efield(1) = efield(1)*dt
1722  call self%kernel_smoother_0%add_charge( xi(1), &
1723  wi(1)* 0.5_f64*(vi(2)+vnew(2))*dt, self%j_dofs_local(:,2) )
1724  call self%kernel_smoother_0%evaluate &
1725  (xi(1), self%efield_to_val(:,2), efield(2) )
1726  efield(2) = efield(2)*dt
1727  end if
1728  vnew(1:2) = vi(1:2) + qoverm * efield(1:2)
1729 
1730  ! Here yo can change the residual criterion
1731  residuals(1) = residuals(1) + (xnew(1)-self%xnew(i_sp,i_part))**2*abs(wi(1))!max( residuals(1), abs(xnew(1)-self%xnew(i_sp,i_part)) )
1732  residuals(2) = residuals(2) + (vnew(1)-self%vnew(i_sp,1,i_part))**2*abs(wi(1))!max( residuals(2), abs(self%vnew(i_sp,1,i_part)-vnew(1)) )
1733  residuals(3) = residuals(3) + (vnew(2)-self%vnew(i_sp,2,i_part))**2*abs(wi(1))!max( residuals(3), abs(self%vnew(i_sp,2,i_part)-vnew(2)) )
1734  self%xnew(i_sp,i_part) = xnew(1)
1735  self%vnew(i_sp,:,i_part) = vnew(1:2)
1736  end do
1737  end do
1738 
1739  ! Update d_n
1740  self%j_dofs = 0.0_f64
1741  ! MPI to sum up contributions from each processor
1742  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,1), &
1743  self%kernel_smoother_1%n_dofs, mpi_sum, self%j_dofs(:,1) )
1744  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
1745  self%kernel_smoother_1%n_dofs, mpi_sum, self%j_dofs(:,2) )
1746 
1747  call self%filter%apply_inplace( self%j_dofs(:,1) )
1748  call self%filter%apply_inplace( self%j_dofs(:,2) )
1749 
1750  ! Filter
1751  !call filter( self%j_dofs, self%j_dofs_local, self%kernel_smoother_1%n_dofs )
1752  !self%j_dofs = self%j_dofs_local
1753 
1754  ! FFT filter
1755  !write(15, *) self%j_dofs(:,1)
1756  !call sll_s_fft_exec_r2r_1d ( self%filter_fw, self%j_dofs(:,1), self%j_dofs_local(:,1) )
1757  !write(14,*) self%j_dofs_local(:,1)
1758  !self%j_dofs_local(ndh+1-5:ndh+1+5,1) = 0.0_f64
1759  !call sll_s_fft_exec_r2r_1d ( self%filter_bw, self%j_dofs_local(:,1), self%j_dofs(:,1) )
1760  !write(16,*) self%j_dofs(:,1)
1761  !stop
1762 
1763  ! Update the electric field.
1764  efield_dofs = self%efield_dofs
1765  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,1), 1, efield_dofs(:,1) )
1766  call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, efield_dofs(:,2) )
1767 
1768  if ( self%maxwell_solver%strong_ampere .eqv. .true. ) then
1769  residuals(4) = self%maxwell_solver%l2norm_squared( efield_dofs_old(:,1)-&
1770  efield_dofs(:,1), self%maxwell_solver%s_deg_0 ) + &
1771  self%maxwell_solver%l2norm_squared( efield_dofs_old(:,2)-&
1772  efield_dofs(:,2), self%maxwell_solver%s_deg_0 -1 )
1773  else
1774  residuals(4) = self%maxwell_solver%l2norm_squared( efield_dofs_old(:,1)-&
1775  efield_dofs(:,1), self%maxwell_solver%s_deg_0 - 1 ) + &
1776  self%maxwell_solver%l2norm_squared( efield_dofs_old(:,2)-&
1777  efield_dofs(:,2), self%maxwell_solver%s_deg_0 )
1778  end if
1779  ! Here you can change the residual norm
1780  residuals(4) = (sum((efield_dofs_old-efield_dofs)**2))*self%delta_x!maxval( efield_dofs_old - efield_dofs )
1781 
1782 
1783  call sll_o_collective_allreduce( sll_v_world_collective, residuals, 4, mpi_max, residuals_all )
1784 
1785  residuals_all = sqrt(residuals_all)
1786  residual = residuals_all(4)!max(residuals_all(4)), maxval(residuals_all(1:3))*0.01_f64)!residuals_all(4)!maxval(residuals_all)!residuals_all(4)! maxval(residuals_all)
1787  !call filter( efield_dofs, self%j_dofs_local, self%kernel_smoother_1%n_dofs )
1788  !efield_dofs = self%j_dofs_local
1789  ! FFT filter
1790  !write(15, *) self%j_dofs(:,1)
1791  !call sll_s_fft_exec_r2r_1d ( self%filter_fw, efield_dofs(:,1), self%j_dofs_local(:,1) )
1792  !write(14,*) self%j_dofs_local(:,1)
1793  !self%j_dofs_local(ndh+1-5:ndh+1+5,1) = 0.0_f64
1794  !call sll_s_fft_exec_r2r_1d ( self%filter_bw, self%j_dofs_local(:,1), efield_dofs(:,1) )
1795  !write(16,*) self%j_dofs(:,1)
1796  !stop
1797 
1798  efield_dofs_old = efield_dofs
1799  end do
1800 
1801 !!$ ! NEW VERSION WITH SWAPP TO OTHER ITERATION
1802 !!$ if ( residual > self%iter_tolerance ) then
1803 !!$ print*, 'Warning: Iteration no.', self%iter_counter+1 ,'did not converge.', residuals_all, niter
1804 !!$ self%n_failed = self%n_failed+1
1805 !!$ call advect_e_newtondisgradE(self, dt )
1806 
1807 !!$ self%iter_failed = .true.
1808 !!$ else
1809 !!$
1810 !!$ !print*, maxval(abs(self%efield_dofs-efield_dofs))
1811 !!$ self%efield_dofs = efield_dofs
1812 !!$ do i_sp = 1, self%particle_group%n_species
1813 !!$ do i_part = 1, self%particle_group%group(i_sp)%n_particles
1814 !!$ vnew(1:2) = self%vnew(i_sp,:,i_part)
1815 !!$ xnew(1) = modulo(self%xnew(i_sp,i_part), self%Lx)
1816 !!$ call self%particle_group%group(i_sp)%set_v( i_part, vnew )
1817 !!$ call self%particle_group%group(i_sp)%set_x( i_part, xnew )
1818 !!$ end do
1819 !!$ end do
1820 !!$
1821 !!$ self%iter_counter = self%iter_counter + 1
1822 !!$ self%niter(self%iter_counter) = niter
1823 !!$ end if
1824 !!$ ! OLD VERSION TAKING FAILED ITERATION
1825  if ( residual > self%iter_tolerance ) then
1826  print*, 'Warning: Iteration no.', self%iter_counter+1 ,'did not converge.', residuals_all, niter
1827  self%n_failed = self%n_failed+1
1828  end if
1829 
1830  !print*, maxval(abs(self%efield_dofs-efield_dofs))
1831  self%efield_dofs = efield_dofs
1832  call self%filter%apply( self%efield_dofs(:,1), self%efield_filter(:,1) )
1833  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
1834  do i_sp = 1, self%particle_group%n_species
1835  do i_part = 1, self%particle_group%group(i_sp)%n_particles
1836  vnew(1:2) = self%vnew(i_sp,:,i_part)
1837  xnew(1) = modulo(self%xnew(i_sp,i_part), self%Lx)
1838  call self%particle_group%group(i_sp)%set_v( i_part, vnew )
1839  call self%particle_group%group(i_sp)%set_x( i_part, xnew )
1840  end do
1841  end do
1842 
1843  self%iter_counter = self%iter_counter + 1
1844  self%niter(self%iter_counter) = niter
1845 
1846  end subroutine advect_e_start_disgrade_pic_vm_1d2v_helper
1847 
1848 
1851  subroutine disgrade_for_start ( self, dt, efield_dofs )
1852  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
1853  sll_real64, intent(in) :: dt
1854  sll_real64, intent(out) :: efield_dofs(:,:)
1855  ! local variables
1856  sll_int32 :: i_part, i_sp
1857  sll_real64 :: vi(3), xi(3), wi(1)
1858  sll_real64 :: qoverm, factor
1859  sll_real64 :: efield(2)
1860  sll_real64 :: rhs0(self%n_dofs0)
1861  sll_real64 :: rhs1(self%n_dofs1)
1862 
1863  ! Set to zero
1864  self%j_dofs_local = 0.0_f64
1865  self%particle_mass_1_local = 0.0_f64
1866  self%particle_mass_0_local = 0.0_f64
1867 
1868 
1869  ! First particle loop
1870  do i_sp=1,self%particle_group%n_species
1871  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
1872  factor = dt**2*0.25_f64*self%betar(2)*qoverm
1873  do i_part = 1,self%particle_group%group(i_sp)%n_particles
1874  vi = self%particle_group%group(i_sp)%get_v(i_part)
1875  xi = self%particle_group%group(i_sp)%get_x(i_part)
1876  xi(1) = xi(1) + 0.5_f64 * dt * vi(1)
1877  xi(1) = modulo(xi(1), self%Lx)
1878 
1879 
1880  ! Get charge for accumulation of j
1881  wi = self%particle_group%group(i_sp)%get_charge(i_part)
1882 
1883  ! Accumulate the particle mass matrix diagonals
1884  call self%kernel_smoother_1%add_particle_mass_full( xi(1), &
1885  wi(1) * factor, &
1886  self%particle_mass_1_local )
1887  call self%kernel_smoother_0%add_particle_mass_full( xi(1), &
1888  wi(1) * factor, &
1889  self%particle_mass_0_local )
1890  ! Accumulate jx
1891  call self%kernel_smoother_1%add_charge( xi(1), wi(1)*vi(1), &
1892  self%j_dofs_local(:,1) )
1893  ! Accumulate jy
1894  call self%kernel_smoother_0%add_charge( xi(1), wi(1)*vi(2), &
1895  self%j_dofs_local(:,2) )
1896  ! Evaulate E_x at particle position and propagate v a half step
1897  call self%kernel_smoother_1%evaluate &
1898  (xi(1), self%efield_dofs(:,1), efield(1) )
1899  ! Evaulate E_y at particle position and propagate v a half step
1900  call self%kernel_smoother_0%evaluate &
1901  (xi(1), self%efield_dofs(:,2), efield(2))
1902  vi(1:2) = vi(1:2) + dt* 0.5_f64* qoverm * efield
1903  self%vnew(i_sp,:,i_part) = vi(1:2)
1904  self%xnew(i_sp,i_part) = xi(1)
1905 
1906  end do
1907  end do
1908 
1909  ! Update d_n
1910  self%j_dofs = 0.0_f64
1911  ! MPI to sum up contributions from each processor
1912  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(1:self%n_dofs1,1), &
1913  self%n_dofs1, mpi_sum, self%j_dofs(1:self%n_dofs1,1) )
1914  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
1915  self%n_dofs0, mpi_sum, self%j_dofs(:,2) )
1916 
1917  self%particle_mass_1%particle_mass = 0.0_f64
1918  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_1_local, &
1919  self%size_particle_mass_1, mpi_sum, self%particle_mass_1%particle_mass)
1920  self%particle_mass_0%particle_mass = 0.0_f64
1921  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_0_local, &
1922  self%size_particle_mass_0, mpi_sum, self%particle_mass_0%particle_mass)
1923 
1924 
1925  ! Compute rhs
1926  self%linear_operator_1%sign = -self%force_sign
1927  call self%linear_operator_1%dot(self%efield_dofs(1:self%n_dofs1,1) , rhs1 )
1928  rhs1 = rhs1 - dt * self%force_sign * self%j_dofs(1:self%n_dofs1,1)
1929  ! Solve the Schur complement
1930  self%linear_operator_1%sign = self%force_sign
1931  call self%linear_solver_1%set_guess(self%efield_dofs(1:self%n_dofs1,1))
1932  call self%linear_solver_1%solve( rhs1, efield_dofs(1:self%n_dofs1,1) )
1933 
1934  ! Compute rhs
1935  self%linear_operator_0%sign = -self%force_sign
1936  call self%linear_operator_0%dot(self%efield_dofs(:,2) , rhs0 )
1937  rhs0 = rhs0 - self%force_sign * dt * self%j_dofs(:,2)
1938 
1939  ! Solve the Schur complement
1940  self%linear_operator_0%sign = self%force_sign
1941  call self%linear_solver_0%set_guess(self%efield_dofs(:,2))
1942  call self%linear_solver_0%solve( rhs0, efield_dofs(:,2) )
1943 
1944 
1945  ! Second particle loop (second half step of particle propagation)
1946  do i_sp=1,self%particle_group%n_species
1947  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
1948  do i_part = 1, self%particle_group%group(i_sp)%n_particles
1949 
1950  vi(1:2) = self%vnew(i_sp,:,i_part)
1951  xi(1) = self%xnew(i_sp,i_part)
1952  ! Evaulate E_x at particle position and propagate v a half step
1953  call self%kernel_smoother_1%evaluate &
1954  (xi(1), efield_dofs(:,1), efield(1))
1955  ! Evaulate E_y at particle position and propagate v a half step
1956  call self%kernel_smoother_0%evaluate &
1957  (xi(1), efield_dofs(:,2), efield(2))
1958  vi(1:2) = vi(1:2) + dt* 0.5_f64* qoverm * efield
1959  self%vnew(i_sp,:,i_part) = vi(1:2)
1960  self%xnew(i_sp,i_part) = modulo(xi(1) + 0.5_f64*dt*vi(1), self%Lx)
1961  end do
1962  end do
1963 
1964 
1965  end subroutine disgrade_for_start
1966 
1967 
1968 
1970  subroutine reinit_fields( self )
1971  class(sll_t_time_propagator_pic_vm_1d2v_helper), intent(inout) :: self
1972 
1973  call self%filter%apply( self%efield_dofs(:,1), self%efield_filter(:,1) )
1974  call self%filter%apply( self%efield_dofs(:,2), self%efield_filter(:,2) )
1975  call self%filter%apply( self%bfield_dofs, self%bfield_filter )
1976 
1977  end subroutine reinit_fields
1978 
1979 
1980 end module sll_m_time_propagator_pic_vm_1d2v_helper
sll_m_collective::sll_o_collective_allreduce
Combines values from all processes and distributes the result back to all processes.
Definition: sll_m_collective.F90:312
sll_m_collective
Parallelizing facility.
Definition: sll_m_collective.F90:128
sll_m_collective::sll_f_get_collective_rank
integer(kind=i32) function, public sll_f_get_collective_rank(col)
Determines the rank of the calling process in the communicator.
Definition: sll_m_collective.F90:533
sll_m_collective::sll_v_world_collective
type(sll_t_collective_t), pointer, public sll_v_world_collective
Control of subset assignment. Processes with the same color are in the same new communicator.
Definition: sll_m_collective.F90:255
sll_m_control_variate
Control variate.
Definition: sll_m_control_variate.F90:6
sll_m_filter_base_1d
filter base class
Definition: sll_m_filter_base_1d.F90:6
sll_m_linear_operator_particle_mass_1d
Definition: sll_m_linear_operator_particle_mass_1d.F90:1
sll_m_linear_operator_particle_mass_cl_1d
Definition: sll_m_linear_operator_particle_mass_cl_1d.F90:1
sll_m_linear_operator_particle_mass_smooth_1d
Definition: sll_m_linear_operator_particle_mass_smooth_1d.F90:1
sll_m_linear_operator_schur_ev_1d
Definition: sll_m_linear_operator_schur_ev_1d.F90:1
sll_m_linear_operator_schur_phiv_1d
Definition: sll_m_linear_operator_schur_phiv_1d.F90:1
sll_m_linear_solver_cg
module for conjugate gradient method in pure form
Definition: sll_m_linear_solver_cg.F90:11
sll_m_maxwell_1d_base
Module interface to solve Maxwell's equations in 1D.
Definition: sll_m_maxwell_1d_base.F90:8
sll_m_particle_group_base
Definition: sll_m_particle_group_base.F90:1
sll_m_particle_mass_1d_base
Definition: sll_m_particle_mass_1d_base.F90:1
sll_m_particle_mesh_coupling_base_1d
Base class for kernel smoothers for accumulation and field evaluation in PIC.
Definition: sll_m_particle_mesh_coupling_base_1d.F90:5
sll_m_particle_mesh_coupling_spline_1d
Kernel smoother for 1d with splines of arbitrary degree placed on a uniform mesh.
Definition: sll_m_particle_mesh_coupling_spline_1d.F90:5
sll_m_particle_mesh_coupling_spline_smooth_1d
Kernel smoother for 2d with splines of arbitrary degree placed on a uniform mesh.
Definition: sll_m_particle_mesh_coupling_spline_smooth_1d.F90:5
sll_m_particle_mesh_coupling_spline_strong_1d
Kernel smoother for 1d with splines of arbitrary degree placed on a uniform mesh. This version is for...
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:6
sll_m_time_propagator_base
Base class for Hamiltonian splittings.
Definition: sll_m_time_propagator_base.F90:6
sll_m_time_propagator_pic_vm_1d2v_helper
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1
sll_m_time_propagator_pic_vm_1d2v_helper::subcycle_xv
subroutine subcycle_xv(self, dt, qoverm, efield_dofs, wi, position, velocity, sub_iter_counter)
Helper function for subcycle (using Picard iteration)
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:938
sll_m_time_propagator_pic_vm_1d2v_helper::advect_e_sub_pic_vm_1d2v_helper
subroutine advect_e_sub_pic_vm_1d2v_helper(self, dt)
Operator for e,x-part with subcycling.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:829
sll_m_time_propagator_pic_vm_1d2v_helper::advect_e_start_disgrade_pic_vm_1d2v_helper
subroutine advect_e_start_disgrade_pic_vm_1d2v_helper(self, dt)
Operator for first variant without subcycling (Picard iteration, started by DISGRADE)
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1649
sll_m_time_propagator_pic_vm_1d2v_helper::initialize_file_pic_vm_1d2v_helper
subroutine initialize_file_pic_vm_1d2v_helper(self, maxwell_solver, kernel_smoother_0, kernel_smoother_1, particle_group, phi_dofs, efield_dofs, bfield_dofs, x_min, Lx, filter, filename, build_particle_mass, boundary_particles, force_sign, control_variate, i_weight, betar, jmean)
Constructor.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1327
sll_m_time_propagator_pic_vm_1d2v_helper::disgrade_for_start
subroutine disgrade_for_start(self, dt, efield_dofs)
DISGRADE as first guess DisgradE as first guess.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1852
sll_m_time_propagator_pic_vm_1d2v_helper::initialize_pic_vm_1d2v_helper
subroutine initialize_pic_vm_1d2v_helper(self, maxwell_solver, kernel_smoother_0, kernel_smoother_1, particle_group, phi_dofs, efield_dofs, bfield_dofs, x_min, Lx, filter, build_particle_mass, boundary_particles, solver_tolerance, iter_tolerance, max_iter, force_sign, control_variate, i_weight, betar, jmean)
Constructor.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1109
sll_m_time_propagator_pic_vm_1d2v_helper::advect_vb_pic_vm_1d2v_helper
subroutine advect_vb_pic_vm_1d2v_helper(self, dt)
advect_vb: Equations to be solved $V_{1,2}^{n+1}= (\cos(B_3)&\sin(B) \ -\sin(B) &\cos(B) ) V_{1,...
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:262
sll_m_time_propagator_pic_vm_1d2v_helper::advect_e_pic_vm_1d2v_helper
subroutine advect_e_pic_vm_1d2v_helper(self, dt)
Solution with Schur complement: $ S_{+}=M_1+\frac{\Delta t^2 q^2}{4 m} (\mathbb{\Lambda}^1)^T \mathbb...
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:326
sll_m_time_propagator_pic_vm_1d2v_helper::advect_ex_pic_vm_1d2v_helper
subroutine advect_ex_pic_vm_1d2v_helper(self, dt)
Operator for first variant without subcycling (Picard iteration, started by DISGRADE)
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:484
sll_m_time_propagator_pic_vm_1d2v_helper::reinit_fields
subroutine reinit_fields(self)
Computes the filtered dofs.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1971
sll_m_time_propagator_pic_vm_1d2v_helper::compute_particle_boundary_current_evaluate
subroutine compute_particle_boundary_current_evaluate(self, xi, xnew, vi, vbar, wi, qoverm, dt)
Helper function for advect_ex.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:730
sll_m_time_propagator_pic_vm_1d2v_helper::advect_x_pic_vm_1d2v_helper
subroutine advect_x_pic_vm_1d2v_helper(self, dt)
Advection of x part separately.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:190
sll_m_time_propagator_pic_vm_1d2v_helper::delete_pic_vm_1d2v_helper
subroutine delete_pic_vm_1d2v_helper(self)
Destructor.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:1594
sll_m_time_propagator_pic_vm_1d2v_helper::compute_particle_boundary
subroutine compute_particle_boundary(self, xold, xnew, vi)
Helper function for advect_x.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:221
sll_m_time_propagator_pic_vm_1d2v_helper::advect_eb_pic_vm_1d2v_helper
subroutine advect_eb_pic_vm_1d2v_helper(self, dt)
advect_eb: Equations to be solved Solution with Schur complement: $ S=M_1+\frac{\Delta t^2}{4} D^\top...
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:308
sll_m_time_propagator_pic_vm_3d3v_cl_helper
Particle pusher based on antisymmetric splitting with AVF for 3d3v Vlasov-Maxwell with coordinate tra...
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:5
sll_m_time_propagator_pic_vm_3d3v_cl_helper::sll_p_boundary_particles_periodic
integer(kind=i32), parameter, public sll_p_boundary_particles_periodic
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:38
sll_m_time_propagator_pic_vm_3d3v_cl_helper::sll_p_boundary_particles_reflection
integer(kind=i32), parameter, public sll_p_boundary_particles_reflection
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:40
sll_m_time_propagator_pic_vm_3d3v_cl_helper::sll_p_boundary_particles_absorption
integer(kind=i32), parameter, public sll_p_boundary_particles_absorption
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:41
sll_m_time_propagator_pic_vm_3d3v_cl_helper::sll_p_boundary_particles_singular
integer(kind=i32), parameter, public sll_p_boundary_particles_singular
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:39
control_variate
real(kind=f64) function, dimension(size(particle, 2)) control_variate(particle)
Definition: sll_pif_general.F90:314
sll_m_control_variate::sll_t_control_variates
Definition: sll_m_control_variate.F90:40
sll_m_filter_base_1d::sll_c_filter_base_1d
Definition: sll_m_filter_base_1d.F90:19
sll_m_linear_operator_particle_mass_1d::sll_t_linear_operator_particle_mass_1d
Definition: sll_m_linear_operator_particle_mass_1d.F90:14
sll_m_linear_operator_particle_mass_cl_1d::sll_t_linear_operator_particle_mass_cl_1d
Definition: sll_m_linear_operator_particle_mass_cl_1d.F90:13
sll_m_linear_operator_particle_mass_smooth_1d::sll_t_linear_operator_particle_mass_smooth_1d
Definition: sll_m_linear_operator_particle_mass_smooth_1d.F90:14
sll_m_linear_operator_schur_ev_1d::sll_t_linear_operator_schur_ev_1d
Definition: sll_m_linear_operator_schur_ev_1d.F90:18
sll_m_linear_operator_schur_phiv_1d::sll_t_linear_operator_schur_phiv_1d
Definition: sll_m_linear_operator_schur_phiv_1d.F90:19
sll_m_linear_solver_cg::sll_t_linear_solver_cg
class for the cg linear solver
Definition: sll_m_linear_solver_cg.F90:38
sll_m_maxwell_1d_base::sll_c_maxwell_1d_base
Definition: sll_m_maxwell_1d_base.F90:25
sll_m_particle_group_base::sll_t_particle_array
Definition: sll_m_particle_group_base.F90:100
sll_m_particle_mass_1d_base::sll_c_particle_mass_1d_base
Definition: sll_m_particle_mass_1d_base.F90:13
sll_m_particle_mesh_coupling_base_1d::sll_c_particle_mesh_coupling_1d
Basic type of a kernel smoother used for PIC simulations.
Definition: sll_m_particle_mesh_coupling_base_1d.F90:25
sll_m_particle_mesh_coupling_spline_1d::sll_t_particle_mesh_coupling_spline_1d
Spline kernel smoother in1d.
Definition: sll_m_particle_mesh_coupling_spline_1d.F90:46
sll_m_particle_mesh_coupling_spline_smooth_1d::sll_t_particle_mesh_coupling_spline_smooth_1d
Spline kernel smoother in1d.
Definition: sll_m_particle_mesh_coupling_spline_smooth_1d.F90:46
sll_m_particle_mesh_coupling_spline_strong_1d::sll_t_particle_mesh_coupling_spline_strong_1d
Basic type of a kernel smoother used for PIC simulations.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:38
sll_m_time_propagator_base::sll_c_time_propagator_base
Type for Hamiltonian splittings.
Definition: sll_m_time_propagator_base.F90:22
sll_m_time_propagator_pic_vm_1d2v_helper::sll_t_time_propagator_pic_vm_1d2v_helper
Hamiltonian splitting type for Vlasov-Maxwell 1d2v.
Definition: sll_m_time_propagator_pic_vm_1d2v_helper.F90:81
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1