Fedosov Model Integration into HemoCell

examples/CMakeLists.txt

@@ -6,6 +6,7 @@ add_subdirectory("curvedflow_with_preinlet")
 add_subdirectory("flowaroundsphere")
 add_subdirectory("microcontraction")
 add_subdirectory("oneCellShear")
+add_subdirectory("oneCellShear_fedosov")
 add_subdirectory("parachuting")
 add_subdirectory("parallelplanes")
 add_subdirectory("pipeflow")

examples/oneCellShear_fedosov/CMakeLists.txt

@@ -0,0 +1,14 @@
+# executable will have the same name as its directory
+get_filename_component(EXEC_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+# write the resulting executable in the _current_ directory
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+
+# add executable from source files in current directory assuming
+# to add more source files, manually register them using `add_executable`
+get_filename_component(MAIN ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+add_executable(${EXEC_NAME} "${CMAKE_CURRENT_SOURCE_DIR}/${MAIN}.cpp")
+
+# link the executable to `hemocell` and `HDF5` dependencies
+target_link_libraries(${EXEC_NAME} ${PROJECT_NAME})
+target_link_libraries(${EXEC_NAME} ${HDF5_C_HL_LIBRARIES} ${HDF5_LIBRARIES})

examples/oneCellShear_fedosov/RBC.pos

@@ -0,0 +1,2 @@
+1
+9.5 9.5 4.5 90 0 0

examples/oneCellShear_fedosov/RBC.xml

@@ -0,0 +1,20 @@
+<?xml version="1.0" ?>
+<hemocell>
+<MaterialModel>
+    <comment>Parameters for the Fedosov RBC constitutive model.</comment>
+    <name>RBC</name>
+    <eta_m> 0.0 </eta_m> <!-- Membrane viscosity. [5e-10 Ns/m]-->
+    <viscosityRatio>1.0</viscosityRatio> <!-- ratio between interior and exterior viscosity -->
+    <enableInteriorViscosity>0</enableInteriorViscosity>
+    <minNumTriangles> 640 </minNumTriangles> <!--Minimun numbers of triangles per cell. Not always exact. [642]-->
+    <radius> 3.91e-6 </radius> <!-- Radius of the RBC in [ 3.96 um] -->
+    <Volume> 81 </Volume> <!-- Volume of the RBC in µm³ -->
+    <x0> 0.4545 </x0> <!-- Ratio of l_0 / l_max -->
+    <m> 2 </m> <!-- POW function exp coefficient -->
+    <mu0> 6.3e-6 </mu0> <!-- Shear modulus [N/m] -->
+    <ka> 0.29 </ka>  <!-- Global Area Constraint [N/m] -->
+    <kd> 0.01 </kd>  <!-- Local Area Constraint [N/m] -->
+    <kv> 1000.0 </kv>  <!-- Volume constraint [N/m^2] -->
+    <kb> 67.0 </kb>  <!-- * KbT, Bending coefficient [J] -->
+</MaterialModel>
+</hemocell>

examples/oneCellShear_fedosov/clean.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+echo "Warning! This will remove all files produced during the CMake build, except the binary!"
+read -p "Are you sure? [y/n]" -n 1 -r
+echo    # (optional) move to a new line
+if [[ $REPLY =~ ^[Yy]$ ]]
+then
+	echo '** Clearing palabos and hemocell compiled libraries...'
+	find ../../build/hemocell/. -type f ! -name 'CMakeLists.txt' -delete
+	echo '** Clearing local build...'
+	rm -r ./build
+fi
+

examples/oneCellShear_fedosov/compile.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+trap "exit" INT
+
+# This script invokes the compilation of the example present in the current
+# directory.
+
+echo "=========== Building =========="
+date
+
+example=${PWD}
+
+if [ ! -d "../../build" ]; then
+  echo "* Running CMake..."
+  mkdir ../../build
+  cd ../../build || exit 1
+  cmake ..
+  cd "$example" || exit 1
+fi
+
+echo "* Compiling..."
+cd ../../build || exit 1
+cmake --build . --target "${example##*/}"
+cd "$example" || exit 1
+
+date
+echo "=========== Done ==========="

examples/oneCellShear_fedosov/config.xml

@@ -0,0 +1,29 @@
+<?xml version="1.0" ?>
+<hemocell>
+
+<parameters>
+    <warmup> 0 </warmup> <!-- Number of LBM iterations to prepare fluid field. -->
+</parameters>
+
+
+<ibm>
+    <radius> 3.91e-6 </radius> <!-- Radius of the particle in [m] (dx) [3.3e-6, 3.91e-6, XX and 4.284 for shapes [0,1,2,3] respectively -->
+</ibm>
+
+<domain>
+    <shearrate> 111.0 </shearrate>   <!--Shear rate for the fluid domain. [s^-1] [25]. -->
+    <rhoP> 1025 </rhoP>   <!--Density of the surrounding fluid, Physical units [kg/m^3]-->
+    <nuP> 1.1e-6 </nuP>   <!-- Kinematic viscosity of the surrounding fluid, physical units [m^2/s]-->
+    <dx> 0.5e-6 </dx> <!--Physical length of 1 Lattice Unit -->
+    <dt> 0.5e-7 </dt> <!-- Time step for the LBM system. A negative value will set Tau=1 and calc. the corresponding time-step. -->
+    <particleEnvelope>25</particleEnvelope>
+    <kBT>4.100531391e-21</kBT> <!-- in SI, m2 kg s-2 (or J) for T=300 -->
+</domain>
+
+<sim>
+    <tmax> 100000 </tmax> <!-- total number of iterations -->
+    <tmeas> 2000 </tmeas> <!-- interval after which data is written --> 
+    <tcheckpoint>500000</tcheckpoint>
+</sim>
+
+</hemocell>

examples/oneCellShear_fedosov/oneCellShear_fedosov.cpp

@@ -0,0 +1,177 @@
+/*
+This file is part of the HemoCell library
+
+HemoCell is developed and maintained by the Computational Science Lab 
+in the University of Amsterdam. Any questions or remarks regarding this library 
+can be sent to: [email protected]
+
+When using the HemoCell library in scientific work please cite the
+corresponding paper: https://doi.org/10.3389/fphys.2017.00563
+
+The HemoCell library is free software: you can redistribute it and/or
+modify it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+The library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include "hemocell.h"
+#include "rbcFedosovModel.h"
+#include "helper/hemocellInit.hh"
+#include "helper/cellInfo.h"
+
+#include "palabos3D.h"
+#include "palabos3D.hh"
+
+using namespace hemo;
+
+int main(int argc, char* argv[])
+{   
+	if(argc < 2)
+	{
+			cout << "Usage: " << argv[0] << " <configuration.xml>" << endl;
+			return -1;
+	}
+
+	HemoCell hemocell(argv[1], argc, argv);
+	Config * cfg = hemocell.cfg;
+
+
+
+
+// ----------------- Read in config file & calc. LBM parameters ---------------------------
+	pcout << "(OneCellShear) (Parameters) calculating shear flow parameters" << endl;
+	plint nz = 10.0*(1e-6/(*cfg)["domain"]["dx"].read<T>());
+  plint nx = 2*nz;
+  plint ny = 2*nz;
+  param::lbm_shear_parameters((*cfg),ny);
+  param::printParameters();
+
+	// ------------------------ Init lattice --------------------------------
+
+	pcout << "(CellStretch) Initializing lattice: " << nx <<"x" << ny <<"x" << nz << " [lu]" << std::endl;
+
+	plint extendedEnvelopeWidth = 2;  // Because we might use ibmKernel with with 2.
+
+	hemocell.lattice = new MultiBlockLattice3D<T,DESCRIPTOR>(
+			defaultMultiBlockPolicy3D().getMultiBlockManagement(nx, ny, nz, extendedEnvelopeWidth),
+			defaultMultiBlockPolicy3D().getBlockCommunicator(),
+			defaultMultiBlockPolicy3D().getCombinedStatistics(),
+			defaultMultiBlockPolicy3D().getMultiCellAccess<T, DESCRIPTOR>(),
+			new GuoExternalForceBGKdynamics<T, DESCRIPTOR>(1.0/param::tau));
+
+	pcout << "(OneCellShear) Re corresponds to u_max = " << (param::re * param::nu_p)/(hemocell.lattice->getBoundingBox().getNy()*param::dx) << " [m/s]" << endl;
+	// -------------------------- Define boundary conditions ---------------------
+
+	OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
+			= createLocalBoundaryCondition3D<T,DESCRIPTOR>();
+
+	hemocell.lattice->toggleInternalStatistics(false);
+
+	iniLatticeSquareCouette_Y(*hemocell.lattice, nx, ny, nz, *boundaryCondition, param::shearrate_lbm);
+
+	hemocell.lattice->initialize();
+  hemocell.outputInSiUnits = true;
+
+	// ----------------------- Init cell models --------------------------
+
+	hemocell.initializeCellfield();
+	hemocell.addCellType<RbcFedosovModel>("RBC", RBC_FROM_SPHERE);
+	vector<int> outputs = {OUTPUT_POSITION,OUTPUT_TRIANGLES,OUTPUT_FORCE,OUTPUT_FORCE_VOLUME,OUTPUT_FORCE_BENDING,OUTPUT_FORCE_LINK,OUTPUT_FORCE_AREA, OUTPUT_FORCE_VISC}; 
+	hemocell.setOutputs("RBC", outputs);
+
+	outputs = {OUTPUT_VELOCITY};
+	hemocell.setFluidOutputs(outputs);
+
+  // -------------------------- Setting periodicity after initializing lattice and cell fields ---------------------
+  hemocell.setSystemPeriodicity(0, true);
+  hemocell.setSystemPeriodicity(2, true);
+
+// ---------------------- Initialise particle positions if it is not a checkpointed run ---------------
+
+	//loading the cellfield
+  if (not cfg->checkpointed) {
+    hemocell.loadParticles();
+    hemocell.writeOutput();
+  } else {
+    pcout << "(OneCellShear) CHECKPOINT found!" << endl;
+    hemocell.loadCheckPoint();
+  }
+
+
+  if (hemocell.iter == 0) { 
+    pcout << "(OneCellShear) fresh start: warming up cell-free fluid domain for "  << (*cfg)["parameters"]["warmup"].read<plint>() << " iterations..." << endl; 
+    for (plint itrt = 0; itrt < (*cfg)["parameters"]["warmup"].read<plint>(); ++itrt) {  
+      hemocell.lattice->collideAndStream();  
+    } 
+  }
+
+  pcout << "(OneCellShea) Shear rate: " << (*cfg)["domain"]["shearrate"].read<T>() << " s^-1." << endl;
+
+  unsigned int tmax = (*cfg)["sim"]["tmax"].read<unsigned int>();
+  unsigned int tmeas = (*cfg)["sim"]["tmeas"].read<unsigned int>();
+  unsigned int tcheckpoint = (*cfg)["sim"]["tcheckpoint"].read<unsigned int>();
+
+  // Get undeformed cell values
+  CellInformationFunctionals::calculateCellVolume(&hemocell);
+  CellInformationFunctionals::calculateCellArea(&hemocell);
+  T volume_eq = (CellInformationFunctionals::info_per_cell[0].volume)/pow(1e-6/param::dx,3);
+  T surface_eq = (CellInformationFunctionals::info_per_cell[0].area)/pow(1e-6/param::dx,2);
+
+  T D0 = 2.0 * (*cfg)["ibm"]["radius"].read<T>() * 1e6;
+
+  // Creating output log file
+  plb_ofstream fOut;
+  if(cfg->checkpointed)
+    fOut.open("stretch.log", std::ofstream::app);
+  else
+    fOut.open("stretch.log");
+
+
+  while (hemocell.iter < tmax ) {
+
+    hemocell.iterate();
+
+    if (hemocell.iter % tmeas == 0) {
+      hemocell.writeOutput();
+
+      // Fill up the static info structure with desired data
+      CellInformationFunctionals::calculateCellVolume(&hemocell);
+      CellInformationFunctionals::calculateCellArea(&hemocell);
+      CellInformationFunctionals::calculateCellPosition(&hemocell);
+      CellInformationFunctionals::calculateCellStretch(&hemocell);
+      CellInformationFunctionals::calculateCellBoundingBox(&hemocell);
+
+      T volume = (CellInformationFunctionals::info_per_cell[0].volume)/pow(1e-6/param::dx,3);
+      T surface = (CellInformationFunctionals::info_per_cell[0].area)/pow(1e-6/param::dx,2);
+      hemo::Array<T,3> position = CellInformationFunctionals::info_per_cell[0].position/(1e-6/param::dx);
+      hemo::Array<T,6> bbox = CellInformationFunctionals::info_per_cell[0].bbox/(1e-6/param::dx);
+      T largest_diam = (CellInformationFunctionals::info_per_cell[0].stretch)/(1e-6/param::dx);
+      T rel_D2 = (largest_diam/D0)*(largest_diam/D0);
+      T def_idx = (rel_D2 - 1.0) / (rel_D2 + 1.0) * 100.0;
+
+      pcout << "\t Cell center at: {" <<position[0]<<","<<position[1]<<","<<position[2] << "} µm" << endl;  
+      pcout << "\t Diameters: {" << bbox[1]-bbox[0] <<", " << bbox[3]-bbox[2] <<", " << bbox[5]-bbox[4] <<"}  µm" << endl;
+      pcout << "\t Surface: " << surface << " µm^2" << " (" << surface / surface_eq * 100.0 << "%)" << "  Volume: " << volume << " µm^3" << " (" << volume / volume_eq * 100.0 << "%)"<< endl;
+      pcout << "\t Largest diameter: " << largest_diam << " µm." << endl;
+      pcout << "\t Deformation index: " << def_idx << " [%]" << endl;
+
+      fOut << hemocell.iter << " " << bbox[1]-bbox[0] << " " << bbox[3]-bbox[2] << " " << bbox[5]-bbox[4] << " " << volume / volume_eq * 100.0 << " " << surface / surface_eq * 100.0 << " " << largest_diam << " " << def_idx << endl;
+
+      CellInformationFunctionals::clear_list();
+    }
+    if (hemocell.iter % tcheckpoint == 0) {
+      hemocell.saveCheckPoint();
+    }
+  }
+
+  fOut.close();
+  pcout << "(OneCellShear) Simulation finished :)" << std::endl;
+  return 0;
+}

examples/oneCellShear_fedosov/shear.gpl

@@ -0,0 +1,24 @@
+gamma = 111
+
+set multiplot layout 2,1 rowsfirst
+
+
+set xlabel "Iterations"
+set ylabel "Largest diameter [um]"
+set title "Shearrate of ".gamma." s^-1"
+set yrange [7.8:8.5]
+plot 'stretch.log' u 1:7 w l t "D"
+
+
+
+set arrow 1 from graph 0,first 100 to graph 1,first 100 nohead
+set arrow 2 from graph 0,first 111 to graph 1,first 111 nohead
+set yrange [95:115]
+set xlabel "Iterations"
+set ylabel "Percentage"
+unset title
+plot 'stretch.log' u 1:5 w l t "Volume", '' u 1:6 w l t "Surface"
+
+unset multiplot
+
+# while(1) { replot; pause 15 }

helper/hemocellInit.hh

@@ -91,4 +91,25 @@ void iniLatticeSquareCouette(plb::MultiBlockLattice3D<T,Descriptor>& lattice,
     lattice.initialize();
 }
 
+template<typename T, template<class U> class Descriptor>
+void iniLatticeSquareCouette_Y(plb::MultiBlockLattice3D<T,Descriptor>& lattice,
+                 plint nx, plint ny, plint nz,
+                 plb::OnLatticeBoundaryCondition3D<T,Descriptor>& boundaryCondition, T shearRate)
+{
+    plb::Box3D top = plb::Box3D(0, nx-1, ny-1, ny-1, 0, nz-1);
+    plb::Box3D bottom = plb::Box3D(0, nx-1, 0, 0, 0, nz-1);
+
+    boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top );
+    boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, bottom );
+
+    T vHalf = (ny-1)*shearRate*0.5;
+    setBoundaryVelocity(lattice, top, plb::Array<T,3>(vHalf,0.0,0.0));
+    setBoundaryVelocity(lattice, bottom, plb::Array<T,3>(-vHalf,0.0,0.0));
+
+    setExternalVector( lattice, lattice.getBoundingBox(),
+            Descriptor<T>::ExternalField::forceBeginsAt, plb::Array<T,Descriptor<T>::d>(0.0,0.0,0.0));
+
+    lattice.initialize();
+}
+
 #endif  // HEMOCELLINIT_HH

mechanics/cellMechanics.h

@@ -76,6 +76,36 @@ class CellMechanics {
   T calculate_etaM(Config & cfg ){
     return cfg["MaterialModel"]["eta_m"].read<T>() * param::dx / param::dt / param::df;
   };
+
+  T calculate_kd(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T kd =  cfg["MaterialModel"]["kd"].read<T>();
+    return kd * param::dx / param::df;    // converting to LBM units
+  };
+
+  T calculate_ka(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T ka =  cfg["MaterialModel"]["ka"].read<T>();
+    return ka * param::dx / param::df;    // converting to LBM units
+  };
+
+  T calculate_kv(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T kv =  cfg["MaterialModel"]["kv"].read<T>();
+    return kv * param::dx * param::dx / param::df;  // converting to LBM units
+  };
+
+  T calculate_kb(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T kb =  cfg["MaterialModel"]["kb"].read<T>();
+    return kb * param::kBT_lbm;  // converting to LBM units
+  };
+
+  T calculate_x0(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T x0 =  cfg["MaterialModel"]["x0"].read<T>();
+    return x0;
+  };
+
+  T calculate_m(Config & cfg, plb::MeshMetrics<T> & meshmetric){
+    T m =  cfg["MaterialModel"]["m"].read<T>();
+    return m;
+  };
 };
 }
 #endif