Cannot validate the stretchCell case using HemoCell

[mojtabaAmiraslanpour]

I started to use the Hemocell example to validate Figure 4 of this paper. However, I ended up with this chart:

As you can see further I increase the value of force, my values in D_max would not match paper data. I have not changed any parameters in this case and I wonder why that happens. I also continued the simulation for 650K iterations but it didn't change that much. Do you have any comments about this?

Thanks, Mojtaba

[MaxvdKolk]

Hi Mojtaba,

The current definition of the RBC mechanics are slightly stiffer compared to the originally published values. This provides improved stability when running simulations with (high) shear flows, while maintaining similar force-displacement curves matching with the validated data set.

The various stiffnesses of the RBC are defined in mechanics/rbcHighOrderModel.cpp. To get closer to the original curves, you could modify the bending stiffness from the current 0.055 (line 157), to 2.467 (or more precisely (pi/2)^2).

Note that the current set of parameters are not unique, so various combinations of stiffnesses will show force-displacement curves matching the validation data set. In the next release will add additional documentation to clarify these constants and the resulting force-displacement curves.

[mojtabaAmiraslanpour]

Hi Mojtaba,

The current definition of the RBC mechanics are slightly stiffer compared to the originally published values. This provides improved stability when running simulations with (high) shear flows, while maintaining similar force-displacement curves matching with the validated data set.

The various stiffnesses of the RBC are defined in mechanics/rbcHighOrderModel.cpp. To get closer to the original curves, you could modify the bending stiffness from the current 0.055 (line 157), to 2.467 (or more precisely (pi/2)^2).

Note that the current set of parameters are not unique, so various combinations of stiffnesses will show force-displacement curves matching the validation data set. In the next release will add additional documentation to clarify these constants and the resulting force-displacement curves.

Thank you @MaxvdKolk for your response!

1. Stiffness of RCB: I tried that and the results were closer to the chart:

2. Resolution: As another approach, I tried increasing the resolution of the LB field and my results were closer for the F = 37 and 187 pN cases (not for the 87 pN case, do not know why) but I am not sure if I have done it correctly. When the simulation starts I got some Sanity check warnings that tell me values are wrong.

(HemoCell) (SanityCheck) WARNING: Fluid dx is not 5e-7 but 2.5e-07 This is unvalidated!

I also tried increasing the resolution for RCB itself so both get increased proportionally but got more sanity warnings! So is there any right way to do it?

3. Constant Volume and Surface: Also I have contacted "Gábor Závodszky" and he told me that there are some additional parameters that enforce constant volume, or surface area, etc. It has been done for numerical stability and stuff. So do you know where I can change them? and what are the correct values for them? Actually, I see that while running surface and volume of the RCB are changing. So I cannot understand this argument.

4. Do you have any more suggestions? Any more parameters to change?

Thanks! Mojtaba

[MaxvdKolk]

Hi Mojtaba,

This looks very close :+1:

I also tried increasing the resolution for RCB itself so both get increased proportionally but got more sanity warnings! So is there any right way to do it?

It might be that Gabor can provide more insight on the exact settings of these resolutions. As far as I know, the sanity checks will report warnings when using a different resolution dx than 5e-7, as (most of) the validations have been evaluated using that specific resolution. So for different resolutions, you might need to carefully run validation tests to see if the desired behaviour is still maintained throughout the simulations.

3. Constant Volume and Surface

The parameters are described in section 2.2 of the publication and are implemented within the cell mechanics model (see rbcHigherOrderModel.cpp). In release 2.3 those constants are given directly in the implementation, for instance:

• maximum surface area change: line 79
• maximum volumetric change: line 104
• maximum bending force: line 157
• maximum stretch length: line 181

As the combination of the validated parameters are not necessarily unique, you can find different combinations of these parameters that result in force-displacement curves that match the experimental data more/less closely. You might find different combinations of these parameters that fit better with the application in mind. For instance, the change to the bending stiffness (as compared to the originally published values) still provide validated behaviour, while being much more (numerically) stable when undergoing high shear flows.

