11! Copyright (c), The Regents of the University of California
22! Terms of use are as specified in LICENSE.tx
3+
4+ #include " assert_macros.h"
5+
36submodule(distribution_m) distribution_s
7+ use assert_m
48 use intrinsic_array_m, only : intrinsic_array_t
5- use do_concurrent_m, only : do_concurrent_sampled_speeds, do_concurrent_my_velocities
6- use assert_m, only : assert
79 implicit none
810
911contains
@@ -18,54 +20,57 @@ pure function monotonically_increasing(f) result(monotonic)
1820 module procedure construct
1921 integer i
2022
21- call assert(all (sample_distribution(:,2 )>= 0.D0 ), " distribution_t%construct: sample_distribution>=0." , &
22- intrinsic_array_t(sample_distribution))
23+ call_assert_diagnose(all (sample_distribution(:,2 )>= 0.D0 ), " distribution_t%construct: sample_distribution>=0." , intrinsic_array_t(sample_distribution))
2324
2425 associate(nintervals = > size (sample_distribution,1 ))
2526 distribution% vel_ = [(sample_distribution(i,1 ), i = 1 , nintervals)] ! Assign speeds to each distribution bin
2627 distribution% cumulative_distribution_ = [0.D0 , [(sum (sample_distribution(1 :i,2 )), i= 1 , nintervals)]]
27-
28- call assert(monotonically_increasing(distribution% cumulative_distribution_), &
29- " distribution_t: monotonically_increasing(distribution%cumulative_distribution_)" , &
30- intrinsic_array_t(distribution% cumulative_distribution_))
3128 end associate
3229
30+ call_assert_diagnose(monotonically_increasing(distribution% cumulative_distribution_), " distribution_t: monotonically increasing cumulative distribution" , intrinsic_array_t(distribution% cumulative_distribution_))
31+
3332 end procedure construct
3433
3534 module procedure cumulative_distribution
36- call assert(allocated (self% cumulative_distribution_), &
37- " distribution_t%cumulative_distribution: allocated(cumulative_distribution_)" )
35+ call_assert(allocated (self% cumulative_distribution_))
3836 my_cumulative_distribution = self% cumulative_distribution_
3937 end procedure
4038
4139 module procedure velocities
4240
4341 double precision , allocatable :: sampled_speeds(:,:), dir(:,:,:)
44- integer step
42+ integer cell, step, k
4543
46- call assert(allocated (self% cumulative_distribution_), &
47- " distribution_t%cumulative_distribution: allocated(cumulative_distribution_)" )
48- call assert(allocated (self% vel_), " distribution_t%cumulative_distribution: allocated(vel_)" )
44+ call_assert(allocated (self% cumulative_distribution_))
45+ call_assert(allocated (self% vel_))
4946
50- ! Sample from the distribution
51- call do_concurrent_sampled_speeds(speeds, self% vel_, self% cumulative_distribution(), sampled_speeds)
52-
53- associate(nsteps = > size (speeds,2 ))
47+ ! Sample from the distribution
48+ associate(ncells = > size (speeds,1 ), nsteps = > size (speeds,2 ))
49+ allocate (sampled_speeds(ncells,nsteps))
50+ do concurrent(cell = 1 :ncells, step = 1 :nsteps)
51+ k = findloc(speeds(cell,step) >= self% cumulative_distribution(), value= .false. , dim= 1 )- 1
52+ sampled_speeds(cell,step) = self% vel_(k)
53+ end do
54+
55+ ! Create unit vectors
56+ dir = directions(:,1 :nsteps,:)
5457
55- ! Create unit vectors
56- dir = directions(:,1 :nsteps,:)
58+ associate(dir_mag = > sqrt (dir(:,:,1 )** 2 + dir(:,:,2 )** 2 + dir(:,:,3 )** 2 ))
59+ associate(dir_mag_ = > merge (dir_mag, epsilon (dir_mag), dir_mag/= 0 .))
60+ dir(:,:,1 ) = dir(:,:,1 )/ dir_mag_
61+ dir(:,:,2 ) = dir(:,:,2 )/ dir_mag_
62+ dir(:,:,3 ) = dir(:,:,3 )/ dir_mag_
63+ end associate
64+ end associate
5765
58- associate(dir_mag = > sqrt (dir(:,:,1 )** 2 + dir(:,:,2 )** 2 + dir(:,:,3 )** 2 ))
59- associate(dir_mag_ = > merge (dir_mag, epsilon (dir_mag), dir_mag/= 0 .))
60- dir(:,:,1 ) = dir(:,:,1 )/ dir_mag_
61- dir(:,:,2 ) = dir(:,:,2 )/ dir_mag_
62- dir(:,:,3 ) = dir(:,:,3 )/ dir_mag_
63- end associate
64- end associate
65-
66- call do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities)
67-
68- end associate
66+ allocate (my_velocities, mold= dir)
67+
68+ do concurrent(step= 1 :nsteps)
69+ my_velocities(:,step,1 ) = sampled_speeds(:,step)* dir(:,step,1 )
70+ my_velocities(:,step,2 ) = sampled_speeds(:,step)* dir(:,step,2 )
71+ my_velocities(:,step,3 ) = sampled_speeds(:,step)* dir(:,step,3 )
72+ end do
73+ end associate
6974
7075 end procedure
7176
0 commit comments