Add AB2 and RK2; Randomized Cell Placement Support for Different Cell Types

common/ModTimeInt.F90

@@ -22,7 +22,7 @@ module ModTimeInt
 
   private
 
-  public :: TimeInt, &
+  public :: TimeInt_Euler, &
     TimeInt_AxiSymm, &
     TimeIntModVC, &
     OneTimeInt, &
@@ -32,12 +32,21 @@ module ModTimeInt
     TimeInt_Init, &
     TimeInt_Finalize, &
     Compute_Rbc_Vel, &
-    Compute_Rbc_Vel_SHB
+    Compute_Rbc_Vel_SHB, &
+    TimeInt_AB2, &
+    TimeInt_RK2
 
   private :: FilterRbcs, &
     ReboxRbcs, &
     PeriodicWall, &
-    AdjustBkgVel
+    AdjustBkgVel, &
+    t_rbcv
+
+
+  ! To store cell velocities for higher-order time integrator method
+  type t_rbcv
+    real(WP), dimension(:,:,:), pointer :: v
+  end type t_rbcv
 
 contains
 
@@ -94,11 +103,179 @@ subroutine TimeInt_Finalize
 
   end subroutine TimeInt_Finalize
 
+
+
+
+!***********************************************************************
+! Time Integrate using generic Runge-Kutta 2nd Order method
+! Implementation from Chapter 5.1 of "Numerical Methods for Ordinary Differential Systems, the Initial Value Problem" by JD Lambert
+  subroutine TimeInt_RK2(alpha)
+    real(WP) :: alpha
+
+    integer :: lt
+    integer :: irbc, iwall
+    type(t_rbc),pointer :: rbc
+    type(t_wall),pointer :: wall
+    real(WP) :: clockBgn, clockEnd
+    integer :: ierr
+
+    real(WP) :: b1, b2, c1, c2
+    type(t_rbcv), dimension(:), allocatable :: k1 ! temporary velocities storage
+
+    ! Set up RK2 constants from alpha
+    c1 = 0
+    c2 = alpha
+    b2 = 1 / (2 * alpha)
+    b1 = 1 - b2
+
+    !set up time
+    time = time0
+
+    !allocate for k1
+    allocate(k1(nrbc))
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      allocate(k1(irbc)%v, MOLD=rbc%v)
+    end do
+
+    !evaluate with RK2 generalized
+    do lt = Nt0+1, Nt
+      clockBgn = MPI_Wtime()
+
+      !compute velocity (k1)
+      call Compute_Rbc_Vel
+      call NoSlipWall
+
+      !save current velocity into k1 tmp-var
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        k1(irbc)%v = rbc%v
+      end do
+
+      !step rbcs forward so that cells are at (x_n + c2*h*U_k)
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        rbc%x = rbc%x + c2*Ts*(k1(irbc)%v)
+      end do
+
+      !calculate velocity (k2)
+      call Compute_Rbc_Vel
+      call NoSlipWall
+
+      !revert cell positions back to x_n
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        rbc%x = rbc%x - c2*Ts*(k1(irbc)%v)
+      end do
+
+      !Use RK2 formula to find x_n+1
+      !note that k1 is stored in the k1 tmp-var, and k2 is in each rbc%v
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        rbc%x = rbc%x + Ts * (b1*(k1(irbc)%v) + b2*rbc%v)
+      end do
+
+      call ReboxRbcs
+
+      !update time
+      time = time + Ts
+
+      !output results
+      clockEnd = MPI_Wtime()
+      call WriteAll(lt, time)
+      if (rootWorld) then
+        write(*, '(A,I9,A,F15.5,A,F12.2)')  &
+        'lt = ', lt, '  T = ', time, ' time cost = ', clockEnd - clockBgn
+        write(*, *)
+      end if
+    end do
+
+    !deallocate k1 tmp var
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      deallocate(k1(irbc)%v)
+    end do
+    deallocate(k1)
+
+  end subroutine TimeInt_RK2
+
+!***********************************************************************
+! Time Integrate using Adams-Bashforth 2nd Order
+  subroutine TimeInt_AB2
+
+    integer :: lt
+    integer :: irbc, iwall
+    type(t_rbc),pointer :: rbc
+    type(t_wall),pointer ::wall
+    real(WP) :: clockBgn, clockEnd
+    integer :: ierr
+    !Velocity at k-1
+    type(t_rbcv), dimension(:), allocatable :: v_1
+
+    !evaluate first timestep with Euler-Forward
+    call Compute_Rbc_Vel
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      rbc%x = rbc%x + Ts * rbc%v
+    end do
+
+
+    !allocate space for each cell's previous velocity in v_1
+    allocate(v_1(nrbc))
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      allocate(v_1(irbc)%v, MOLD=rbc%v)
+    end do
+
+    !evaluate X_2...n with AB2
+    do lt = Nt0+2, Nt
+      clockBgn = MPI_WTime() !Start timing
+
+      !save previous velocity (U_n-1)
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        v_1(irbc)%v = rbc%v
+      end do
+
+      !compute current velocity (U_n)
+      call Compute_Rbc_Vel
+      call NoSlipWall
+
+      !Evolve RBCs with AB2:
+      do irbc = 1, nrbc
+        rbc => rbcs(irbc)
+        rbc%x = rbc%x + Ts*(3*rbc%v - v_1(irbc)%v)/2
+      end do
+
+      call ReboxRbcs
+
+      !update time
+      time = time + Ts
+      clockEnd = MPI_WTime()
+
+      !Output results
+      call WriteAll(lt, time)
+      if (rootWorld) then
+        write(*, '(A,I9,A,F15.5,A,F12.2)')  &
+        'lt = ', lt, '  T = ', time, ' time cost = ', clockEnd - clockBgn
+        write(*, *)
+      end if
+    end do
+
+    !deallocate tmp-var
+    do irbc = 1, nrbc
+      rbc => rbcs(irbc)
+      deallocate(v_1(irbc)%v)
+    end do
+    deallocate(v_1)
+  end subroutine TimeInt_AB2
+
+
 !**********************************************************************
 ! Time integrate using 
 ! Note:
 !  Time steps from Nt0+1 to Nt
