forked from garrelt/C2-Ray3Dm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
evolve8.F90
1898 lines (1564 loc) · 63 KB
/
evolve8.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!>
!! \brief This module contains routines for calculating the ionization and temperature evolution of the entire grid (3D).
!!
!! Module for Capreole / C2-Ray (f90)
!!
!! \b Author: Garrelt Mellema
!!
!! \b Date:
!!
!! \b Version: 3D, no OpenMP
module evolve
! This module contains routines having to do with the calculation of
! the ionization evolution of the entire grid (3D).
! This version has been adapted for efficiency in order to be able
! to calculate large meshes.
! - evolve3D : step through grid
! - evolve2D : step through z-plane
! - evolve0D : take care of one grid point
! - evolve0D_global: update entire grid
! Needs:
! doric : ionization calculation for one point + photo-ionization rates
! tped : temperature,pressure,electron density calculation
use precision, only: dp
use my_mpi ! supplies all the MPI and OpenMP definitions and variables
use file_admin, only: logf,iterdump
use sizes, only: Ndim, mesh
use grid, only: x,y,z,vol,dr
use material, only: ndens, xh, temper, get_temperature_point
use sourceprops, only: NumSrc, srcpos, NormFlux !SrcSeries
use radiation, only: NumFreqBnd
use photonstatistics, only: state_before, calculate_photon_statistics, &
photon_loss, LLS_loss, report_photonstatistics, state_after, total_rates, &
total_ionizations, update_grandtotal_photonstatistics
use c2ray_parameters, only: convergence_fraction, subboxsize, max_subbox, S_star_nominal
implicit none
private
public :: evolve3D, phih_grid, evolve_ini
!> Periodic boundary conditions, has to be true for this version
logical,parameter :: periodic_bc = .true.
!> Minimum number of MPI processes for using the master-slave setup
integer, parameter :: min_numproc_master_slave=10
! Grid variables
!> H Photo-ionization rate on the entire grid
real(kind=dp),dimension(:,:,:),allocatable :: phih_grid
!> Time-averaged H ionization fraction
real(kind=dp),dimension(:,:,:,:),allocatable :: xh_av
!real(kind=dp),dimension(:,:,:),allocatable :: xh_av
!> Intermediate result for H ionization fraction
real(kind=dp),dimension(:,:,:,:),allocatable :: xh_intermed
!real(kind=dp),dimension(:,:,:),allocatable :: xh_intermed
!> H0 Column density (outgoing)
real(kind=dp),dimension(:,:,:),allocatable :: coldensh_out
!> Buffer for MPI communication
real(kind=dp),dimension(:,:,:),allocatable :: buffer
!> Photon loss from the grid
real(kind=dp) :: photon_loss_all(1:NumFreqBnd)
!> Photon loss from one source
real(kind=dp),dimension(:,:),allocatable :: photon_loss_src_thread
real(kind=dp) :: photon_loss_src(1:NumFreqBnd)
! mesh positions of end points for RT
integer,dimension(Ndim) :: lastpos_l !< mesh position of left end point for RT
integer,dimension(Ndim) :: lastpos_r !< mesh position of right end point for RT
integer,dimension(Ndim) :: last_l !< mesh position of left end point for RT
integer,dimension(Ndim) :: last_r !< mesh position of right end point for RT
integer :: sum_nbox !< sum of all nboxes (on one processor)
integer :: sum_nbox_all !< sum of all nboxes (on all processors)
integer :: tn !< thread number
contains
! =======================================================================
!> Allocate the arrays needed for evolve
subroutine evolve_ini ()
allocate(phih_grid(mesh(1),mesh(2),mesh(3)))
allocate(xh_av(mesh(1),mesh(2),mesh(3),0:1))
!allocate(xh_av(mesh(1),mesh(2),mesh(3)))
allocate(xh_intermed(mesh(1),mesh(2),mesh(3),0:1))
!allocate(xh_intermed(mesh(1),mesh(2),mesh(3)))
allocate(coldensh_out(mesh(1),mesh(2),mesh(3)))
allocate(buffer(mesh(1),mesh(2),mesh(3)))
allocate(photon_loss_src_thread(1:NumFreqBnd,nthreads))
end subroutine evolve_ini
! ===========================================================================
!> Evolve the entire grid over a time step dt
subroutine evolve3D (dt,restart)
! Calculates the evolution of the hydrogen ionization state
! Author: Garrelt Mellema
! Date: 28-Feb-2008 (21-Aug-2006 (f77/OMP: 13-Jun-2005))
! Version: Multiple sources / Using average fractions to converge
! loop over sources
! History:
! 11-Jun-2004 (GM) : grid arrays now passed via common (in grid.h)
! and material arrays also (in material.h).
! 11-Jun-2004 (GM) : adapted for multiple sources.
! 3-Jan-2005 (GM) : reintegrated with updated Ifront3D
! 20-May-2005 (GM) : split original eveolve0D into two routines
! 13-Jun-2005 (HM) : OpenMP version : Hugh Merz
! 21-Aug-2006 (GM) : MPI parallelization over the sources (static model).
! 28-Feb-2008 (GM) : Added master-slave model for distributing
! over the processes. The program decides which
! model to use.
! The time step
real(kind=dp),intent(in) :: dt !< time step
integer,intent(in) :: restart !< restart flag
! Loop variables
integer :: niter ! iteration counter
! Wall clock counting
! 8 bytes to beat the maxcount
integer(kind=8) :: wallclock1
integer(kind=8) :: wallclock2
integer(kind=8) :: countspersec
! Flag variable (passed back from evolve0D_global)
integer :: conv_flag
! Minimum number of cells which are allowed to be non-converged
integer :: conv_criterion
#ifdef MPI
integer :: mympierror
#endif
! End of declarations
! Initialize wall clock counter (for dumps)
call system_clock(wallclock1)
! Initial state (for photon statistics)
call state_before (xh)
! initialize average and intermediate results to initial values
if (restart == 0) then
xh_av(:,:,:,:)=xh(:,:,:,:)
xh_intermed(:,:,:,:)=xh(:,:,:,:)
!xh_av(:,:,:)=xh(:,:,:,1)
!xh_intermed(:,:,:)=xh(:,:,:,1)
niter=0 ! iteration starts at zero
conv_flag=mesh(1)*mesh(2)*mesh(3) ! initialize non-convergence
else
! Reload xh_av,xh_intermed,photon_loss,niter
call start_from_dump(niter)
call global_pass (conv_flag,dt)
endif
! Set the conv_criterion, if there are few sources we should make
! sure that things are converged around these sources.
conv_criterion=min(int(convergence_fraction*mesh(1)*mesh(2)*mesh(3)), &
(NumSrc-1)/3)
! Iterate to reach convergence for multiple sources
do
! Update xh if converged and exit
if (conv_flag <= conv_criterion) then
xh(:,:,:,:)=xh_intermed(:,:,:,:)
!xh(:,:,:)=xh_intermed(:,:,:)
exit
else
if (niter > 100) then
! Complain about slow convergence
if (rank == 0) write(logf,*) 'Multiple sources not converging'
exit
endif
endif
! Iteration loop counter
niter=niter+1
call pass_all_sources (niter,dt)
! Report subbox statistics
if (rank == 0) &
write(logf,*) "Average number of subboxes: ",sum_nbox_all/NumSrc
if (rank == 0) then
call system_clock(wallclock2,countspersec)
! Write iteration dump if more than 15 minutes have passed.
! system_clock starts counting at 0 when it reaches
! a max value. To catch this, test also for negative
! values of wallclock2-wallclock1
write(logf,*) "Time and limit are: ", &
wallclock2-wallclock1, 15.0*60.0*countspersec
if (wallclock2-wallclock1 > 15*60*countspersec .or. &
wallclock2-wallclock1 < 0 ) then
call write_iteration_dump(niter)
wallclock1=wallclock2
endif
endif
call global_pass (conv_flag,dt)
enddo
! Calculate photon statistics
call calculate_photon_statistics (dt,xh,xh_av)
call report_photonstatistics (dt)
call update_grandtotal_photonstatistics (dt)
end subroutine evolve3D
! ===========================================================================
subroutine write_iteration_dump (niter)
integer,intent(in) :: niter ! iteration counter
integer :: ndump=0
character(len=20) :: iterfile
ndump=ndump+1
if (mod(ndump,2) == 0) then
iterfile="iterdump2.bin"
else
iterfile="iterdump1.bin"
endif
open(unit=iterdump,file=iterfile,form="unformatted", &
status="unknown")
write(iterdump) niter
write(iterdump) photon_loss_all
write(iterdump) phih_grid
write(iterdump) xh_av
write(iterdump) xh_intermed
close(iterdump)
end subroutine write_iteration_dump
! ===========================================================================
subroutine start_from_dump(niter)
integer,intent(out) :: niter ! iteration counter
#ifdef MPI
integer :: mympierror
#endif
if (rank == 0) then
open(unit=iterdump,file="iterdump.bin",form="unformatted", &
status="old")
read(iterdump) niter
read(iterdump) photon_loss_all
read(iterdump) phih_grid
read(iterdump) xh_av
read(iterdump) xh_intermed
close(iterdump)
write(logf,*) "Read iteration ",niter," from dump file"
write(logf,*) 'photon loss counter: ',photon_loss_all
write(logf,*) "Intermediate result for mean ionization fraction: ", &
sum(xh_intermed(:,:,:,1))/real(mesh(1)*mesh(2)*mesh(3))
!sum(xh_intermed(:,:,:))/real(mesh(1)*mesh(2)*mesh(3))
endif
#ifdef MPI
! Distribute the input parameters to the other nodes
call MPI_BCAST(niter,1, &
MPI_INTEGER,0,MPI_COMM_NEW,mympierror)
call MPI_BCAST(photon_loss_all,NumFreqBnd, &
MPI_DOUBLE_PRECISION,0,MPI_COMM_NEW,mympierror)
call MPI_BCAST(phih_grid,mesh(1)*mesh(2)*mesh(3), &
MPI_DOUBLE_PRECISION,0,MPI_COMM_NEW,mympierror)
call MPI_BCAST(xh_av,mesh(1)*mesh(2)*mesh(3)*2, &
!call MPI_BCAST(xh_av,mesh(1)*mesh(2)*mesh(3), &
MPI_DOUBLE_PRECISION,0,MPI_COMM_NEW,mympierror)
call MPI_BCAST(xh_intermed,mesh(1)*mesh(2)*mesh(3)*2, &
!call MPI_BCAST(xh_intermed,mesh(1)*mesh(2)*mesh(3), &
MPI_DOUBLE_PRECISION,0,&
MPI_COMM_NEW,mympierror)
#endif
end subroutine start_from_dump
! ===========================================================================
subroutine pass_all_sources(niter,dt)
! For random permutation of sources:
use m_ctrper, only: ctrper
integer,intent(in) :: niter ! iteration counter
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
real(kind=dp) :: LLS_loss_all
#ifdef MPI
integer :: mympierror
#endif
if (rank == 0) write(logf,*) 'Doing all sources '
! reset global rates to zero for this iteration
phih_grid(:,:,:)=0.0
! reset photon loss counters
photon_loss(:) = 0.0
LLS_loss = 0.0 ! make this a NumFreqBnd vector if needed later (GM/101129)
! Reset sum of subboxes counter
sum_nbox=0
! Make a randomized list of sources :: call in serial
! disabled / GM110512
!if ( rank == 0 ) call ctrper (SrcSeries(1:NumSrc),1.0)
#ifdef MPI
! Distribute the source list to the other nodes
! disabled / GM110512
!call MPI_BCAST(SrcSeries,NumSrc,MPI_INTEGER,0,MPI_COMM_NEW,mympierror)
#endif
! Ray trace the whole grid for all sources.
! We can do this in two ways, depending on
! the number of processors. For many processors
! the master-slave setup should be more efficient.
if (npr > min_numproc_master_slave) then
call do_grid_master_slave (dt,niter)
else
call do_grid_static (dt,niter)
endif
#ifdef MPI
! accumulate (sum) the MPI distributed photon losses
call MPI_ALLREDUCE(photon_loss, photon_loss_all, NumFreqBnd, &
MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_NEW, mympierror)
! accumulate (sum) the MPI distributed photon losses
call MPI_ALLREDUCE(LLS_loss, LLS_loss_all, 1, &
MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_NEW, mympierror)
! Put LLS_loss_all back in the LLS variable
LLS_loss = LLS_loss_all
! accumulate (sum) the MPI distributed phih_grid
call MPI_ALLREDUCE(phih_grid, buffer, mesh(1)*mesh(2)*mesh(3), &
MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_NEW, mympierror)
! Overwrite the processor local values with the accumulated value
phih_grid(:,:,:)=buffer(:,:,:)
! accumulate (sum) the MPI distributed sum of number of boxes
call MPI_ALLREDUCE(sum_nbox, sum_nbox_all, 1, &
MPI_INTEGER, MPI_SUM, MPI_COMM_NEW, mympierror)
! Only on the first iteration does evolve2D (evolve0D) change the
! ionization fractions
if (niter == 1) then
! accumulate (max) MPI distributed xh_av
call MPI_ALLREDUCE(xh_av(:,:,:,1), buffer, mesh(1)*mesh(2)*mesh(3), &
!call MPI_ALLREDUCE(xh_av(:,:,:), buffer, mesh(1)*mesh(2)*mesh(3), &
MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_NEW, mympierror)
! Overwrite the processor local values with the accumulated value
xh_av(:,:,:,1) = buffer(:,:,:)
xh_av(:,:,:,0) = max(0.0_dp,min(1.0_dp,1.0-xh_av(:,:,:,1)))
!xh_av(:,:,:) = buffer(:,:,:)
! accumulate (max) MPI distributed xh_intermed
call MPI_ALLREDUCE(xh_intermed(:,:,:,1), buffer, &
!call MPI_ALLREDUCE(xh_intermed(:,:,:), buffer, &
mesh(1)*mesh(2)*mesh(3), MPI_DOUBLE_PRECISION, MPI_MAX, &
MPI_COMM_NEW, mympierror)
! Overwrite the processor local values with the accumulated value
xh_intermed(:,:,:,1)=buffer(:,:,:)
xh_intermed(:,:,:,0)=max(0.0_dp,min(1.0_dp,1.0-xh_intermed(:,:,:,1)))
!xh_intermed(:,:,:)=buffer(:,:,:)
endif
#else
photon_loss_all(:)=photon_loss(:)
sum_nbox_all=sum_nbox
#endif
end subroutine pass_all_sources
! ===========================================================================
subroutine global_pass (conv_flag,dt)
! Flag variable (passed back from evolve0D_global)
integer,intent(out) :: conv_flag
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
integer :: i,j,k ! mesh position
! Mesh position of the cell being treated
integer,dimension(Ndim) :: pos
! Report photon losses over grid boundary
if (rank == 0) write(logf,*) 'photon loss counter: ',photon_loss_all(:)
! Turn total photon loss into a mean per cell (used in evolve0d_global)
! GM/110225: Possibly this should be
! photon_loss_all(:)/(real(mesh(1))*real(mesh(2))
! Since this is the correct answer for an optically thin volume
! Compromise: photon_loss_all(:)/(real(mesh(1))**3)*sum(xh_av(:,:,:,1))
! then if the box is fully ionized the optically thin 1/N^2 is used
! and if the box is more neutral something close to 1/N^3 is used.
! Not sure
photon_loss(:)=photon_loss_all(:)/(real(mesh(1))*real(mesh(2))*real(mesh(3)))
! Report minimum value of xh_av(0) to check for zeros
if (rank == 0) then
write(logf,*) "min xh_av(0): ",minval(xh_av(:,:,:,0))
!write(logf,*) "min xh_av(0): ",minval(1.0_dp-xh_av(:,:,:))
endif
! Apply total photo-ionization rates from all sources (phih_grid)
conv_flag=0 ! will be used to check for convergence
! Loop through the entire mesh
if (rank == 0) write(logf,*) 'Doing global '
do k=1,mesh(3)
do j=1,mesh(2)
do i=1,mesh(1)
pos=(/ i,j,k /)
call evolve0D_global(dt,pos,conv_flag)
enddo
enddo
enddo
! Report on convergence and intermediate result
if (rank == 0) then
write(logf,*) "Number of non-converged points: ",conv_flag
write(logf,*) "Intermediate result for mean H ionization fraction: ", &
sum(xh_intermed(:,:,:,1))/real(mesh(1)*mesh(2)*mesh(3))
!sum(xh_intermed(:,:,:))/real(mesh(1)*mesh(2)*mesh(3))
endif
! Report on photon conservation
call calculate_photon_statistics (dt,xh_intermed,xh_av)
call report_photonstatistics (dt)
end subroutine global_pass
! ===========================================================================
!> Ray tracing the entire grid for all the sources using the
!! master-slave model for distributing the sources over the
!! MPI processes.
subroutine do_grid_master_slave (dt,niter)
! Ray tracing the entire grid for all the sources using the
! master-slave model for distributing the sources over the
! MPI processes.
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
integer,intent(in) :: niter !< interation counter, passed on to evolve0D
if (rank == 0) then
call do_grid_master ()
else
call do_grid_slave (dt,niter)
endif
end subroutine do_grid_master_slave
! ===========================================================================
!> The master task in the master-slave setup for distributing
!! the ray-tracing over the sources over the MPI processes.
subroutine do_grid_master ()
! The master task in the master-slave setup for distributing
! the ray-tracing over the sources over the MPI processes.
integer :: ns1
integer :: sources_done,whomax,who,answer
! counter for master-slave process
integer,dimension(:),allocatable :: counts
#ifdef MPI
integer :: mympierror
#endif
#ifdef MPI
! Source Loop - Master Slave with rank=0 as Master
sources_done = 0
ns1 = 0
! Allocate counter for master-slave process
if (.not.(allocated(counts))) allocate(counts(0:npr-1))
! send tasks to slaves
whomax = min(NumSrc,npr-1)
do who=1,whomax
if (ns1 <= NumSrc) then
ns1=ns1+1
call MPI_Send (ns1, 1, MPI_INTEGER, who, 1, &
MPI_COMM_NEW, mympierror)
endif
enddo
do while (sources_done < NumSrc)
! wait for an answer from a slave.
call MPI_Recv (answer, & ! address of receive buffer
1, & ! number of items to receive
MPI_INTEGER, & ! type of data
MPI_ANY_SOURCE, & ! can receive from any other
1, & ! tag
MPI_COMM_NEW, & ! communicator
mympi_status, & ! status
mympierror)
who = mympi_status(MPI_SOURCE) ! find out who sent us the answer
sources_done=sources_done+1 ! and the number of sources done
! put the slave on work again,
! but only if not all tasks have been sent.
! we use the value of num to detect this */
if (ns1 < NumSrc) then
ns1=ns1+1
call MPI_Send (ns1, 1, MPI_INTEGER, &
who, &
1, &
MPI_COMM_NEW, &
mympierror)
endif
enddo
! Now master sends a message to the slaves to signify that they
! should end the calculations. We use a special tag for that:
do who = 1,npr-1
call MPI_Send (0, 1, MPI_INTEGER, &
who, &
2, & ! tag
MPI_COMM_NEW, &
mympierror)
! the slave will send to master the number of calculations
! that have been performed.
! We put this number in the counts array.
call MPI_Recv (counts(who), & ! address of receive buffer
1, & ! number of items to receive
MPI_INTEGER, & ! type of data
who, & ! receive from process who
7, & ! tag
MPI_COMM_NEW, & ! communicator
mympi_status, & ! status
mympierror)
enddo
write(logf,*) 'Mean number of sources per processor: ', &
real(NumSrc)/real(npr-1)
write(logf,*) 'Counted mean number of sources per processor: ', &
real(sum(counts(1:npr-1)))/real(npr-1)
write(logf,*) 'Minimum and maximum number of sources ', &
'per processor: ', &
minval(counts(1:npr-1)),maxval(counts(1:npr-1))
flush(logf)
#endif
end subroutine do_grid_master
! ===========================================================================
!> The slave task in the master-slave setup for distributing
!! the ray-tracing over the sources over the MPI processes.
subroutine do_grid_slave(dt,niter)
! The slave task in the master-slave setup for distributing
! the ray-tracing over the sources over the MPI processes.
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
integer,intent(in) :: niter !< interation counter, passed on to evolve0D
integer :: local_count
integer :: ns1
#ifdef MPI
integer :: mympierror
#endif
#ifdef MPI
local_count=0
call MPI_Recv (ns1, & ! address of receive buffer
1, & ! number of items to receive
MPI_INTEGER, & ! type of data
0, & ! can receive from master only
MPI_ANY_TAG, & ! can expect two values, so
! we use the wildcard MPI_ANY_TAG
! here
MPI_COMM_NEW, & ! communicator
mympi_status, & ! status
mympierror)
! if tag equals 2, then skip the calculations
if (mympi_status(MPI_TAG) /= 2) then
do
#ifdef MPILOG
! Report
write(logf,*) 'Processor ',rank,' received: ',ns1
write(logf,*) ' that is source ',ns1 !SrcSeries(ns1)
write(logf,*) ' at:',srcpos(:,ns1)
flush(logf)
#endif
! Do the source at hand
call do_source(dt,ns1,niter)
! Update local counter
local_count=local_count+1
#ifdef MPILOG
! Report ionization fractions
write(logf,*) sum(xh_intermed(:,:,:,1))/ &
!write(logf,*) sum(xh_intermed(:,:,:))/ &
real(mesh(1)*mesh(2)*mesh(3))
write(logf,*) sum(xh_av(:,:,:,1))/real(mesh(1)*mesh(2)*mesh(3))
!write(logf,*) sum(xh_av(:,:,:))/real(mesh(1)*mesh(2)*mesh(3))
write(logf,*) local_count
#endif
! Send 'answer'
call MPI_Send (local_count, 1, & ! sending one int
MPI_INTEGER, 0, & ! to master
1, & ! tag
MPI_COMM_NEW, & ! communicator
mympierror)
! Receive new source number
call MPI_Recv (ns1, & ! address of receive buffer
1, & ! number of items to receive
MPI_INTEGER, & ! type of data
0, & ! can receive from master only
MPI_ANY_TAG, & ! can expect two values, so
! we use the wildcard MPI_ANY_TAG
! here
MPI_COMM_NEW, & ! communicator
mympi_status, & ! status
mympierror)
! leave this loop if tag equals 2
if (mympi_status(MPI_TAG) == 2) then
#ifdef MPILOG
write(logf,*) 'Stop received'
flush(logf)
#endif
exit
endif
enddo
endif
! this is the point that is reached when a task is received with
! tag = 2
! send the number of calculations to master and return
#ifdef MPILOG
! Report
write(logf,*) 'Processor ',rank,' did ',local_count,' sources'
flush(logf)
#endif
call MPI_Send (local_count, &
1, &
MPI_INTEGER, & ! sending one int
0, & ! to master
7, & ! tag
MPI_COMM_NEW,& ! communicator
mympierror)
#endif
end subroutine do_grid_slave
! ===========================================================================
!> Does the ray-tracing over the sources by distributing
!! the sources evenly over the available MPI processes-
subroutine do_grid_static (dt,niter)
! Does the ray-tracing over the sources by distributing
! the sources evenly over the available MPI processes-
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
integer,intent(in) :: niter !< interation counter, passed on to evolve0D
integer :: ns1
! Source Loop - distributed for the MPI nodes
do ns1=1+rank,NumSrc,npr
call do_source(dt,ns1,niter)
enddo
end subroutine do_grid_static
! ===========================================================================
!> Does the ray-tracing over the entire 3D grid for one source.
!! The number of this source in the current list is ns1.
subroutine do_source(dt,ns1,niter)
! Does the ray-tracing over the entire 3D grid for one source.
! The number of this source in the current list is ns1.
real(kind=dp),intent(in) :: dt !< time step, passed on to evolve0D
integer, intent(in) :: ns1 !< number of the source being done
integer,intent(in) :: niter !< interation counter, passed on to evolve0D
integer :: naxis,nplane,nquadrant
integer :: ns
integer :: k
integer :: nbox
integer ::nnt
! Mesh position of the cell being treated
integer,dimension(Ndim) :: rtpos
! Pick up source number from the source list
!ns=SrcSeries(ns1)
ns=ns1
! reset column densities for new source point
! coldensh_out is unique for each source point
coldensh_out(:,:,:)=0.0
! Find the mesh position for the end points of the loop
! We trace until we reach max_subbox (set in c2ray_parameters)
! or the end of the grid. In the periodic case the end of the
! grid is always mesh/2 away from the source. If the grid is
! even-sized we trave mesh/2 cells to the left and mesh/2-1
! cell to the right. If it is odd, it is mesh/2 in either direction.
! The mod(mesh,2) takes care of handling this.
if (periodic_bc) then
lastpos_r(:)=srcpos(:,ns)+min(max_subbox,mesh(:)/2-1+mod(mesh(:),2))
lastpos_l(:)=srcpos(:,ns)-min(max_subbox,mesh(:)/2)
else
lastpos_r(:)=min(srcpos(:,ns)+max_subbox,mesh(:))
lastpos_l(:)=max(srcpos(:,ns)-max_subbox,1)
endif
! Loop through grid in the order required by
! short characteristics
! Transfer is done in a set of cubes of increasing size.
! If the HII region is small we do not waste time calculating
! column densities of parts of the grid where no radiation
! penetrates. To test whether the current subbox is large
! enough we use the photon_loss_src. If this is non-zero,
! photons are leaving this subbox and we need to do another
! one. We also stop once we have done the whole grid.
nbox=0 ! subbox counter
photon_loss_src(:)=NormFlux(ns)*S_star_nominal !-1.0 ! to pass the first while test
last_r(:)=srcpos(:,ns) ! to pass the first while test
last_l(:)=srcpos(:,ns) ! to pass the first while test
! Loop through boxes of increasing size
! NOTE: make this limit on the photon_loss a fraction of
! a source flux loss_fraction*NormFlux(ns)*S_star_nominal)
do while (all(photon_loss_src(:) > 1e-6*NormFlux(ns)*S_star_nominal) &
!do while (all(photon_loss_src(:) /= 0.0) &
.and. last_r(3) < lastpos_r(3) &
.and. last_l(3) > lastpos_l(3))
nbox=nbox+1 ! increase subbox counter
photon_loss_src(:) = 0.0 ! reset photon_loss_src to zero
photon_loss_src_thread(:,:) = 0.0 ! reset photon_loss_src to zero
last_r(:)=min(srcpos(:,ns)+subboxsize*nbox,lastpos_r(:))
last_l(:)=max(srcpos(:,ns)-subboxsize*nbox,lastpos_l(:))
! OpenMP: if we have multiple OpenMP threads (nthreads > 1) we
! parallelize over the threads by doing independent parts of
! the mesh.
if (nthreads > 1) then ! OpenMP parallelization
! First do source point (on first pass)
if (nbox == 1) then
rtpos(:)=srcpos(:,ns)
call evolve0D(dt,rtpos,ns,niter)
endif
! do independent areas of the mesh in parallel using OpenMP
!$omp parallel default(shared) private(tn)
!!!reduction(+:photon_loss_src)
! Find out your thread number
#ifdef MY_OPENMP
tn=omp_get_thread_num()+1
#else
tn=1
#endif
! Then do the the axes
!$omp do schedule(dynamic,1)
do naxis=1,6
call evolve1D_axis(dt,ns,niter,naxis)
enddo
!$omp end do
! Then the source planes
!$omp do schedule (dynamic,1)
do nplane=1,12
call evolve2D_plane(dt,ns,niter,nplane)
end do
!$omp end do
! Then the quadrants
!$omp do schedule (dynamic,1)
do nquadrant=1,8
call evolve3D_quadrant(dt,ns,niter,nquadrant)
end do
!$omp end do
!$omp end parallel
! Collect photon losses for each thread
do nnt=1,nthreads
photon_loss_src(:)=photon_loss_src(:) + &
photon_loss_src_thread(:,nnt)
enddo
else ! No OpenMP parallelization
! 1. transfer in the upper part of the grid
! (srcpos(3)-plane and above)
do k=srcpos(3,ns),last_r(3)
rtpos(3)=k
call evolve2D(dt,rtpos,ns,niter)
end do
! 2. transfer in the lower part of the grid (below srcpos(3))
do k=srcpos(3,ns)-1,last_l(3),-1
rtpos(3)=k
call evolve2D(dt,rtpos,ns,niter)
end do
! No OpenMP threads so we use position 1
photon_loss_src(:)=photon_loss_src_thread(:,1)
endif
enddo
! Record the final photon loss, this is the photon loss that leaves
! the grid.
photon_loss(:)=photon_loss(:) + photon_loss_src(:)
! Sum the total number of subboxes used for reporting later
sum_nbox=sum_nbox+nbox
end subroutine do_source
! ===========================================================================
!> Traverse a z-plane (z=rtpos(3)) by sweeping in the x and y
!! directions.
subroutine evolve2D(dt,rtpos,ns,niter)
! Traverse a z-plane (z=rtpos(3)) by sweeping in the x and y
! directions.
real(kind=dp),intent(in) :: dt !! passed on to evolve0D
integer,dimension(Ndim),intent(inout) :: rtpos !< mesh position, pos(3) is
!! intent(in)
integer,intent(in) :: ns !< current source
integer,intent(in) :: niter !< passed on to evolve0D
integer :: i,j ! mesh positions
! sweep in `positive' j direction
do j=srcpos(2,ns),last_r(2)
rtpos(2)=j
do i=srcpos(1,ns),last_r(1)
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) ! `positive' i
end do
do i=srcpos(1,ns)-1,last_l(1),-1
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) ! `negative' i
end do
end do
! sweep in `negative' j direction
do j=srcpos(2,ns)-1,last_l(2),-1
rtpos(2)=j
do i=srcpos(1,ns),last_r(1)
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) ! `positive' i
end do
do i=srcpos(1,ns)-1,last_l(1),-1
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) ! `negative' i
end do
end do
end subroutine evolve2D
! ===========================================================================
! Ray tracing for the axes going through the source point
! should be called after having done the source point
subroutine evolve1D_axis(dt,ns,niter,naxis)
real(kind=dp),intent(in) :: dt ! passed on to evolve0D
integer,intent(in) :: ns ! current source
integer,intent(in) :: niter ! passed on to evolve0D
integer,intent(in) :: naxis ! axis to do
integer :: i,j,k
integer,dimension(Ndim) :: rtpos ! mesh position
select case (naxis)
case(1)
! sweep in +i direction
rtpos(2:3)=srcpos(2:3,ns)
do i=srcpos(1,ns)+1,last_r(1)
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) !# `positive' i
enddo
case(2)
! sweep in -i direction
rtpos(2:3)=srcpos(2:3,ns)
do i=srcpos(1,ns)-1,last_l(1),-1
rtpos(1)=i
call evolve0D(dt,rtpos,ns,niter) !# `negative' i
end do
case(3)
! sweep in +j direction
rtpos(1)=srcpos(1,ns)
rtpos(3)=srcpos(3,ns)
do j=srcpos(2,ns)+1,last_r(2)
rtpos(2)=j
call evolve0D(dt,rtpos,ns,niter) !# `positive' j
end do
case(4)
! sweep in -j direction
rtpos(1)=srcpos(1,ns)
rtpos(3)=srcpos(3,ns)
do j=srcpos(2,ns)-1,last_l(2),-1
rtpos(2)=j
call evolve0D(dt,rtpos,ns,niter) !# `negative' j
end do
case(5)
! sweep in +k direction
rtpos(1:2)=srcpos(1:2,ns)
do k=srcpos(3,ns)+1,last_r(3)
rtpos(3)=k
call evolve0D(dt,rtpos,ns,niter) !# `positive' k
end do
case(6)
! sweep in -k direction
rtpos(1:2)=srcpos(1:2,ns)
do k=srcpos(3,ns)-1,last_l(3),-1
rtpos(3)=k
call evolve0D(dt,rtpos,ns,niter) !# `negative' k
end do
end select
end subroutine evolve1D_axis
! ===========================================================================
!> Ray tracing for planes containing the source point
!! should be called after evolve1D_axis
subroutine evolve2D_plane(dt,ns,niter,nplane)
! find column density for the axes going through the source point
! should be called after having done the source point
real(kind=dp),intent(in) :: dt ! passed on to evolve0D
integer,intent(in) :: ns ! current source
integer,intent(in) :: niter ! passed on to evolve0D
integer,intent(in) :: nplane ! plane to do
integer :: i,j,k
integer,dimension(Ndim) :: rtpos ! mesh position
select case (nplane)
case(1)
! sweep in +i,+j direction
rtpos(3)=srcpos(3,ns)