flow.h 8.97 KB
``````#ifndef STIM_FLOW_H
#define STIM_FLOW_H

#include <vector>
#include <algorithm>

//STIM include
#include <stim/math/vec3.h>
#include <stim/parser/arguments.h>
#include <stim/biomodels/network.h>

#ifdef __CUDACC__
#include <cublas_v2.h>
#include <stim/cuda/cudatools/error.h>
#endif

namespace stim {
template <typename A, typename B, typename C>
struct triple {
A first;
B second;
C third;
};

template <typename T>
struct bridge {
std::vector<unsigned> v;	// vertices' indices
std::vector<typename stim::vec3<T> > V;	// vertices' coordinates
T l;		// length
T deltaP;	// pressure drop
T Q;		// volume flow rate
};

template <typename T>
class flow {

private:

// calculate the cofactor of elemen[row][col]
void get_minor(T** src, T** dest, int row, int col, int order) {

// index of element to be copied
int rowCount = 0;
int colCount = 0;

for (int i = 0; i < order; i++) {
if (i != row) {
colCount = 0;
for (int j = 0; j < order; j++) {
// when j is not the element
if (j != col) {
dest[rowCount][colCount] = src[i][j];
colCount++;
}
}
rowCount++;
}
}
}

// calculate the det()
T determinant(T** mat, int order) {

// degenate case when n = 1
if (order == 1)
return mat[0][0];

T det = 0.0;		// determinant value

// allocate the cofactor matrix
T** minor = (T**)malloc((order - 1) * sizeof(T*));
for (int i = 0; i < order - 1; i++)
minor[i] = (T*)malloc((order - 1) * sizeof(T));

for (int i = 0; i < order; i++) {

// get minor of element(0, i)
get_minor(mat, minor, 0, i, order);

// recursion
det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1);
}

// release memory
for (int i = 0; i < order - 1; i++)
free(minor[i]);
free(minor);

return det;
}

public:
T** C;																	// Conductance
std::vector<typename stim::triple<unsigned, unsigned, float> > Q;		// volume flow rate
std::vector<T> QQ;														// Q' vector
std::vector<T> P;														// initial pressure
std::vector<T> pressure;												// final pressure

//std::vector<typename stim::triple<unsigned, unsigned, T> > V;		 // velocity
//std::vector<typename stim::triple<unsigned, unsigned, T> > Q;		 // volume flow rate
//std::vector<typename stim::triple<unsigned, unsigned, T> > deltaP; // pressure drop

flow() {}				// default constructor

void init(unsigned n_e, unsigned n_v) {

C = new T*[n_v]();
for (unsigned i = 0; i < n_v; i++) {
C[i] = new T[n_v]();
}

QQ.resize(n_v);
P.resize(n_v);
pressure.resize(n_v);

Q.resize(n_e);
}

void reset(unsigned n_v) {

for (unsigned i = 0; i < n_v; i++) {
for (unsigned j = 0; j < n_v; j++) {
C[i][j] = 0;
}
}
}

void clear(unsigned n_v) {

for (unsigned i = 0; i < n_v; i++)
delete[] C[i];
delete[] C;
}

/// Calculate the inverse of A and store the result in C
void inversion(T** A, int order, T* C) {

#ifdef __CUDACC__

// convert from double pointer to single pointer, make it flat
T* Aflat = (T*)malloc(order * order * sizeof(T));
for (unsigned i = 0; i < order; i++)
for (unsigned j = 0; j < order; j++)
Aflat[i * order + j] = A[i][j];

// create device pointer
T* d_Aflat;		// flat original matrix
T* d_Cflat;	// flat inverse matrix
T** d_A;		// put the flat original matrix into another array of pointer
T** d_C;
int *d_P;
int *d_INFO;

// allocate memory on device
HANDLE_ERROR(cudaMalloc((void**)&d_Aflat, order * order * sizeof(T)));
HANDLE_ERROR(cudaMalloc((void**)&d_Cflat, order * order * sizeof(T)));
HANDLE_ERROR(cudaMalloc((void**)&d_A, sizeof(T*)));
HANDLE_ERROR(cudaMalloc((void**)&d_C, sizeof(T*)));
HANDLE_ERROR(cudaMalloc((void**)&d_P, order * 1 * sizeof(int)));
HANDLE_ERROR(cudaMalloc((void**)&d_INFO, 1 * sizeof(int)));

// copy matrix from host to device
HANDLE_ERROR(cudaMemcpy(d_Aflat, Aflat, order * order * sizeof(T), cudaMemcpyHostToDevice));

// copy matrix from device to device
HANDLE_ERROR(cudaMemcpy(d_A, &d_Aflat, sizeof(T*), cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(d_C, &d_Cflat, sizeof(T*), cudaMemcpyHostToDevice));

// calculate the inverse of matrix based on cuBLAS
cublasHandle_t handle;
CUBLAS_HANDLE_ERROR(cublasCreate_v2(&handle));	// create cuBLAS handle object

CUBLAS_HANDLE_ERROR(cublasSgetrfBatched(handle, order, d_A, order, d_P, d_INFO, 1));

int INFO = 0;
HANDLE_ERROR(cudaMemcpy(&INFO, d_INFO, sizeof(int), cudaMemcpyDeviceToHost));
if (INFO == order)
{
std::cout << "Factorization Failed : Matrix is singular." << std::endl;
exit(1);
}

CUBLAS_HANDLE_ERROR(cublasSgetriBatched(handle, order, (const T **)d_A, order, d_P, d_C, order, d_INFO, 1));

CUBLAS_HANDLE_ERROR(cublasDestroy_v2(handle));

// copy inverse matrix from device to device
HANDLE_ERROR(cudaMemcpy(&d_Cflat, d_C, sizeof(T*), cudaMemcpyDeviceToHost));

// copy inverse matrix from device to host
HANDLE_ERROR(cudaMemcpy(C, d_Cflat, order * order * sizeof(T), cudaMemcpyDeviceToHost));

// clear up
free(Aflat);
HANDLE_ERROR(cudaFree(d_Aflat));
HANDLE_ERROR(cudaFree(d_Cflat));
HANDLE_ERROR(cudaFree(d_A));
HANDLE_ERROR(cudaFree(d_C));
HANDLE_ERROR(cudaFree(d_P));
HANDLE_ERROR(cudaFree(d_INFO));

#else
// get the determinant of a
double det = 1.0 / determinant(A, order);

// memory allocation
T* tmp = (T*)malloc((order - 1)*(order - 1) * sizeof(T));
T** minor = (T**)malloc((order - 1) * sizeof(T*));
for (int i = 0; i < order - 1; i++)
minor[i] = tmp + (i * (order - 1));

for (int j = 0; j < order; j++) {
for (int i = 0; i < order; i++) {
// get the co-factor (matrix) of A(j,i)
get_minor(A, minor, j, i, order);
C[i][j] = det * determinant(minor, order - 1);
if ((i + j) % 2 == 1)
C[i][j] = -C[i][j];
}
}

// release memory
free(tmp);
free(minor);
#endif
}
};
}

