DISPATCH
particle_solver_mod.f90
1 !===============================================================================
2 !> Kick-Drift-Kick (KDK) N-body solver
3 !===============================================================================
5  USE io_unit_mod
6  USE trace_mod
7  USE particle_mod
9  USE dll_mod
10  implicit none
11  private
12  type, public, extends(particle_list_t):: particle_solver_t
13  contains
14  procedure:: init
15  procedure:: update
16  procedure:: courant_time
17  end type
18  real(8), save:: courant=0.5
19  type(particle_solver_t), target, public:: particle_solver
20 CONTAINS
21 
22 !===============================================================================
23 !> Initialize a particle list
24 !===============================================================================
25 SUBROUTINE init (self, name)
26  class(particle_solver_t):: self
27  character(len=*), optional:: name
28  logical, save:: first_time=.true.
29  namelist /kdk_params/ courant
30  !-----------------------------------------------------------------------------
31  call self%particle_list_t%init
32  if (first_time) then
33  !$omp critical (input_cr)
34  if (first_time) then
35  first_time = .false.
36  rewind(io_unit%input)
37  read (io_unit%input, kdk_params)
38  write (io_unit%output, kdk_params)
39  end if
40  !$omp end critical (input_cr)
41  end if
42 END SUBROUTINE init
43 
44 !===============================================================================
45 !> KDK update
46 !===============================================================================
47 SUBROUTINE update (self, p)
48  class(particle_solver_t):: self
49  class(particle_t), pointer:: p
50  real(8):: dt, a(3)
51  integer:: new, nt
52  !.............................................................................
53  new = 1 + mod(p%it,self%nt)
54  associate(v=>p%v(:,new), r=>p%r(:,new))
55  dt = self%courant_time (p) ! estimate dt
56  a = self%force(p)/p%mass ! acceleration
57  v = v + 0.5_8*a*dt ! kick 1
58  r = r + v*dt ! drift
59  a = self%force(p)/p%mass ! acceleration
60  v = v + 0.5_8*a*dt ! kick 2
61  end associate
62  p%time = p%time+dt
63  p%t(new) = p%time
64  p%it = new
65  nt = self%nt
66  p%iit(1:nt-1) = p%iit(2:nt)
67  p%iit(nt) = new
68 END SUBROUTINE update
69 
70 !===============================================================================
71 !> Courant conditions: Estimate, as inexpensively as possible, an approximately
72 !> reflexive time step. If an estimate of the acceleration results in it not
73 !> being important for the timestep size, then choose a fraction of the cell
74 !> travel time (whatever that means in this context)
75 !===============================================================================
76 FUNCTION courant_time (self, p) RESULT (dt)
77  class(particle_solver_t):: self
78  class(particle_t), pointer:: p
79  real(8):: dt
80  !.............................................................................
81  real(8):: a(3), ac, vc, dt1, dt2
82  integer:: new, i0, i1
83  !-----------------------------------------------------------------------------
84  new = 1 + mod(p%it,self%nt)
85  i1 = p%iit(self%nt-1)
86  i0 = p%iit(self%nt-2)
87  associate(v=>p%v(:,i1))
88  a = (p%v(:,i1)-p%v(:,i0))/(p%t(i1)-p%t(i0))
89  ac = 1.0/sqrt(sum(a**2)+tiny(1.0_8))
90  vc = 1.0/sqrt(sum(v**2)+tiny(1.0_8))
91  dt1 = p%ds*vc*courant
92  dt2 = p%ds*ac/vc
93  if (dt2 < dt1) then
94  a = self%force(p)/p%mass
95  ac = 1.0/sqrt(sum(a**2)+tiny(1.0_8))
96  dt2 = p%ds*ac/vc*courant
97  end if
98  dt = min(dt1,dt2)
99  end associate
100 END FUNCTION courant_time
101 
102 END MODULE particle_solver_mod
Particle data type, extends a dll_node data type, so it can become part of a particle_list_t data typ...
Definition: particle_mod.f90:5
Doubly linked list (DLL), carrying anything, as simply as possible.
Definition: dll_mod.f90:4
Kick-Drift-Kick (KDK) N-body solver.
Particle list, extending a doubly-linked list. Each particle maintains arrays with previous positions...