Note: we are currently preparing release 2.4 which features more extensive documentation and which refactors these constants. In the new version, these will be defined explicitly under config/constant_defaults.h as

// Note: these are already expressed as their squared value
constexpr T MaxCellVolumetricChange = 0.01;
constexpr T MaxCellSurfaceAreaChange = 0.09;
constexpr T MaxCellBendingAngle = 0.0555;
constexpr T MaxCellPersistenceLength = 9.0;

which will make it easier for you to change these numbers in one location and observe the effects on the force-displacement curves. Hopefully this should be available somewhere next week.

Actually, I see that while running surface and volume of the RCB are changing. So I cannot understand this argument

Yes, there is some change of surface and volume possible of the RBC. For instance, you can take a look at the stretch response of the cell (section 2.2, the first bullet point and equation). The constant there tau_l corresponds to one of the constants defined in the code. Once the stretch reaches a certain amount, the force term will steeply increase and prevent further extension of the cell. You can think of these as barrier-like functions preventing too much change in either stretch, bending, surface area, or volume. Some change is possible, but bounded by these constants.

Do you have any more suggestions? Any more parameters to change?

You could investigate different values of the used time step (in this example dt = 1e-7 I believe), although I am not sure how bit of an influence that makes for this specific validation case. There are also parameters to set how often the particle model and fluid models are updated (stepMaterialEvery, stepParticleEvery), however for this validation case updating every iteration (default behaviour) should be fine.

Also make sure that you observe a steady state (similar to the second figure in your first post) when playing around with the different parameters. it might be that you require more/less when trying different sets of constants.

[mojtabaAmiraslanpour]

Thank you @MaxvdKolk to answer my questions thoroughly.

1. I mentioned that tau_bend = pi/6 in the paper but you suggested putting (pi/2)^2. Are there any reasons? I expect to get the same results in the paper by setting the same values and that's why I am not trying my own combinations of parameters.

2. The other thing I did not see a mesh dependency study in the paper (or maybe I missed it?). So I am not sure if my resolution is the same as the one in the paper.

Best

[MaxvdKolk]

Hi Mojtaba,

Agreed, not completely sure why the publication mentions tau_bend = pi/6. The code contained pi/2. Following the expression for the minimal representable curvature r_min = L_0/2 * sin(tau_bend / 2) with using here L_0 = 0.5 micron and tau_bend = pi/2, we find r_min = 0.176 micron. This is similar to the minimal curvature mentioned in the paper (0.18 micron). I did not notice the pi/6 mentioned in that paragraph.. and I have to apologise for not knowing if either I am misinterpreting these equations and the derivations or that there is a typo in the publication..

Note that the constants are given directly as their squared values. In all these force equations they will only appear squared, so we only store tau_bend^2 rather than tau_bend itself.

The other thing I did not see a mesh dependency study in the paper (or maybe I missed it?). So I am not sure if my resolution is the same as the one in the paper.

As far as I am aware, the study has been done for a specific configuration, using dx = 0.5 micron for the lattice grids. Using other resolutions should also work, although you might to reevaluate these curves to be completely sure. The default configurations of the stretching cell example match with the values used during the numerical experiments for the validation.

[mojtabaAmiraslanpour]

