Added diffusion and angle distribution

src/matcha/distribution_m.f90

@@ -8,17 +8,18 @@ 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
 
@@ -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
 

src/matcha/distribution_s.F90

@@ -2,7 +2,8 @@
 ! 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_sampled_speeds, do_concurrent_my_velocities
+  use do_concurrent_m, only : do_concurrent_sample, do_concurrent_my_velocities  
   use assert_m, only : assert
   implicit none
 
@@ -30,43 +31,156 @@ pure function monotonically_increasing(f) result(monotonic)
         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
+    double precision dot,adot,magv
+
 
     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)
-
-     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)
-
+     !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)
+
+     !print*,'sampled_angles = ',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)
+!       
+!     end associate
+
+       associate(ncells => size(speeds,1), nsteps => size(speeds,2))
+
+       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
+
+               !magv = dsqrt(my_velocities(cell,step,1)**2 + &
+               !             my_velocities(cell,step,2)**2 + &
+               !             my_velocities(cell,step,3)**2)
+               !dot = (my_velocities(cell,step,1)*x3 + &
+               !       my_velocities(cell,step,2)*y3 + &
+               !       my_velocities(cell,step,3)*z3)/magv
+               !adot = acos(dot)
+               !print*,'adot = ',adot,sampled_angles(cell,step)
+
+            end if
+
+         end do
+       end do
+
      end associate
 
+
+
   end procedure
 
 end submodule distribution_s

src/matcha/do_concurrent_m.f90

@@ -3,10 +3,18 @@ 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_sampled_speeds, do_concurrent_my_velocities, do_concurrent_k
+  public :: do_concurrent_angles, do_concurrent_speeds
   public :: do_concurrent_output_distribution
 
   interface
+
+    pure module subroutine do_concurrent_sample(rvalues, val, cumulative_distribution, sampled_values) bind(C)
+      implicit none
+      real(c_double), intent(in) :: rvalues(:,:), val(:), cumulative_distribution(:)
+      real(c_double), intent(out), allocatable :: sampled_values(:,:)
+    end subroutine
+
 
     pure module subroutine do_concurrent_sampled_speeds(speeds, vel, cumulative_distribution, sampled_speeds) bind(C)
       implicit none
@@ -40,7 +48,14 @@ module subroutine do_concurrent_speeds(history, speeds) bind(C)
       type(t_cell_collection_bind_C_t), intent(in) :: history(:)
       real(c_double), intent(out), allocatable :: speeds(:)
     end subroutine
-
+
+    module subroutine do_concurrent_angles(history, angles) bind(C)
+      implicit none
+      type(t_cell_collection_bind_C_t), intent(in) :: history(:)
+      real(c_double), intent(out), allocatable :: angles(:)
+    end subroutine
+
+
   end interface
 
 end module do_concurrent_m

src/matcha/do_concurrent_s.f90

@@ -3,7 +3,22 @@
   implicit none
 
 contains
-
+
+ module procedure do_concurrent_sample
+
+    integer cell, step
+
+    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(rvalues(cell,step) >= cumulative_distribution, value=.false., dim=1)-1)
+          sampled_values(cell,step) = val(k)
+        end associate
+      end do
+    end associate
+
+  end procedure
+
   module procedure do_concurrent_sampled_speeds
 
     integer cell, step
@@ -68,14 +83,14 @@
 
     real(c_double), pointer :: positions(:,:)
     real(c_double), allocatable :: x(:,:,:)
-  
+
     associate(npositions => size(history), ncells => history(1)%positions_shape(1))
 
       allocate(x(npositions,ncells,nspacedims))
 
       do i=1,npositions
-         call c_f_pointer(history(i)%positions_ptr, positions, history(1)%positions_shape)
-         x(i,:,:) = positions
+        call c_f_pointer(history(i)%positions_ptr, positions, history(1)%positions_shape)
+        x(i,:,:) = positions
       end do
 
       associate(t => history%time)
@@ -89,8 +104,46 @@
           end associate
         end do
       end associate
+
     end associate
 
   end procedure
+
+  module procedure do_concurrent_angles
+
+    integer i, j, k
+    integer, parameter :: nspacedims=3
+
+    real(c_double), pointer :: positions(:,:)
+    real(c_double), allocatable :: x(:,:,:)
+    real(c_double) v1mag,v2mag,eps
+
+    eps = 1.d-16
+    associate(npositions => size(history), ncells => history(1)%positions_shape(1))
+
+      allocate(x(npositions,ncells,nspacedims))
+
+      do i=1,npositions
+        call c_f_pointer(history(i)%positions_ptr, positions, history(1)%positions_shape)
+        x(i,:,:) = positions
+      end do
+
+      allocate(angles(ncells*(npositions-2)) )      
+      do concurrent(i = 1:npositions-2, j = 1:ncells)
+        associate( &
+          v1 => x(i+1,j,:) - x(i,j,:), &
+          v2 => x(i+2,j,:) - x(i+1,j,:), &
+          ij => i + (j-1)*(npositions-2) &
+         )
+          v1mag = sqrt(sum([(v1(k)**2, k=1,nspacedims)]))
+          v2mag = sqrt(sum([(v2(k)**2, k=1,nspacedims)]))
+          angles(ij) = acos( (v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3))/(v1mag*v2mag+eps) )
+        end associate
+      end do
+
+    end associate
+
+  end procedure
+
 
 end submodule do_concurrent_s