Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Sickle Cell Model and Related Tools to RBC3D Simulation #1

Merged
merged 22 commits into from
Aug 1, 2023

Conversation

smullangi3
Copy link
Contributor

@smullangi3 smullangi3 commented Jul 29, 2023

This pull-request adds an input file SickleCell.dat, as well as functions in ModIO.F90 for reading the input file into a cell for simulation.

ModRBC.F90 modifications include:

  • another possible sickle cell model (commented). Unfortunately, this model was ill-conditioned for the simulation
  • The "RBC_MakeUncurvedSickleSpheroid" subroutine, which creates an uncurved prolate spheroid. To produce the "SickleCell.dat" Sickle Cell model, this uncurved prolate spheroid was then inserted into a flow: vbkg = 8 for 5000 timesteps of size 0.00014, and with mechanical properties equivalent to a healthy Biconcave RBC (celltype 1 in case/tube.F90).

Tube.F90 also contains modifications to include the Sickle cell as celltype #3, where mechanical properties such as shear and bending modulus are set to 3x and 4x higher than a conventional healthy RBC respectively.

As seen in Initcond.F90, we can create a Sickle Cell in a simulation by replacing:
call Rbc_Create(rbc, nlat0, dealias)
call RBC_MakeBiconcave(rbc, radEqv, xc)
with:
call ImportReadRbc('Input/SickleCell.dat', rbc, xc)
wherever we want a Sickle Cell instead of a Biconcave RBC.

Currently, Initcond sets up a simulation of 8 cells in single file inside a tube, where cells alternate between healthy and sickle.

This code has been tested through simulations, comprising of:

  • 1 sickle cell
  • Many sickle cells
  • Combinations of sickle and other cells
    In all cases, errors don't occur and the simulation behaves as expected (as long as the timestep size isn't too large to cause divergence - a suitable size is set in the below input file case_different_celltypes/Input/tube.in).

(code commits are a bit all-over, i recommend probably doing a squash merge once this goes in)

1. ! refRad
3 ! nCellTypes
1 1 1 ! viscRat
1. 1. 1. ! refRad
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this doesn't seem right, shoudn't it be 1 number? what about viscRat? Regardless, please create multiple "cases" that serve different example purposes.

Copy link
Contributor Author

@smullangi3 smullangi3 Jul 30, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ModConf.F90 line 245:

read(funit, *) nCellTypes; print *, 'nCellTypes = ',nCelltypes read(funit, *) viscRatIN(1:nCellTypes) print *, 'viscRat = ', viscRatIN(1:nCellTypes) read(funit, *) refRadIN(1:nCellTypes) print *, 'refRad = ', refRadIN(1:nCellTypes)

so, it looks like you need as many viscRats and refRads as the number of cell types.

Yes, I will move this sickle+rbc simulation to a new case folder.

10000 ! Nt
0.0014 ! Ts
10000 ! Nt
0.00014 ! Ts
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this very small time step is an example of why this example case should be separate from a health rbc one, so we need more than one example

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, will create new case-folder for rev2 for the changes.

common/ModIO.F90 Outdated
close(cell_unit)
end subroutine WriteRBCPlain

subroutine ReadRBCPlain(fn, rbc, xc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

by plain do you mean healthy?

Copy link
Contributor Author

@smullangi3 smullangi3 Jul 30, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, no this was bad naming on my part. I just mean that:
WriteRBCPlain takes a cell and a filename, and writes just the cell's plain information (w/o any transformations). Then, ReadRBCPlain can reread in this cell model for a different simulation.

For next rev, I will add better comments and naming convention

call Rbc_Create(rbc, nlat0, dealias)
call RBC_MakeBiconcave(rbc, radEqv, xc)
else
call ReadRBCPlain('Input/SickleCell.dat', rbc, xc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so plain is sickle? also, what it seems like the input data depends on the values of nlat and nlon, or no? if so, what does one do if they want a sickle cell with different nlat and nlon?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, in this case, we have a file of a Sickle Cell (previously exported) which is imported using the ReadRBCPlain subroutine. I will fix naming and explain better via comments in next revision.

The imported cell will be the same exact shape as the cell which was previously exported to the defined file. So, if we want a Sickle Cell with different nlat/nlon, we would have to recreate the sickle cell with those different measurements and then re-export it.

case/tube.F90 Outdated
call RBC_ComputeGeometry(rbcRef)

rbcRef => rbcRefs(2)
call RBC_Create(rbcRef, nlat0)
call RBC_MakeSphere(rbcRef, radEqv)
call RBC_ComputeGeometry(rbcRef)

rbcRef => rbcRefs(3)
call ReadRBCPlain('Input/SickleCell.dat', rbcRef)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so plain is sickle cell? what happens if the number of points on the sickle cell nlat and nlon doesn't match the values the user inputs?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the sickle cell's nlat/nlon in this case will just be the same as the nlat/nlon of the exported cell-shape. So, even if the user changes the input nlat0, the sickle cell disregards that and uses the model's data from the Input/SickleCell.dat file.

@@ -175,6 +177,115 @@ subroutine RBC_MakeSphere(cell, r, xc)
end if

end subroutine RBC_MakeSphere

subroutine RBC_MakeSickleGraham(cell, r, xc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please make this more descriptive than Graham, also put a reference to the paper where these ideas come from in the comments. so far the code in this PR has no comments (!!)

Copy link
Contributor Author

@smullangi3 smullangi3 Jul 30, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, will change to better meaningful names and add explanations in the comments for rev2

@@ -0,0 +1,4 @@
12 24
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this looks like nlat and nlon and a dealiasing factor.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep, this was the nlat and nlon (12, 24) of the cell in the file. So, when I import the cell, it will always be with these dimensions (defined in the file when exported).

Makefile.in Outdated
@@ -1,5 +1,5 @@
# Directories (change this to where you put RBC3D!)
WORK_DIR = /storage/home/hcoda1/6/sbryngelson3/p-sbryngelson3-0/test-rbc/RBC3D
WORK_DIR = /storage/home/hcoda1/6/smullangi3/RBC3D
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove this from PR

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will do in rev2

8 ! PBspln_Ewld

3 ! nCellTypes
1 1 3 ! viscRat (Sickle Cells experience ~3 to 4x viscosity as a healthy RBC) Reference: https://doi.org/10.1016/j.actbio.2012.07.011
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

make this 1. unclear if the simulation is stable for viscRat /= 1.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will change back to all 1s in rev3

common/ModIO.F90 Show resolved Hide resolved
! Unfortunately this shape doesn't work very well since points are ill-conditioned
! and diverge for this type of simulation, so we do not proceed with this sickle-cell model

! subroutine RBC_MakeSickle(cell, rad, xc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

uncomment this and then add an write('this configuration is not confirmed working yet' and an mpi_abort() call. do not leave stale code in the codebase.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will do this in rev3

! end if
! end do

! end function RBC_SolveSickleRho
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

uncomment this function

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure - also this function uses a different module for evaluating quartic roots. So I'll uncomment this and include that new ModPolyRoots module in rev3

@sbryngelson sbryngelson merged commit ea82c01 into comp-physics:master Aug 1, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants