Visualization Tutorial

.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.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

docs/visualization.md

@@ -0,0 +1,91 @@
+# Visualizing a Simulation
+
+For this visualization tutorial, I ran `examples/randomized_case` with a smaller tube length (15), radius (4), and hematocrit (.18), but you can use any example case file. Then, I downloaded the files onto my local computer. From my local version of Paraview, I clicked open (top left button) and loaded `wall000000000.dat` and `x000000000.dat*` via group import. My window now looks like this after turning the wall opacity down:
+
+![alt text](images/image-6.png)
+
+Now, you can click the `Play` button at the top to see the cell output at each frame. This is at frame 21 of the simulation I ran. However, we can make the visualization look nicer with a few extra steps.
+
+![alt text](images/image-7.png)
+
+## Smooth Surfaces
+
+The wall and cell surfaces can be smoother if we use the generate surface normals filter. You can apply a filter by pressing `option + space` or `ctrl + space` on your keyboard and then typing in the name of the filter. To fully apply a filter, you have to select `Apply` in `Properties`.
+
+1. Select `wall000000000.dat` in the pipeline browser
+2. Select `Apply` in `Properties`
+3. Apply `Extract Surface` filter
+4. Select `Apply` in `Properties` again
+5. Apply `Generate Surface Normals` filter
+6. Repeat steps 2-6 with `x000000000.dat*` selected
+
+Now, your pipeline browser and render view should look like this:
+
+![alt text](images/image-8.png)
+
+## Colors
+
+We can change the color of the cells to red by following these steps:
+
+1. Open the properties panel for `GenerateSurfaceNormals2`
+2. Select `Coloring = Solid Color`
+3. Select `Edit` and using hex value `#980808` or a similar red color
+
+We can also change the background color to white although I do like paraview blue:
+
+1. Open propties panel and go to `Background`
+2. Deselect `Use Color Palette for Background`
+3. Select white for background color
+
+Now, our simulation looks like this:
+
+![alt text](images/image-9.png)
+
+## Ray-tracing and Lighting
+
+Going to the `Lighting` section of the properties panel and turning `Specular` to a higher value will make cells shiny and result in a nicer visualization.
+
+We can also use a few different rendering options to improve it further:
+
+1. Open properties panel
+2. Go to `Render Passes` section
+3. Click `Use Tone Mapping` and `Use Ambient Occlusion`
+4. Click use `Camera Parallel Projection` if you want a flatter look, but I'm keeping it off for this angle
+
+![alt text](images/image-11.png)
+
+To get started with ray-tracing the simulation, you can:
+
+1. Open properties panel and go to the `Ray Traced Rendering` section
+2. Select `Enable Ray Tracing`
+3. Select `OSPRay pathtracer` for the backend
+4. Set `Samples Per Pixel` to 5 or higher to get rid of graininess in the rendered image
+5. Set `Background Mode` to `Both` to keep whatever background you had previously
+6. Select `GenerateSurfaceNormals1` to get the wall
+7. Go to `Ray Tracing` section
+8. Select `Glass_thin` under the `Material` dropdown or another suitable material
+ * Note that you can add a material to the cells too if it doesn't change their color. `Value Indexed` works.
+
+Now, our simulation looks like this:
+
+![alt text](images/image-12.png)
+
+I can use a more advanced lighting interpolation called `PBR` instead of `Gouraud` to make the cells look different:
+
+![alt text](images/image-13.png)
+
+## Adding a Reflective Surface (Box)
+
+We can add a box underneath our tube so the cells have something to reflect off of. You can add this by typing `options + space + Box`, and then setting the box dimensions based off of your tube size. These are the dimensions I used:
+
+![alt text](images/image-14.png)
+
+I also changed the material of the box to something reflective, specifically `Metal_Lead_mirror`, but there might be a better material. 
+
+There are lots more Paraview options that you can change, but these are the basics! Ray-tracing takes up a lot of resources, so you might want to follow the remote visualization instructions [here](https://github.com/comp-physics/Scientific-Visualization?tab=readme-ov-file) if you're getting screenshots for a full simulation video. The `OSPRay raycaster` with `Shadows` turned on is a less intensive ray-tracer but still provides some visual interest.
+
+![alt text](images/image-16.png)
+
+# Making a Video
+
+In the previous steps, we edited a snapshot of one timestep of the simulation. To create a video showing the simulation data through all the timesteps, you can follow the steps in [this repo](https://github.com/comp-physics/Scientific-Visualization/blob/master/Tutorials/creating-an-annimation.md). It has instructions on how to save all the screenshots at each timestep into a folder and then clip them together with FFmpeg. Note that it may take several hours to save all the screenshots with ray-tracing, so it's probably better to not ray-trace the simulation if you just want to see the general flow.

examples/cutout/Input/tube.in

@@ -0,0 +1,30 @@
+0.44 !0.04 ! alpha_Ewld !changed for big tube
+1.E-3 ! eps_Ewld
+8 ! PBspln_Ewld
+
+3 ! nCellTypes
+1 1 1 ! viscRat
+1. 1. 1. ! refRad
+.false. ! Deflate
+
+0.0 !
+0.0 !
+0.0 !
+
+30000 ! Nt
+0.0014 ! Ts
+
+25 ! cell_out
+-10000 ! wall_out
+-10 ! pGrad_out
+-10 ! flow_out
+-1 ! ftotout
+100 ! restart_out
+
+'D/restart.LATEST.dat'
+
+0.02 ! epsDist
+10. ! forceCoef
+10. ! viscRatThresh
+.false. ! rigidsep
+0. 0. 0. 0. ! fmags