Added grid and gradient for regular communication between images

src/matcha/grid_m.f90

@@ -0,0 +1,49 @@
+! Copyright (c), The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module gridr_m
+  implicit none
+
+  private
+  public :: gridr_t
+
+  type gridr_t
+    private
+  contains
+    procedure :: gridparameters
+    procedure :: gradient
+  end type  
+
+  interface gridr_t
+
+   pure module function construct() result(gridr)
+     implicit none
+     type(gridr_t) gridr
+   end function
+
+  end interface
+
+  interface
+
+    pure module function gridparameters(self,gb,ge,ng) result(gridp)
+     implicit none
+     class(gridr_t), intent(in) :: self
+     integer, intent(in) :: ng(:)
+     double precision, intent(in) :: gb,ge
+     double precision gridp(7)
+    end function gridparameters     
+
+    module function gradient(self, my_num_cells, ng, dx, gb, tconc, x) result(gx)      
+      implicit none
+      class(gridr_t), intent(in) :: self
+      integer, intent(in) :: my_num_cells
+      integer, intent(in) :: ng(:)
+      double precision, intent(in) :: dx(:)
+      double precision, intent(in) :: gb,tconc
+      double precision, intent(in) :: x(:,:)
+      double precision, allocatable :: gx(:,:)
+    end function gradient
+
+  end interface
+
+end module gridr_m

src/matcha/grid_s.F90

@@ -0,0 +1,70 @@
+! Copyright (c), The Regents of the University of California
+! Terms of use are as specified in LICENSE.tx
+submodule(gridr_m) gridr_s
+  use intrinsic_array_m, only : intrinsic_array_t
+
+#ifdef USE_CAFFEINE
+   use caffeine_assert_m, only : assert
+#else
+   use assert_m, only : assert
+#endif
+
+  implicit none
+
+contains
+
+  module procedure construct
+  end procedure construct
+
+  module procedure gridparameters
+
+    gridp(1) = (ge - gb)/dble(ng(1))
+    gridp(2) = (ge - gb)/dble(ng(2))
+    gridp(3) = (ge - gb)/dble(ng(3))
+
+    gridp(4) = 1.d0/(2.d0*gridp(1))
+    gridp(5) = 1.d0/(2.d0*gridp(2))
+    gridp(6) = 1.d0/(2.d0*gridp(3))
+
+    gridp(7) = 1.d0/(gridp(1)*gridp(2)*gridp(3))
+
+  end procedure gridparameters
+
+
+  module procedure gradient
+    integer i,j,k,ii,jj,kk
+    double precision, allocatable :: concentration(:,:,:)
+    double precision, allocatable :: grad(:,:,:,:)
+
+    allocate(concentration(ng(1),ng(2),ng(3)))
+    allocate(grad(ng(1),ng(2),ng(3),3))    
+
+    concentration = 0.d0
+    grad = 0.d0
+
+    do i = 1,my_num_cells
+       ii = int((x(i,1) - gb)/dx(1)) + 1
+       jj = int((x(i,2) - gb)/dx(2)) + 1
+       kk = int((x(i,3) - gb)/dx(3)) + 1
+       concentration(ii,jj,kk) = concentration(ii,jj,kk) + tconc
+    end do
+
+    call co_sum(concentration)
+
+    do concurrent(i = 2:ng(1)-1, j = 2:ng(2)-1, k = 2:ng(3)-1)  
+       grad(i,j,k,1) = dx(4)*(concentration(i+1,j,k)-concentration(i-1,j,k))
+       grad(i,j,k,2) = dx(5)*(concentration(i,j+1,k)-concentration(i,j-1,k))
+       grad(i,j,k,3) = dx(6)*(concentration(i,j,k+1)-concentration(i,j,k-1))
+    end do
+
+    allocate(gx(my_num_cells,3))
+    do concurrent(i = 1:my_num_cells)
+       ii = int((x(i,1) - gb)/dx(1)) + 1
+       jj = int((x(i,2) - gb)/dx(2)) + 1
+       kk = int((x(i,3) - gb)/dx(3)) + 1
+       gx(i,:) = grad(ii,jj,kk,:)
+    end do
+
+  end procedure gradient
+
+end submodule gridr_s

src/matcha/input_m.f90

@@ -8,16 +8,21 @@ module input_m
 
   type input_t
     private
-    integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4
-    double precision :: time_step_ = 0.1D0
-    !double precision, allocatable :: sample_distribution_(:,:)
-    !allocate(sample_distribution_(num_intervals_,2))
+    integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4, ngrid_ = 50
+
+    double precision :: time_step_ = 0.1D-4
+    double precision :: grid_begin_ = -500.d0, grid_end_ = 600.d0, cytokine_ = 1.d-6, gfac_ = 1.d-9
   contains
     procedure :: num_cells
     procedure :: num_positions
     procedure :: num_dimensions
     procedure :: num_intervals