Good points @MaxvdKolk! I was not paying attention to the squared values of the constants! I corrected them (used tau_bend = pi/2), used dx=0.5e-6, and ran for a higher number of iterations (420'000), and here is what I got:

Overall better results than the other two.

1. The one thing that I cannot understand now is the initial volume value. It is given as an input inside the RCB.xml file but it actually is imported from the STL file, right? I put 71 um^3 for the volume value (which is mentioned in the paper but was 90 um^3 by default) but I get a value of 81 um^3 in the simulation log file.

2. Also, the membrane viscosity eta_m is set to be zero. is this correct?

Best, Mojtaba

[gzavo]

Hi @mojtabaAmiraslanpour, for your last two questions: 1) It is also possible to import cell shapes from STL, but in all the examples RBCs and PLTs are generated at runtime by the following algorithm: a) we generate a sphere with homogeneous surface vertex distribution, b) we subdivide the sphere to get to the requested vertex count, c) we project the surface vertices to the shape of the cell (e.g., biconcave for RBCs). This allows us to change the number of vertices from the config file, which is not possible using an STL. The volume is a tricky questions, because we have two types of it: (1) the mesh volume, which is used for instance at stretching tests, (2) the larger effective volume which is used for cell collision interaction. If you generate the cell, and not load it from STL, the volume depends on the radius you give in the xml. E.g., r = 3.91e-6 => V_mesh = 81.116 um3; V_effective = 97.687 um3 The effective volume happens due to the fluid-structure coupling. The immersed boundary method uses interpolation kernels (of size 0.5 um by default in HemoCell), that increases the outer radius of the cell by 0.25 um to an effective radius of 4.16e-6 um (which leads to the larger volume). Technically it means that two cells already interact when their IBM kernels start to overlap, which happens before their actual meshes touch. So when calculating bulk blood flow properties (e.g., viscosity) the effective volume yields the correct predictions. On the other hand if you are measuring the deformation of a single cell, you want the mesh radius instead of the effective interaction radius. 2) Yes, it is turned off by default, because it causes very little difference in most cases, however it forces significantly smaller time-steps. If you need very accurate deformation velocities you can turn it on (i.e. give it a non-zero value). Also note that it is different from cell viscosity.

Hope these helped!

[mojtabaAmiraslanpour]

Thank you @gzavo! It did help a lot. So if we take a look at the RCB.xml file:

...
    <radius> 3.91e-6 </radius> <!-- Radius of the RBC in [ 3.96 um] -->
    <Volume> 71 </Volume> <!-- Volume of the RBC in µm³ -->
...

We see there are both radius and volume defined here. What you say is that only the radius is playing a role in the current examples and not the volume. So Why do a have volume as a parameter in here after all?

I am just asking because in the paper it has been mentioned that the volume of the RCB is considered to be 71 um.

[gzavo]

You're right @mojtabaAmiraslanpour, I forgot to answer that part. It is used for some density calculation in the output. The reason is that the mesh volume is not the same as the effective volume, so you can supply the desired volume value.

[mojtabaAmiraslanpour]

Hello all, Right now I am trying to validate oneCellShear and pipeFlow cases which are mentioned in the paper. For the pipeflow case, I see that we need to put Re as an input so it can calculate the plasma velocity. I can't find any Re mentioned in the paper. The only thing I noted is the caption of figure 8 which says they used an average velocity of 1.5 cm/s which gives a Re=1.74 for D=128 um case. However, we know that velocity is not the same in all cases. So where can I find that?

P.S: Does it depend on the Re number, at all?

Thanks!

[mojtabaAmiraslanpour]

For the case of oneCellShear example, I got something like this for DI-iter chart: (the shear rate is set to be 111 s^-1)

So it gives an average DI value of about 13% which is not in agreement with figure 5 in the paper. When I increase tau_bening from 0.055 to 2.467 the cell is distorted and not reaching a steady or periodic state. Any comments on that?

[mojtabaAmiraslanpour]

@gzavo @MaxvdKolk Any comments on these?

The other issue I have is using the provided initial_states in the pipeFlow case. I see there is only one file called cells.pos in the folders, while we need two RBC.pos and PLT.pos files for our simulation. How should I get these two files separately?

[MaxvdKolk]

Hi @mojtabaAmiraslanpour,

Sorry for the delay, I was out of office for a couple days. Regarding the oneCellShear I would need to take a closer look to understand what might be going wrong... I will try to get back to you as soon as possible, hopefully before end of the week.

