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