1!! Copyright (C) 2019 N. Tancogne-Dejean
2!! Copyright (C) 2020 M. Oliveira
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
22module classical_particle_oct_m
23  use clock_oct_m
24  use global_oct_m
25  use interaction_abst_oct_m
26  use interaction_gravity_oct_m
27  use ghost_interaction_oct_m
28  use io_oct_m
29  use iso_c_binding
30  use messages_oct_m
31  use mpi_oct_m
32  use namespace_oct_m
33  use parser_oct_m
34  use profiling_oct_m
35  use propagator_abst_oct_m
36  use quantity_oct_m
37  use space_oct_m
38  use system_abst_oct_m
39  use write_iter_oct_m
40
41  implicit none
42
43  private
44  public ::                 &
45    classical_particle_t,   &
46    classical_particle_init
47
48  type, extends(system_abst_t) :: classical_particle_t
49    private
50
51    FLOAT :: mass
52    FLOAT :: pos(1:MAX_DIM)
53    FLOAT :: vel(1:MAX_DIM)
54    FLOAT :: acc(1:MAX_DIM)
55    FLOAT, allocatable :: prev_acc(:,:) !< A storage of the prior times.
56    FLOAT :: save_pos(1:MAX_DIM)   !< A storage for the SCF loops
57    FLOAT :: save_vel(1:MAX_DIM)   !< A storage for the SCF loops
58    FLOAT :: tot_force(1:MAX_DIM)
59    FLOAT :: prev_tot_force(1:MAX_DIM) !< Used for the SCF convergence criterium
60
61    type(c_ptr) :: output_handle
62  contains
63    procedure :: add_interaction_partner => classical_particle_add_interaction_partner
64    procedure :: has_interaction => classical_particle_has_interaction
65    procedure :: initial_conditions => classical_particle_initial_conditions
66    procedure :: do_td_operation => classical_particle_do_td
67    procedure :: iteration_info => classical_particle_iteration_info
68    procedure :: output_start => classical_particle_output_start
69    procedure :: output_write => classical_particle_output_write
70    procedure :: output_finish => classical_particle_output_finish
71    procedure :: is_tolerance_reached => classical_particle_is_tolerance_reached
72    procedure :: store_current_status => classical_particle_store_current_status
73    procedure :: update_quantity => classical_particle_update_quantity
74    procedure :: update_exposed_quantity => classical_particle_update_exposed_quantity
75    procedure :: copy_quantities_to_interaction => classical_particle_copy_quantities_to_interaction
76    procedure :: update_interactions_start => classical_particle_update_interactions_start
77    procedure :: update_interactions_finish => classical_particle_update_interactions_finish
78    final :: classical_particle_finalize
79  end type classical_particle_t
80
81  interface classical_particle_t
82    procedure classical_particle_constructor
83  end interface classical_particle_t
84
85contains
86
87  ! ---------------------------------------------------------
88  !> The factory routine (or constructor) allocates a pointer of the
89  !! corresponding type and then calls the init routine which is a type-bound
90  !! procedure of the corresponding type. With this design, also derived
91  !! classes can use the init routine of the parent class.
92  function classical_particle_constructor(namespace) result(sys)
93    class(classical_particle_t), pointer    :: sys
94    type(namespace_t),           intent(in) :: namespace
95
96    PUSH_SUB(classical_particle_constructor)
97
98    SAFE_ALLOCATE(sys)
99
100    call classical_particle_init(sys, namespace)
101
102    POP_SUB(classical_particle_constructor)
103  end function classical_particle_constructor
104
105  ! ---------------------------------------------------------
106  !> The init routine is a module level procedure
107  !! This has the advantage that different classes can have different
108  !! signatures for the initialization routines because they are not
109  !! type-bound and thus also not inherited.
110  ! ---------------------------------------------------------
111  subroutine classical_particle_init(this, namespace)
112    class(classical_particle_t), intent(inout) :: this
113    type(namespace_t),           intent(in)    :: namespace
114
115    PUSH_SUB(classical_particle_init)
116
117    this%namespace = namespace
118
119    call messages_print_stress(stdout, "Classical Particle", namespace=namespace)
120
121    call space_init(this%space, namespace)
122
123    !%Variable ClassicalParticleMass
124    !%Type float
125    !%Section ClassicalParticles
126    !%Description
127    !% Mass of classical particle in Kg.
128    !%End
129    call parse_variable(namespace, 'ClassicalParticleMass', M_ONE, this%mass)
130    call messages_print_var_value(stdout, 'ClassicalParticleMass', this%mass)
131
132    this%quantities(POSITION)%required = .true.
133    this%quantities(VELOCITY)%required = .true.
134    this%quantities(POSITION)%protected = .true.
135    this%quantities(VELOCITY)%protected = .true.
136
137    call messages_print_stress(stdout, namespace=namespace)
138
139    POP_SUB(classical_particle_init)
140  end subroutine classical_particle_init
141
142  ! ---------------------------------------------------------
143  subroutine classical_particle_add_interaction_partner(this, partner)
144    class(classical_particle_t), target, intent(inout) :: this
145    class(system_abst_t),            intent(inout) :: partner
146
147    class(ghost_interaction_t), pointer :: ghost
148    class(interaction_gravity_t), pointer :: gravity
149    type(interaction_gravity_t) :: gravity_t
150
151    PUSH_SUB(classical_particle_add_interaction_partner)
152
153    if (partner%has_interaction(gravity_t)) then
154      gravity => interaction_gravity_t(this%space%dim, partner)
155      this%quantities(POSITION)%required = .true.
156      this%quantities(MASS)%required = .true.
157      gravity%system_mass => this%mass
158      gravity%system_pos  => this%pos
159      call this%interactions%add(gravity)
160    else
161      ghost => ghost_interaction_t(partner)
162      call this%interactions%add(ghost)
163    end if
164
165    POP_SUB(classical_particle_add_interaction_partner)
166  end subroutine classical_particle_add_interaction_partner
167
168  ! ---------------------------------------------------------
169  logical function classical_particle_has_interaction(this, interaction)
170    class(classical_particle_t),   intent(in) :: this
171    class(interaction_abst_t), intent(in) :: interaction
172
173    PUSH_SUB(classical_particle_has_interaction)
174
175    select type (interaction)
176    type is (interaction_gravity_t)
177      classical_particle_has_interaction = .true.
178    class default
179      classical_particle_has_interaction = .false.
180    end select
181
182    POP_SUB(classical_particle_has_interaction)
183  end function classical_particle_has_interaction
184
185  ! ---------------------------------------------------------
186  subroutine classical_particle_initial_conditions(this, from_scratch)
187    class(classical_particle_t), intent(inout) :: this
188    logical,                 intent(in)    :: from_scratch
189
190    integer :: n_rows, idir
191    type(block_t) :: blk
192
193    PUSH_SUB(classical_particle_initial_conditions)
194
195    if (.not. from_scratch) then
196      message(1) = "Classical particle propagation is currently only allowed from scratch"
197      call messages_fatal(1, namespace=this%namespace)
198    end if
199
200    !%Variable ClassicalParticleInitialPosition
201    !%Type block
202    !%Section ClassicalParticles
203    !%Description
204    !% Initial position of classical particle, in Km.
205    !%End
206    this%pos = M_ZERO
207    if (parse_block(this%namespace, 'ClassicalParticleInitialPosition', blk) == 0) then
208      n_rows = parse_block_n(blk)
209      if (n_rows > 1) call  messages_input_error(this%namespace, 'ClassicalParticleInitialPosition')
210
211      do idir = 1, this%space%dim
212        call parse_block_float(blk, 0, idir - 1, this%pos(idir))
213      end do
214      call parse_block_end(blk)
215    end if
216    call messages_print_var_value(stdout, 'ClassicalParticleInitialPosition', this%pos(1:this%space%dim))
217
218    !%Variable ClassicalParticleInitialVelocity
219    !%Type block
220    !%Section ClassicalParticles
221    !%Description
222    !% Initial velocity of classical particle in Km/s.
223    !%End
224    this%vel = M_ZERO
225    if (parse_block(this%namespace, 'ClassicalParticleInitialVelocity', blk) == 0) then
226      n_rows = parse_block_n(blk)
227      if (n_rows > 1) call  messages_input_error(this%namespace, 'ClassicalParticleInitialVelocity')
228      do idir = 1, this%space%dim
229        call parse_block_float(blk, 0, idir - 1, this%vel(idir))
230      end do
231      call parse_block_end(blk)
232    end if
233    call messages_print_var_value(stdout, 'ClassicalParticleInitialVelocity', this%vel(1:this%space%dim))
234
235    POP_SUB(classical_particle_initial_conditions)
236  end subroutine classical_particle_initial_conditions
237
238  ! ---------------------------------------------------------
239  subroutine classical_particle_do_td(this, operation)
240    class(classical_particle_t), intent(inout) :: this
241    integer,                 intent(in)    :: operation
242
243    integer :: ii
244
245    PUSH_SUB(classical_particle_do_td)
246
247    select case(operation)
248    case (VERLET_START)
249      SAFE_ALLOCATE(this%prev_acc(1:this%space%dim, 1))
250      this%acc(1:this%space%dim) = this%tot_force(1:this%space%dim) / this%mass
251
252    case (VERLET_FINISH, BEEMAN_FINISH)
253      SAFE_DEALLOCATE_A(this%prev_acc)
254
255    case (VERLET_UPDATE_POS)
256      this%pos(1:this%space%dim) = this%pos(1:this%space%dim) + this%prop%dt * this%vel(1:this%space%dim) &
257                                 + M_HALF * this%prop%dt**2 * this%acc(1:this%space%dim)
258
259      call this%quantities(POSITION)%clock%increment()
260
261    case (VERLET_COMPUTE_ACC)
262      do ii = size(this%prev_acc, dim=2) - 1, 1, -1
263        this%prev_acc(1:this%space%dim, ii + 1) = this%prev_acc(1:this%space%dim, ii)
264      end do
265      this%prev_acc(1:this%space%dim, 1) = this%acc(1:this%space%dim)
266      this%acc(1:this%space%dim) = this%tot_force(1:this%space%dim) / this%mass
267
268    case (VERLET_COMPUTE_VEL)
269      this%vel(1:this%space%dim) = this%vel(1:this%space%dim) &
270        + M_HALF * this%prop%dt * (this%prev_acc(1:this%space%dim, 1) + this%acc(1:this%space%dim))
271
272      call this%quantities(VELOCITY)%clock%increment()
273
274
275    case (BEEMAN_START)
276      SAFE_ALLOCATE(this%prev_acc(1:this%space%dim, 2))
277      this%acc(1:this%space%dim) = this%tot_force(1:this%space%dim) / this%mass
278      this%prev_acc(1:this%space%dim, 1) = this%acc(1:this%space%dim)
279
280    case (BEEMAN_PREDICT_POS)
281      this%pos(1:this%space%dim) = this%pos(1:this%space%dim) + this%prop%dt * this%vel(1:this%space%dim) &
282                                 + M_ONE/CNST(6.0) * this%prop%dt**2  &
283                                 * (M_FOUR*this%acc(1:this%space%dim) - this%prev_acc(1:this%space%dim, 1))
284
285      if (.not. this%prop%predictor_corrector) then
286        call this%quantities(POSITION)%clock%increment()
287      end if
288
289    case (BEEMAN_PREDICT_VEL)
290      this%vel(1:this%space%dim) = this%vel(1:this%space%dim)  &
291        + M_ONE/CNST(6.0) * this%prop%dt * ( M_TWO * this%acc(1:this%space%dim) + &
292        CNST(5.0) * this%prev_acc(1:this%space%dim, 1) - this%prev_acc(1:this%space%dim, 2))
293
294      call this%quantities(VELOCITY)%clock%increment()
295
296    case (BEEMAN_CORRECT_POS)
297      this%pos(1:this%space%dim) = this%save_pos(1:this%space%dim) + this%prop%dt * this%save_vel(1:this%space%dim) &
298                                 + M_ONE/CNST(6.0) * this%prop%dt**2  &
299                                 * (this%acc(1:this%space%dim) + M_TWO * this%prev_acc(1:this%space%dim, 1))
300
301      ! We set it to the propagation time to avoid double increment
302      call this%quantities(POSITION)%clock%set_time(this%prop%clock)
303
304    case (BEEMAN_CORRECT_VEL)
305      this%vel(1:this%space%dim) = this%save_vel(1:this%space%dim) &
306                                 + M_HALF * this%prop%dt * (this%acc(1:this%space%dim) + this%prev_acc(1:this%space%dim, 1))
307
308      ! We set it to the propagation time to avoid double increment
309      call this%quantities(VELOCITY)%clock%set_time(this%prop%clock)
310
311    case default
312      message(1) = "Unsupported TD operation."
313      call messages_fatal(1, namespace=this%namespace)
314    end select
315
316   POP_SUB(classical_particle_do_td)
317  end subroutine classical_particle_do_td
318
319  ! ---------------------------------------------------------
320  logical function classical_particle_is_tolerance_reached(this, tol) result(converged)
321    class(classical_particle_t),   intent(in)    :: this
322    FLOAT,                     intent(in)    :: tol
323
324    PUSH_SUB(classical_particle_is_tolerance_reached)
325
326    ! Here we put the criterion that acceleration change is below the tolerance
327    converged = .false.
328    if ( (sum((this%prev_tot_force(1:this%space%dim) - this%tot_force(1:this%space%dim))**2)/ this%mass) < tol**2) then
329      converged = .true.
330    end if
331
332    if (debug%info) then
333      write(message(1), '(a, e12.6, a, e12.6)') "Debug: -- Change in acceleration  ", &
334        sqrt(sum((this%prev_tot_force(1:this%space%dim) - this%tot_force(1:this%space%dim))**2))/this%mass, " and tolerance ", tol
335      call messages_info(1)
336    end if
337
338    POP_SUB(classical_particle_is_tolerance_reached)
339  end function classical_particle_is_tolerance_reached
340
341  ! ---------------------------------------------------------
342  subroutine classical_particle_store_current_status(this)
343    class(classical_particle_t),   intent(inout)    :: this
344
345    PUSH_SUB(classical_particle_store_current_status)
346
347    this%save_pos(1:this%space%dim) = this%pos(1:this%space%dim)
348    this%save_vel(1:this%space%dim) = this%vel(1:this%space%dim)
349
350    POP_SUB(classical_particle_store_current_status)
351  end subroutine classical_particle_store_current_status
352
353  ! ---------------------------------------------------------
354  subroutine classical_particle_iteration_info(this)
355    class(classical_particle_t), intent(in) :: this
356
357    integer :: idir
358    character(len=20) :: fmt
359
360    PUSH_SUB(classical_particle_iteration_info)
361
362    write(message(1),'(2X,A,1X,A)') "Classical particle:", trim(this%namespace%get())
363
364    write(fmt,'("(4X,A,1X,",I2,"e14.6)")') this%space%dim
365    write(message(2),fmt) "Coordinates: ", (this%pos(idir), idir = 1, this%space%dim)
366    write(message(3),fmt) "Velocity:    ", (this%vel(idir), idir = 1, this%space%dim)
367    write(message(4),fmt) "Acceleration:", (this%acc(idir), idir = 1, this%space%dim)
368    write(message(5),fmt) "Force:       ", (this%tot_force(idir), idir = 1, this%space%dim)
369    write(message(6),'(4x,A,I8.7)')  'Clock tick:      ', this%clock%get_tick()
370    write(message(7),'(4x,A,e14.6)') 'Simulation time: ', this%clock%get_sim_time()
371    call messages_info(7)
372
373    POP_SUB(classical_particle_iteration_info)
374  end subroutine classical_particle_iteration_info
375
376  ! ---------------------------------------------------------
377  subroutine classical_particle_output_start(this)
378    class(classical_particle_t), intent(inout) :: this
379
380    PUSH_SUB(classical_particle_output_start)
381
382    ! Create output handle
383    call io_mkdir('td.general', this%namespace)
384    if (mpi_grp_is_root(mpi_world)) then
385      call write_iter_init(this%output_handle, 0, this%prop%dt, trim(io_workpath("td.general/coordinates", this%namespace)))
386    end if
387
388    ! Output info for first iteration
389    call this%output_write(0)
390
391    POP_SUB(classical_particle_output_start)
392  end subroutine classical_particle_output_start
393
394  ! ---------------------------------------------------------
395  subroutine classical_particle_output_finish(this)
396    class(classical_particle_t), intent(inout) :: this
397
398    PUSH_SUB(classical_particle_output_finish)
399
400    if (mpi_grp_is_root(mpi_world)) then
401      call write_iter_end(this%output_handle)
402    end if
403
404    POP_SUB(classical_particle_output_finish)
405  end subroutine classical_particle_output_finish
406
407  ! ---------------------------------------------------------
408  subroutine classical_particle_output_write(this, iter)
409    class(classical_particle_t), intent(inout) :: this
410    integer,                 intent(in)    :: iter
411
412    integer :: idir
413    character(len=50) :: aux
414
415    if(.not.mpi_grp_is_root(mpi_world)) return ! only first node outputs
416
417    PUSH_SUB(classical_particle_output_write)
418
419    if(iter == 0) then
420      ! header
421      call write_iter_clear(this%output_handle)
422      call write_iter_string(this%output_handle,'################################################################################')
423      call write_iter_nl(this%output_handle)
424      call write_iter_string(this%output_handle,'# HEADER')
425      call write_iter_nl(this%output_handle)
426
427      ! first line: column names
428      call write_iter_header_start(this%output_handle)
429
430      do idir = 1, this%space%dim
431        write(aux, '(a2,i3,a1)') 'x(', idir, ')'
432        call write_iter_header(this%output_handle, aux)
433      end do
434      do idir = 1, this%space%dim
435        write(aux, '(a2,i3,a1)') 'v(', idir, ')'
436        call write_iter_header(this%output_handle, aux)
437      end do
438      do idir = 1, this%space%dim
439        write(aux, '(a2,i3,a1)') 'f(', idir, ')'
440        call write_iter_header(this%output_handle, aux)
441      end do
442      call write_iter_nl(this%output_handle)
443
444      ! second line: units
445      call write_iter_string(this%output_handle, '#[Iter n.]')
446      call write_iter_header(this%output_handle, '[seconds]')
447      call write_iter_string(this%output_handle, 'Positions in Km, Velocities in Km/s')
448      call write_iter_nl(this%output_handle)
449
450      call write_iter_string(this%output_handle,'################################################################################')
451      call write_iter_nl(this%output_handle)
452    end if
453
454    call write_iter_start(this%output_handle)
455    call write_iter_double(this%output_handle, this%pos, this%space%dim)
456    call write_iter_double(this%output_handle, this%vel, this%space%dim)
457    call write_iter_double(this%output_handle, this%tot_force, this%space%dim)
458    call write_iter_nl(this%output_handle)
459
460    POP_SUB(classical_particle_output_write)
461  end subroutine classical_particle_output_write
462
463  ! ---------------------------------------------------------
464  subroutine classical_particle_update_quantity(this, iq, requested_time)
465    class(classical_particle_t),   intent(inout) :: this
466    integer,                   intent(in)    :: iq
467    class(clock_t),            intent(in)    :: requested_time
468
469    PUSH_SUB(classical_particle_update_quantity)
470
471    ! We are not allowed to update protected quantities!
472    ASSERT(.not. this%quantities(iq)%protected)
473
474    select case (iq)
475    case (MASS)
476      ! The classical particle has a mass, but it is not necessary to update it, as it does not change with time.
477      call this%quantities(iq)%clock%set_time(requested_time)
478    case default
479      message(1) = "Incompatible quantity."
480      call messages_fatal(1)
481    end select
482
483    POP_SUB(classical_particle_update_quantity)
484  end subroutine classical_particle_update_quantity
485
486  ! ---------------------------------------------------------
487  subroutine classical_particle_update_exposed_quantity(this, iq, requested_time)
488    class(classical_particle_t),   intent(inout) :: this
489    integer,                   intent(in)    :: iq
490    class(clock_t),            intent(in)    :: requested_time
491
492    PUSH_SUB(classical_particle_update_exposed_quantity)
493
494    ! We are not allowed to update protected quantities!
495    ASSERT(.not. this%quantities(iq)%protected)
496
497    select case (iq)
498    case (MASS)
499      ! The classical particle has a mass, but it does not require any update, as it does not change with time.
500      call this%quantities(iq)%clock%set_time(this%clock)
501    case default
502      message(1) = "Incompatible quantity."
503      call messages_fatal(1)
504    end select
505
506    POP_SUB(classical_particle_update_exposed_quantity)
507  end subroutine classical_particle_update_exposed_quantity
508
509  ! ---------------------------------------------------------
510  subroutine classical_particle_copy_quantities_to_interaction(this, interaction)
511    class(classical_particle_t),          intent(inout) :: this
512    class(interaction_abst_t),        intent(inout) :: interaction
513
514    PUSH_SUB(classical_particle_copy_quantities_to_interaction)
515
516    select type (interaction)
517    type is (interaction_gravity_t)
518      interaction%partner_mass = this%mass
519      interaction%partner_pos = this%pos
520    type is (ghost_interaction_t)
521      ! Nothing to copy
522    class default
523      message(1) = "Unsupported interaction."
524      call messages_fatal(1)
525    end select
526
527    POP_SUB(classical_particle_copy_quantities_to_interaction)
528  end subroutine classical_particle_copy_quantities_to_interaction
529
530  ! ---------------------------------------------------------
531  subroutine classical_particle_update_interactions_start(this)
532    class(classical_particle_t), intent(inout) :: this
533
534    PUSH_SUB(classical_particle_update_interactions_start)
535
536    ! Store previous force, as it is used as SCF criterium
537    this%prev_tot_force(1:this%space%dim) = this%tot_force(1:this%space%dim)
538
539    POP_SUB(classical_particle_update_interactions_start)
540  end subroutine classical_particle_update_interactions_start
541
542  ! ---------------------------------------------------------
543  subroutine classical_particle_update_interactions_finish(this)
544    class(classical_particle_t), intent(inout) :: this
545
546    type(interaction_iterator_t) :: iter
547
548    PUSH_SUB(classical_particle_update_interactions_finish)
549
550    ! Compute the total force acting on the classical particle
551    this%tot_force(1:this%space%dim) = M_ZERO
552    call iter%start(this%interactions)
553    do while (iter%has_next())
554      select type (interaction => iter%get_next_interaction())
555      type is (interaction_gravity_t)
556        this%tot_force(1:this%space%dim) = this%tot_force(1:this%space%dim) + interaction%force(1:this%space%dim)
557      type is (ghost_interaction_t)
558        ! Nothing to do
559      class default
560        message(1) = "Unknown interaction by the classical particle " + this%namespace%get()
561        call messages_fatal(1)
562      end select
563    end do
564
565    POP_SUB(classical_particle_update_interactions_finish)
566  end subroutine classical_particle_update_interactions_finish
567
568  ! ---------------------------------------------------------
569  subroutine classical_particle_finalize(this)
570    type(classical_particle_t), intent(inout) :: this
571
572    type(interaction_iterator_t) :: iter
573    class(interaction_abst_t), pointer :: interaction
574
575    PUSH_SUB(classical_particle_finalize)
576
577    deallocate(this%prop)
578
579    call iter%start(this%interactions)
580    do while (iter%has_next())
581      interaction => iter%get_next_interaction()
582      SAFE_DEALLOCATE_P(interaction)
583    end do
584
585    POP_SUB(classical_particle_finalize)
586  end subroutine classical_particle_finalize
587
588end module classical_particle_oct_m
589
590!! Local Variables:
591!! mode: f90
592!! coding: utf-8
593!! End:
594