ICM
This module implements the algorithm described in [1] for PBC using cuFFT to solve the stokes operator. Fluid properties are stored in a staggered grid for improved translational invariance [2]. ICM solves the incompressible fluctuating Navier-Stokes equation and couples the resulting fluid with immersed particles particles with the Immerse Boundary Method (IBM).
Currently the case with excess mass 0 is encoded.
Use it as any other module:
Hydro::ICM::Parameters par;
par.temperature = 1.0;
par.density = 1.0;
par.viscosity = 1.0;
//The hydrodynamic radius is given by ~0.91*Lbox/cellDimension in ICM, so the hydrodynamic radius given here cannot be enforced exactly, rather the most approximate one will be computed by selecting an appropiate cellDimension (which will be an FFT wise number).
par.hydrodynamicRadius = 1.0;
//Instead of the hydrodynamic radius the cell dimensions can be provided directly. Giving rh~0.91·L/cellDim
//par.cells=make_int3(64,64,64); //cells.z == 1 means 2D
//In any case a message will be issued with the used hydrodynamic radius.
par.dt = 0.01;
par.box = Box(64);
//par.sumThermalDrift = false; //Default is true, controls if the thermal drift is taken into account or not.
auto bdhi = make_shared<BDHI::ICM>(pd, sys, par);
//add any interactor
bdhi->forwardTime();
The function getSelfMobility() will return the short time self mobility of the particles. The function getHydrodynamicRadius() will return the used hydrodynamic radius. You should expect a variation of +- 1% between this value and the one you measure from simulations.
eq. 36 in [1] can be rewritten as:
Unperturbed fluid velocity at the next time step (unperturbed means that the excess mass of the particles has not been taken into account yet, this effect is introduced later as a correction):
Where:
And g are the forces acting on the fluid.
P represents the projection into the divergence-free space (comes from removing the pressure term from the Navier-Stokes formulation)
With PBC P can be written in Fourier space as:
D: Divergence operator
G: Gradient operator
L: Laplacian operator
S: Spreading operator (transmits particle properties to fluid)
W: Symmetric noise tensor, Gaussian numbers with mean 0 and std 1
On the other hand the particles and fluid are coupled with:
Where J is the interpolation operator defined as
Where dV is the volume of a grid cell.
Only the case where the excess mass is 0 is encoded, so in this case the code takes and updates the particle positions accordingly.
The particles are updated following the mid point (predictor-corrector) scheme developed in [1].
About J and S:
Particles to fluid:
Fluid to particle:
Special care is needed due to each coordinate of the velocity being defined in a different grid (staggered grid).
Where δ is a spreading kernel (smeared delta).
The default kernel used is the 3-point Peskin kernel, see IBM_kernels.cuh.
But others can be selected with arbitrary support to better describe or satisfy certain conditions.
More information is available at [4].
The thermal drift is computed using random finite differences as explained in [3].
[1] Inertial coupling method for particles in an incompressible fluctuating fluid. Florencio Balboa, Rafael Delgado-Buscalioni, Boyce E. Griffith and Aleksandar Donev. 2014 https://doi.org/10.1016/j.cma.2013.10.029.
[2] Staggered Schemes for Fluctuating Hydrodynamics. Florencio Balboa, et.al. 2012.
[3] Brownian Dynamics without Green's Functions. Steve Delong, Florencio Balboa, et. al. 2014
[4] The Immerse Boundary Method. Charles S. Peskin 2002.