Visualization Tutorial (#49)

.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