-  subroutine TimeInt
+  subroutine TimeInt_Euler
 
     integer :: lt
     integer :: irbc, iwall
@@ -192,7 +369,7 @@ subroutine TimeInt
       end if
     end do ! lt
 
-  end subroutine TimeInt
+  end subroutine TimeInt_Euler
 
   subroutine TimeInt_AxiSymm
 

examples/case/tube.F90

@@ -20,7 +20,7 @@ program cells_in_tube
 
   call InitAll
 
-  call TimeInt
+  call TimeInt_Euler
 
   call FinalizeAll
 

examples/case_diff_celltypes/tube.F90

@@ -20,7 +20,7 @@ program cells_in_tube
 
   call InitAll
 
-  call TimeInt
+  call TimeInt_Euler
 
   call FinalizeAll
 

examples/randomized_case/Input/tube.in

@@ -2,17 +2,17 @@
 1.E-3                   ! eps_Ewld
 8                       ! PBspln_Ewld
 
-1                       ! nCellTypes
-1                       ! viscRat
-1.                      ! refRad
+3                       ! nCellTypes
+1 1 1                   ! viscRat
+1. 1. 1.                ! refRad
 .false.                 ! Deflate
 
 0.0     !
 0.0     !
 0.0     !
 
 10000                    ! Nt
-0.0014                   ! Ts
+0.00014                  ! Ts
 
 100                      ! cell_out
 -10000                   ! wall_out

examples/randomized_case/initcond_random.F90

@@ -9,7 +9,7 @@ program randomized_cell_gen
     use ModIO
     use ModBasicMath
 
-    integer,parameter :: ranseed = 484
+    integer,parameter :: ranseed = 161269
 
     !initial condition setup parameters
     real(WP), parameter :: hematocrit = 0.20
@@ -31,7 +31,7 @@ program randomized_cell_gen
     !calculate number of cells for the defined hematocrit, assuming all blood cells are healthy RBCs for volume
     !hematocrit = 4 * nrbc / (tube_radius^2 * tube_length)
     nrbcMax = (tubelen * tuber**2 * hematocrit) / 4
