Create 2/3d lbm solver for saturated/unsaturated granular material with GPU implementation

# Create 2/3d lbm solver for saturated/unsaturated granular material with GPU implementation
## Summary
The base class of 2/3 lbm solver for saturated/unsaturated granular material is created and implemented with GPU architecture. Particle distribution are initialized and reserved on the device, enabling parallelization of collision and streaming phases on GPU.
## Motivation
Use it to compute the multiphase LBM problem. More features can be added to the base class and promisingly expedited on GPU. Consequently, a SoA data structure will be carefully established and employed, yielding the identical statistics with the CPU-based model much faster.
## Design Detail
1. Algorithm design
The computation is mainly based on parallelism of the algorithm and forcing schemes for the multiphase LBM modelling, including the basic momentum exchanges between different phases as in conventional LBM modelling, and the phase separation realized by adding a Shan-Chen interaction force term.
2. Parallelized computation on GPU: Functions that are able to be parallelized consist of computing density/pressure field, calculating SC interation force, updating equilibrium distribution depending on the fluid velocity, and collision/streaming operators.
3. LBM base design
To speed up the iterations on GPU, memory allocation and synchronization between the host and device should be considered in the current model, requiring changes made on the certain types of variables and adding specific member functions to realize the data transfer.
3.1 New variable for GPU-based class
The variables enumerated below supplement the ones in the previous model.
Scalar type
Tq, Tdim, block_size, grid_size, (having trouble with this part),length etc.
Vector type
Note that we use flattened array to represent the distribution collison_f and streaming_f, density rho, pressure p, force force, and other significant variables which are involved in the collision and streaming. Though vector type is now supported by cuda thrust, we found no extra advantages of using such type over the flattened array.
Map type
cluster_to_solids_map, solid_to_clusters_map (not yet added to the model but planned for future)
Other type
lattice_info(templatized base class for weights and velocities)
3.2 Current functions need to be modified
Since the unified memory is accessible from CPU or GPU, the arrays can be initialized on the host and there is no need to call curand_uniform to generate random numbers on the device.
```c++
//! Reserve unified memory for arrays and initialize distribution
template <>
lbm::LBM::initialization(){
size_t nbytes_real = lx_ * ly_ * lz_ * sizeof(real);
size_t nbytes_bool = lx_ * ly_ * lz_ * sizeof(bool);
size_t nbytes_size_t = lx_ * ly_ * lz_ * sizeof(size_t);

cudaErrorCheck(cudaMallocManaged((void**)&rho, nbytes_real));
cudaErrorCheck(cudaMallocManaged((void**)&pressure, nbytes_real));
cudaErrorCheck(cudaMallocManaged((void**)&force, 3 * nbytes_real));
cudaErrorCheck(cudaMallocManaged((void**)&collision_f, Q * nbytes_real));
cudaErrorCheck(cudaMallocManaged((void**)&streaming_f, Q * nbytes_real));

for (size_t z = 0; z < lz_; z++) {
for (size_t y = 0; y < ly_; y++) {
for (size_t x = 0; x < lx_; x++) {
idx=x+y*length+z*length*length;
rho[idx]=rho0 * (1.0 + IniPerturbRate *
((real)rand() / (real)RAND_MAX - 0.5));
for (size_t i = 0; i < Tq; i++) f[idx*i] = lattice_info.weights[i] * rho[idx];
}
}
}
}
```

3.3 New functions for iteration (need discussion here)
__device__ void lbm::LBM::do_streaming();
__device__ void lbm::LBM::do_collision();
__global__ void iteration();

4. Solver design as a global function or other types(need discussion here)
Now that a class named LBM in the namespace of lbm is designed and consists of all necessary arrays which are well initialized and allocated on the device, the next step is how to transfer the iteration to the device. Here are two approaches for discussion:
(1) add two device functions, do_collision and do_streaming as member of LBM and pass the entire LBM class to the global function iteration
(2) define a global function with pointers to arrays passed in, for example, __global__ collision(* collision_f, …) ; then, define a member function lbm::LBM::collision(this->collision_f) . I personally support this approach but the use of this-> is not to my knowledge.
## Drawbacks(will be completed after our first successful or failed attempts)