From 237f2f3e014942200dd9303345cd6b11e93ae5cc Mon Sep 17 00:00:00 2001 From: davytorres Date: Mon, 24 Oct 2022 00:32:53 -0600 Subject: [PATCH] Added angle distribution --- src/matcha/distribution_m.f90 | 24 ++++-- src/matcha/distribution_s.F90 | 129 ++++++++++++++++++++++++++++----- src/matcha/do_concurrent_m.f90 | 8 +- src/matcha/do_concurrent_s.f90 | 10 +-- src/matcha/input_m.f90 | 19 ++++- src/matcha/input_s.f90 | 46 ++++++++++++ src/matcha_s.F90 | 29 ++++---- 7 files changed, 215 insertions(+), 50 deletions(-) diff --git a/src/matcha/distribution_m.f90 b/src/matcha/distribution_m.f90 index 5c132fe..b7e6041 100644 --- a/src/matcha/distribution_m.f90 +++ b/src/matcha/distribution_m.f90 @@ -8,20 +8,21 @@ module distribution_m type distribution_t private - double precision, allocatable, dimension(:) :: vel_, cumulative_distribution_ + double precision, allocatable, dimension(:) :: vel_,cumulative_distribution_,angle_,cumulative_angle_distribution_ contains procedure :: cumulative_distribution + procedure :: cumulative_angle_distribution procedure :: velocities end type interface distribution_t - pure module function construct(sample_distribution) result(distribution) + pure module function construct(sample_distribution,sample_angle_distribution) result(distribution) implicit none - double precision, intent(in) :: sample_distribution(:,:) + double precision, intent(in) :: sample_distribution(:,:),sample_angle_distribution(:,:) type(distribution_t) distribution - end function - + end function + end interface interface @@ -30,13 +31,20 @@ pure module function cumulative_distribution(self) result(my_cumulative_distribu implicit none class(distribution_t), intent(in) :: self double precision, allocatable :: my_cumulative_distribution(:) - end function + end function cumulative_distribution + + pure module function cumulative_angle_distribution(self) result(my_cumulative_angle_distribution) + implicit none + class(distribution_t), intent(in) :: self + double precision, allocatable :: my_cumulative_angle_distribution(:) + end function cumulative_angle_distribution + - pure module function velocities(self, speeds, directions) result(my_velocities) + pure module function velocities(self, speeds, angles, directions) result(my_velocities) !! Return the t_cell_collection_t object's velocity vectors implicit none class(distribution_t), intent(in) :: self - double precision, intent(in) :: speeds(:,:), directions(:,:,:) + double precision, intent(in) :: speeds(:,:),angles(:,:),directions(:,:,:) double precision, allocatable :: my_velocities(:,:,:) end function velocities diff --git a/src/matcha/distribution_s.F90 b/src/matcha/distribution_s.F90 index 93b25d5..c63e78e 100644 --- a/src/matcha/distribution_s.F90 +++ b/src/matcha/distribution_s.F90 @@ -2,7 +2,7 @@ ! Terms of use are as specified in LICENSE.tx submodule(distribution_m) distribution_s use intrinsic_array_m, only : intrinsic_array_t - use do_concurrent_m, only : do_concurrent_sampled_speeds, do_concurrent_my_velocities + use do_concurrent_m, only : do_concurrent_sample, do_concurrent_my_velocities #ifdef USE_CAFFEINE use caffeine_assert_m, only : assert @@ -35,44 +35,137 @@ pure function monotonically_increasing(f) result(monotonic) "distribution_t: monotonically_increasing(distribution%cumulative_distribution_)", & intrinsic_array_t(distribution%cumulative_distribution_)) end associate + + call assert(all(sample_angle_distribution(:,2)>=0.D0), "distribution_t%construct: sample_distribution>=0.", & + intrinsic_array_t(sample_angle_distribution)) + + associate(nintervals => size(sample_angle_distribution,1)) + distribution%angle_ = [(sample_angle_distribution(i,1), i =1, nintervals)] ! Assign speeds to each distribution bin + distribution%cumulative_angle_distribution_ = [0.D0, [(sum(sample_angle_distribution(1:i,2)), i=1, nintervals)]] + + call assert(monotonically_increasing(distribution%cumulative_angle_distribution_), & + "distribution_t: monotonically_increasing(distribution%cumulative_distribution_)", & + intrinsic_array_t(distribution%cumulative_angle_distribution_)) + end associate + end procedure construct module procedure cumulative_distribution call assert(allocated(self%cumulative_distribution_), & "distribution_t%cumulative_distribution: allocated(cumulative_distribution_)") my_cumulative_distribution = self%cumulative_distribution_ + end procedure + + module procedure cumulative_angle_distribution + call assert(allocated(self%cumulative_angle_distribution_), & + "distribution_t%cumulative_angle_distribution: allocated(cumulative_angle_distribution_)") + my_cumulative_angle_distribution = self%cumulative_angle_distribution_ end procedure + module procedure velocities - double precision, allocatable :: sampled_speeds(:,:), dir(:,:,:) - integer step + double precision, allocatable :: sampled_speeds(:,:), sampled_angles(:,:), dir(:,:,:) + integer cell,step + double precision dir_mag_ + + double precision pi,twopi + double precision a,b,c + double precision x1,y1,z1,x2,y2,z2,x3,y3,z3,vec1mag + double precision vxp,vyp,vzp,eps,vec3mag + double precision random_phi call assert(allocated(self%cumulative_distribution_), & "distribution_t%cumulative_distribution: allocated(cumulative_distribution_)") call assert(allocated(self%vel_), "distribution_t%cumulative_distribution: allocated(vel_)") + call assert(allocated(self%cumulative_angle_distribution_), & + "distribution_t%cumulative_angle_distribution: allocated(cumulative_angle_distribution_)") + call assert(allocated(self%angle_), "distribution_t%cumulative_angle_distribution: allocated(angle_)") + ! Sample from the distribution - call do_concurrent_sampled_speeds(speeds, self%vel_, self%cumulative_distribution(), sampled_speeds) + call do_concurrent_sample(speeds, self%vel_, self%cumulative_distribution(), sampled_speeds) + call do_concurrent_sample(angles, self%angle_, self%cumulative_angle_distribution(), sampled_angles) - associate(nsteps => size(speeds,2)) - - ! Create unit vectors - dir = directions(:,1:nsteps,:) - - associate(dir_mag => sqrt(dir(:,:,1)**2 +dir(:,:,2)**2 + dir(:,:,3)**2)) - associate(dir_mag_ => merge(dir_mag, epsilon(dir_mag), dir_mag/=0.)) - dir(:,:,1) = dir(:,:,1)/dir_mag_ - dir(:,:,2) = dir(:,:,2)/dir_mag_ - dir(:,:,3) = dir(:,:,3)/dir_mag_ - end associate - end associate - - call do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities) + associate(ncells => size(speeds,1), nsteps => size(speeds,2)) + +! ! Create unit vectors +! dir = directions(:,1:nsteps,:) +! +! associate(dir_mag => sqrt(dir(:,:,1)**2 +dir(:,:,2)**2 + dir(:,:,3)**2)) +! associate(dir_mag_ => merge(dir_mag, epsilon(dir_mag), dir_mag/=0.)) +! dir(:,:,1) = dir(:,:,1)/dir_mag_ +! dir(:,:,2) = dir(:,:,2)/dir_mag_ +! dir(:,:,3) = dir(:,:,3)/dir_mag_ +! end associate +! end associate +! call do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities) + + allocate(my_velocities(ncells,nsteps,3)) + + pi = acos(-1.d0) + twopi = 2.d0*pi + eps = 1.d-12 + +! do cell = 1,ncells + do concurrent(cell = 1:ncells) + do step = 1,nsteps + random_phi = 2.d0*pi*directions(cell,step,1) + vxp = sampled_speeds(cell,step)*sin(sampled_angles(cell,step))*cos(random_phi) + vyp = sampled_speeds(cell,step)*sin(sampled_angles(cell,step))*sin(random_phi) + vzp = sampled_speeds(cell,step)*cos(sampled_angles(cell,step)) + if (step .eq. 1) then + my_velocities(cell,step,1) = vxp + my_velocities(cell,step,2) = vyp + my_velocities(cell,step,3) = vzp + else + x3 = my_velocities(cell,step-1,1) + y3 = my_velocities(cell,step-1,2) + z3 = my_velocities(cell,step-1,3) + vec3mag = sqrt(x3**2 + y3**2 + z3**2)+eps + x3 = x3/vec3mag + y3 = y3/vec3mag + z3 = z3/vec3mag + + if (abs(z3) .gt. eps) then + a = 0.d0 + b = 1.d0 + c = 0.d0 + else + a = 0.d0 + b = 0.d0 + c = 1.d0 + end if + + !curl(x3,y3,z3,a,b,c,x1,y1,z1) + x1 = y3*c - z3*b + y1 = z3*a - x3*c + z1 = x3*b - y3*a + vec1mag = sqrt(x1**2 + y1**2 + z1**2) + eps + x1 = x1/vec1mag + y1 = y1/vec1mag + z1 = z1/vec1mag + + ! curl(x3,y3,z3,x1,y1,z1,x2,y2,z2) + x2 = y3*z1 - z3*y1 + y2 = z3*x1 - x3*z1 + z2 = x3*y1 - y3*x1 + + + my_velocities(cell,step,1) = vxp*x1 + vyp*x2 + vzp*x3 + my_velocities(cell,step,2) = vxp*y1 + vyp*y2 + vzp*y3 + my_velocities(cell,step,3) = vxp*z1 + vyp*z2 + vzp*z3 + + end if + + end do + end do + end associate end procedure end submodule distribution_s + diff --git a/src/matcha/do_concurrent_m.f90 b/src/matcha/do_concurrent_m.f90 index 21869c1..6c9fdfe 100644 --- a/src/matcha/do_concurrent_m.f90 +++ b/src/matcha/do_concurrent_m.f90 @@ -3,15 +3,15 @@ module do_concurrent_m use t_cell_collection_m, only : t_cell_collection_t, t_cell_collection_bind_C_t implicit none private - public :: do_concurrent_sampled_speeds, do_concurrent_my_velocities, do_concurrent_k, do_concurrent_speeds + public :: do_concurrent_sample, do_concurrent_my_velocities, do_concurrent_k, do_concurrent_speeds public :: do_concurrent_output_distribution interface - pure module subroutine do_concurrent_sampled_speeds(speeds, vel, cumulative_distribution, sampled_speeds) bind(C) + pure module subroutine do_concurrent_sample(rvalues, val, cumulative_distribution, sampled_values) bind(C) implicit none - real(c_double), intent(in) :: speeds(:,:), vel(:), cumulative_distribution(:) - real(c_double), intent(out), allocatable :: sampled_speeds(:,:) + real(c_double), intent(in) :: rvalues(:,:), val(:), cumulative_distribution(:) + real(c_double), intent(out), allocatable :: sampled_values(:,:) end subroutine pure module subroutine do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities) bind(C) diff --git a/src/matcha/do_concurrent_s.f90 b/src/matcha/do_concurrent_s.f90 index 80c37e7..5fb3224 100644 --- a/src/matcha/do_concurrent_s.f90 +++ b/src/matcha/do_concurrent_s.f90 @@ -4,15 +4,15 @@ contains - module procedure do_concurrent_sampled_speeds + module procedure do_concurrent_sample integer cell, step - associate(ncells => size(speeds,1), nsteps => size(speeds,2)) - allocate(sampled_speeds(ncells,nsteps)) + associate(ncells => size(rvalues,1), nsteps => size(rvalues,2)) + allocate(sampled_values(ncells,nsteps)) do concurrent(cell = 1:ncells, step = 1:nsteps) - associate(k => findloc(speeds(cell,step) >= cumulative_distribution, value=.false., dim=1)-1) - sampled_speeds(cell,step) = vel(k) + associate(k => findloc(rvalues(cell,step) >= cumulative_distribution, value=.false., dim=1)-1) + sampled_values(cell,step) = val(k) end associate end do end associate diff --git a/src/matcha/input_m.f90 b/src/matcha/input_m.f90 index 9d4a7c9..d8233e6 100644 --- a/src/matcha/input_m.f90 +++ b/src/matcha/input_m.f90 @@ -8,7 +8,7 @@ module input_m type input_t private - integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4 + integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4, num_angle_intervals_ = 4 double precision :: time_step_ = 0.1D0 !double precision, allocatable :: sample_distribution_(:,:) !allocate(sample_distribution_(num_intervals_,2)) @@ -17,8 +17,10 @@ module input_m procedure :: num_positions procedure :: num_dimensions procedure :: num_intervals + procedure :: num_angle_intervals procedure :: time_step procedure :: sample_distribution + procedure :: sample_angle_distribution end type interface @@ -46,6 +48,12 @@ pure module function num_intervals(self) result(n) class(input_t), intent(in) :: self integer n end function + + pure module function num_angle_intervals(self) result(n) + implicit none + class(input_t), intent(in) :: self + integer n + end function pure module function time_step(self) result(dt) implicit none @@ -57,7 +65,14 @@ pure module function sample_distribution(self) result(empirical_distribution) implicit none class(input_t), intent(in) :: self double precision, allocatable :: empirical_distribution(:,:) - end function sample_distribution + end function sample_distribution + + pure module function sample_angle_distribution(self) result(empirical_angle_distribution) + implicit none + class(input_t), intent(in) :: self + double precision, allocatable :: empirical_angle_distribution(:,:) + end function sample_angle_distribution + end interface diff --git a/src/matcha/input_s.f90 b/src/matcha/input_s.f90 index 3492f24..7115145 100644 --- a/src/matcha/input_s.f90 +++ b/src/matcha/input_s.f90 @@ -21,6 +21,10 @@ n = self%num_intervals_ end procedure + module procedure num_angle_intervals + n = self%num_angle_intervals_ + end procedure + module procedure time_step dt = self%time_step_ end procedure @@ -64,5 +68,47 @@ end do end procedure + + module procedure sample_angle_distribution + integer i,nintervals + double precision angle_lower,angle_upper + double precision range,pi,dangle,sumy + double precision angle_lower_bin + double precision angle_upper_bin + double precision angle_middle_bin + double precision pi + double precision, allocatable :: angles(:),probability(:) + nintervals = self%num_angle_intervals_ + + pi = acos(-1.d0) + angle_lower = 0.d0 + angle_upper = pi + range = angle_upper - angle_lower + pi = acos(-1.d0) + dangle = range/dble(nintervals) + allocate(angles(nintervals),probability(nintervals)) +! Create normal distribution + sumy = 0.d0 + do i = 1,nintervals + angle_lower_bin = angle_lower + dble(i-1)*dangle + angle_upper_bin = angle_lower + dble(i)*dangle + angle_middle_bin = 0.5d0*(angle_lower_bin + angle_upper_bin) + angles(i) = angle_middle_bin + probability(i) = 1.d0/dble(nintervals) ! Use uniform distribution + sumy = sumy + probability(i) + end do + + do i = 1,nintervals + probability(i) = probability(i)/sumy + end do + + allocate(empirical_angle_distribution(nintervals,2)) + + do i = 1,nintervals + empirical_angle_distribution(i,1) = angles(i) + empirical_angle_distribution(i,2) = probability(i) + end do + + end procedure end submodule input_s diff --git a/src/matcha_s.F90 b/src/matcha_s.F90 index 27d1c4f..ae6de52 100644 --- a/src/matcha_s.F90 +++ b/src/matcha_s.F90 @@ -21,14 +21,16 @@ ndim => input%num_dimensions(), & nintervals => input%num_intervals(), & dt => input%time_step(), & - empirical_distribution => input%sample_distribution() & + empirical_distribution => input%sample_distribution(), & + empirical_angle_distribution => input%sample_angle_distribution() & ) block double precision, parameter :: scale = 100.D0 double precision, allocatable :: random_positions(:,:), random_4vectors(:,:,:) type(distribution_t) distribution - integer, parameter :: nveldim = 4 + type(distribution_t) distribution_angle + integer, parameter :: nveldim = 5 integer step type(data_partition_t) data_partition @@ -45,17 +47,18 @@ associate(nsteps => npositions -1) allocate(random_4vectors(my_num_cells,nsteps,nveldim)) call random_number(random_4vectors) - distribution = distribution_t(empirical_distribution) - - associate(random_speeds => random_4vectors(:,:,1), random_directions => random_4vectors(:,:,2:4)) - associate(v => distribution%velocities(random_speeds, random_directions)) - allocate(history(nsteps)) - history(1) = t_cell_collection_t(scale*random_positions, time=0.D0) - do step = 2, nsteps - associate(x => history(step-1)%positions(), t => history(step-1)%time()) - history(step) = t_cell_collection_t(x + v(:,step-1,:)*dt, t + dt) - end associate - end do + distribution = distribution_t(empirical_distribution,empirical_angle_distribution) + associate(random_speeds=>random_4vectors(:,:,1),random_directions=>random_4vectors(:,:,3:5)) + associate(random_angles => random_4vectors(:,:,2)) + associate(v => distribution%velocities(random_speeds,random_angles,random_directions)) + allocate(history(nsteps)) + history(1) = t_cell_collection_t(scale*random_positions, time=0.D0) + do step = 2, nsteps + associate(x => history(step-1)%positions(), t => history(step-1)%time()) + history(step) = t_cell_collection_t(x + v(:,step-1,:)*dt, t + dt) + end associate + end do + end associate end associate end associate end associate