Visualization Tutorial (#49)

Author: suzm10

.github/workflows/mac.yml

@@ -9,6 +9,7 @@ jobs:
   self:
     name: Mac Runner
     runs-on: macos-latest
+    if: github.repository == 'Comp-Physics/RBC3D'
     continue-on-error: true
     steps:
       - name: Checkout

.gitignore

@@ -16,3 +16,5 @@ packages*
 traction
 build
 *.sbatch
+*/*/*/restart*
+examples/*/*.py

README.md

@@ -30,11 +30,11 @@ brew install gcc mpich gfortran pkg-config wget cmake
 ./rbc.sh install-mac
 ```
 
-and then set these environment variables in your `~/.zshrc` or `~/.bashrc`. Note that `$HOME` will need to be replaced with the folder you cloned RBC3D in.
+and then from the RBC3D root directory, run these commands but replace `.zshrc` with where you store environment variables:
 
 ```shell
-export PETSC_DIR=$HOME/RBC3D/packages/petsc-3.19.6
-export PETSC_ARCH=arch-darwin-c-opt
+rootdir=`pwd`
+echo -e "export PETSC_DIR=$rootdir/packages/petsc-3.21.3 \nexport PETSC_ARCH=arch-darwin-c-opt" >> ~/.zshrc
 ```
 
 Then to execute and run a case, you can:
@@ -56,8 +56,7 @@ mpiexec -n 1 ../../build/case/minit
 mpiexec -n 2 ../../build/case/mtube
 ```
 
-To run a case with more cells and nodes, you should use a supercomputing cluster. Instructions on how to build RBC3D on a cluster are [available here](https://github.com/comp-physics/RBC3D/blob/master/install/readme.md).
-
+To run a case with more cells and nodes, you should use a supercomputing cluster. Instructions on how to build RBC3D on a cluster are [available here](https://github.com/comp-physics/RBC3D/blob/master/docs/clusters.md).
 
 ### Papers that use RBC3D
 

common/ModConf.F90

@@ -249,11 +249,11 @@ subroutine ReadConfig(fn)
       print *, 'refRad = ', refRadIN(1:nCellTypes)
 
       read (funit, *) Deflate; print *, 'Deflate = ', Deflate
-      !print *, 'before error'
+      ! print *, 'before error'
       read (funit, *) pGradTar(1); print *, 'pGradTar = ', pGradTar(1)
       read (funit, *) pGradTar(2); print *, 'pGradTar = ', pGradTar(2)
       read (funit, *) pGradTar(3); print *, 'pGradTar = ', pGradTar(3)
-      !print *, 'after error'
+      ! print *, 'after error'
       read (funit, *) Nt; print *, 'Nt = ', Nt
       read (funit, *) Ts; print *, 'Ts = ', Ts
 

common/ModDataStruct.F90

@@ -92,7 +92,7 @@ module ModDataStruct
 !
 !======================================================================
 ! Cell type
-!  celltype --   1  red cell ; 2 leukocyte ; 3 platelet
+!  celltype --   1 red cell ; 2 leukocyte ; 3 platelet
 !
 !======================================================================
 ! Solid body modes

common/ModIO.F90

@@ -553,8 +553,6 @@ subroutine ReadWallMesh(fn, wall)
     integer, allocatable :: connect(:, :)
     character(*), parameter :: func_name = "ReadWallMesh"
 
-    print *, "fn: ", trim(fn)
-
     ! Check whether the file exists
     ierr = NF90_OPEN(trim(fn), NF90_NOWRITE, ncid)
     if (ierr .ne. 0) then
@@ -806,20 +804,13 @@ subroutine ReadRestart(fn)
     ! global parameters
     if (rootWorld) then
       write (*, *) 'Read restart file ', TRIM(fn)
-      print *, 'error 0'
       open (restart_unit, file=trim(fn), form='unformatted', action='read')
 
-      print *, 'error 1'
       read (restart_unit) Lb
       iLb = 1./Lb
-      print *, 'error 2'
-      print *, 'z'
       read (restart_unit) Nt0
-      print *, 'a'
       read (restart_unit) time0
-      print *, 'b'
       read (restart_unit) vBkg
-      print *, '0'
       write (*, *) 'Lb = ', Lb
       write (*, *) 'Nt0 = ', Nt0
       write (*, *) 'time0 = ', time0

common/ModPostProcess.F90

@@ -103,11 +103,11 @@ function CellFlowRate() result(flowrate)
       zc = 0.5*(minval(rbc%x(:, :, 3)) + maxval(rbc%x(:, :, 3)))
 
       do ilon = 1, rbc%nlon
-      do ilat = 1, rbc%nlat
-        ds = rbc%detJ(ilat, ilon)*rbc%w(ilat)
-        vn = dot_product(rbc%g(ilat, ilon, :), rbc%a3(ilat, ilon, :))
-        flowrate = flowrate + (rbc%x(ilat, ilon, 3) - zc)*vn*ds
-      end do ! ilat
+        do ilat = 1, rbc%nlat
+          ds = rbc%detJ(ilat, ilon)*rbc%w(ilat)
+          vn = dot_product(rbc%g(ilat, ilon, :), rbc%a3(ilat, ilon, :))
+          flowrate = flowrate + (rbc%x(ilat, ilon, 3) - zc)*vn*ds
+        end do ! ilat
       end do ! ilon
     end do ! irbc
 
@@ -134,12 +134,12 @@ subroutine ComputeParticleStress(tau)
       end do ! ii
 
       do ilon = 1, rbc%nlon
-      do ilat = 1, rbc%nlat
-        x = rbc%x(ilat, ilon, :) - xc
-        f = rbc%f(ilat, ilon, :)
-        dS = rbc%w(ilat)*rbc%detJ(ilat, ilon)
-        forall (ii=1:3, jj=1:3) tau(ii, jj) = tau(ii, jj) - dS*f(ii)*x(jj)
-      end do ! ilat
+        do ilat = 1, rbc%nlat
+          x = rbc%x(ilat, ilon, :) - xc
+          f = rbc%f(ilat, ilon, :)
+          dS = rbc%w(ilat)*rbc%detJ(ilat, ilon)
+          forall (ii=1:3, jj=1:3) tau(ii, jj) = tau(ii, jj) - dS*f(ii)*x(jj)
+        end do ! ilat
       end do ! ilon
     end do ! irbc
 

common/ModRbc.F90

@@ -423,38 +423,38 @@ subroutine RBC_ComputeGeometry(cell)
 
     ! Surface tangential and normals
     do ilon = 1, nlon
-    do ilat = 1, nlat
-      a1 = cell%a1(ilat, ilon, :)
-      a2 = cell%a2(ilat, ilon, :)
+      do ilat = 1, nlat
+        a1 = cell%a1(ilat, ilon, :)
+        a2 = cell%a2(ilat, ilon, :)
 
-      a(1, 1) = dot_product(a1, a1)
-      a(1, 2) = dot_product(a1, a2)
-      a(2, 1) = a(1, 2)
-      a(2, 2) = dot_product(a2, a2)
+        a(1, 1) = dot_product(a1, a1)
+        a(1, 2) = dot_product(a1, a2)
+        a(2, 1) = a(1, 2)
+        a(2, 2) = dot_product(a2, a2)
 
-      detA = a(1, 1)*a(2, 2) - a(1, 2)*a(2, 1)
-      iDetA = 1./detA
+        detA = a(1, 1)*a(2, 2) - a(1, 2)*a(2, 1)
+        iDetA = 1./detA
 
-      a_rcp(1, 1) = iDeta*a(2, 2)
-      a_rcp(1, 2) = -iDeta*a(1, 2)
-      a_rcp(2, 1) = -iDeta*a(2, 1)
-      a_rcp(2, 2) = iDeta*a(1, 1)
+        a_rcp(1, 1) = iDeta*a(2, 2)
+        a_rcp(1, 2) = -iDeta*a(1, 2)
+        a_rcp(2, 1) = -iDeta*a(2, 1)
+        a_rcp(2, 2) = iDeta*a(1, 1)
 
-      a1_rcp = a_rcp(1, 1)*a1 + a_rcp(1, 2)*a2
-      a2_rcp = a_rcp(2, 1)*a1 + a_rcp(2, 2)*a2
+        a1_rcp = a_rcp(1, 1)*a1 + a_rcp(1, 2)*a2
+        a2_rcp = a_rcp(2, 1)*a1 + a_rcp(2, 2)*a2
 
-      ! Store results
-      cell%a1_rcp(ilat, ilon, :) = a1_rcp
-      cell%a2_rcp(ilat, ilon, :) = a2_rcp
+        ! Store results
+        cell%a1_rcp(ilat, ilon, :) = a1_rcp
+        cell%a2_rcp(ilat, ilon, :) = a2_rcp
 
-      cell%a(ilat, ilon, :, :) = a
-      cell%a_rcp(ilat, ilon, :, :) = a_rcp
+        cell%a(ilat, ilon, :, :) = a
+        cell%a_rcp(ilat, ilon, :, :) = a_rcp
 
-      cell%detj(ilat, ilon) = sqrt(detA)/sin(cell%th(ilat))
+        cell%detj(ilat, ilon) = sqrt(detA)/sin(cell%th(ilat))
 
-      a3 = CrossProd(a1, a2)
-      cell%a3(ilat, ilon, :) = a3/NORM2(a3)
-    end do ! ilat
+        a3 = CrossProd(a1, a2)
+        cell%a3(ilat, ilon, :) = a3/norm2(a3)
+      end do ! ilat
     end do ! ilon
 
     ! Compute curvature tensor
@@ -920,22 +920,22 @@ subroutine Shell_ElasStrs(cell, cellRef, tau)
     nlon = cell%nlon
 
     do ilat = 1, nlat
-    do ilon = 1, nlon
-      ! V2 is the stretch tensor
-      V2 = matmul(cell%a(ilat, ilon, :, :), cellRef%a_rcp(ilat, ilon, :, :))
+      do ilon = 1, nlon
+        ! V2 is the stretch tensor
+        V2 = matmul(cell%a(ilat, ilon, :, :), cellRef%a_rcp(ilat, ilon, :, :))
 
-      lbd1 = V2(1, 1) + V2(2, 2) - 2.0
-      lbd2 = V2(1, 1)*V2(2, 2) - V2(1, 2)*V2(2, 1) - 1.0
-      JS = sqrt(V2(1, 1)*V2(2, 2) - V2(1, 2)*V2(2, 1))
+        lbd1 = V2(1, 1) + V2(2, 2) - 2.0
+        lbd2 = V2(1, 1)*V2(2, 2) - V2(1, 2)*V2(2, 1) - 1.0
+        JS = sqrt(V2(1, 1)*V2(2, 2) - V2(1, 2)*V2(2, 1))
 
-      ! Take the covariant form of V2
-      V2 = cellRef%a_rcp(ilat, ilon, :, :)
+        ! Take the covariant form of V2
+        V2 = cellRef%a_rcp(ilat, ilon, :, :)
 
-      tau_tmp = 0.5*ES/JS*(lbd1 + 1.0)*V2 &
-                + 0.5*JS*(ED*lbd2 - ES)*cell%a_rcp(ilat, ilon, :, :)
+        tau_tmp = 0.5*ES/JS*(lbd1 + 1.0)*V2 &
+                  + 0.5*JS*(ED*lbd2 - ES)*cell%a_rcp(ilat, ilon, :, :)
 
-      tau(ilat, ilon, :, :) = tau_tmp
-    end do ! ilon
+        tau(ilat, ilon, :, :) = tau_tmp
+      end do ! ilon
     end do ! ilat
 
   end subroutine Shell_ElasStrs
@@ -957,12 +957,12 @@ subroutine Shell_Bend(cell, cellRef, bnd)
     nlon = cell%nlon
 
     do ilat = 1, nlat
-    do ilon = 1, nlon
-      b = matmul(cell%a_rcp(ilat, ilon, :, :), cell%b(ilat, ilon, :, :))
-      bref = matmul(cellRef%a_rcp(ilat, ilon, :, :), cellRef%b(ilat, ilon, :, :))
+      do ilon = 1, nlon
+        b = matmul(cell%a_rcp(ilat, ilon, :, :), cell%b(ilat, ilon, :, :))
+        bref = matmul(cellRef%a_rcp(ilat, ilon, :, :), cellRef%b(ilat, ilon, :, :))
 
-      bnd(ilat, ilon, :, :) = -cell%EB*matmul(b - bref, cell%a_rcp(ilat, ilon, :, :))
-    end do ! ilon
+        bnd(ilat, ilon, :, :) = -cell%EB*matmul(b - bref, cell%a_rcp(ilat, ilon, :, :))
+      end do ! ilon
     end do ! ilat
 
   end subroutine Shell_Bend

common/ModVelSolver.F90

@@ -78,8 +78,7 @@ subroutine Solve_RBC_Vel(v)
       call MatShellSetOperation(mat_lhs, MATOP_MULT, MyMatMult, ierr)
 
       call KSPCreate(PETSC_COMM_SELF, ksp_lhs, ierr)
-      ! call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr)
-      call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, ierr)
+      call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, ierr) ! call KSPSetOperators(ksp_lhs, mat_lhs, mat_lhs, SAME_NONZERO_PATTERN, ierr)
       call KSPSetType(ksp_lhs, KSPGMRES, ierr)
       call KSPGetPC(ksp_lhs, pc_lhs, ierr)
       call KSPSetInitialGuessNonzero(ksp_lhs, PETSC_TRUE, ierr)  ! TRIAL

install/clusters.md docs/clusters.mdinstall/clusters.md renamed to docs/clusters.md

@@ -36,13 +36,17 @@ module load gcc/12.3.0 mvapich2/2.3.7-1 netcdf-c hdf5/1.14.1-2-mva2 intel-oneapi
 ./rbc.sh install-ice
 ```
 
-Before you can run cmake, you must set these environment variables. You can place them in your `~/.bashrc`. If you didn't place `RBC3D` in your `$HOME` directory, then replace it with where you placed `RBC3D`.
+## Environment Variables
+
+Before you can run cmake, you must set `PETSC_DIR` and `PETSC_ARCH` environment variables. You can place them in your `~/.bashrc`. This path depends on where you placed RBC3D. To get the path to where you placed it you can run this from your RBC3D root directory:
 
 ```shell
-export PETSC_DIR=$HOME/RBC3D/packages/petsc-3.19.6
-export PETSC_ARCH=arch-linux-c-opt
+rootdir=`pwd`
+echo -e "export PETSC_DIR=$rootdir/packages/petsc-3.21.3 \nexport PETSC_ARCH=arch-linux-c-opt" >> ~/.bashrc
 ```
 
+## Running an Example Case
+
 Then to execute and run a case, you can:
 ```shell
 mkdir build