Thanks for pointing out the issue regarding the pipe flow cell positions. I was not aware of this and will make sure we resolve this before we provide the next release. Below my thoughts on how to resolve it at the moment.

---

For the initialisation of the pipe flow problems, I think that using a single cell.pos is old behaviour, maybe @gzavo can confirm this? Looking at the structure of the cell.pos, my guess to for an intermediate solution would be the following: I think you can split the cell.pos files into two parts named RBC.pos and PLT.pos. You can do so based on the first two lines in cell.pos that indicate the number of RBC and PLT respectively. To clarify:

# cell.pos
nRBC
nPLT
pos_rbc_1
...
pos_rbc_nRBC
pos_plt_1
...
pos_plt_nPLT

# RBC.pos
nRBC
pos_rbc_1
...
pos_rbc_nRBC

# PLT.pos
nPLT
pos_plt_1
...
pos_plt_nPLT

With this, you should be able to use the RBC and PLT positions that were used for the validation paper. Just make sure your domain is large enough (use the same resolution as the paper) such that all cell positions (both RBC and PLT) fit within the domain. Otherwise, the initialisation of the cells will ignore those cells that fall outside of the simulation domain.

[mojtabaAmiraslanpour]

Thank you @MaxvdKolk,

1. So I guessed right for the cell.pos file. I did the same in my simulations.
2. Right now I am working on the D32H20 case, So I have put 64 for refDirN. Is that right? (since each dx is equal to 0.5e-6)
3. Can you also give me your comments on the Re number value? I am getting different results for different Re numbers.
4. Also as for the last question, can you give me any insights on how you have computed Relative Apparent Viscosity? I see it is the division of two velocities in the code. I didn't find any definitions in relevant papers.

Thanks

[gzavo]

Hi @mojtabaAmiraslanpour,