-
+    
     !set periodic boundary box based on tube shape
     Lb(1) = tuber * 2 + 0.5
     Lb(2) = tuber * 2 + 0.5
@@ -62,7 +62,7 @@ program randomized_cell_gen
     do i = 1,nrbcMax
 
         clockBgn = MPI_Wtime()
-        call place_cell
+        call place_cell(1)    
         clockEnd = MPI_Wtime()
 
         nrbc = nrbc + 1
@@ -75,10 +75,15 @@ program randomized_cell_gen
     !Write Cells Out to Cell TecPlot File
     write(fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat'
     call WriteManyRBCs(fn, nrbc, rbcs )
-    write(fn, FMT=fn_FMT) 'D/', '1x', 0, '.dat'
-    call WriteManyRBCsByType(fn, nrbc, rbcs, 1)
-    write(fn, FMT=fn_FMT) 'D/', '3x', 0, '.dat'
-    call WriteManyRBCsByType(fn, nrbc, rbcs, 3)
+
+    !Uncomment these lines to write cells to separate files by their celltype
+    !write(fn, FMT=fn_FMT) 'D/', '1x', 0, '.dat'
+    !call WriteManyRBCsByType(fn, nrbc, rbcs, 1)
+    !write(fn, FMT=fn_FMT) 'D/', '2x', 0, '.dat'
+    !call WriteManyRBCsByType(fn, nrbc, rbcs, 1)
+    !write(fn, FMT=fn_FMT) 'D/', '3x', 0, '.dat'
+    !call WriteManyRBCsByType(fn, nrbc, rbcs, 3)
+
     !Write Walls Out to Wall TecPlot File
     write(fn, FMT=fn_FMT) 'D/', 'wall', 0, '.dat'
     call WriteManyWalls(fn, nwall, walls )    
@@ -128,7 +133,7 @@ subroutine recenter
 end subroutine recenter
 
 !find an open spot in the simulation to place a cell 
-subroutine place_cell
+subroutine place_cell(celltype)
 
     type(t_Rbc) :: newcell
     type(t_Rbc), pointer :: cell 
@@ -138,12 +143,24 @@ subroutine place_cell
     integer :: nlat0
     logical :: place_success
 
+    integer :: celltype
+
     nlat0 = 12
 
     !create template newcell
-    call Rbc_Create(newcell, nlat0, 3)
-    call Rbc_MakeBiconcave(newcell, 1.)
-    newcell%celltype = 1
+    select case(celltype)
+        case(1)
+            call RBC_Create(newcell, nlat0)
+            call RBC_MakeBiconcave(newcell, 1.)
+        case(2)
+            call RBC_Create(newcell, nlat0)
+            call RBC_MakeLeukocyte(newcell, 1.)
+        case(3)
+            call ImportReadRBC('Input/SickleCell.dat', newcell)
+        case default
+            stop "bad cellcase"
+    end select
+    newcell%celltype = celltype
 
     place_success = .false.
 
@@ -239,18 +256,18 @@ logical function check_cell_collision(cell1, cell2)
     check_cell_collision = .false.
 
     !each point in cell1
-    do i = 1, cell1%nlat
-    do j = 1, cell1%nlon
+    do i = 1, cell1%nlat-1
+    do j = 1, cell1%nlon-1
 
     !define 2 diagonal points c1p & c2p for collision detection
     c1p = cell1%x(i, j, : )
-    c2p = cell1%x(modulo(i + 1, cell1%nlat), modulo(j + 1, cell1%nlon), : )
+    c2p = cell1%x(i + 1, j + 1, : )
 
     !each point in cell2
     do i2 = 1, cell2%nlat
     do j2 = 1, cell2%nlon
 
-        sp = cell2%x(i, j, : )
+        sp = cell2%x(i2, j2, : )
 
         !check if the cell2 point lies inside the cube delineated by c1p & c2p
         if ( &
@@ -270,4 +287,4 @@ logical function check_cell_collision(cell1, cell2)
 end function check_cell_collision
 
 
-end program randomized_cell_gen
+end program randomized_cell_gen