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_3d3v_trafo_helper.F90
Go to the documentation of this file.
1 
3 module sll_m_time_propagator_pic_vm_3d3v_trafo_helper
4  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5 #include "sll_assert.h"
6 #include "sll_errors.h"
7 #include "sll_memory.h"
8 #include "sll_working_precision.h"
9 
10  use sll_m_collective, only: &
11  sll_f_get_collective_rank, &
12  sll_o_collective_allreduce, &
13  sll_v_world_collective
14 
15  use sll_m_control_variate, only: &
16  sll_t_control_variates
17 
18  use sll_m_time_propagator_pic_vm_3d3v_cl_helper, only: &
19  sll_p_boundary_particles_periodic, &
20  sll_p_boundary_particles_singular, &
21  sll_p_boundary_particles_reflection, &
22  sll_p_boundary_particles_absorption, &
23  sll_s_compute_particle_boundary_simple, &
24  sll_s_compute_matrix_inverse
25 
26  use sll_m_initial_distribution, only: &
27  sll_t_params_cos_gaussian_screwpinch
28 
29  use sll_m_linear_operator_block, only : &
30  sll_t_linear_operator_block
31 
32  use sll_m_particle_mass_3d_base, only: &
33  sll_c_particle_mass_3d_base
34 
35  use sll_m_linear_operator_particle_mass_3d_diag, only : &
36  sll_t_linear_operator_particle_mass_3d_diag
37 
38  use sll_m_linear_operator_particle_mass_cl_3d_diag, only : &
39  sll_t_linear_operator_particle_mass_cl_3d_diag
40 
41  use sll_m_linear_operator_particle_mass_3d_od, only : &
42  sll_t_linear_operator_particle_mass_3d_od
43 
44  use sll_m_linear_operator_particle_mass_cl_3d_od, only : &
45  sll_t_linear_operator_particle_mass_cl_3d_od
46 
47  use sll_m_linear_operator_schur_ev_3d, only : &
48  sll_t_linear_operator_schur_ev_3d
49 
50  use sll_m_linear_operator_schur_phiv_3d, only : &
51  sll_t_linear_operator_schur_phiv_3d
52 
53  use sll_m_linear_solver_cg, only : &
54  sll_t_linear_solver_cg
55 
56  use sll_m_mapping_3d, only: &
57  sll_t_mapping_3d
58 
59  use sll_m_maxwell_3d_base, only: &
60  sll_c_maxwell_3d_base
61 
62  use mpi, only: &
63  mpi_sum, &
64  mpi_max
65 
66  use sll_m_particle_group_base, only: &
67  sll_t_particle_array
68 
69  use sll_m_particle_mesh_coupling_base_3d, only: &
70  sll_c_particle_mesh_coupling_3d
71 
72  use sll_m_preconditioner_fft, only : &
73  sll_t_preconditioner_fft
74 
75  use sll_m_preconditioner_singular, only : &
76  sll_t_preconditioner_singular
77 
78  use sll_m_profile_functions, only: &
79  sll_t_profile_functions
80 
81 
82  implicit none
83 
84  public :: &
85  sll_t_time_propagator_pic_vm_3d3v_trafo_helper
86 
87  private
88  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
90  type :: sll_t_time_propagator_pic_vm_3d3v_trafo_helper
91  class(sll_c_maxwell_3d_base), pointer :: maxwell_solver
92  class(sll_c_particle_mesh_coupling_3d), pointer :: particle_mesh_coupling
93  class(sll_t_particle_array), pointer :: particle_group
94 
95  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_1
96  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_2
97  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_3
98  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_12
99  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_21
100  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_13
101  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_31
102  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_23
103  class(sll_c_particle_mass_3d_base), allocatable :: particle_mass_32
104 
105  type( sll_t_linear_operator_block ) :: particle_mass_op
106  type( sll_t_linear_operator_schur_ev_3d ) :: linear_op_schur_ev
107  type( sll_t_linear_operator_schur_phiv_3d ) :: linear_op_schur_phiv
108  type(sll_t_preconditioner_fft) :: preconditioner_fft
109  type(sll_t_preconditioner_singular) :: preconditioner1
110  type( sll_t_linear_solver_cg ) :: linear_solver_schur_ev
111  type( sll_t_linear_solver_cg ) :: linear_solver_schur_phiv
112 
113  type( sll_t_mapping_3d ), pointer :: map
114 
115  sll_int32 :: spline_degree(3)
116  sll_real64 :: lx(3)
117  sll_real64 :: x_min(3)
118  sll_real64 :: x_max(3)
119  sll_int32 :: n_total0
120  sll_int32 :: n_total1
121  sll_int32 :: nspan_1(3), nspan_2(3)
122 
123  sll_real64 :: betar(2) = 1._f64
124 
125  sll_real64, allocatable :: vec1(:)
126  sll_real64, allocatable :: vec2(:)
127 
128  sll_real64, pointer :: phi_dofs(:)
129  sll_real64, pointer :: efield_dofs(:)
130  sll_real64, pointer :: bfield_dofs(:)
131  sll_real64, allocatable :: j_dofs(:)
132  sll_real64, allocatable :: j_dofs_local(:)
133  sll_real64, allocatable :: particle_mass_1_local(:,:)
134  sll_real64, allocatable :: particle_mass_2_local(:,:)
135  sll_real64, allocatable :: particle_mass_3_local(:,:)
136  sll_real64, allocatable :: particle_mass_12_local(:,:)
137  sll_real64, allocatable :: particle_mass_21_local(:,:)
138  sll_real64, allocatable :: particle_mass_23_local(:,:)
139  sll_real64, allocatable :: particle_mass_32_local(:,:)
140  sll_real64, allocatable :: particle_mass_31_local(:,:)
141  sll_real64, allocatable :: particle_mass_13_local(:,:)
142 
143  sll_real64 :: solver_tolerance = 1d-12
144  sll_real64 :: iter_tolerance = 1d-10
145  sll_int32 :: max_iter = 10
146 
147  sll_int32 :: boundary_particles = 100
148  sll_int32 :: counter_left = 0
149  sll_int32 :: counter_right = 0
150  sll_real64, pointer :: rhob(:) => null()
151 
152  sll_real64 :: force_sign = 1._f64
153  logical :: adiabatic_electrons = .false.
154  logical :: jmean = .false.
155  logical :: boundary = .false.
156 
157  sll_int32 :: n_particles_max
158  sll_real64, allocatable :: xnew(:,:,:)
159  sll_real64, allocatable :: vnew(:,:,:)
160  sll_real64, allocatable :: efield_dofs_new(:)
161  sll_real64, allocatable :: efield_dofs_work(:)
162  sll_real64, allocatable :: phi_dofs_new(:)
163  sll_real64, allocatable :: phi_dofs_work(:)
164 
165  sll_int32 :: n_failed = 0
166  sll_int32 :: iter_counter = 0
167  sll_int32 :: niter(100000)
168 
169  ! For control variate
170  class(sll_t_control_variates), pointer :: control_variate => null()
171  sll_int32 :: i_weight = 1
172 
173  contains
174  procedure :: advect_x => advect_x_pic_vm_3d3v_trafo
175  procedure :: advect_vb => advect_vb_pic_vm_3d3v_trafo
176  procedure :: advect_eb => advect_eb_pic_vm_3d3v_trafo
177  procedure :: advect_e => advect_e_pic_vm_3d3v_trafo
178  procedure :: advect_ex => advect_ex_pic_vm_3d3v_trafo
179 
180  procedure :: init => initialize_pic_vm_3d3v_trafo
181  procedure :: init_from_file => initialize_file_pic_vm_3d3v_trafo
182  procedure :: free => delete_pic_vm_3d3v_trafo
183 
184  end type sll_t_time_propagator_pic_vm_3d3v_trafo_helper
185 
186 contains
187 
188 
189  !---------------------------------------------------------------------------!
192  subroutine advect_x_pic_vm_3d3v_trafo ( self, dt )
193  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
194  sll_real64, intent( in ) :: dt
195  !local variables
196  sll_int32 :: i_part, i, j, i_sp, iter
197  sll_real64 :: xi(3), xnew(3), xh(3), vi(3), vt(3), jmat(3,3), wall(3)
198  sll_real64 :: err
199  iter = 0
200  do i_sp = 1, self%particle_group%n_species
201  do i_part = 1, self%particle_group%group(i_sp)%n_particles
202  ! Read out particle position and velocity
203  xi = self%particle_group%group(i_sp)%get_x(i_part)
204  vi = self%particle_group%group(i_sp)%get_v(i_part)
205 
206  if( self%map%inverse)then
207  xh = self%map%get_x(xi)
208  xh = xh + dt * vi
209  xnew = self%map%get_xi(xh)
210  call sll_s_compute_particle_boundary( self, xi, xnew, vi )
211  else
212  !Predictor-Corrector with loop for corrector step
213  jmat=self%map%jacobian_matrix_inverse(xi)
214  !Transformation of the v coordinates
215  do j = 1, 3
216  vt(j) = jmat(j,1)*vi(1) + jmat(j,2)*vi(2) + jmat(j,3)*vi(3)
217  end do
218  !x^\star=\mod(x^n+dt*DF^{-1}(x^n)v^n,1)
219  xnew = xi + dt * vt
220  i = 0
221  err= maxval(abs(xi - xnew))
222  do while(i < self%max_iter .and. err > self%iter_tolerance)
223  call sll_s_compute_particle_boundary_simple( self%boundary_particles, self%counter_left, self%counter_right, xi, xnew )
224  jmat=self%map%jacobian_matrix_inverse( xnew )
225  do j = 1, 3
226  xh(j) = xi(j) + 0.5_f64 * dt * (vt(j) + jmat(j,1)*vi(1) + jmat(j,2)*vi(2)+ jmat(j,3)*vi(3))
227  end do
228  err = maxval(abs(xh - xnew))
229  xnew = xh
230  i = i + 1
231  end do
232  iter = iter + i - 1
233  call sll_s_compute_particle_boundary( self, xi, xnew, vi )
234  end if
235  call self%particle_group%group(i_sp)%set_x ( i_part, xnew )
236  call self%particle_group%group(i_sp)%set_v ( i_part, vi )
237  ! Update particle weights
238  if (self%particle_group%group(i_sp)%n_weights == 3 ) then
239  wall = self%particle_group%group(i_sp)%get_weights(i_part)
240  wall(3) = self%control_variate%cv(i_sp)%update_df_weight( xnew, vi, 0.0_f64, wall(1), wall(2) )
241  call self%particle_group%group(i_sp)%set_weights( i_part, wall )
242  end if
243  end do
244  end do
245 
246  end subroutine advect_x_pic_vm_3d3v_trafo
247 
248 
250  subroutine sll_s_compute_particle_boundary( self, xold, xnew, vi )
251  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
252  sll_real64, intent( inout ) :: xold(3)
253  sll_real64, intent( inout ) :: xnew(3)
254  sll_real64, intent( inout ) :: vi(3)
255  !local variables
256  sll_real64 :: xmid(3), xt(3), xbar, dx
257  sll_real64 :: jmatrix(3,3)
258 
259  if(xnew(1) < 0._f64 .or. xnew(1) > 1._f64 )then
260  if(xnew(1) < 0._f64 )then
261  xbar = 0._f64
262  self%counter_left = self%counter_left+1
263  else if(xnew(1) > 1._f64)then
264  xbar = 1._f64
265  self%counter_right = self%counter_right+1
266  end if
267  dx = (xbar- xold(1))/(xnew(1)-xold(1))
268  xmid = xold + dx * (xnew-xold)
269  xmid(1) = xbar
270 
271  select case(self%boundary_particles)
272  case(sll_p_boundary_particles_singular)
273  if(xnew(1) < 0._f64 )then
274  xnew(1) = -xnew(1)
275  xnew(2) = xnew(2) + 0.5_f64
276  else if(xnew(1) > 1._f64 )then
277  jmatrix = self%map%jacobian_matrix_inverse_transposed(xmid)
278  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
279  xnew(1) = 2._f64 - xnew(1)
280  xnew(2) = 1._f64 - xnew(2)
281  end if
282  case(sll_p_boundary_particles_reflection)
283  xt = xmid
284  xt(2:3) = modulo(xt(2:3),1._f64)
285  jmatrix = self%map%jacobian_matrix_inverse_transposed(xt)
286  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
287  xnew(1) = 2._f64*xbar-xnew(1)
288  case(sll_p_boundary_particles_absorption)
289  case( sll_p_boundary_particles_periodic)
290  xnew(1) = modulo(xnew(1), 1._f64)
291  case default
292  xnew(1) = modulo(xnew(1), 1._f64)
293  end select
294  else if(xnew(1) < -1._f64 .or. xnew(1) > 2._f64 ) then
295  print*, xnew
296  sll_error("particle boundary", "particle out of bound")
297  end if
298  xnew(2:3) = modulo(xnew(2:3), 1._f64)
299 
300  end subroutine sll_s_compute_particle_boundary
301 
302 
303  !---------------------------------------------------------------------------!
306  subroutine advect_vb_pic_vm_3d3v_trafo ( self, dt )
307  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
308  sll_real64, intent( in ) :: dt
309  !local variables
310  sll_int32 :: i_part, i_sp, j
311  sll_real64 :: qmdt
312  sll_real64 :: vi(3), xi(3), wall(3)
313  sll_real64 :: bfield(3), jmatrix(3,3), rhs(3)
314  sll_real64 :: vt(3), c(3)
315 
316  do i_sp = 1, self%particle_group%n_species
317  qmdt = self%particle_group%group(i_sp)%species%q_over_m()*dt*0.5_f64;
318  do i_part = 1, self%particle_group%group(i_sp)%n_particles
319  vi = self%particle_group%group(i_sp)%get_v(i_part)
320  xi = self%particle_group%group(i_sp)%get_x(i_part)
321 
322  call self%particle_mesh_coupling%evaluate &
323  (xi, [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)-1], self%bfield_dofs(1:self%n_total0), bfield(1))
324  call self%particle_mesh_coupling%evaluate &
325  (xi, [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)-1],self%bfield_dofs(self%n_total0+1:self%n_total0+self%n_total1), bfield(2))
326  call self%particle_mesh_coupling%evaluate &
327  (xi, [self%spline_degree(1)-1, self%spline_degree(2)-1, self%spline_degree(3)],self%bfield_dofs(self%n_total0+self%n_total1+1:self%n_total0+2*self%n_total1), bfield(3))
328 
329  jmatrix=self%map%jacobian_matrix_inverse(xi)
330 
331  !VT = DF^{-1} * vi
332  do j=1,3
333  vt(j)=jmatrix(j,1)*vi(1)+jmatrix(j,2)*vi(2)+jmatrix(j,3)*vi(3)
334  end do
335  !c = VT x bfield
336  c(1)=bfield(3)*vt(2)-bfield(2)*vt(3)
337  c(2)=bfield(1)*vt(3)-bfield(3)*vt(1)
338  c(3)=bfield(2)*vt(1)-bfield(1)*vt(2)
339  !rhs = vi + sign * DF^{-T} * c
340  do j=1,3
341  rhs(j)= vi(j) + qmdt*(jmatrix(1,j)*c(1)+jmatrix(2,j)*c(2)+jmatrix(3,j)*c(3))
342  end do
343 
344  call sll_s_compute_matrix_inverse(rhs, vi, bfield, jmatrix, qmdt)
345 
346  call self%particle_group%group(i_sp)%set_v( i_part, vi )
347  ! Update particle weights
348  if (self%particle_group%group(i_sp)%n_weights == 3 ) then
349  wall = self%particle_group%group(i_sp)%get_weights(i_part)
350  wall(3) = self%control_variate%cv(i_sp)%update_df_weight( xi, vi, 0.0_f64, wall(1), wall(2) )
351  call self%particle_group%group(i_sp)%set_weights( i_part, wall )
352  end if
353  end do
354  end do
355 
356  end subroutine advect_vb_pic_vm_3d3v_trafo
357 
358 
359  !---------------------------------------------------------------------------!
364  subroutine advect_eb_pic_vm_3d3v_trafo ( self, dt )
365  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
366  sll_real64, intent( in ) :: dt
367 
368  call self%maxwell_solver%compute_curl_part( dt, self%efield_dofs, self%bfield_dofs, self%betar(1) )
369 
370  end subroutine advect_eb_pic_vm_3d3v_trafo
371 
372 
373  !---------------------------------------------------------------------------!
378  subroutine advect_e_pic_vm_3d3v_trafo( self, dt )
379  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
380  sll_real64, intent( in ) :: dt
381  ! local variables
382  sll_int32 :: i_part, j, i_sp
383  sll_real64 :: vi(3), vt(3), xi(3), wi(1), wall(3)
384  sll_real64 :: metric(3,3), jmat(3,3)
385  sll_real64 :: factor, qoverm
386  sll_real64 :: efield(3), ephys(3)
387  sll_real64 :: rhs(self%n_total1+2*self%n_total0)
388 
389  ! Set to zero
390  self%j_dofs_local = 0.0_f64
391  self%j_dofs = 0.0_f64
392  self%particle_mass_1_local = 0.0_f64
393  self%particle_mass_2_local = 0.0_f64
394  self%particle_mass_3_local = 0.0_f64
395  if(self%map%flag2d)then
396  self%particle_mass_12_local = 0.0_f64
397  self%particle_mass_21_local = 0.0_f64
398  if(self%map%flag3d)then
399  self%particle_mass_23_local = 0.0_f64
400  self%particle_mass_32_local = 0.0_f64
401  self%particle_mass_31_local = 0.0_f64
402  self%particle_mass_13_local = 0.0_f64
403  end if
404  end if
405 
406  self%particle_mass_1%particle_mass = 0.0_f64
407  self%particle_mass_2%particle_mass = 0.0_f64
408  self%particle_mass_3%particle_mass = 0.0_f64
409  if(self%map%flag2d)then
410  self%particle_mass_12%particle_mass = 0.0_f64
411  self%particle_mass_21%particle_mass = 0.0_f64
412  if(self%map%flag3d)then
413  self%particle_mass_23%particle_mass = 0.0_f64
414  self%particle_mass_32%particle_mass = 0.0_f64
415  self%particle_mass_13%particle_mass = 0.0_f64
416  self%particle_mass_31%particle_mass = 0.0_f64
417  end if
418  end if
419 
420  ! First particle loop
421  do i_sp = 1, self%particle_group%n_species
422  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
423  if( self%adiabatic_electrons ) then
424  factor = dt**2 * 0.25_f64* qoverm
425  else
426  factor = dt**2 * self%betar(2) * 0.25_f64* qoverm
427  end if
428  do i_part = 1, self%particle_group%group(i_sp)%n_particles
429  vi = self%particle_group%group(i_sp)%get_v(i_part)
430  xi = self%particle_group%group(i_sp)%get_x(i_part)
431 
432  ! Get charge for accumulation of j
433  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
434 
435  metric= self%map%metric_inverse( xi )
436  jmat=self%map%jacobian_matrix_inverse( xi )
437  do j=1,3
438  vt(j)=vi(1)*jmat(j,1)+vi(2)*jmat(j,2)+vi(3)*jmat(j,3)
439  end do
440 
441  ! Accumulate j
442  call self%particle_mesh_coupling%add_charge( xi, wi(1)*vt(1), &
443  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
444  self%j_dofs_local(1:self%n_total1) )
445  call self%particle_mesh_coupling%add_charge( xi, wi(1)*vt(2), &
446  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
447  self%j_dofs_local(1+self%n_total1:self%n_total1+self%n_total0) )
448  call self%particle_mesh_coupling%add_charge( xi, wi(1)*vt(3), &
449  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
450  self%j_dofs_local(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0) )
451 
452 
453  ! Accumulate the particle mass matrix
454  call self%particle_mesh_coupling%add_particle_mass( xi, wi(1)*metric(1,1) * factor, &
455  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
456  self%particle_mass_1_local )
457  call self%particle_mesh_coupling%add_particle_mass( xi, wi(1)*metric(2,2) * factor, &
458  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
459  self%particle_mass_2_local )
460  call self%particle_mesh_coupling%add_particle_mass( xi, wi(1)*metric(3,3) * factor, &
461  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
462  self%particle_mass_3_local )
463  if(self%map%flag2d)then
464  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(1,2) * factor, &
465  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
466  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)],&
467  self%particle_mass_12_local )
468  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(2,1) * factor, &
469  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)],&
470  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
471  self%particle_mass_21_local )
472  if(self%map%flag3d)then
473  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(2,3) * factor, &
474  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
475  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
476  self%particle_mass_23_local )
477  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(3,2) * factor, &
478  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
479  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
480  self%particle_mass_32_local )
481  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(3,1) * factor, &
482  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
483  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
484  self%particle_mass_31_local )
485  call self%particle_mesh_coupling%add_particle_mass_od( xi, wi(1)*metric(1,3) * factor, &
486  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
487  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
488  self%particle_mass_13_local )
489  end if
490  end if
491 
492 
493  ! Evaulate E at particle position and propagate v a half step
494  call self%particle_mesh_coupling%evaluate &
495  (xi, [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
496  self%efield_dofs(1:self%n_total1), efield(1))
497  call self%particle_mesh_coupling%evaluate &
498  (xi, [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
499  self%efield_dofs(1+self%n_total1:self%n_total1+self%n_total0), efield(2))
500  call self%particle_mesh_coupling%evaluate &
501  (xi, [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
502  self%efield_dofs(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0), efield(3))
503 
504  ! physical efield, jmatrix is not transposed
505  do j=1, 3
506  ephys(j) = jmat(1,j)* efield(1)+jmat(2,j)* efield(2)+jmat(3,j)* efield(3)
507  end do
508 
509  ! velocity update
510  vi = vi + dt* 0.5_f64* qoverm * ephys
511 
512  call self%particle_group%group(i_sp)%set_v( i_part, vi )
513  ! Update particle weights
514  if (self%particle_group%group(i_sp)%n_weights == 3 ) then
515  wall = self%particle_group%group(i_sp)%get_weights(i_part)
516  wall(3) = self%control_variate%cv(i_sp)%update_df_weight( xi, vi, 0.0_f64, wall(1), wall(2) )
517  call self%particle_group%group(i_sp)%set_weights( i_part, wall )
518  end if
519  end do
520  end do
521 
522  ! MPI to sum up contributions from each processor
523  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local, &
524  self%n_total1+self%n_total0*2, mpi_sum, self%j_dofs )
525 
526  if ( self%jmean ) then
527  self%j_dofs(1:self%n_total1) = self%j_dofs(1:self%n_total1) - sum(self%j_dofs(1:self%n_total1))/real(self%n_total1, f64)
528  self%j_dofs(1+self%n_total1:self%n_total1+self%n_total0) = self%j_dofs(1+self%n_total1:self%n_total1+self%n_total0) - sum(self%j_dofs(1+self%n_total1:self%n_total1+self%n_total0))/real(self%n_total0, f64)
529  self%j_dofs(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0) = self%j_dofs(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0) - sum(self%j_dofs(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0))/real(self%n_total0, f64)
530  end if
531 
532  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_1_local, &
533  self%n_total1*self%nspan_1(1) , mpi_sum, self%particle_mass_1%particle_mass)
534  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_2_local, &
535  self%n_total0*self%nspan_1(2) , mpi_sum, self%particle_mass_2%particle_mass)
536  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_3_local, &
537  self%n_total0*self%nspan_1(3) , mpi_sum, self%particle_mass_3%particle_mass)
538 
539  if(self%map%flag2d)then
540  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_12_local, &
541  self%n_total1*self%nspan_2(1) , mpi_sum, self%particle_mass_12%particle_mass)
542  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_21_local, &
543  self%n_total0*self%nspan_2(1) , mpi_sum, self%particle_mass_21%particle_mass)
544  if(self%map%flag3d)then
545  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_23_local, &
546  self%n_total0*self%nspan_2(2) , mpi_sum, self%particle_mass_23%particle_mass)
547  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_31_local, &
548  self%n_total0*self%nspan_2(3) , mpi_sum, self%particle_mass_31%particle_mass)
549  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_32_local, &
550  self%n_total0*self%nspan_2(2) , mpi_sum, self%particle_mass_32%particle_mass)
551  call sll_o_collective_allreduce( sll_v_world_collective, self%particle_mass_13_local, &
552  self%n_total1*self%nspan_2(3) , mpi_sum, self%particle_mass_13%particle_mass)
553  end if
554  end if
555 
556  call self%particle_mass_op%set( 1, 1, self%particle_mass_1 )
557  call self%particle_mass_op%set( 2, 2, self%particle_mass_2 )
558  call self%particle_mass_op%set( 3, 3, self%particle_mass_3 )
559  if(self%map%flag2d)then
560  call self%particle_mass_op%set( 1, 2, self%particle_mass_12 )
561  call self%particle_mass_op%set( 2, 1, self%particle_mass_21 )
562  if(self%map%flag3d)then
563  call self%particle_mass_op%set( 2, 3, self%particle_mass_23 )
564  call self%particle_mass_op%set( 3, 2, self%particle_mass_32 )
565  call self%particle_mass_op%set( 3, 1, self%particle_mass_31 )
566  call self%particle_mass_op%set( 1, 3, self%particle_mass_13 )
567  end if
568  end if
569 
570  if( self%adiabatic_electrons ) then
571  ! Compute rhs
572  self%linear_op_schur_phiv%sign = -1._f64
573  call self%linear_op_schur_phiv%dot(self%phi_dofs, rhs(1: self%n_total0 ))
574  call self%maxwell_solver%multiply_gt( self%j_dofs, rhs(self%n_total0+1: 2*self%n_total0 ) )
575  rhs(1: self%n_total0 ) = rhs(1: self%n_total0 ) + dt * rhs(self%n_total0+1: 2*self%n_total0 )
576 
577  ! Solve the Schur complement
578  self%linear_op_schur_phiv%sign = 1._f64
579  call self%linear_solver_schur_phiv%set_guess(self%phi_dofs)
580  call self%linear_solver_schur_phiv%solve( rhs(1: self%n_total0 ), self%phi_dofs )
581  if( self%boundary ) then !perfect boundary conditions
582  do j = 1, self%maxwell_solver%n_dofs(2)*self%maxwell_solver%n_dofs(3)
583  self%phi_dofs(1+(j-1)*self%maxwell_solver%n_dofs(1)) = 0._f64
584  self%phi_dofs(j*self%maxwell_solver%n_dofs(1)) = 0._f64
585  end do
586  end if
587  call self%maxwell_solver%multiply_g( self%phi_dofs, self%efield_dofs )
588  self%efield_dofs = -self%efield_dofs
589  else
590  ! Compute rhs
591  self%linear_op_schur_ev%sign = -self%force_sign
592  call self%linear_op_schur_ev%dot(self%efield_dofs, rhs)
593  rhs = rhs - self%force_sign *dt * self%betar(2) * self%j_dofs
594 
595  ! Solve the Schur complement
596  self%linear_op_schur_ev%sign = self%force_sign
597  call self%linear_solver_schur_ev%set_guess(self%efield_dofs)
598  call self%linear_solver_schur_ev%solve(rhs, self%efield_dofs)
599  end if
600 
601  if( self%boundary ) then !perfect boundary conditions
602  do j = 1, self%maxwell_solver%n_dofs(2)*self%maxwell_solver%n_dofs(3)
603  self%efield_dofs(self%n_total1+1+(j-1)*self%maxwell_solver%n_dofs(1)) = 0._f64
604  self%efield_dofs(self%n_total1+j*self%maxwell_solver%n_dofs(1)) = 0._f64
605  self%efield_dofs(self%n_total1+self%n_total0+1+(j-1)*self%maxwell_solver%n_dofs(1)) = 0._f64
606  self%efield_dofs(self%n_total1+self%n_total0+j*self%maxwell_solver%n_dofs(1)) = 0._f64
607  end do
608  end if
609  ! Second particle loop (second half step of particle propagation)
610  do i_sp = 1, self%particle_group%n_species
611  qoverm = self%particle_group%group(i_sp)%species%q_over_m();
612  do i_part = 1, self%particle_group%group(i_sp)%n_particles
613  vi = self%particle_group%group(i_sp)%get_v(i_part)
614  xi = self%particle_group%group(i_sp)%get_x(i_part)
615 
616  ! Evaulate E at particle position and propagate v a half step
617  call self%particle_mesh_coupling%evaluate &
618  (xi, [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
619  self%efield_dofs(1:self%n_total1), efield(1))
620  call self%particle_mesh_coupling%evaluate &
621  (xi, [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
622  self%efield_dofs(1+self%n_total1:self%n_total1+self%n_total0), efield(2))
623  call self%particle_mesh_coupling%evaluate &
624  (xi, [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
625  self%efield_dofs(1+self%n_total1+self%n_total0:self%n_total1+2*self%n_total0), efield(3))
626 
627  jmat = self%map%jacobian_matrix_inverse_transposed(xi)
628  do j=1, 3
629  ephys(j) = jmat(j,1)* efield(1)+jmat(j,2)* efield(2)+jmat(j,3)* efield(3)
630  end do
631 
632  ! velocity update
633  vi = vi + dt* 0.5_f64* qoverm* ephys
634 
635  call self%particle_group%group(i_sp)%set_v( i_part, vi )
636 
637  ! Update particle weights
638  if (self%particle_group%group(i_sp)%n_weights == 3 ) then
639  wall = self%particle_group%group(i_sp)%get_weights(i_part)
640  wall(3) = self%control_variate%cv(i_sp)%update_df_weight( xi, vi, 0.0_f64, wall(1), wall(2) )
641  call self%particle_group%group(i_sp)%set_weights( i_part, wall )
642  end if
643  end do
644  end do
645 
646  end subroutine advect_e_pic_vm_3d3v_trafo
647 
648 
649  !---------------------------------------------------------------------------!
654  subroutine advect_ex_pic_vm_3d3v_trafo( self, dt )
655  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
656  sll_real64, intent( in ) :: dt
657  ! local variables
658  sll_int32 :: i_part, i_sp, j
659  sll_real64 :: vi(3), vh(3), xi(3), xs(3), wi(1), mi(1), xnew(3), vbar(3), vnew(3)
660  sll_real64 :: sign, wall(3)
661  sll_real64 :: jmat(3,3), jmatrix(3,3)
662  sll_int32 :: niter
663  sll_real64 :: residual(8), residual_local(8)
664  sll_real64 :: xmid(3), xbar, dx
665  sll_real64 :: efield(3)
666 
667  self%efield_dofs_new = self%efield_dofs
668  self%phi_dofs_new = self%phi_dofs
669  do i_sp=1,self%particle_group%n_species
670  do i_part = 1,self%particle_group%group(i_sp)%n_particles
671  vi = self%particle_group%group(i_sp)%get_v(i_part)
672  xi = self%particle_group%group(i_sp)%get_x(i_part)
673 
674  ! Get charge for accumulation of j
675  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
676 
677  if( self%map%inverse) then
678  xs = self%map%get_x(xi)
679  xs = xs + dt * vbar
680  xnew = self%map%get_xi(xs)
681  else
682  jmat = self%map%jacobian_matrix_inverse_transposed( xi )
683  do j = 1, 3
684  vh(j) = jmat(1,j)*vi(1) + jmat(2,j)*vi(2) + jmat(3,j)*vi(3)
685  end do
686  xnew = xi + dt * vh
687  end if
688 
689  if(xnew(1) < 0._f64 .or. xnew(1) > 1._f64 )then
690  if(xnew(1) < 0._f64 )then
691  xbar = 0._f64
692  else if(xnew(1) > 1._f64)then
693  xbar = 1._f64
694  end if
695  dx = (xbar- xi(1))/(xnew(1)-xi(1))
696  xmid = xi + dx * (xnew-xi)
697  xmid(1) = xbar
698  xnew = xmid
699  else
700  dx = 1._f64
701  end if
702  vh = (xnew-xi) * wi(1)
703  call self%particle_mesh_coupling%add_current_evaluate( xi, xnew, vh, self%efield_dofs, &
704  self%j_dofs_local, efield )
705  do j = 1, 3
706  vnew(j) = vi(j) + dx*sign*(jmat(j,1)*efield(1) + jmat(j,2)*efield(2) + jmat(j,3)*efield(3))
707  end do
708  xnew(2:3) = modulo(xnew(2:3), 1._f64)
709 
710  self%vnew(i_sp,:, i_part) = vnew
711  self%xnew(i_sp,:, i_part) = xnew
712  end do
713  end do
714  self%j_dofs = 0.0_f64
715  ! MPI to sum up contributions from each processor
716  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local, &
717  self%n_total1+2*self%n_total0, mpi_sum, self%j_dofs )
718 
719  if( self%adiabatic_electrons) then
720  call self%maxwell_solver%compute_phi_from_j( self%j_dofs, self%phi_dofs_new, self%efield_dofs_new )
721  else
722  call self%maxwell_solver%compute_E_from_j( self%betar(2)*self%j_dofs, self%efield_dofs_new )
723  end if
724 
725  niter = 0
726  residual = 1.0_f64
727 
728  do while ( (maxval(residual(1:7)) > self%iter_tolerance) .and. niter < self%max_iter )
729  niter = niter+1
730  self%efield_dofs_work = 0.5_f64*( self%efield_dofs + self%efield_dofs_new )
731  self%j_dofs_local = 0.0_f64
732  residual_local = 0._f64
733  ! Particle loop
734  do i_sp = 1, self%particle_group%n_species
735  sign = dt*self%particle_group%group(i_sp)%species%q_over_m();
736  do i_part = 1,self%particle_group%group(i_sp)%n_particles
737  vi = self%particle_group%group(i_sp)%get_v(i_part)
738  xi = self%particle_group%group(i_sp)%get_x(i_part)
739 
740  ! Get charge for accumulation of j
741  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
742  mi = self%particle_group%group(i_sp)%get_mass(i_part, self%i_weight)
743  vbar = 0.5_f64 * (self%vnew(i_sp,:, i_part)+vi)
744  xnew = self%xnew(i_sp,:, i_part)
745 
746  if( self%map%inverse) then
747  xs = self%map%get_x(xi)
748  xs = xs + dt * vbar
749  xnew = self%map%get_xi(xs)
750  else
751  jmat = self%map%jacobian_matrix_inverse_transposed( xi )
752  jmatrix=self%map%jacobian_matrix_inverse_transposed( xnew )
753  do j = 1, 3
754  vh(j) = 0.5_f64 * ((jmatrix(1,j)+jmat(1,j))*vbar(1) + (jmatrix(2,j)+jmat(2,j))*vbar(2) + (jmatrix(3,j)+ jmat(3,j))*vbar(3))
755  end do
756  xnew = xi + dt * vh
757  end if
758 
759  call compute_particle_boundary_current_evaluate_iter( self, xi, xnew, vi, wi, sign )
760 
761  residual_local(1) = residual_local(1) + (xnew(1)-self%xnew(i_sp,1,i_part))**2*abs(wi(1))
762  residual_local(2) = residual_local(2) + (xnew(2)-self%xnew(i_sp,2,i_part))**2*abs(wi(1))
763  residual_local(3) = residual_local(3) + (xnew(3)-self%xnew(i_sp,3,i_part))**2*abs(wi(1))
764 
765  residual_local(4) = residual_local(4) + (vi(1)-self%vnew(i_sp,1,i_part))**2*abs(wi(1))
766  residual_local(5) = residual_local(5) + (vi(2)-self%vnew(i_sp,2,i_part))**2*abs(wi(1))
767  residual_local(6) = residual_local(6) + (vi(3)-self%vnew(i_sp,3,i_part))**2*abs(wi(1))
768 
769  residual_local(7) = residual_local(7) + mi(1)*0.5_f64*( vi(1)**2 + vi(2)**2 + vi(3)**2-self%vnew(i_sp,1,i_part)**2 -self%vnew(i_sp,2,i_part)**2 -self%vnew(i_sp,3,i_part)**2 )
770 
771  self%xnew(i_sp, :, i_part) = xnew
772  self%vnew(i_sp, :, i_part) = vi
773  end do
774  end do
775 
776  self%j_dofs = 0.0_f64
777  ! MPI to sum up contributions from each processor
778  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local, &
779  self%n_total1+2*self%n_total0, mpi_sum, self%j_dofs )
780 
781  if( self%adiabatic_electrons) then
782  self%phi_dofs_work = self%phi_dofs
783  call self%maxwell_solver%compute_phi_from_j( self%j_dofs, self%phi_dofs_work, self%efield_dofs_work )
784 
785  ! Compute residuals based on phi
786  residual_local(8) = residual_local(8) + 0.5_f64*(self%maxwell_solver%inner_product( self%phi_dofs_work, self%phi_dofs_work, 0 ) - self%maxwell_solver%inner_product( self%phi_dofs_new, self%phi_dofs_new, 0 ) )
787 
788  residual_local(7) = (sum((self%phi_dofs_work-self%phi_dofs_new)**2))*product(self%particle_mesh_coupling%delta_x)
789  else
790  self%efield_dofs_work = self%efield_dofs
791  call self%maxwell_solver%compute_E_from_j( self%betar(2)*self%j_dofs, self%efield_dofs_work )
792 
793  ! Compute residuals based on e
794  residual_local(8) = residual_local(8) + 0.5_f64*(self%maxwell_solver%inner_product( self%efield_dofs_work, self%efield_dofs_work, 1 ) - self%maxwell_solver%inner_product( self%efield_dofs_new, self%efield_dofs_new, 1 ) )
795 
796  residual_local(7) = (sum((self%efield_dofs_work-self%efield_dofs_new)**2))*product(self%particle_mesh_coupling%delta_x)
797  end if
798 
799  call sll_o_collective_allreduce( sll_v_world_collective, residual_local, 8, mpi_max, residual )
800  residual(1:7) = sqrt(residual(1:7))
801  residual(8) = abs(residual(8))
802 
803  self%phi_dofs_new = self%phi_dofs_work
804  self%efield_dofs_new = self%efield_dofs_work
805  end do
806 
807  if ( maxval(residual(1:7)) > self%iter_tolerance ) then
808  print*, 'Warning: Iteration no.', self%iter_counter+1 ,'did not converge.', residual, niter
809  self%n_failed = self%n_failed+1
810  end if
811 
812 
813  self%efield_dofs_work = 0.5_f64*( self%efield_dofs + self%efield_dofs_new )
814  self%j_dofs_local = 0.0_f64
815 
816  ! Particle loop
817  do i_sp=1,self%particle_group%n_species
818  sign = dt*self%particle_group%group(i_sp)%species%q_over_m();
819  do i_part = 1,self%particle_group%group(i_sp)%n_particles
820  vi = self%particle_group%group(i_sp)%get_v(i_part)
821  xi = self%particle_group%group(i_sp)%get_x(i_part)
822 
823  ! Get charge for accumulation of j
824  wi = self%particle_group%group(i_sp)%get_charge(i_part, self%i_weight)
825  vbar = 0.5_f64 * (self%vnew(i_sp,:, i_part)+vi)
826  xnew = self%xnew(i_sp,:, i_part)
827 
828  if( self%map%inverse) then
829  xs = self%map%get_x(xi)
830  xs = xs + dt * vbar
831  xnew = self%map%get_xi(xs)
832  else
833  jmat = self%map%jacobian_matrix_inverse_transposed( xi )
834  jmatrix=self%map%jacobian_matrix_inverse_transposed( xnew )
835  do j = 1, 3
836  vh(j) = 0.5_f64 * ((jmatrix(1,j)+jmat(1,j))*vbar(1) + (jmatrix(2,j)+jmat(2,j))*vbar(2) + (jmatrix(3,j)+ jmat(3,j))*vbar(3))
837  end do
838  xnew = xi + dt * vh
839  end if
840  call compute_particle_boundary_current_evaluate( self, xi, xnew, vi, wi, sign )
841 
842  call self%particle_group%group(i_sp)%set_v( i_part, vi )
843  call self%particle_group%group(i_sp)%set_x( i_part, xnew )
844  ! Update particle weights
845  if (self%particle_group%group(i_sp)%n_weights == 3 ) then
846  wall = self%particle_group%group(i_sp)%get_weights(i_part)
847  wall(3) = self%control_variate%cv(i_sp)%update_df_weight( xnew, vi, 0.0_f64, wall(1), wall(2) )
848  call self%particle_group%group(i_sp)%set_weights( i_part, wall )
849  end if
850  end do
851  end do
852 
853  self%j_dofs = 0.0_f64
854  ! MPI to sum up contributions from each processor
855  call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local, &
856  self%n_total1+2*self%n_total0, mpi_sum, self%j_dofs )
857 
858  if( self%adiabatic_electrons) then
859  call self%maxwell_solver%compute_phi_from_j( self%j_dofs, self%phi_dofs, self%efield_dofs )
860  else
861  call self%maxwell_solver%compute_E_from_j( self%betar(2)*self%j_dofs, self%efield_dofs )
862  end if
863 
864  self%iter_counter = self%iter_counter + 1
865  self%niter(self%iter_counter) = niter
866 
867  end subroutine advect_ex_pic_vm_3d3v_trafo
868 
869 
871  subroutine compute_particle_boundary_current_evaluate_iter( self, xi, xnew, vi, wi, sign )
872  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
873  sll_real64, intent( in ) :: xi(3)
874  sll_real64, intent( inout ) :: xnew(3)
875  sll_real64, intent( inout ) :: vi(3)
876  sll_real64, intent( in ) :: wi(1)
877  sll_real64, intent( in ) :: sign
878  !local variables
879  sll_real64 :: xmid(3), xt(3), vh(3), xbar, dx
880  sll_real64 :: jmatrix(3,3), jmat(3,3), efield(3)
881  sll_int32 :: j
882 
883  jmat = self%map%jacobian_matrix_inverse_transposed( xi )
884  if(xnew(1) < -1._f64 .or. xnew(1) > 2._f64)then
885  print*, xnew
886  sll_error("particle boundary", "particle out of bound")
887  else if(xnew(1) < 0._f64 .or. xnew(1) > 1._f64 )then
888  if(xnew(1) < 0._f64 )then
889  xbar = 0._f64
890  else if(xnew(1) > 1._f64)then
891  xbar = 1._f64
892  end if
893  dx = (xbar- xi(1))/(xnew(1)-xi(1))
894  xmid = xi + dx * (xnew-xi)
895  xmid(1) = xbar
896 
897  vh = (xmid-xi) * wi(1)
898  call self%particle_mesh_coupling%add_current_evaluate( xi, xmid, vh, self%efield_dofs_work, &
899  self%j_dofs_local, efield )
900 
901  xt = xmid
902  xt(2:3) = modulo(xt(2:3), 1._f64)
903  jmatrix = self%map%jacobian_matrix_inverse_transposed(xt)
904  do j = 1, 3
905  vi(j) = vi(j) + dx*sign*0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
906  end do
907  select case(self%boundary_particles)
908  case(sll_p_boundary_particles_singular)
909  call self%particle_mesh_coupling%add_charge(xmid, wi(1), self%spline_degree, self%rhob)
910  xmid(2) = xbar + (1._f64-2._f64*xbar)*xmid(2) + 0.5_f64-0.5_f64*xbar
911  xt(2:3) = modulo(xt(2:3), 1._f64)
912  jmatrix = self%map%jacobian_matrix_inverse_transposed(xt)
913  call self%particle_mesh_coupling%add_charge(xmid, -wi(1), self%spline_degree, self%rhob)
914  xnew(1) = 2._f64*xbar-xnew(1)
915  xnew(2) = xbar + (1._f64-2._f64*xbar)*xnew(2) + 0.5_f64-0.5_f64*xbar
916  if(xnew(1) > 1._f64 )then
917  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
918  end if
919  case(sll_p_boundary_particles_reflection)
920  xnew(1) = 2._f64*xbar-xnew(1)
921  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
922  case(sll_p_boundary_particles_absorption)
923  call self%particle_mesh_coupling%add_charge(xmid, wi(1), self%spline_degree, self%rhob)
924  case( sll_p_boundary_particles_periodic)
925  xnew(1) = modulo(xnew(1), 1._f64)
926  xmid(1) = 1._f64-xbar
927  case default
928  xnew(1) = modulo(xnew(1), 1._f64)
929  xmid(1) = 1._f64-xbar
930  end select
931  vh = (xnew - xmid)*wi(1)
932  call self%particle_mesh_coupling%add_current_evaluate( xmid, xnew, vh, self%efield_dofs_work, &
933  self%j_dofs_local, efield )
934  xnew(2:3) = modulo(xnew(2:3), 1._f64)
935  jmat = self%map%jacobian_matrix_inverse_transposed(xnew)
936  do j = 1, 3
937  vi(j) = vi(j) + (1._f64-dx) * sign *0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
938  end do
939  if(self%boundary_particles == sll_p_boundary_particles_reflection) then
940  do j = 1, 3 !calculation of velocity in logical coordinates
941  vh(j) = jmat(1,j)*vi(1) + jmat(2,j)*vi(2) + jmat(3,j)*vi(3)
942  end do
943  vh(1) = - vh(1)
944  jmat = self%map%jacobian_matrix(xnew)
945  do j = 1, 3 !calculation of velocity in physical coordinates
946  vi(j) = jmat(j,1)*vh(1) + jmat(j,2)*vh(2) + jmat(j,3)*vh(3)
947  end do
948  end if
949  else
950  vh = (xnew-xi) * wi(1)
951  call self%particle_mesh_coupling%add_current_evaluate( xi, xnew, vh, self%efield_dofs_work, &
952  self%j_dofs_local, efield )
953  xnew(2:3) = modulo(xnew(2:3), 1._f64)
954  jmatrix = self%map%jacobian_matrix_inverse_transposed(xnew)
955  do j = 1, 3
956  vi(j) = vi(j) + sign *0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
957  end do
958  end if
959 
960  end subroutine compute_particle_boundary_current_evaluate_iter
961 
962 
964  subroutine compute_particle_boundary_current_evaluate( self, xi, xnew, vi, wi, sign )
965  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
966  sll_real64, intent( in ) :: xi(3)
967  sll_real64, intent( inout ) :: xnew(3)
968  sll_real64, intent( inout ) :: vi(3)
969  sll_real64, intent( in ) :: wi(1)
970  sll_real64, intent( in ) :: sign
971  !local variables
972  sll_real64 :: xmid(3), xt(3), vh(3), xbar, dx
973  sll_real64 :: jmatrix(3,3), jmat(3,3), efield(3)
974  sll_int32 :: j
975 
976  jmat = self%map%jacobian_matrix_inverse_transposed( xi )
977 
978  if(xnew(1) < -1._f64 .or. xnew(1) > 2._f64)then
979  print*, xnew
980  sll_error("particle boundary", "particle out of bound")
981  else if(xnew(1) < 0._f64 .or. xnew(1) > 1._f64 )then
982  if(xnew(1) < 0._f64 )then
983  xbar = 0._f64
984  self%counter_left = self%counter_left+1
985  else if(xnew(1) > 1._f64)then
986  xbar = 1._f64
987  self%counter_right = self%counter_right+1
988  end if
989  dx = (xbar- xi(1))/(xnew(1)-xi(1))
990  xmid = xi + dx * (xnew-xi)
991  xmid(1) = xbar
992 
993  vh = (xmid-xi) * wi(1)
994  call self%particle_mesh_coupling%add_current_evaluate( xi, xmid, vh, self%efield_dofs_work, &
995  self%j_dofs_local, efield )
996 
997  xt = xmid
998  xt(2:3) = modulo(xt(2:3), 1._f64)
999  jmatrix = self%map%jacobian_matrix_inverse_transposed(xt)
1000  do j = 1, 3
1001  vi(j) = vi(j) + dx*sign *0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
1002  end do
1003  select case(self%boundary_particles)
1004  case(sll_p_boundary_particles_singular)
1005  call self%particle_mesh_coupling%add_charge(xmid, wi(1), self%spline_degree, self%rhob)
1006  xmid(2) = xbar + (1._f64-2._f64*xbar)*xmid(2) + 0.5_f64-0.5_f64*xbar
1007  xt(2:3) = modulo(xt(2:3), 1._f64)
1008  jmatrix = self%map%jacobian_matrix_inverse_transposed(xt)
1009  call self%particle_mesh_coupling%add_charge(xmid, -wi(1), self%spline_degree, self%rhob)
1010  xnew(1) = 2._f64*xbar-xnew(1)
1011  xnew(2) = xbar + (1._f64-2._f64*xbar)*xnew(2) + 0.5_f64-0.5_f64*xbar
1012  if(xnew(1) > 1._f64 )then
1013  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
1014  end if
1015  case(sll_p_boundary_particles_reflection)
1016  xnew(1) = 2._f64*xbar-xnew(1)
1017  vi = vi-2._f64*(vi(1)*jmatrix(1,1)+vi(2)*jmatrix(2,1)+vi(3)*jmatrix(3,1))*jmatrix(:,1)/sum(jmatrix(:,1)**2)
1018  case(sll_p_boundary_particles_absorption)
1019  call self%particle_mesh_coupling%add_charge(xmid, wi(1), self%spline_degree, self%rhob)
1020  case( sll_p_boundary_particles_periodic)
1021  xnew(1) = modulo(xnew(1), 1._f64)
1022  xmid(1) = 1._f64-xbar
1023  case default
1024  xnew(1) = modulo(xnew(1), 1._f64)
1025  xmid(1) = 1._f64-xbar
1026  end select
1027  vh = (xnew - xmid)*wi(1)
1028  call self%particle_mesh_coupling%add_current_evaluate( xmid, xnew, vh, self%efield_dofs_work, &
1029  self%j_dofs_local, efield )
1030  xnew(2:3) = modulo(xnew(2:3), 1._f64)
1031  jmat = self%map%jacobian_matrix_inverse_transposed(xnew)
1032  do j = 1, 3
1033  vi(j) = vi(j) + (1._f64-dx)*sign *0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
1034  end do
1035  else
1036  vh = (xnew-xi) * wi(1)
1037  call self%particle_mesh_coupling%add_current_evaluate( xi, xnew, vh, self%efield_dofs_work, &
1038  self%j_dofs_local, efield )
1039  xnew(2:3) = modulo(xnew(2:3), 1._f64)
1040  jmatrix = self%map%jacobian_matrix_inverse_transposed(xnew)
1041  do j = 1, 3
1042  vi(j) = vi(j) + sign *0.5_f64*((jmatrix(j,1)+jmat(j,1))*efield(1) + (jmatrix(j,2)+jmat(j,2))*efield(2) + (jmatrix(j,3)+ jmat(j,3))*efield(3))
1043  end do
1044  end if
1045 
1046  end subroutine compute_particle_boundary_current_evaluate
1047 
1048 
1049  !---------------------------------------------------------------------------!
1051  subroutine initialize_pic_vm_3d3v_trafo(&
1052  self, &
1053  maxwell_solver, &
1054  particle_mesh_coupling, &
1055  particle_group, &
1056  phi_dofs, &
1057  efield_dofs, &
1058  bfield_dofs, &
1059  x_min, &
1060  Lx, &
1061  map, &
1062  boundary_particles, &
1063  solver_tolerance, &
1064  iter_tolerance, max_iter, &
1065  betar,&
1066  force_sign, &
1067  rhob, &
1068  control_variate, &
1069  jmean)
1070  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( out ) :: self
1071  class(sll_c_maxwell_3d_base), target, intent( in ) :: maxwell_solver
1072  class(sll_c_particle_mesh_coupling_3d), target, intent(in) :: particle_mesh_coupling
1073  class(sll_t_particle_array), target, intent( in ) :: particle_group
1074  sll_real64, target, intent( in ) :: phi_dofs(:)
1075  sll_real64, target, intent( in ) :: efield_dofs(:)
1076  sll_real64, target, intent( in ) :: bfield_dofs(:)
1077  sll_real64, intent( in ) :: x_min(3)
1078  sll_real64, intent( in ) :: lx(3)
1079  type(sll_t_mapping_3d), target, intent( inout ) :: map
1080  sll_int32, optional, intent( in ) :: boundary_particles
1081  sll_real64, optional, intent( in ) :: solver_tolerance
1082  sll_real64, optional, intent( in ) :: iter_tolerance
1083  sll_int32, optional, intent( in ) :: max_iter
1084  sll_real64, optional, intent( in ) :: betar(2)
1085  sll_real64, optional, intent(in) :: force_sign
1086  sll_real64, optional, target, intent( in ) :: rhob(:)
1087  class(sll_t_control_variates), optional, target, intent(in) :: control_variate
1088  logical, optional, intent(in) :: jmean
1089  !local variables
1090  sll_int32 :: ierr, j
1091 
1092  if (present(solver_tolerance) ) then
1093  self%solver_tolerance = solver_tolerance
1094  end if
1095 
1096  if (present(iter_tolerance) ) then
1097  self%iter_tolerance = iter_tolerance
1098  end if
1099 
1100  if (present(max_iter) ) then
1101  self%max_iter = max_iter
1102  end if
1103 
1104  if( present(force_sign) )then
1105  self%force_sign = force_sign
1106  end if
1107 
1108  if (self%force_sign == 1._f64) then
1109  if( particle_group%group(1)%species%q > 0._f64) self%adiabatic_electrons = .true.
1110  end if
1111 
1112  if (present(jmean)) then
1113  self%jmean = jmean
1114  end if
1115 
1116  if (present(boundary_particles) ) then
1117  self%boundary_particles = boundary_particles
1118  end if
1119 
1120  if( present(rhob) )then
1121  self%rhob => rhob
1122  end if
1123 
1124 
1125  self%maxwell_solver => maxwell_solver
1126  self%particle_mesh_coupling => particle_mesh_coupling
1127  self%particle_group => particle_group
1128  self%phi_dofs => phi_dofs
1129  self%efield_dofs => efield_dofs
1130  self%bfield_dofs => bfield_dofs
1131  self%map => map
1132  self%n_total0 = self%maxwell_solver%n_total0
1133  self%n_total1 = self%maxwell_solver%n_total1
1134  self%spline_degree = self%particle_mesh_coupling%spline_degree
1135  self%x_min = self%map%get_x([0._f64,0._f64,0._f64])
1136  self%Lx = self%map%Lx
1137  self%x_max = x_min + lx
1138 
1139 
1140  self%nspan_1(1) = (2*self%spline_degree(1)-1)*(2*self%spline_degree(2)+1)*(2*self%spline_degree(3)+1)
1141  self%nspan_1(2) = (2*self%spline_degree(1)+1)*(2*self%spline_degree(2)-1)*(2*self%spline_degree(3)+1)
1142  self%nspan_1(3) = (2*self%spline_degree(1)+1)*(2*self%spline_degree(2)+1)*(2*self%spline_degree(3)-1)
1143 
1144  self%nspan_2(1) = (2*self%spline_degree(1))*(2*self%spline_degree(2))*(2*self%spline_degree(3)+1)
1145  self%nspan_2(2) = (2*self%spline_degree(1)+1)*(2*self%spline_degree(2))*(2*self%spline_degree(3))
1146  self%nspan_2(3) = (2*self%spline_degree(1))*(2*self%spline_degree(2)+1)*(2*self%spline_degree(3))
1147 
1148  sll_allocate( self%vec1(1:self%n_total1+self%n_total0*2), ierr )
1149  sll_allocate( self%vec2(1:self%n_total1+self%n_total0*2), ierr )
1150  self%vec1 = 0._f64
1151  self%vec2 = 0._f64
1152  sll_allocate( self%j_dofs(1:self%n_total1+self%n_total0*2), ierr )
1153  sll_allocate( self%j_dofs_local(1:self%n_total1+self%n_total0*2), ierr )
1154  sll_allocate( self%particle_mass_1_local(1:self%nspan_1(1), 1:self%n_total1), ierr )
1155  sll_allocate( self%particle_mass_2_local(1:self%nspan_1(2), 1:self%n_total0), ierr )
1156  sll_allocate( self%particle_mass_3_local(1:self%nspan_1(3), 1:self%n_total0), ierr )
1157  sll_allocate( self%particle_mass_12_local(1:self%nspan_2(1), 1:self%n_total1), ierr )
1158  sll_allocate( self%particle_mass_23_local(1:self%nspan_2(2), 1:self%n_total0), ierr )
1159  sll_allocate( self%particle_mass_31_local(1:self%nspan_2(3), 1:self%n_total0), ierr )
1160  sll_allocate( self%particle_mass_21_local(1:self%nspan_2(1), 1:self%n_total0), ierr )
1161  sll_allocate( self%particle_mass_32_local(1:self%nspan_2(2), 1:self%n_total0), ierr )
1162  sll_allocate( self%particle_mass_13_local(1:self%nspan_2(3), 1:self%n_total1), ierr )
1163 
1164  if( self%particle_mesh_coupling%n_cells(1)+self%spline_degree(1) == self%maxwell_solver%n_dofs(1) ) then
1165  self%boundary = .true.
1166  allocate( sll_t_linear_operator_particle_mass_cl_3d_diag :: self%particle_mass_1 )
1167  allocate( sll_t_linear_operator_particle_mass_cl_3d_diag :: self%particle_mass_2 )
1168  allocate( sll_t_linear_operator_particle_mass_cl_3d_diag :: self%particle_mass_3 )
1169  if(self%map%flag2d)then
1170  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_12 )
1171  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_21 )
1172  if(self%map%flag3d)then
1173  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_13 )
1174  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_31 )
1175  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_23 )
1176  allocate( sll_t_linear_operator_particle_mass_cl_3d_od :: self%particle_mass_32 )
1177  end if
1178  end if
1179  else
1180  allocate( sll_t_linear_operator_particle_mass_3d_diag :: self%particle_mass_1 )
1181  allocate( sll_t_linear_operator_particle_mass_3d_diag :: self%particle_mass_2 )
1182  allocate( sll_t_linear_operator_particle_mass_3d_diag :: self%particle_mass_3 )
1183  if(self%map%flag2d)then
1184  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_12 )
1185  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_21 )
1186  if(self%map%flag3d)then
1187  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_13 )
1188  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_31 )
1189  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_23 )
1190  allocate( sll_t_linear_operator_particle_mass_3d_od :: self%particle_mass_32 )
1191  end if
1192  end if
1193  end if
1194 
1195  call self%particle_mass_1%create( self%particle_mesh_coupling%n_cells, [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)] )
1196  call self%particle_mass_2%create( self%particle_mesh_coupling%n_cells, [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)] )
1197  call self%particle_mass_3%create( self%particle_mesh_coupling%n_cells, [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1] )
1198 
1199  if(self%map%flag2d)then
1200  call self%particle_mass_12%create( self%particle_mesh_coupling%n_cells, &
1201  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
1202  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)] )
1203  call self%particle_mass_21%create( self%particle_mesh_coupling%n_cells, &
1204  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
1205  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)] )
1206  if(self%map%flag3d)then
1207  call self%particle_mass_23%create( self%particle_mesh_coupling%n_cells, &
1208  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)], &
1209  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1] )
1210  call self%particle_mass_32%create( self%particle_mesh_coupling%n_cells, &
1211  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
1212  [self%spline_degree(1), self%spline_degree(2)-1, self%spline_degree(3)] )
1213  call self%particle_mass_31%create( self%particle_mesh_coupling%n_cells, &
1214  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1], &
1215  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)] )
1216  call self%particle_mass_13%create( self%particle_mesh_coupling%n_cells, &
1217  [self%spline_degree(1)-1, self%spline_degree(2), self%spline_degree(3)], &
1218  [self%spline_degree(1), self%spline_degree(2), self%spline_degree(3)-1] )
1219  end if
1220  end if
1221  call self%particle_mass_op%create( 3, 3 )
1222 
1223  if( self%adiabatic_electrons ) then
1224  call self%linear_op_schur_phiv%create( self%maxwell_solver, self%particle_mass_op )
1225  call self%linear_solver_schur_phiv%create( self%linear_op_schur_phiv )
1226  self%linear_solver_schur_phiv%atol = self%solver_tolerance
1227  !self%linear_solver_schur_phiv%verbose = .true.
1228  else
1229  call self%linear_op_schur_ev%create( self%maxwell_solver, self%particle_mass_op )
1230  call self%preconditioner_fft%init( self%maxwell_solver%Lx, self%maxwell_solver%n_cells, self%maxwell_solver%s_deg_0, self%boundary )
1231 
1232  if( self%boundary ) then
1233  self%vec1 = 1._f64
1234  call self%maxwell_solver%multiply_mass( [1], self%vec1, self%vec2 )
1235  do j = 1, self%maxwell_solver%n_dofs(2)*self%maxwell_solver%n_dofs(3)
1236  self%vec2(self%n_total1+1+(j-1)*self%maxwell_solver%n_dofs(1)) = 1._f64
1237  self%vec2(self%n_total1+j*self%maxwell_solver%n_dofs(1)) = 1._f64
1238  self%vec2(self%n_total1+self%n_total0+1+(j-1)*self%maxwell_solver%n_dofs(1)) = 1._f64
1239  self%vec2(self%n_total1+self%n_total0+j*self%maxwell_solver%n_dofs(1)) = 1._f64
1240  end do
1241  self%vec2 = 1._f64/sqrt(abs(self%vec2))
1242  call self%preconditioner1%create( self%preconditioner_fft%inverse_mass1_3d, self%vec2, 2*self%n_total0+self%n_total1 )
1243  call self%linear_solver_schur_ev%create( self%linear_op_schur_ev, self%preconditioner1)
1244  else
1245  call self%linear_solver_schur_ev%create( self%linear_op_schur_ev, self%preconditioner_fft%inverse_mass1_3d )
1246  end if
1247  self%linear_solver_schur_ev%atol = self%solver_tolerance/maxval(self%Lx)
1248  !self%linear_solver_schur_ev%verbose = .true.
1249  end if
1250 
1251  if (present(betar)) then
1252  self%betar = betar
1253  else
1254  self%betar = 1.0_f64
1255  end if
1256 
1257  self%n_particles_max = self%particle_group%group(1)%n_particles
1258  do j = 2,self%particle_group%n_species
1259  self%n_particles_max = max(self%n_particles_max, self%particle_group%group(j)%n_particles )
1260  end do
1261 
1262  sll_allocate( self%xnew(self%particle_group%n_species,3,self%n_particles_max), ierr )
1263  sll_allocate( self%vnew(self%particle_group%n_species,3,self%n_particles_max), ierr )
1264  sll_allocate( self%efield_dofs_new(self%n_total1+2*self%n_total0), ierr )
1265  sll_allocate( self%efield_dofs_work(self%n_total1+2*self%n_total0), ierr )
1266  sll_allocate( self%phi_dofs_work(self%n_total0), ierr )
1267  sll_allocate( self%phi_dofs_new(self%n_total0), ierr )
1268  self%xnew = 0._f64
1269  self%vnew = 0._f64
1270  self%efield_dofs_new = 0._f64
1271  self%phi_dofs_new = 0._f64
1272  self%phi_dofs_work = 0._f64
1273 
1274  if (present(control_variate)) then
1275  allocate(self%control_variate )
1276  allocate(self%control_variate%cv(self%particle_group%n_species) )
1277  self%control_variate => control_variate
1278  self%i_weight = self%particle_group%group(1)%n_weights
1279  end if
1280 
1281  end subroutine initialize_pic_vm_3d3v_trafo
1282 
1283 
1285  subroutine initialize_file_pic_vm_3d3v_trafo(&
1286  self, &
1287  maxwell_solver, &
1288  particle_mesh_coupling, &
1289  particle_group, &
1290  phi_dofs, &
1291  efield_dofs, &
1292  bfield_dofs, &
1293  x_min, &
1294  Lx, &
1295  map, &
1296  filename, &
1297  boundary_particles, &
1298  betar, &
1299  force_sign, &
1300  rhob, &
1301  control_variate, &
1302  jmean)
1303  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( out ) :: self
1304  class(sll_c_maxwell_3d_base), target, intent( in ) :: maxwell_solver
1305  class(sll_c_particle_mesh_coupling_3d), target, intent(in) :: particle_mesh_coupling
1306  class(sll_t_particle_array), target, intent( in ) :: particle_group
1307  sll_real64, target, intent( in ) :: phi_dofs(:)
1308  sll_real64, target, intent( in ) :: efield_dofs(:)
1309  sll_real64, target, intent( in ) :: bfield_dofs(:)
1310  sll_real64, intent( in ) :: x_min(3)
1311  sll_real64, intent( in ) :: lx(3)
1312  type(sll_t_mapping_3d), target, intent( inout ) :: map
1313  character(len=*), intent( in ) :: filename
1314  sll_int32, optional, intent( in ) :: boundary_particles
1315  sll_real64, optional, intent( in ) :: betar(2)
1316  sll_real64, optional, intent(in) :: force_sign
1317  sll_real64, optional, target, intent( in ) :: rhob(:)
1318  class(sll_t_control_variates), optional, target, intent(in) :: control_variate
1319  logical, optional, intent(in) :: jmean
1320  !local variables
1321  character(len=256) :: file_prefix
1322  sll_int32 :: input_file, rank
1323  sll_int32 :: io_stat, io_stat0, io_stat1, file_id
1324  sll_real64 :: maxwell_tolerance, iter_tolerance, betar_set(2), force_sign_set
1325  sll_int32 :: max_iter, boundary_particles_set
1326  logical :: jmean_set
1327 
1328  namelist /output/ file_prefix
1329  namelist /time_solver/ maxwell_tolerance
1330  namelist /time_iterate/ iter_tolerance, max_iter
1331 
1332  rank = sll_f_get_collective_rank(sll_v_world_collective)
1333 
1334  if( present(boundary_particles) )then
1335  boundary_particles_set = boundary_particles
1336  else
1337  boundary_particles_set = 100
1338  end if
1339 
1340  if( present(betar) )then
1341  betar_set = betar
1342  else
1343  betar_set = 1._f64
1344  end if
1345 
1346  if( present(force_sign) )then
1347  force_sign_set = force_sign
1348  else
1349  force_sign_set = 1._f64
1350  end if
1351 
1352  if (present(jmean)) then
1353  jmean_set = jmean
1354  else
1355  jmean_set = .false.
1356  end if
1357 
1358  if( present( control_variate) ) then
1359  ! Read in solver tolerance
1360  open(newunit = input_file, file=filename, status='old',iostat=io_stat)
1361  if (io_stat /= 0) then
1362  if (rank == 0 ) then
1363  print*, 'sll_m_time_propagator_pic_vm_3d3v_trafo_helper: Input file does not exist. Set default tolerance.'
1364  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1365  write(file_id, *) 'solver tolerance:', 1d-12
1366  write(file_id, *) 'iter_tolerance:', 1d-12
1367  write(file_id, *) 'max_iter:', 10
1368  close(file_id)
1369  end if
1370  call self%init( maxwell_solver, &
1371  particle_mesh_coupling, &
1372  particle_group, &
1373  phi_dofs, &
1374  efield_dofs, &
1375  bfield_dofs, &
1376  x_min, &
1377  lx, &
1378  map, &
1379  boundary_particles_set, &
1380  betar=betar_set, &
1381  force_sign=force_sign_set, &
1382  rhob = rhob, &
1383  control_variate = control_variate,&
1384  jmean=jmean_set)
1385  else
1386  read(input_file, output,iostat=io_stat0)
1387  read(input_file, time_solver,iostat=io_stat)
1388  read(input_file, time_iterate,iostat=io_stat1)
1389  if (io_stat /= 0 .and. io_stat1 /= 0) then
1390  if (rank == 0 ) then
1391  print*, 'sll_m_time_propagator_pic_vm_3d3v_trafo_helper: Input parameter does not exist. Set default tolerance.'
1392  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1393  write(file_id, *) 'solver tolerance:', 1d-12
1394  write(file_id, *) 'iter_tolerance:', 1d-12
1395  write(file_id, *) 'max_iter:', 10
1396  close(file_id)
1397  end if
1398  call self%init( maxwell_solver, &
1399  particle_mesh_coupling, &
1400  particle_group, &
1401  phi_dofs, &
1402  efield_dofs, &
1403  bfield_dofs, &
1404  x_min, &
1405  lx, &
1406  map, &
1407  boundary_particles_set, &
1408  betar=betar_set, &
1409  force_sign=force_sign_set, &
1410  rhob = rhob, &
1411  control_variate = control_variate,&
1412  jmean=jmean_set)
1413  else if (io_stat == 0 .and. io_stat1 /= 0) then
1414  if (rank == 0 ) then
1415  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1416  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1417  write(file_id, *) 'iter_tolerance:', 1d-12
1418  write(file_id, *) 'max_iter:', 10
1419  close(file_id)
1420  end if
1421  call self%init( maxwell_solver, &
1422  particle_mesh_coupling, &
1423  particle_group, &
1424  phi_dofs, &
1425  efield_dofs, &
1426  bfield_dofs, &
1427  x_min, &
1428  lx, &
1429  map, &
1430  boundary_particles_set, &
1431  maxwell_tolerance, &
1432  betar=betar_set, &
1433  force_sign=force_sign_set, &
1434  rhob = rhob, &
1435  control_variate = control_variate)
1436  else if (io_stat /= 0 .and. io_stat1 == 0) then
1437  if (rank == 0 ) then
1438  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1439  write(file_id, *) 'solver tolerance:', 1d-12
1440  write(file_id, *) 'iter_tolerance:', iter_tolerance
1441  write(file_id, *) 'max_iter:', max_iter
1442  close(file_id)
1443  end if
1444  call self%init( maxwell_solver, &
1445  particle_mesh_coupling, &
1446  particle_group, &
1447  phi_dofs, &
1448  efield_dofs, &
1449  bfield_dofs, &
1450  x_min, &
1451  lx, &
1452  map, &
1453  boundary_particles_set, &
1454  iter_tolerance=iter_tolerance, max_iter=max_iter, &
1455  betar=betar_set, &
1456  force_sign=force_sign_set, &
1457  rhob = rhob, &
1458  control_variate = control_variate,&
1459  jmean=jmean_set)
1460  else
1461  if (rank == 0 ) then
1462  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1463  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1464  write(file_id, *) 'iter_tolerance:', iter_tolerance
1465  write(file_id, *) 'max_iter:', max_iter
1466  close(file_id)
1467  end if
1468  call self%init( maxwell_solver, &
1469  particle_mesh_coupling, &
1470  particle_group, &
1471  phi_dofs, &
1472  efield_dofs, &
1473  bfield_dofs, &
1474  x_min, &
1475  lx, &
1476  map, &
1477  boundary_particles_set, &
1478  maxwell_tolerance, &
1479  iter_tolerance, max_iter, &
1480  betar_set, &
1481  force_sign=force_sign_set, &
1482  rhob = rhob, &
1483  control_variate = control_variate,&
1484  jmean=jmean_set)
1485  end if
1486  close(input_file)
1487  end if
1488  else
1489  ! Read in solver tolerance
1490  open(newunit = input_file, file=filename, status='old',iostat=io_stat)
1491  if (io_stat /= 0) then
1492  if (rank == 0 ) then
1493  print*, 'sll_m_time_propagator_pic_vm_3d3v_trafo_helper: Input file does not exist. Set default tolerance.'
1494  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1495  write(file_id, *) 'solver tolerance:', 1d-12
1496  write(file_id, *) 'iter_tolerance:', 1d-12
1497  write(file_id, *) 'max_iter:', 10
1498  close(file_id)
1499  end if
1500  call self%init( maxwell_solver, &
1501  particle_mesh_coupling, &
1502  particle_group, &
1503  phi_dofs, &
1504  efield_dofs, &
1505  bfield_dofs, &
1506  x_min, &
1507  lx, &
1508  map, &
1509  boundary_particles_set, &
1510  betar=betar_set, &
1511  force_sign=force_sign_set,&
1512  rhob = rhob, &
1513  jmean=jmean_set)
1514  else
1515  read(input_file, output,iostat=io_stat0)
1516  read(input_file, time_solver,iostat=io_stat)
1517  read(input_file, time_iterate,iostat=io_stat1)
1518  if (io_stat /= 0 .and. io_stat1 /= 0) then
1519  if (rank == 0 ) then
1520  print*, 'sll_m_time_propagator_pic_vm_3d3v_trafo_helper: Input parameter does not exist. Set default tolerance.'
1521  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1522  write(file_id, *) 'solver tolerance:', 1d-12
1523  write(file_id, *) 'iter_tolerance:', 1d-12
1524  write(file_id, *) 'max_iter:', 10
1525  close(file_id)
1526  end if
1527  call self%init( maxwell_solver, &
1528  particle_mesh_coupling, &
1529  particle_group, &
1530  phi_dofs, &
1531  efield_dofs, &
1532  bfield_dofs, &
1533  x_min, &
1534  lx, &
1535  map, &
1536  boundary_particles_set, &
1537  betar=betar_set, &
1538  force_sign=force_sign_set,&
1539  rhob = rhob, &
1540  jmean=jmean_set)
1541  else if (io_stat == 0 .and. io_stat1 /= 0) then
1542  if (rank == 0 ) then
1543  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1544  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1545  write(file_id, *) 'iter_tolerance:', 1d-12
1546  write(file_id, *) 'max_iter:', 10
1547  close(file_id)
1548  end if
1549  call self%init( maxwell_solver, &
1550  particle_mesh_coupling, &
1551  particle_group, &
1552  phi_dofs, &
1553  efield_dofs, &
1554  bfield_dofs, &
1555  x_min, &
1556  lx, &
1557  map, &
1558  boundary_particles_set, &
1559  maxwell_tolerance, &
1560  betar=betar_set, &
1561  force_sign=force_sign_set,&
1562  rhob = rhob, &
1563  jmean=jmean_set)
1564  else if (io_stat /= 0 .and. io_stat1 == 0) then
1565  if (rank == 0 ) then
1566  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1567  write(file_id, *) 'solver tolerance:', 1d-12
1568  write(file_id, *) 'iter_tolerance:', iter_tolerance
1569  write(file_id, *) 'max_iter:', max_iter
1570  close(file_id)
1571  end if
1572  call self%init( maxwell_solver, &
1573  particle_mesh_coupling, &
1574  particle_group, &
1575  phi_dofs, &
1576  efield_dofs, &
1577  bfield_dofs, &
1578  x_min, &
1579  lx, &
1580  map, &
1581  boundary_particles_set, &
1582  iter_tolerance=iter_tolerance, max_iter=max_iter, &
1583  betar=betar_set, &
1584  force_sign=force_sign_set,&
1585  rhob = rhob, &
1586  jmean=jmean_set)
1587  else
1588  if (rank == 0 ) then
1589  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1590  write(file_id, *) 'solver tolerance:', maxwell_tolerance
1591  write(file_id, *) 'iter_tolerance:', iter_tolerance
1592  write(file_id, *) 'max_iter:', max_iter
1593  close(file_id)
1594  end if
1595  call self%init( maxwell_solver, &
1596  particle_mesh_coupling, &
1597  particle_group, &
1598  phi_dofs, &
1599  efield_dofs, &
1600  bfield_dofs, &
1601  x_min, &
1602  lx, &
1603  map, &
1604  boundary_particles_set, &
1605  maxwell_tolerance, &
1606  iter_tolerance, max_iter, &
1607  betar_set, &
1608  force_sign=force_sign_set,&
1609  rhob = rhob, &
1610  jmean=jmean_set)
1611  end if
1612  close(input_file)
1613  end if
1614  end if
1615 
1616 
1617  end subroutine initialize_file_pic_vm_3d3v_trafo
1618 
1619 
1620  !---------------------------------------------------------------------------!
1622  subroutine delete_pic_vm_3d3v_trafo(self)
1623  class(sll_t_time_propagator_pic_vm_3d3v_trafo_helper), intent( inout ) :: self
1624  !local variables
1625  sll_int32 :: file_id, j
1626 
1627  if(self%boundary)then
1628  print*, 'left boundary', self%counter_left
1629  print*, 'right boundary', self%counter_right
1630  end if
1631  if(self%adiabatic_electrons) then
1632  call self%linear_solver_schur_phiv%free()
1633  call self%linear_op_schur_phiv%free()
1634  else
1635  call self%preconditioner_fft%free()
1636  call self%linear_solver_schur_ev%free()
1637  call self%linear_op_schur_ev%free()
1638  end if
1639  call self%particle_mass_op%free()
1640  call self%particle_mass_1%free()
1641  call self%particle_mass_2%free()
1642  call self%particle_mass_3%free()
1643  deallocate( self%particle_mass_1 )
1644  deallocate( self%particle_mass_2 )
1645  deallocate( self%particle_mass_3 )
1646  if(self%map%flag2d)then
1647  call self%particle_mass_12%free()
1648  call self%particle_mass_21%free()
1649  deallocate( self%particle_mass_12 )
1650  deallocate( self%particle_mass_21 )
1651  if(self%map%flag3d)then
1652  call self%particle_mass_23%free()
1653  call self%particle_mass_32%free()
1654  call self%particle_mass_31%free()
1655  call self%particle_mass_13%free()
1656  deallocate( self%particle_mass_13 )
1657  deallocate( self%particle_mass_31 )
1658  deallocate( self%particle_mass_23 )
1659  deallocate( self%particle_mass_32 )
1660  end if
1661  end if
1662 
1663  deallocate( self%j_dofs )
1664  deallocate( self%j_dofs_local )
1665  deallocate( self%particle_mass_1_local )
1666  deallocate( self%particle_mass_2_local )
1667  deallocate( self%particle_mass_3_local )
1668  deallocate( self%particle_mass_12_local )
1669  deallocate( self%particle_mass_23_local )
1670  deallocate( self%particle_mass_31_local )
1671 
1672  self%maxwell_solver => null()
1673  self%particle_mesh_coupling => null()
1674  self%particle_group => null()
1675  self%efield_dofs => null()
1676  self%bfield_dofs => null()
1677  self%map => null()
1678 
1679  print*, 'No. of failed iterations:', self%n_failed
1680  open(newunit=file_id, file="outer-iters.dat", action='write')
1681  do j=1,self%iter_counter
1682  write(file_id,*) self%niter(j)
1683  end do
1684 
1685  end subroutine delete_pic_vm_3d3v_trafo
1686 
1687 
1688 end module sll_m_time_propagator_pic_vm_3d3v_trafo_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_initial_distribution
Parameters to define common initial distributions.
Definition: sll_m_initial_distribution.F90:5
sll_m_linear_operator_block
module for a block linear operator
Definition: sll_m_linear_operator_block.F90:11
sll_m_linear_operator_particle_mass_3d_diag
Definition: sll_m_linear_operator_particle_mass_3d_diag.F90:1
sll_m_linear_operator_particle_mass_3d_od
Definition: sll_m_linear_operator_particle_mass_3d_od.F90:1
sll_m_linear_operator_particle_mass_cl_3d_diag
Definition: sll_m_linear_operator_particle_mass_cl_3d_diag.F90:1
sll_m_linear_operator_particle_mass_cl_3d_od
Definition: sll_m_linear_operator_particle_mass_cl_3d_od.F90:1
sll_m_linear_operator_schur_ev_3d
Definition: sll_m_linear_operator_schur_ev_3d.F90:1
sll_m_linear_operator_schur_phiv_3d
Definition: sll_m_linear_operator_schur_phiv_3d.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_mapping_3d
Module interfaces for coordinate transformation.
Definition: sll_m_mapping_3d.F90:13
sll_m_maxwell_3d_base
Module interface to solve Maxwell's equations in 3D.
Definition: sll_m_maxwell_3d_base.F90:7
sll_m_particle_group_base
Definition: sll_m_particle_group_base.F90:1
sll_m_particle_mass_3d_base
Definition: sll_m_particle_mass_3d_base.F90:1
sll_m_particle_mesh_coupling_base_3d
Base class for kernel smoothers for accumulation and field evaluation in PIC.
Definition: sll_m_particle_mesh_coupling_base_3d.F90:5
sll_m_preconditioner_fft
Module interface to solve Maxwell's equations with coordinate transformation in 3D The linear systems...
Definition: sll_m_preconditioner_fft.F90:10
sll_m_preconditioner_singular
Module interface to solve Maxwell's equations with coordinate transformation in 3D The linear systems...
Definition: sll_m_preconditioner_singular.F90:10
sll_m_profile_functions
functions for initial profile of the particle distribution function
Definition: sll_m_profile_functions.F90:4
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_s_compute_particle_boundary_simple
subroutine, public sll_s_compute_particle_boundary_simple(boundary_particles, counter_left, counter_right, xold, xnew)
compute new position
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:345
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_s_compute_matrix_inverse
subroutine, public sll_s_compute_matrix_inverse(x, y, bf, jm, sign)
invert matrix
Definition: sll_m_time_propagator_pic_vm_3d3v_cl_helper.F90:68
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
sll_m_time_propagator_pic_vm_3d3v_trafo_helper
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:3
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::delete_pic_vm_3d3v_trafo
subroutine delete_pic_vm_3d3v_trafo(self)
Destructor.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:1623
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::initialize_file_pic_vm_3d3v_trafo
subroutine initialize_file_pic_vm_3d3v_trafo(self, maxwell_solver, particle_mesh_coupling, particle_group, phi_dofs, efield_dofs, bfield_dofs, x_min, Lx, map, filename, boundary_particles, betar, force_sign, rhob, control_variate, jmean)
Constructor.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:1303
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::advect_ex_pic_vm_3d3v_trafo
subroutine advect_ex_pic_vm_3d3v_trafo(self, dt)
advect_ex: Equations to be solved $\frac{\Xi^{n+1}-\Xi^n}{\Delta t}=\frac{DF^{-1}(\Xi^{n+1})+DF^{-1}(...
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:655
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::advect_vb_pic_vm_3d3v_trafo
subroutine advect_vb_pic_vm_3d3v_trafo(self, dt)
advect_vb: Equations to be solved $(\mathbb{I}-\Delta \frac{\Delta t q}{2 m} DF^{-\top} \mathbb{B}(\X...
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:307
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::compute_particle_boundary_current_evaluate_iter
subroutine compute_particle_boundary_current_evaluate_iter(self, xi, xnew, vi, wi, sign)
Helper function for advect_ex.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:872
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::advect_e_pic_vm_3d3v_trafo
subroutine advect_e_pic_vm_3d3v_trafo(self, dt)
advect_e: Equations to be solved Solution with Schur complement: $ S_{+}=M_1+\frac{\Delta t^2 q^2}{4 ...
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:379
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::sll_s_compute_particle_boundary
subroutine sll_s_compute_particle_boundary(self, xold, xnew, vi)
Helper function for advect_x.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:251
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::advect_x_pic_vm_3d3v_trafo
subroutine advect_x_pic_vm_3d3v_trafo(self, dt)
Finalization.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:193
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::compute_particle_boundary_current_evaluate
subroutine compute_particle_boundary_current_evaluate(self, xi, xnew, vi, wi, sign)
Helper function for advect_ex.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:965
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::advect_eb_pic_vm_3d3v_trafo
subroutine advect_eb_pic_vm_3d3v_trafo(self, dt)
advect_eb: Equations to be solved Solution with Schur complement: $ S=M_1+\frac{\Delta t^2}{4} C^\top...
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:365
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::initialize_pic_vm_3d3v_trafo
subroutine initialize_pic_vm_3d3v_trafo(self, maxwell_solver, particle_mesh_coupling, particle_group, phi_dofs, efield_dofs, bfield_dofs, x_min, Lx, map, boundary_particles, solver_tolerance, iter_tolerance, max_iter, betar, force_sign, rhob, control_variate, jmean)
Constructor.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:1070
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_initial_distribution::sll_t_params_cos_gaussian_screwpinch
Definition: sll_m_initial_distribution.F90:100
sll_m_linear_operator_block::sll_t_linear_operator_block
class for a linear operator_block
Definition: sll_m_linear_operator_block.F90:34
sll_m_linear_operator_particle_mass_3d_diag::sll_t_linear_operator_particle_mass_3d_diag
Definition: sll_m_linear_operator_particle_mass_3d_diag.F90:13
sll_m_linear_operator_particle_mass_3d_od::sll_t_linear_operator_particle_mass_3d_od
Definition: sll_m_linear_operator_particle_mass_3d_od.F90:14
sll_m_linear_operator_particle_mass_cl_3d_diag::sll_t_linear_operator_particle_mass_cl_3d_diag
Definition: sll_m_linear_operator_particle_mass_cl_3d_diag.F90:13
sll_m_linear_operator_particle_mass_cl_3d_od::sll_t_linear_operator_particle_mass_cl_3d_od
Definition: sll_m_linear_operator_particle_mass_cl_3d_od.F90:14
sll_m_linear_operator_schur_ev_3d::sll_t_linear_operator_schur_ev_3d
Definition: sll_m_linear_operator_schur_ev_3d.F90:16
sll_m_linear_operator_schur_phiv_3d::sll_t_linear_operator_schur_phiv_3d
Definition: sll_m_linear_operator_schur_phiv_3d.F90:16
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_mapping_3d::sll_t_mapping_3d
type collecting functions for analytical coordinate mapping
Definition: sll_m_mapping_3d.F90:68
sll_m_maxwell_3d_base::sll_c_maxwell_3d_base
Definition: sll_m_maxwell_3d_base.F90:27
sll_m_particle_group_base::sll_t_particle_array
Definition: sll_m_particle_group_base.F90:100
sll_m_particle_mass_3d_base::sll_c_particle_mass_3d_base
Definition: sll_m_particle_mass_3d_base.F90:13
sll_m_particle_mesh_coupling_base_3d::sll_c_particle_mesh_coupling_3d
Basic type of a kernel smoother used for PIC simulations.
Definition: sll_m_particle_mesh_coupling_base_3d.F90:26
sll_m_preconditioner_fft::sll_t_preconditioner_fft
Definition: sll_m_preconditioner_fft.F90:35
sll_m_preconditioner_singular::sll_t_preconditioner_singular
Definition: sll_m_preconditioner_singular.F90:30
sll_m_profile_functions::sll_t_profile_functions
Definition: sll_m_profile_functions.F90:19
sll_m_time_propagator_pic_vm_3d3v_trafo_helper::sll_t_time_propagator_pic_vm_3d3v_trafo_helper
Helper for implicit time propagator for 3d3v Vlasov-Maxwell with coordinate transformation.
Definition: sll_m_time_propagator_pic_vm_3d3v_trafo_helper.F90:90
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1