Blame view

stim/biomodels/flow.h 8.97 KB
28b882ce   Jiaming Guo   new flow header
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
  #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 r;		// radii
  		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++;
ea2056c0   cherub   added flow.h
58
  				}
ea2056c0   cherub   added flow.h
59
60
  			}
  		}
ea2056c0   cherub   added flow.h
61
  
28b882ce   Jiaming Guo   new flow header
62
63
  		// calculate the det()
  		T determinant(T** mat, int order) {
ea2056c0   cherub   added flow.h
64
  
28b882ce   Jiaming Guo   new flow header
65
66
67
  			// degenate case when n = 1
  			if (order == 1)
  				return mat[0][0];
ea2056c0   cherub   added flow.h
68
  
28b882ce   Jiaming Guo   new flow header
69
  			T det = 0.0;		// determinant value
ea2056c0   cherub   added flow.h
70
  
28b882ce   Jiaming Guo   new flow header
71
72
73
74
  			// 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));
ea2056c0   cherub   added flow.h
75
  
ea2056c0   cherub   added flow.h
76
  
28b882ce   Jiaming Guo   new flow header
77
  			for (int i = 0; i < order; i++) {
ea2056c0   cherub   added flow.h
78
  
28b882ce   Jiaming Guo   new flow header
79
80
  				// get minor of element(0, i)
  				get_minor(mat, minor, 0, i, order);
ea2056c0   cherub   added flow.h
81
  
28b882ce   Jiaming Guo   new flow header
82
83
84
  				// recursion
  				det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1);
  			}
ea2056c0   cherub   added flow.h
85
  
28b882ce   Jiaming Guo   new flow header
86
87
88
89
  			// release memory
  			for (int i = 0; i < order - 1; i++)
  				free(minor[i]);
  			free(minor);
ea2056c0   cherub   added flow.h
90
  
28b882ce   Jiaming Guo   new flow header
91
  			return det;
ea2056c0   cherub   added flow.h
92
  		}
ea2056c0   cherub   added flow.h
93
  
28b882ce   Jiaming Guo   new flow header
94
  	public:
9db76bdb   Jiaming Guo   change flow object
95
96
97
98
99
  		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
28b882ce   Jiaming Guo   new flow header
100
  
9db76bdb   Jiaming Guo   change flow object
101
102
103
  		//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
ea2056c0   cherub   added flow.h
104
  
9db76bdb   Jiaming Guo   change flow object
105
  		flow() {}				// default constructor
ea2056c0   cherub   added flow.h
106
  
9db76bdb   Jiaming Guo   change flow object
107
108
109
110
111
  		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]();
28b882ce   Jiaming Guo   new flow header
112
  			}
ea2056c0   cherub   added flow.h
113
  
9db76bdb   Jiaming Guo   change flow object
114
115
116
  			QQ.resize(n_v);
  			P.resize(n_v);
  			pressure.resize(n_v);
ea2056c0   cherub   added flow.h
117
  
9db76bdb   Jiaming Guo   change flow object
118
  			Q.resize(n_e);
28b882ce   Jiaming Guo   new flow header
119
  		}
ea2056c0   cherub   added flow.h
120
  
b5c53037   Jiaming Guo   fixed memory leak
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
  		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;
  		}
  
a19e3c80   David Mayerich   matrix edits for ...
137
138
  		/// Calculate the inverse of A and store the result in C
  		void inversion(T** A, int order, T* C) {
28b882ce   Jiaming Guo   new flow header
139
140
141
142
143
144
145
146
  
  #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];
28b882ce   Jiaming Guo   new flow header
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
  
  			// 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));
ea2056c0   cherub   added flow.h
170
  
28b882ce   Jiaming Guo   new flow header
171
172
173
  			// calculate the inverse of matrix based on cuBLAS
  			cublasHandle_t handle;		
  			CUBLAS_HANDLE_ERROR(cublasCreate_v2(&handle));	// create cuBLAS handle object
ea2056c0   cherub   added flow.h
174
  
28b882ce   Jiaming Guo   new flow header
175
  			CUBLAS_HANDLE_ERROR(cublasSgetrfBatched(handle, order, d_A, order, d_P, d_INFO, 1));
ea2056c0   cherub   added flow.h
176
  
28b882ce   Jiaming Guo   new flow header
177
178
179
  			int INFO = 0;
  			HANDLE_ERROR(cudaMemcpy(&INFO, d_INFO, sizeof(int), cudaMemcpyDeviceToHost));
  			if (INFO == order)
ea2056c0   cherub   added flow.h
180
  			{
28b882ce   Jiaming Guo   new flow header
181
182
183
184
185
186
187
188
189
190
191
192
193
  				std::cout << "Factorization Failed : Matrix is singular." << std::endl;
  				cudaDeviceReset();
  				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
a19e3c80   David Mayerich   matrix edits for ...
194
  			HANDLE_ERROR(cudaMemcpy(C, d_Cflat, order * order * sizeof(T), cudaMemcpyDeviceToHost));
28b882ce   Jiaming Guo   new flow header
195
  
28b882ce   Jiaming Guo   new flow header
196
197
  			// clear up
  			free(Aflat);
28b882ce   Jiaming Guo   new flow header
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
  			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];
ea2056c0   cherub   added flow.h
222
223
  				}
  			}
28b882ce   Jiaming Guo   new flow header
224
225
226
227
228
  
  			// release memory
  			free(tmp);
  			free(minor);
  #endif
ea2056c0   cherub   added flow.h
229
  		}
28b882ce   Jiaming Guo   new flow header
230
  	};
ea2056c0   cherub   added flow.h
231
232
  }
  
9db76bdb   Jiaming Guo   change flow object
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
  #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]);
  //}