+    procedure :: ngrid
     procedure :: time_step
+    procedure :: grid_begin
+    procedure :: grid_end
+    procedure :: cytokine
+    procedure :: gfac
     procedure :: sample_distribution
   end type
 
@@ -46,13 +51,44 @@ pure module function num_intervals(self) result(n)
       class(input_t), intent(in) :: self
       integer n
     end function 
+
+    pure module function ngrid(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
       class(input_t), intent(in) :: self
       double precision dt
     end function time_step
+
+    pure module function grid_begin(self) result(gg)
+      implicit none
+      class(input_t), intent(in) :: self
+      double precision gg
+    end function grid_begin
+
+    pure module function grid_end(self) result(gg)
+      implicit none
+      class(input_t), intent(in) :: self
+      double precision gg
+    end function grid_end
+
+    pure module function cytokine(self) result(gg)
+      implicit none
+      class(input_t), intent(in) :: self
+      double precision gg
+    end function cytokine
+
+    pure module function gfac(self) result(gg)
+      implicit none
+      class(input_t), intent(in) :: self
+      double precision gg
+    end function gfac
 
+
     pure module function sample_distribution(self) result(empirical_distribution)
       implicit none
       class(input_t), intent(in) :: self

src/matcha/input_s.f90

@@ -21,9 +21,31 @@
       n = self%num_intervals_
   end procedure
 
+  module procedure ngrid
+      n = self%ngrid_
+  end procedure
+
   module procedure time_step
       dt = self%time_step_
   end procedure
+
+  module procedure grid_begin
+      gg = self%grid_begin_
+  end procedure
+
+  module procedure grid_end
+      gg = self%grid_end_
+  end procedure
+
+  module procedure cytokine
+      gg = self%cytokine_
+  end procedure
+
+  module procedure gfac
+      gg = self%gfac_
+  end procedure        
+
+
 
   module procedure sample_distribution
      integer i,nintervals

src/matcha_s.F90

@@ -3,8 +3,8 @@
 submodule(matcha_m) matcha_s
   use t_cell_collection_m, only : t_cell_collection_t
   use distribution_m, only : distribution_t
+  use gridr_m, only : gridr_t
   use data_partition_m, only : data_partition_t
-
 #ifdef USE_CAFFEINE
    use caffeine_m, only : this_image => caf_this_image
 #endif
@@ -14,31 +14,49 @@
 contains
 
   module procedure matcha
-
+  
     associate( &
       ncells => input%num_cells(), &
       npositions => input%num_positions(), &
       ndim => input%num_dimensions(), &
       nintervals => input%num_intervals(), &
+      num_grid => input%ngrid(), &
       dt => input%time_step(), &
+      gb => input%grid_begin(), &
+      ge => input%grid_end(), &
+      cytokinev => input%cytokine(), &
+      gfacv => input%gfac(), &
       empirical_distribution => input%sample_distribution() &
-    )
-
+      )
+      
       block
         double precision, parameter :: scale = 100.D0
         double precision, allocatable :: random_positions(:,:), random_4vectors(:,:,:)
+        double precision, allocatable :: gdx(:)
+        double precision gxs
         type(distribution_t) distribution
+        type(gridr_t) gridr
+        integer, allocatable :: ng(:)
         integer, parameter :: nveldim = 4
+        integer i,j,k,ii,jj,kk
         integer step
+        double precision tconc
         type(data_partition_t) data_partition
+
+        allocate(ng(ndim))
+        allocate(gdx(2*ndim+1))
+        ng(:) = num_grid
+        gridr = gridr_t()
+        gdx = gridr%gridparameters(gb,ge,ng)
+        tconc = cytokinev*gfacv*gdx(7)
+
 
         call data_partition%define_partitions(cardinality=ncells)
 
         associate(me => this_image())
           associate(my_num_cells => data_partition%last(me) - data_partition%first(me) + 1)
 
             call random_init(repeatable=.true., image_distinct=.true.)
-
             allocate(random_positions(my_num_cells,ndim))
             call random_number(random_positions)  
 
@@ -52,8 +70,11 @@
                   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)
+                     associate(x => history(step-1)%positions(), t => history(step-1)%time())
+                       associate(gx => gridr%gradient(my_num_cells,ng,gdx,gb,tconc,x))  
+                         history(step) = t_cell_collection_t(x + (v(:,step-1,:) + gx(:,:))*dt, t + dt)
+                       end associate
+!                       history(step) = t_cell_collection_t(x + v(:,step-1,:)*dt, t + dt)
                     end associate
                   end do
                 end associate