#endif

//// calculate the flow rate of 3D model(circle cross section)
//void calculate_flow_rate(unsigned e, T r) {
//	stim::triple<unsigned, unsigned, T> tmp_Q;
//	tmp_Q.first = V[e].first;			// copy the vertices information
//	tmp_Q.second = V[e].second;
//	tmp_Q.third = V[e].third * stim::PI * pow(r, 2);	// UNITS: uL/s
//	Q.push_back(tmp_Q);					// push back the volume flow rate information for every edge
//}

//// calculate the flow rate of 2D model(rectangular cross section)
//void calculate_flow_rate(unsigned e, T r, T h) {
//	stim::triple<unsigned, unsigned, T> tmp_Q;
//	tmp_Q.first = V[e].first;			// copy the vertices information
//	tmp_Q.second = V[e].second;
//	tmp_Q.third = V[e].third * h * r;					// UNITS: uL/s = mm^3/s
//	Q.push_back(tmp_Q);					// push back the volume flow rate information for every edge
//}

//// calculate the pressure drop of 3D model(circle cross section)
//void calculate_deltaP(unsigned e, T u, T l, T r) {
//	stim::triple<unsigned, unsigned, T> tmp_deltaP;
//	tmp_deltaP.first = V[e].first;		// copy the vertices information
//	tmp_deltaP.second = V[e].second;
//	tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4));		// UNITS: g/mm/s^2 = Pa
//	deltaP.push_back(tmp_deltaP);		// push back the volume flow rate information for every edge
//}

//// calculate the pressure drop of 2D model(rectangular cross section)
//void calculate_deltaP(unsigned e, T u, T l, T r, T h) {
//	stim::triple<unsigned, unsigned, T> tmp_deltaP;
//	tmp_deltaP.first = V[e].first;		// copy the vertices information
//	tmp_deltaP.second = V[e].second;
//	tmp_deltaP.third = (12 * u * l * Q[e].third) / (h * pow(r, 3));	// UNITS: g/mm/s^2 = Pa
//	deltaP.push_back(tmp_deltaP);		// push back the volume flow rate information for every edge
//}

//// better way to do this???
//// find the maximum and minimum pressure positions
//void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) {
//	std::vector<T> P(n_v, FLT_MAX);
//	// set one to reference
//	P[Q[0].first] = 0.0;
//	unsigned first = 0;
//	unsigned second = 0;
//	// calculate all the relative pressure in brute force manner
//	for (unsigned e = 0; e < n_e; e++) {
//		// assuming the obj file stores in a straight order, in other words, like swc file
//		first = Q[e].first;
//		second = Q[e].second;
//		if (P[first] != FLT_MAX)		// if pressure at start vertex is known
//			P[second] = P[first] - deltaP[e].third;
//		else if (P[second] != FLT_MAX)	// if pressure at end vertex is known
//			P[first] = P[second] + deltaP[e].third;
//	}

//	// find the maximum and minimum pressure position
//	auto m1 = std::max_element(P.begin(), P.end());		// temporarily max number
//	auto m2 = std::min_element(P.begin(), P.end());		// temporarily min number

//	max = std::distance(P.begin(), m1);
//	min = std::distance(P.begin(), m2);

//	T tmp_m = *m2;
//	// Now set the lowest pressure port to reference pressure(0.0 Pa)
//	for (unsigned i = 0; i < n_v; i++)
//		P[i] -= tmp_m;

//	for (unsigned i = 0; i < n_v; i++)
//		pressure.push_back(P[i]);
//}``````