1. Max is right, we used to use a single file for cell positions several years ago, when we only had RBCs and PLTs. The old files will be removed from the package, I much rather suggest using the packCell tool to generate the necessary cell distributions. Note the difference between tank hematocrit and tube hematocrit.
2. yes
3. Blood is a shear thinning fluid, a different Re in the same geometry means different shear, hence different viscosity. As I recall we used Re=1.0 in all cases of pipeflow reported in that paper. For the experiments this value is sadly not reported in the literature, so Re=1.0 is an educated guess based on the technology available in the early 90s, when most of the reference measurements (e.g., Pries '92) were made.
4. Relative apparent viscosity is the actual bulk viscosity of blood in the whole geometry compared to plasma viscosity. The two velocities are the actual velocity with the given driving body-force and the velocity predicted by Pousseuille law based on pure plasma viscosity and the same body-force.

[mojtabaAmiraslanpour]

Thank you @gzavo

4. Regarding this, I want to ask if you have considered the cytoplasm's viscosity too? Because a considerable amount of blood flow is occupied by cytoplasm which is inside the cells: So how have you considered this? In this picture, I see we have a continuous phase between inside and outside of the cells which in general is not the case. even if we consider the same viscosity for plasma and cytoplasm, again we would have some discontinuity, right?

[gzavo]

Hi @mojtabaAmiraslanpour , I suggest looking into this paper for answers on cytoplasmic viscosity.

[MaxvdKolk]

@mojtabaAmiraslanpour Regarding the oneCellShear example, you will see improved results when increasing the size of the simulation domain. The domain used in the example is rather small, mostly for illustrative purposes. This should reduce any effects of the boundary conditions on the deformations observed in the RBC.

Note, when increasing the domain size you should modify the initial position of the cell accordingly (in the .pos file) to ensure the red blood cell remains centred in the domain.

[mojtabaAmiraslanpour]

@gzavo Thank you for the paper!

1. Is this method already implemented in the Hemocell? If yes, where can I set the viscosity ratio (Λ) inside the code to see the effects?

2. @MaxvdKolk I am already running the cases with larger domains. Can this case run in parallel (by changing the code)? It actually takes so much time to solve on a single core. Besides, I am not changing the domain size in x and z-directions due to the periodic BC. Is that correct?

3. Also, I want to know if you have used Lees-Edwards BC in the oneCellShear case since it actually eliminates wall effects and we can run it with a small domain as well. Right? If this case is not using that, are there any examples that use it?

Best

[mojtabaAmiraslanpour]

I am trying to find the average DI value for the pipeFlow case. To do that I need to find the largest diameter of each cell, computing the DI for each cell, and then making an average. I use CellInformationFunctionals class to find the stretch value for cell 0:

CellInformationFunctionals::calculateCellStretch(&hemocell);
T largest_diam = 0.;
largest_diam = (CellInformationFunctionals::info_per_cell[0].stretch)/(1e-6/param::dx);

It actually works and shows the largest diameter for cell 0. But when I am trying to do it for other cid's it gives me 0 value. Am I doing something wrong here?

[MaxvdKolk]

I am already running the cases with larger domains. Can this case run in parallel (by changing the code)? It actually takes so much time to solve on a single core. Besides, I am not changing the domain size in x and z-directions due to the periodic BC. Is that correct?

Yes, as far as I am aware this case can be evaluated in parallel without difficulties. You can run the case normally with, for example, mpirun -n $nproc ./oneCellShear config.xml and set $nproc to the number of cores you desire (it will not perform better with more cores than hardware cores available in your system). Keeping the other dimensions the same size, should be OK. The boundary effects that influence the deformation index, mostly arise from the surfaces at which the shear velocity is prescribed.

Also, I want to know if you have used Lees-Edwards BC in the oneCellShear case since it actually eliminates wall effects and we can run it with a small domain as well. Right? If this case is not using that, are there any examples that use it?

The oneCellShear example does not use Lees-Edwards type boundary conditions. The example simply prescribes the velocity field on the opposing surfaces between which the shear flow is desired. The prescribed velocity magnitude is then based on the user-specified shear rate and the dimensions of the domain. It should be possible to utilise Lees-Edwards boundary conditions for this example. However, we do not yet provide a detailed example of the such boundary conditions as they are still an experimental feature and undergoing ongoing development. You are, of course, free to explore the Lees-Edwards functionality in the source code and try to use its functionalities. However, my suggestion would be to simply use the prescribed velocity on a (slightly) enlarged domain until we provide a documented example describing the Lees-Edwards functionality in detail.

It actually works and shows the largest diameter for cell 0. But when I am trying to do it for other cid's it gives me 0 value. Am I doing something wrong here?

This might be that you are accessing cell IDs that are not present in the actual simulations. This could be caused by using a *.pos file in which some cells are truncated, i.e. they are removed during initialisation as they overlap with the outer boundary of the simulation domain. In that case, the cell IDs might not be a complete range of integers and could skip one or more values.

You could, for instance, loop over all key,value pairs present in the info_per_cell map and access the stretch value for the individual cells to compute their deformation indices, something like:

while (hemocell.iter < tmax) {
    hemocell.iterate()

    ...

    CellInformationFunctionals::calculateCellStretch(&hemocell);

    // Loop over (ID, cell_info) pairs in the `info_per_cell` map
    for (auto const& cell_info: CellInformationFunctionals::info_per_cell) {
        auto diameter = (cell_info.second.stretch)/(1e-6/param::dx);
        // Note: D0 is a given constant, see the `oneCellShear` example.
        auto relative_diameter = (diameter/D0)*(diameter/D0);                                                                                      
        auto DI = (relative_diameter - 1.0)/(relative_diameter + 1.0)/100.0;
        hlog << "\n\tDeformation index:" << DI;
    }

    CellInformationFunctionals::clear_list();
}

You might need to do some additional post-processing on the obtained DI values, depending on your the values you desire to compute (average, max/min, etc).

[mojtabaAmiraslanpour]

@MaxvdKolk Thanks for your comprehensive explanations.

1. One thing that I came up with is that it seems DI-shear diagrams differ in these two papers:

Závodszky, Gábor, et al. "Cellular level in-silico modeling of blood rheology with an improved material model for red blood cells." Frontiers in physiology 8 (2017): 563.

De Haan, Mike, et al. "Numerical investigation of the effects of red blood cell cytoplasmic viscosity contrasts on single cell and bulk transport behaviour." Applied Sciences 8.9 (2018): 1616.

Here is the comparison:

I was running multiple simulations to get a DI of 7% for a shear rate of 111 1/s but mentioned that I got the correct value of 8% according to the 2018 paper. So now I wonder which one is correct?

2. The other thing I could not understand nor solve is this behavior in Couette BC:

The bottom face has a zero velocity value. In other words, the velocity is prescribed on the second cell from the bottom. However, I see it is correctly defined in the code:

    plb::Box3D top = plb::Box3D(0, nx-1, 0, ny-1, nz-1, nz-1);
    plb::Box3D bottom = plb::Box3D(0, nx-1, 0, ny-1, 0, 0);

This asymmetry causes a drift velocity for RBC in the +x direction (to the left in the image) as you may mention in the first image. Can you explain this?

Best, Mojtaba

[mojtabaAmiraslanpour]

Hello again,

1. I am trying to reproduce the results in Figure 12 of the Haan. et al paper. I use packCells to produce the initial cell positions for RBCs and PLTs using this command: packCells 140 70 70 -h 0.2 -r Obviously, I am giving x=140 um, y=70 um, z=70 um (since x is the length of the pipe which is double diameter value).

Then I am running the case for Re = 1.3 and not adding any codes to include for the interior viscosity. What I get is this after 40 ms of running:

I get a value of ~1.18 which is not even close to 1.25 reported in the paper. Am I doing something wrong here?

2. Is there any options to only produce the XMF file for the last time step?

[MaxvdKolk]

Hi @mojtabaAmiraslanpour,

So now I wonder which one is correct?

Unfortunately, I cannot comment in great detail on the differences you observe between these publications, as I am not that familiar with the details they discuss. One thing to note, there are multiple combinations of cell parameters that still result in validated behaviour, i.e. fit well between the error bars of the experimental results. So I am not sure if there is a single answer of the right value for the DI for a given shear flow condition. Additionally, the resolution---both of the LBM mesh and the RBC---can influence the reported DI values.

Can you explain this?

That bottom edge is a (known) post-processing artefact. It is related to how the domains are visualised with overlapping envelopes when running in parallel. If you place the RBC in the centre of the domain, it should be OK. Note though, that numerically the cell might not be perfectly centred, so drift is likely to remain present (maybe less pronounced).

I get a value of ~1.18 which is not even close to 1.25 reported in the paper. Am I doing something wrong here?

Not sure if you are doing something wrong, again, I am not too familiar with these papers and their configurations. Something to note regarding the cell packing, is that you might see less cells in the domain as compared to the original packing. For instance, cells that cross the simulation boundary at initialisation time are dropped (to prevent partial cells in the simulation). So, you might be running the simulation with less RBCs and PLTs than you expect.

To investigate, you might compare the RBC and PLT counts, i.e. the first line in these files (head -n 1 RBC.pos) to the reported cell counts in the logs emitted by HemoCell. I suppose that a different number of RBC will influence the used haematocrit in the domain, and thereby influence the apparent viscosity of the system.

Note though, that you are within 4% to 5% of the reported values (1.18 vs 1.25) where the crosses in your graph seem to align with the 1.2 tick mark on the right side (steady state), which might be caused by different resolutions or cell counts in the system.

Is there any options to only produce the XMF file for the last time step?

You can align the tmeas interval at which the output is written with the maximum iteration count.

[mojtabaAmiraslanpour]

Thanks,

1. The thing is when the cell drifts, it eventually reaches the boundary in x directions and is omitted. So then I would get a DI value of -100%. It happens so fast in high shear rates and won't let me see the steady-state DI value. Of course, I can scale the domain in the x-direction to delay it but it would increase the computational effort. Apparently, the cell cannot be re-entered into the domain when it is out.

2. I am trying to run oneCellShear case with interior viscosity ON and try to compile it with the relevant flag. But I get this error when compiling:

   Manually-specified variables were not used by the project:
   
   INTERIOR_VISCOSITY

   I also enabled INTERIOR_VISCOSITY in the constant_defaults.h file, but didn't help.

   #ifndef INTERIOR_VISCOSITY
   #define INTERIOR_VISCOSITY
   #endif

   Are there any workarounds for this?

[MaxvdKolk]

The thing is when the cell drifts, it eventually reaches the boundary in x directions and is omitted.

Yes, eventually the cell might drift outside of the simulation domain and will indeed be deleted if you do not opt for periodicity. You can of course, as you mentioned, increase the domain in x-direction to such an extend that the cell remains within the domain for the full simulation. However, this does become computationally inefficient when there is significant velocity in x-direction. Alternatively, you might choose to consider periodic boundary conditions along the x-direction. Doing so will ensure the cell is reinserted in the simulation domain when leaving through the periodic boundary. To enable this you can include the following (c.f. examples/parachuting.cpp)

HemoCell hemocell(argv[1], argc, argv);
// ...
hemocell.setSystemPeriodicity(0, true);
// ...

Are there any workarounds for this?

No workarounds should be required to enable interior viscosity in your (custom) example cases. If you want to include these effects, you can do so by enabling this feature in the CMakeLists.txt definition for your example. By default (assuming you started from examples/template/CMakeLists.txt), this links the example's .cpp file with the default hemocell ("${PROJECT_NAME}") library, which you might want to change to ${PROJECT_NAME}_interior_viscosity". Then, CMake will ensure your example is linked with the corresponding hemocell library including interior viscosity.

You might want to compare the following configuration for an example without and with interior viscosity:

examples/parachuting/CMakeLists.txt
examples/cellCollision_interior_viscosity/CMakeLists.txt

and the related section of our documentation compiling-examples-with-special-features.

[mojtabaAmiraslanpour]

@MaxvdKolk Just to understand, I see there is a line of code inside the cellCollision_interior_viscosity case that reads:

hemocell.setInteriorViscosityTimeScaleSeperation((*cfg)["sim"]["interiorViscosity"].read<int>(),(*cfg)["sim"]["interiorViscosityEntireGrid"].read<int>());

So don't I need to add this line to my custom example case? And also can you explain a bit about the two parameters being read here. i.e interiorViscosity and interiorViscosityEntireGrid?

I see that we set viscosityRatio in the RBC.xml file.

[mojtabaAmiraslanpour]

I am having a problem setting the Hematocrit value and getting the correct RBC.pos and PLT.pos files. I compare my results to the available .pos files in pipeflow/initial_states. I use the packCellstool to reproduce the D64_Ht21 case. So here is what I do:

packcells 128 64 64 -h 0.21 -r

And what I get is: RBC: 1135

But I should get 1832. I found that inside the packCellscode, the volume of a single RBC is put to be 97 µm^3:

double rbcVolNominal = 97.0; // This is set to match the model used in HemoCell!

I don't know what it means by the comment.

Now, from the oneCellShear case I find the reported RBC volume (at initial state) to be 81.1 µm^3. So when I use that I get 1357 which I am still far away from what I should get.

By running the pipeFlowcase with the initial_state .pos files, I get fairly good results when compared to the paper. That is why I am trying to reproduce them using the packCellstool.

[MaxvdKolk]

So don't I need to add this line to my custom example case? And also can you explain a bit about the two parameters being read here. i.e interiorViscosity and interiorViscosityEntireGrid?

I am not exactly aware on the details regarding interior viscosity in HemoCell, as I have not worked much with this so far. It seems indeed that you need to specify that line, to set the interiorViscosityTimescale accordingly. Otherwise, you might run into a warning in the logs, mentioning that you have forgotten to doing so. If I understand it correctly, this time scale separation indicates at what interval of the iterations the material models regarding the interior viscosity are evaluated and/or updated. Maybe @gzavo can give you more information on the exact consequences of the values set here, and the differences between interiorViscosity and interiorViscosityEntireGrid.

But I should get 1832. So when I use that I get 1357 which I am still far away from what I should get.

Not necessarily. First of all, the packings in those initial states have been generated with older packing tools, so there might be difference there, or due to the random initialisation that might differ between systems. Furthermore, as mentioned in a previous comment by @gzavo, there is a difference between "tank" and "tube" hematocrit. You are now packing the cells within a domain of 128 64 64 micron, i.e. a box-like domain, while the pipe flow simulations are only considered in a tube within that domain. Thus, part of the cells that are generated are dropped from the simulation, as they are packed outside the inner part of the tube (near the corners of the box). This means, that you should pack the box-like domain with a higher hematocrit target, to ensure the hematocrit within the inner part of the tube match with your desired hematocrit values.

You might want to compare the cell count of the RBCs in the box-like packing, with the cells reported by HemoCell for the tube-like domain. There should be a difference, due to the dropped cells. So, by increasing -h 0.21 to account for the tube-like domain, you should be able to get closer to your target cell count.

[mojtabaAmiraslanpour]

@MaxvdKolk Thank you,

Regarding the hematocrit value, actually it is about the removed particles on the pipe surface (cut-through particles) and not the particles near the corners (since hematocrit is an intensive property and is not dependent on the size of the domain). So I was working to get the real-time hematocrit value in the solver and started to write the code. However, I mentioned that the number of cells reported by getNumberOfCellsFromType(&hemocell, "RBC"); is not the same as the range that I am looping on: for (auto const& cell_info: CellInformationFunctionals::info_per_cell) { ... } This confused me. The idea is to compute the volume of the whole RBC's and divide it by the volume of the domain. This is actually how the packCells tool works, as I validated it using a simple excel worksheet.

Finally, I can easily set the -h value in packCells and see the result as solver runs.

[mojtabaAmiraslanpour]

Just found out that when I run the solver with one processor, I get the correct number of cells for both functions. I don't know how to achieve this when running in parallel.

[MaxvdKolk]

Regarding the hematocrit value, actually it is about the removed particles on the pipe surface (cut-through particles) and not the particles near the corners

Sure, I was hinting at the remaining domain outside of the pipe surface, towards the corners of the square cross-section.

The total cell counts should be OK when running with any number of processors, as those are gathered and reduced from all domains into a single, total cell count. However, the info_per_cell class itself is not threadsafe, as mentioned in the comments near the class defined in helper/cellInfo.h. So, when iterating on multiple threads simultaneously, you might see unexpected behaviour. If you want to perform specific operations across atomic blocks (threads) you might want to implement a specific Functional as the parallelisation is then handled by the (underlying) library.

[mojtabaAmiraslanpour]

@MaxvdKolk I am struggling to implement the new Functionalto handle the parallelization. Can you help me find a simple way to find the average volume of all RBCs when I am running in parallel? I do not want to go deep into the code at this stage.

I am just doing this to find the exact hematocrit value at each time step since I am not able to validate the pipeflowcase results. I am thinking the hematocrit value is not correct. All other parameters are the same as the paper.

[MaxvdKolk]

Would it be sufficient to use the total cell count to determine the hematocrit value for the pipe flow problem? For instance use the cell count from CellInformationFunctionals::getTotalNumberOfCells(&hemocell) and determine an estimated hematocrit value considering the pipe's volume?

Or do you aim to evaluate the hematocrit using the deformed cells' volumes?

[mojtabaAmiraslanpour]

Or do you aim to evaluate the hematocrit using the deformed cells' volumes?

I have estimated the hematocrit value by the method you mentioned and I still don't get the correct value. I actually want to include deformed cell volumes as well.