NEURON
rxd_intracellular.cpp
Go to the documentation of this file.
1 #include <../../nrnconf.h>
2 #include <stdio.h>
3 #include <assert.h>
4 #include <string.h>
5 #include "grids.h"
6 #include "rxd.h"
7 #include <nrnwrap_Python.h>
8 #include <unistd.h>
9 
10 #include <cmath>
11 
12 extern int NUM_THREADS;
13 extern TaskQueue* AllTasks;
14 extern double* states;
15 const int ICS_PREFETCH = 3;
16 
17 /*
18  * Sets the data to be used by the grids for 1D/3D hybrid models
19  */
20 extern "C" void set_hybrid_data(int64_t* num_1d_indices_per_grid,
21  int64_t* num_3d_indices_per_grid,
22  int64_t* hybrid_indices1d,
23  int64_t* hybrid_indices3d,
24  int64_t* num_3d_indices_per_1d_seg,
25  int64_t* hybrid_grid_ids,
26  double* rates,
27  double* volumes1d,
28  double* volumes3d,
29  double* dxs) {
30  Grid_node* grid;
31  int i, j, k, id;
32  int grid_id_check = 0;
33 
34  int index_ctr_1d = 0;
35  int index_ctr_3d = 0;
36 
37  int num_grid_3d_indices;
38  int num_grid_1d_indices;
39 
40  double dx;
41 
42  // loop over grids
43  for (id = 0, grid = Parallel_grids[0]; grid != NULL; grid = grid->next, id++) {
44  // if the grid we are on is the next grid in the hybrid grids
45  if (id == hybrid_grid_ids[grid_id_check]) {
46  num_grid_1d_indices = num_1d_indices_per_grid[grid_id_check];
47  num_grid_3d_indices = num_3d_indices_per_grid[grid_id_check];
48  // Allocate memory for hybrid data
49  grid->hybrid = true;
50  grid->hybrid_data->indices1d = (long*) malloc(sizeof(long) * num_grid_1d_indices);
51  grid->hybrid_data->num_3d_indices_per_1d_seg = (long*) malloc(sizeof(long) *
52  num_grid_1d_indices);
53  grid->hybrid_data->volumes1d = (double*) malloc(sizeof(double) * num_grid_1d_indices);
54 
55 
56  grid->hybrid_data->indices3d = (long*) malloc(sizeof(long) * num_grid_3d_indices);
57  grid->hybrid_data->rates = (double*) malloc(sizeof(double) * num_grid_3d_indices);
58  grid->hybrid_data->volumes3d = (double*) malloc(sizeof(double) * num_grid_3d_indices);
59 
60  dx = *dxs++;
61 
62  // Assign grid data
63  grid->hybrid_data->num_1d_indices = num_grid_1d_indices;
64  for (i = 0, k = 0; i < num_grid_1d_indices; i++, index_ctr_1d++) {
65  grid->hybrid_data->indices1d[i] = hybrid_indices1d[index_ctr_1d];
66  grid->hybrid_data->num_3d_indices_per_1d_seg[i] =
67  num_3d_indices_per_1d_seg[index_ctr_1d];
68  grid->hybrid_data->volumes1d[i] = volumes1d[index_ctr_1d];
69 
70  for (j = 0; j < num_3d_indices_per_1d_seg[index_ctr_1d]; j++, index_ctr_3d++, k++) {
71  grid->hybrid_data->indices3d[k] = hybrid_indices3d[index_ctr_3d];
72  grid->hybrid_data->rates[k] = rates[index_ctr_3d];
73  grid->hybrid_data->volumes3d[k] = volumes3d[index_ctr_3d];
74  ((ICS_Grid_node*) grid)->_ics_alphas[hybrid_indices3d[index_ctr_3d]] =
75  volumes3d[index_ctr_3d] / dx;
76  }
77  }
78  grid_id_check++;
79  }
80  }
81 }
82 
83 /* solve Ax=b where A is a diagonally dominant tridiagonal matrix
84  * N - length of the matrix
85  * l_diag - pointer to the lower diagonal (length N-1)
86  * diag - pointer to diagonal (length N)
87  * u_diag - pointer to the upper diagonal (length N-1)
88  * B - pointer to the RHS (length N)
89  * The solution (x) will be stored in B.
90  * l_diag, diag and u_diag are not changed.
91  * c - scratch pad, preallocated memory for (N-1) doubles
92  */
93 static int solve_dd_tridiag(int N,
94  const double* l_diag,
95  const double* diag,
96  const double* u_diag,
97  double* b,
98  double* c) {
99  int i;
100 
101  c[0] = u_diag[0] / diag[0];
102  b[0] = b[0] / diag[0];
103 
104  for (i = 1; i < N - 1; i++) {
105  c[i] = u_diag[i] / (diag[i] - l_diag[i - 1] * c[i - 1]);
106  b[i] = (b[i] - l_diag[i - 1] * b[i - 1]) / (diag[i] - l_diag[i - 1] * c[i - 1]);
107  }
108  b[N - 1] = (b[N - 1] - l_diag[N - 2] * b[N - 2]) / (diag[N - 1] - l_diag[N - 2] * c[N - 2]);
109 
110 
111  /*back substitution*/
112  for (i = N - 2; i >= 0; i--) {
113  b[i] = b[i] - c[i] * b[i + 1];
114  }
115  return 0;
116 }
117 
118 // Homogeneous diffusion coefficient
119 void ics_find_deltas(long start,
120  long stop,
121  long node_start,
122  double* delta,
123  long* line_defs,
124  long* ordered_nodes,
125  double* states,
126  double dc,
127  double* alphas) {
128  long ordered_index = node_start;
129  long current_index = -1;
130  long next_index = -1;
131  double prev_state;
132  double current_state;
133  double next_state;
134  double prev_alpha;
135  double current_alpha;
136  double next_alpha;
137  for (int i = start; i < stop - 1; i += 2) {
138  long line_start = ordered_nodes[ordered_index];
139  int line_length = line_defs[i + 1];
140  if (line_length > 1) {
141  current_index = line_start;
142  ordered_index++;
143  next_index = ordered_nodes[ordered_index];
144  current_state = states[current_index];
145  next_state = states[next_index];
146  current_alpha = alphas[current_index];
147  next_alpha = alphas[next_index];
148  delta[current_index] = dc * next_alpha * current_alpha * (next_state - current_state) /
149  (next_alpha + current_alpha);
150  ordered_index++;
151  for (int j = 1; j < line_length - 1; j++) {
152  prev_state = current_state;
153  current_index = next_index;
154  next_index = ordered_nodes[ordered_index];
155  current_state = next_state;
156  next_state = states[next_index];
157  prev_alpha = current_alpha;
158  current_alpha = next_alpha;
159  next_alpha = alphas[next_index];
160  delta[current_index] =
161  dc * ((next_alpha * current_alpha / (next_alpha + current_alpha)) *
162  (next_state - current_state) -
163  (prev_alpha * current_alpha / (prev_alpha + current_alpha)) *
164  (current_state - prev_state));
165  ordered_index++;
166  }
167  // Here next_state is actually current_state and current_state is prev_state
168  delta[next_index] = dc * (next_alpha * current_alpha) * (current_state - next_state) /
169  (next_alpha + current_alpha);
170  } else {
171  ordered_index++;
172  delta[line_start] = 0.0;
173  }
174  }
175 }
176 
177 // Inhomogeneous diffusion coefficient
178 void ics_find_deltas(long start,
179  long stop,
180  long node_start,
181  double* delta,
182  long* line_defs,
183  long* ordered_nodes,
184  double* states,
185  double* dc,
186  double* alphas) {
187  long ordered_index = node_start;
188  long current_index = -1;
189  long next_index = -1;
190  double prev_state;
191  double current_state;
192  double next_state;
193  double prev_alpha;
194  double current_alpha;
195  double next_alpha;
196  for (int i = start; i < stop - 1; i += 2) {
197  long line_start = ordered_nodes[ordered_index];
198  int line_length = line_defs[i + 1];
199  if (line_length > 1) {
200  current_index = line_start;
201  ordered_index++;
202  next_index = ordered_nodes[ordered_index];
203  current_state = states[current_index];
204  next_state = states[next_index];
205  current_alpha = alphas[current_index];
206  next_alpha = alphas[next_index];
207  delta[current_index] = dc[next_index] * next_alpha * current_alpha *
208  (next_state - current_state) / (next_alpha + current_alpha);
209  ordered_index++;
210  for (int j = 1; j < line_length - 1; j++) {
211  prev_state = current_state;
212  current_index = next_index;
213  next_index = ordered_nodes[ordered_index];
214  current_state = next_state;
215  next_state = states[next_index];
216  prev_alpha = current_alpha;
217  current_alpha = next_alpha;
218  next_alpha = alphas[next_index];
219  delta[current_index] =
220  dc[next_index] * (next_alpha * current_alpha * (next_state - current_state) /
221  (next_alpha + current_alpha)) -
222  dc[current_index] * (prev_alpha * current_alpha * (current_state - prev_state) /
223  (prev_alpha + current_alpha));
224  ordered_index++;
225  }
226  // Here next_state is actually current_state and current_state is prev_state
227  delta[next_index] = dc[next_index] * (next_alpha * current_alpha) *
228  (current_state - next_state) / (next_alpha + current_alpha);
229  } else {
230  ordered_index++;
231  delta[line_start] = 0.0;
232  }
233  }
234 }
235 
236 static void* do_ics_deltas(void* dataptr) {
237  ICSAdiGridData* data = (ICSAdiGridData*) dataptr;
238  ICSAdiDirection* ics_adi_dir = data->ics_adi_dir;
239  ICS_Grid_node* g = data->g;
240  int line_start = data->line_start;
241  int line_stop = data->line_stop;
242  int node_start = data->ordered_start;
243  double* states = g->states;
244  double* deltas = ics_adi_dir->deltas;
245  long* line_defs = ics_adi_dir->ordered_line_defs;
246  long* ordered_nodes = ics_adi_dir->ordered_nodes;
247  double* alphas = g->_ics_alphas;
248 
249  if (ics_adi_dir->dcgrid == NULL)
250  ics_find_deltas(line_start,
251  line_stop,
252  node_start,
253  deltas,
254  line_defs,
255  ordered_nodes,
256  states,
257  ics_adi_dir->dc,
258  alphas);
259  else
260  ics_find_deltas(line_start,
261  line_stop,
262  node_start,
263  deltas,
264  line_defs,
265  ordered_nodes,
266  states,
267  ics_adi_dir->dcgrid,
268  alphas);
269 
270 
271  return NULL;
272 }
273 
274 void run_threaded_deltas(ICS_Grid_node* g, ICSAdiDirection* ics_adi_dir) {
275  int i;
276  for (i = 0; i < NUM_THREADS; i++) {
277  g->ics_tasks[i].line_start = ics_adi_dir->line_start_stop_indices[2 * i];
278  g->ics_tasks[i].line_stop = ics_adi_dir->line_start_stop_indices[(2 * i) + 1];
279  g->ics_tasks[i].ordered_start =
280  ics_adi_dir->ordered_start_stop_indices[(2 * i)]; // Change what I'm storing in
281  // ordered_start_stop_indices so
282  // index is just i
283  g->ics_tasks[i].ics_adi_dir = ics_adi_dir;
284  }
285  /* launch threads */
286  for (i = 0; i < NUM_THREADS - 1; i++) {
287  TaskQueue_add_task(AllTasks, &do_ics_deltas, &(g->ics_tasks[i]), NULL);
288  }
289  /* run one task in the main thread */
290  do_ics_deltas(&(g->ics_tasks[NUM_THREADS - 1]));
291  /* wait for them to finish */
292  TaskQueue_sync(AllTasks);
293 }
294 
295 // Inhomogeneous diffusion coefficient
296 void ics_dg_adi_x_inhom(ICS_Grid_node* g,
297  int line_start,
298  int line_stop,
299  int node_start,
300  double,
301  double* states,
302  double* RHS,
303  double* scratchpad,
304  double* u_diag,
305  double* diag,
306  double* l_diag) {
307  double* delta_x = g->ics_adi_dir_x->deltas;
308  double* delta_y = g->ics_adi_dir_y->deltas;
309  double* delta_z = g->ics_adi_dir_z->deltas;
310  double* states_cur = g->states_cur;
311  double* alphas = g->_ics_alphas;
312  double* dc = g->ics_adi_dir_x->dcgrid;
313  double dx = g->ics_adi_dir_x->d;
314  double dy = g->ics_adi_dir_y->d;
315  double dz = g->ics_adi_dir_z->d;
316  double dt = *dt_ptr;
317  long next_index = -1;
318  long prev_index = -1;
319  double next;
320  double prev;
321 
322  long* x_lines = g->ics_adi_dir_x->ordered_line_defs;
323  long* ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
324 
325  long current_index;
326  long ordered_index = node_start;
327  for (int i = line_start; i < line_stop - 1; i += 2) {
328  long N = x_lines[i + 1];
329  long ordered_index_start = ordered_index;
330  for (int j = 0; j < N; j++) {
331  current_index = ordered_nodes[ordered_index];
332  RHS[j] = (dt / alphas[current_index]) *
333  (delta_x[current_index] / SQ(dx) + 2.0 * delta_y[current_index] / SQ(dy) +
334  2.0 * delta_z[current_index] / SQ(dz)) +
335  states[current_index] + states_cur[current_index];
336  ordered_index++;
337  }
338 
339  ordered_index = ordered_index_start;
340  current_index = ordered_nodes[ordered_index];
341  ordered_index++;
342  next_index = ordered_nodes[ordered_index];
343  next = alphas[next_index] * dc[next_index] / (alphas[next_index] + alphas[current_index]);
344  diag[0] = 1.0 + dt * next / SQ(dx);
345  u_diag[0] = -dt * next / SQ(dx);
346  ordered_index++;
347  for (int c = 1; c < N - 1; c++) {
348  prev_index = current_index;
349  current_index = next_index;
350  next_index = ordered_nodes[ordered_index];
351  prev = alphas[prev_index] * dc[current_index] /
352  (alphas[prev_index] + alphas[current_index]);
353  next = alphas[next_index] * dc[next_index] /
354  (alphas[next_index] + alphas[current_index]);
355  l_diag[c - 1] = -dt * prev / SQ(dx);
356  diag[c] = 1. + dt * (prev + next) / SQ(dx);
357  u_diag[c] = -dt * next / SQ(dx);
358  ordered_index++;
359  }
360  prev = alphas[current_index] * dc[next_index] /
361  (alphas[current_index] + alphas[next_index]);
362  diag[N - 1] = 1.0 + dt * prev / SQ(dx);
363  l_diag[N - 2] = -dt * prev / SQ(dx);
364 
365 
366  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
367 
368  ordered_index = ordered_index_start;
369  for (int k = 0; k < N; k++) {
370  current_index = ordered_nodes[ordered_index];
371  states[current_index] = RHS[k];
372  ordered_index++;
373  }
374  }
375 }
376 
377 void ics_dg_adi_y_inhom(ICS_Grid_node* g,
378  int line_start,
379  int line_stop,
380  int node_start,
381  double,
382  double* states,
383  double* RHS,
384  double* scratchpad,
385  double* u_diag,
386  double* diag,
387  double* l_diag) {
388  double* delta = g->ics_adi_dir_y->deltas;
389  long* lines = g->ics_adi_dir_y->ordered_line_defs;
390  double* alphas = g->_ics_alphas;
391  double* dc = g->ics_adi_dir_y->dcgrid;
392  double dy = g->ics_adi_dir_y->d;
393  double dt = *dt_ptr;
394  long next_index = -1;
395  long prev_index = -1;
396  double next;
397  double prev;
398 
399  long current_index;
400  long* ordered_y_nodes = g->ics_adi_dir_y->ordered_nodes;
401  long ordered_index = node_start;
402 
403  for (int i = line_start; i < line_stop - 1; i += 2) {
404  long N = lines[i + 1];
405  long ordered_index_start = ordered_index;
406 
407  for (int j = 0; j < N; j++) {
408  current_index = ordered_y_nodes[ordered_index];
409  RHS[j] = states[current_index] -
410  dt * delta[current_index] / (alphas[current_index] * SQ(dy));
411  ordered_index++;
412  }
413 
414  ordered_index = ordered_index_start;
415  current_index = ordered_y_nodes[ordered_index];
416  ordered_index++;
417  next_index = ordered_y_nodes[ordered_index];
418  next = alphas[next_index] * dc[next_index] / (alphas[next_index] + alphas[current_index]);
419  diag[0] = 1.0 + dt * next / SQ(dy);
420  u_diag[0] = -dt * next / SQ(dy);
421  ordered_index++;
422  for (int c = 1; c < N - 1; c++) {
423  prev_index = current_index;
424  current_index = next_index;
425  next_index = ordered_y_nodes[ordered_index];
426  prev = alphas[prev_index] * dc[prev_index] /
427  (alphas[prev_index] + alphas[current_index]);
428  next = alphas[next_index] * dc[next_index] /
429  (alphas[next_index] + alphas[current_index]);
430  l_diag[c - 1] = -dt * prev / SQ(dy);
431  diag[c] = 1. + dt * (prev + next) / SQ(dy);
432  u_diag[c] = -dt * next / SQ(dy);
433  ordered_index++;
434  }
435  prev = alphas[current_index] * dc[current_index] /
436  (alphas[current_index] + alphas[next_index]);
437  diag[N - 1] = 1.0 + dt * prev / SQ(dy);
438  l_diag[N - 2] = -dt * prev / SQ(dy);
439 
440  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
441 
442  ordered_index = ordered_index_start;
443  for (int k = 0; k < N; k++) {
444  current_index = ordered_y_nodes[ordered_index];
445  states[current_index] = RHS[k];
446  ordered_index++;
447  }
448  }
449 }
450 void ics_dg_adi_z_inhom(ICS_Grid_node* g,
451  int line_start,
452  int line_stop,
453  int node_start,
454  double,
455  double* states,
456  double* RHS,
457  double* scratchpad,
458  double* u_diag,
459  double* diag,
460  double* l_diag) {
461  double* delta = g->ics_adi_dir_z->deltas;
462  long* lines = g->ics_adi_dir_z->ordered_line_defs;
463  long* ordered_z_nodes = g->ics_adi_dir_z->ordered_nodes;
464  double* alphas = g->_ics_alphas;
465  double* dc = g->ics_adi_dir_z->dcgrid;
466  double dz = g->ics_adi_dir_z->d;
467  double dt = *dt_ptr;
468  long next_index = -1;
469  long prev_index = -1;
470  double next;
471  double prev;
472 
473  long current_index;
474  long ordered_index = node_start;
475  for (int i = line_start; i < line_stop - 1; i += 2) {
476  long N = lines[i + 1];
477  long ordered_index_start = ordered_index;
478 
479  for (int j = 0; j < N; j++) {
480  current_index = ordered_z_nodes[ordered_index];
481  RHS[j] = states[current_index] -
482  dt * delta[current_index] / (SQ(dz) * alphas[current_index]);
483  ordered_index++;
484  }
485 
486  ordered_index = ordered_index_start;
487  current_index = ordered_z_nodes[ordered_index];
488  ordered_index++;
489  next_index = ordered_z_nodes[ordered_index];
490  next = alphas[next_index] * dc[next_index] / (alphas[next_index] + alphas[current_index]);
491  diag[0] = 1.0 + dt * next / SQ(dz);
492  u_diag[0] = -dt * next / SQ(dz);
493  ordered_index++;
494  for (int c = 1; c < N - 1; c++) {
495  prev_index = current_index;
496  current_index = next_index;
497  next_index = ordered_z_nodes[ordered_index];
498  prev = alphas[prev_index] * dc[prev_index] /
499  (alphas[prev_index] + alphas[current_index]);
500  next = alphas[next_index] * dc[next_index] /
501  (alphas[next_index] + alphas[current_index]);
502  l_diag[c - 1] = -dt * prev / SQ(dz);
503  diag[c] = 1. + dt * (prev + next) / SQ(dz);
504  u_diag[c] = -dt * next / SQ(dz);
505  ordered_index++;
506  }
507  prev = alphas[current_index] * dc[current_index] /
508  (alphas[current_index] + alphas[next_index]);
509  diag[N - 1] = 1.0 + dt * prev / SQ(dz);
510  l_diag[N - 2] = -dt * prev / SQ(dz);
511  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
512 
513  ordered_index = ordered_index_start;
514  for (int k = 0; k < N; k++) {
515  current_index = ordered_z_nodes[ordered_index];
516  states[current_index] = RHS[k];
517  ordered_index++;
518  }
519  }
520 }
521 
522 // Homogeneous diffusion coefficient
523 void ics_dg_adi_x(ICS_Grid_node* g,
524  int line_start,
525  int line_stop,
526  int node_start,
527  double,
528  double* states,
529  double* RHS,
530  double* scratchpad,
531  double* u_diag,
532  double* diag,
533  double* l_diag) {
534  double* delta_x = g->ics_adi_dir_x->deltas;
535  double* delta_y = g->ics_adi_dir_y->deltas;
536  double* delta_z = g->ics_adi_dir_z->deltas;
537  double* states_cur = g->states_cur;
538  double* alphas = g->_ics_alphas;
539  double dc = g->ics_adi_dir_x->dc;
540  double dx = g->ics_adi_dir_x->d;
541  double dy = g->ics_adi_dir_y->d;
542  double dz = g->ics_adi_dir_z->d;
543  double dt = *dt_ptr;
544  long next_index = -1;
545  long prev_index = -1;
546  double next;
547  double prev;
548 
549  long* x_lines = g->ics_adi_dir_x->ordered_line_defs;
550  long* ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
551 
552  long current_index;
553  long ordered_index = node_start;
554  for (int i = line_start; i < line_stop - 1; i += 2) {
555  long N = x_lines[i + 1];
556  long ordered_index_start = ordered_index;
557  for (int j = 0; j < N; j++) {
558  current_index = ordered_nodes[ordered_index];
559  RHS[j] = (dt / alphas[current_index]) *
560  (delta_x[current_index] / (SQ(dx)) + 2 * delta_y[current_index] / SQ(dy) +
561  2 * delta_z[current_index] / SQ(dz)) +
562  states[current_index] + states_cur[current_index];
563  ordered_index++;
564  }
565 
566  ordered_index = ordered_index_start;
567  current_index = ordered_nodes[ordered_index];
568  ordered_index++;
569  next_index = ordered_nodes[ordered_index];
570  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
571  diag[0] = 1.0 + dt * next / SQ(dx);
572  u_diag[0] = -dt * next / SQ(dx);
573  ordered_index++;
574  for (int c = 1; c < N - 1; c++) {
575  prev_index = current_index;
576  current_index = next_index;
577  next_index = ordered_nodes[ordered_index];
578  prev = alphas[prev_index] * dc / (alphas[prev_index] + alphas[current_index]);
579  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
580  l_diag[c - 1] = -dt * prev / SQ(dx);
581  diag[c] = 1. + dt * (prev + next) / SQ(dx);
582  u_diag[c] = -dt * next / SQ(dx);
583  ordered_index++;
584  }
585  prev = alphas[current_index] * dc / (alphas[current_index] + alphas[next_index]);
586  diag[N - 1] = 1.0 + dt * prev / SQ(dx);
587  l_diag[N - 2] = -dt * prev / SQ(dx);
588 
589 
590  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
591 
592  ordered_index = ordered_index_start;
593  for (int k = 0; k < N; k++) {
594  current_index = ordered_nodes[ordered_index];
595  states[current_index] = RHS[k];
596  ordered_index++;
597  }
598  }
599 }
600 
601 void ics_dg_adi_y(ICS_Grid_node* g,
602  int line_start,
603  int line_stop,
604  int node_start,
605  double,
606  double* states,
607  double* RHS,
608  double* scratchpad,
609  double* u_diag,
610  double* diag,
611  double* l_diag) {
612  double* delta = g->ics_adi_dir_y->deltas;
613  long* lines = g->ics_adi_dir_y->ordered_line_defs;
614  double* alphas = g->_ics_alphas;
615  double dc = g->ics_adi_dir_y->dc;
616  double dy = g->ics_adi_dir_y->d;
617  double dt = *dt_ptr;
618  long next_index = -1;
619  long prev_index = -1;
620  double next;
621  double prev;
622 
623  long current_index;
624  long* ordered_y_nodes = g->ics_adi_dir_y->ordered_nodes;
625  long ordered_index = node_start;
626 
627  for (int i = line_start; i < line_stop - 1; i += 2) {
628  long N = lines[i + 1];
629  long ordered_index_start = ordered_index;
630 
631  for (int j = 0; j < N; j++) {
632  current_index = ordered_y_nodes[ordered_index];
633  RHS[j] = states[current_index] -
634  dt * delta[current_index] / (alphas[current_index] * SQ(dy));
635  ordered_index++;
636  }
637 
638  ordered_index = ordered_index_start;
639  current_index = ordered_y_nodes[ordered_index];
640  ordered_index++;
641  next_index = ordered_y_nodes[ordered_index];
642  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
643  diag[0] = 1.0 + dt * next / SQ(dy);
644  u_diag[0] = -dt * next / SQ(dy);
645  ordered_index++;
646  for (int c = 1; c < N - 1; c++) {
647  prev_index = current_index;
648  current_index = next_index;
649  next_index = ordered_y_nodes[ordered_index];
650  prev = alphas[prev_index] * dc / (alphas[prev_index] + alphas[current_index]);
651  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
652  l_diag[c - 1] = -dt * prev / SQ(dy);
653  diag[c] = 1. + dt * (prev + next) / SQ(dy);
654  u_diag[c] = -dt * next / SQ(dy);
655  ordered_index++;
656  }
657  prev = alphas[current_index] * dc / (alphas[current_index] + alphas[next_index]);
658  diag[N - 1] = 1.0 + dt * prev / SQ(dy);
659  l_diag[N - 2] = -dt * prev / SQ(dy);
660 
661  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
662 
663  ordered_index = ordered_index_start;
664  for (int k = 0; k < N; k++) {
665  current_index = ordered_y_nodes[ordered_index];
666  states[current_index] = RHS[k];
667  ordered_index++;
668  }
669  }
670 }
671 void ics_dg_adi_z(ICS_Grid_node* g,
672  int line_start,
673  int line_stop,
674  int node_start,
675  double,
676  double* states,
677  double* RHS,
678  double* scratchpad,
679  double* u_diag,
680  double* diag,
681  double* l_diag) {
682  double* delta = g->ics_adi_dir_z->deltas;
683  long* lines = g->ics_adi_dir_z->ordered_line_defs;
684  long* ordered_z_nodes = g->ics_adi_dir_z->ordered_nodes;
685  double* alphas = g->_ics_alphas;
686  double dc = g->ics_adi_dir_z->dc;
687  double dz = g->ics_adi_dir_z->d;
688  double dt = *dt_ptr;
689  long next_index = -1;
690  long prev_index = -1;
691  double next;
692  double prev;
693 
694  long current_index;
695  long ordered_index = node_start;
696  for (int i = line_start; i < line_stop - 1; i += 2) {
697  long N = lines[i + 1];
698  long ordered_index_start = ordered_index;
699 
700  for (int j = 0; j < N; j++) {
701  current_index = ordered_z_nodes[ordered_index];
702  RHS[j] = states[current_index] -
703  dt * delta[current_index] / (SQ(dz) * alphas[current_index]);
704  ordered_index++;
705  }
706 
707  ordered_index = ordered_index_start;
708  current_index = ordered_z_nodes[ordered_index];
709  ordered_index++;
710  next_index = ordered_z_nodes[ordered_index];
711  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
712  diag[0] = 1.0 + dt * next / SQ(dz);
713  u_diag[0] = -dt * next / SQ(dz);
714  ordered_index++;
715  for (int c = 1; c < N - 1; c++) {
716  prev_index = current_index;
717  current_index = next_index;
718  next_index = ordered_z_nodes[ordered_index];
719  prev = alphas[prev_index] * dc / (alphas[prev_index] + alphas[current_index]);
720  next = alphas[next_index] * dc / (alphas[next_index] + alphas[current_index]);
721  l_diag[c - 1] = -dt * prev / SQ(dz);
722  diag[c] = 1. + dt * (prev + next) / SQ(dz);
723  u_diag[c] = -dt * next / SQ(dz);
724  ordered_index++;
725  }
726  prev = alphas[current_index] * dc / (alphas[current_index] + alphas[next_index]);
727  diag[N - 1] = 1.0 + dt * prev / SQ(dz);
728  l_diag[N - 2] = -dt * prev / SQ(dz);
729  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
730 
731  ordered_index = ordered_index_start;
732  for (int k = 0; k < N; k++) {
733  current_index = ordered_z_nodes[ordered_index];
734  states[current_index] = RHS[k];
735  ordered_index++;
736  }
737  }
738 }
739 
740 void* do_ics_dg_adi(void* dataptr) {
741  ICSAdiGridData* data = (ICSAdiGridData*) dataptr;
742  ICSAdiDirection* ics_adi_dir = data->ics_adi_dir;
743  ICS_Grid_node* g = data->g;
744  void (*ics_dg_adi_dir)(ICS_Grid_node * g,
745  int,
746  int,
747  int,
748  double,
749  double*,
750  double*,
751  double*,
752  double*,
753  double*,
754  double*) = ics_adi_dir->ics_dg_adi_dir;
755  int line_start = data->line_start;
756  int line_stop = data->line_stop;
757  int node_start = data->ordered_start;
758  double dt = *dt_ptr;
759  double r = dt / (ics_adi_dir->d * ics_adi_dir->d);
760  ics_dg_adi_dir(g,
761  line_start,
762  line_stop,
763  node_start,
764  r,
765  g->states,
766  data->RHS,
767  data->scratchpad,
768  data->u_diag,
769  data->diag,
770  data->l_diag);
771  return NULL;
772 }
773 
774 void ICS_Grid_node::run_threaded_ics_dg_adi(ICSAdiDirection* ics_adi_dir) {
775  int i;
776  for (i = 0; i < NUM_THREADS; i++) {
777  ics_tasks[i].line_start = ics_adi_dir->line_start_stop_indices[2 * i];
778  ics_tasks[i].line_stop = ics_adi_dir->line_start_stop_indices[(2 * i) + 1];
779  ics_tasks[i].ordered_start =
780  ics_adi_dir->ordered_start_stop_indices[(2 * i)]; // Change what I'm storing in
781  // ordered_start_stop_indices so
782  // index is just i
783  ics_tasks[i].ics_adi_dir = ics_adi_dir;
784  }
785  /* launch threads */
786  for (i = 0; i < NUM_THREADS - 1; i++) {
787  TaskQueue_add_task(AllTasks, &do_ics_dg_adi, &ics_tasks[i], NULL);
788  }
789  /* run one task in the main thread */
790  do_ics_dg_adi(&ics_tasks[NUM_THREADS - 1]);
791  /* wait for them to finish */
792  TaskQueue_sync(AllTasks);
793 }
794 
795 
796 /* ------- Variable Step ICS Code ------- */
797 static inline double flux(const double alphai,
798  const double alphaj,
799  const double statei,
800  const double statej) {
801  return 2.0 * alphai * alphaj * (statei - statej) / ((alphai + alphaj));
802 }
803 
804 static void variable_step_delta(long start,
805  long stop,
806  long node_start,
807  double* ydot,
808  long* line_defs,
809  long* ordered_nodes,
810  double const* const states,
811  double r,
812  double* alphas) {
813  long ordered_index = node_start;
814  long current_index = -1;
815  long next_index = -1;
816  double prev_state;
817  double current_state;
818  double next_state;
819  double prev_alpha;
820  double current_alpha;
821  double next_alpha;
822  for (int i = start; i < stop - 1; i += 2) {
823  long line_start = ordered_nodes[ordered_index];
824  int line_length = line_defs[i + 1];
825  if (line_length > 1) {
826  current_index = line_start;
827  ordered_index++;
828  next_index = ordered_nodes[ordered_index];
829  current_state = states[current_index];
830  next_state = states[next_index];
831  current_alpha = alphas[current_index];
832  next_alpha = alphas[next_index];
833  ydot[current_index] += (r / current_alpha) *
834  flux(next_alpha, current_alpha, next_state, current_state);
835  ordered_index++;
836  for (int j = 1; j < line_length - 1; j++) {
837  prev_state = current_state;
838  current_index = next_index;
839  next_index = ordered_nodes[ordered_index];
840  current_state = next_state;
841  next_state = states[next_index];
842  prev_alpha = current_alpha;
843  current_alpha = next_alpha;
844  next_alpha = alphas[next_index];
845  ydot[current_index] += (r / current_alpha) *
846  (flux(prev_alpha, current_alpha, prev_state, current_state) +
847  flux(next_alpha, current_alpha, next_state, current_state));
848  ordered_index++;
849  }
850  ydot[next_index] += r * flux(current_alpha, next_alpha, current_state, next_state) /
851  next_alpha;
852  } else {
853  ordered_index++;
854  // ydot[line_start] += 0.0;
855  }
856  }
857 }
858 
859 static void variable_step_delta(long start,
860  long stop,
861  long node_start,
862  double* ydot,
863  long* line_defs,
864  long* ordered_nodes,
865  double const* const states,
866  double r,
867  double* dcs,
868  double* alphas) {
869  long ordered_index = node_start;
870  long current_index = -1;
871  long next_index = -1;
872  double prev_state;
873  double current_state;
874  double next_state;
875  double prev_alpha;
876  double current_alpha;
877  double next_alpha;
878  double current_dc, next_dc;
879  for (int i = start; i < stop - 1; i += 2) {
880  long line_start = ordered_nodes[ordered_index];
881  int line_length = line_defs[i + 1];
882  if (line_length > 1) {
883  current_index = line_start;
884  ordered_index++;
885  next_index = ordered_nodes[ordered_index];
886  current_state = states[current_index];
887  next_state = states[next_index];
888  current_alpha = alphas[current_index];
889  next_alpha = alphas[next_index];
890  current_dc = dcs[current_index];
891  next_dc = dcs[next_index];
892  ydot[current_index] += (r / current_alpha) * next_dc *
893  flux(next_alpha, current_alpha, next_state, current_state);
894  ordered_index++;
895  for (int j = 1; j < line_length - 1; j++) {
896  prev_state = current_state;
897  current_index = next_index;
898  next_index = ordered_nodes[ordered_index];
899  current_state = next_state;
900  next_state = states[next_index];
901  prev_alpha = current_alpha;
902  current_alpha = next_alpha;
903  next_alpha = alphas[next_index];
904  current_dc = next_dc;
905  next_dc = dcs[next_index];
906  ydot[current_index] +=
907  (r / current_alpha) *
908  (current_dc * flux(prev_alpha, current_alpha, prev_state, current_state) +
909  next_dc * flux(next_alpha, current_alpha, next_state, current_state));
910  ordered_index++;
911  }
912  ydot[next_index] += r * next_dc *
913  flux(current_alpha, next_alpha, current_state, next_state) /
914  next_alpha;
915  } else {
916  ordered_index++;
917  // ydot[line_start] += 0.0;
918  }
919  }
920 }
921 void _ics_rhs_variable_step_helper(ICS_Grid_node* g, double const* const states, double* ydot) {
922  double rate_x, rate_y, rate_z;
923  // Find rate for each direction
924  double dx = g->ics_adi_dir_x->d, dy = g->ics_adi_dir_y->d, dz = g->ics_adi_dir_z->d;
925 
926  long x_line_start = g->ics_adi_dir_x->line_start_stop_indices[0];
927  long x_line_stop = g->ics_adi_dir_x->line_start_stop_indices[NUM_THREADS * 2 - 1];
928  long x_node_start = g->ics_adi_dir_x->ordered_start_stop_indices[0];
929  long* x_line_defs = g->ics_adi_dir_x->ordered_line_defs;
930  long* x_ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
931 
932  long y_line_start = g->ics_adi_dir_y->line_start_stop_indices[0];
933  long y_line_stop = g->ics_adi_dir_y->line_start_stop_indices[NUM_THREADS * 2 - 1];
934  long y_node_start = g->ics_adi_dir_y->ordered_start_stop_indices[0];
935  long* y_line_defs = g->ics_adi_dir_y->ordered_line_defs;
936  long* y_ordered_nodes = g->ics_adi_dir_y->ordered_nodes;
937 
938  long z_line_start = g->ics_adi_dir_z->line_start_stop_indices[0];
939  long z_line_stop = g->ics_adi_dir_z->line_start_stop_indices[NUM_THREADS * 2 - 1];
940  long z_node_start = g->ics_adi_dir_z->ordered_start_stop_indices[0];
941  long* z_line_defs = g->ics_adi_dir_z->ordered_line_defs;
942  long* z_ordered_nodes = g->ics_adi_dir_z->ordered_nodes;
943 
944  double* alphas = g->_ics_alphas;
945 
946  if (g->ics_adi_dir_x->dcgrid == NULL) {
947  rate_x = g->ics_adi_dir_x->dc / (dx * dx);
948  rate_y = g->ics_adi_dir_y->dc / (dy * dy);
949  rate_z = g->ics_adi_dir_z->dc / (dz * dz);
950  // x contribution
951  variable_step_delta(x_line_start,
952  x_line_stop,
953  x_node_start,
954  ydot,
955  x_line_defs,
956  x_ordered_nodes,
957  states,
958  rate_x,
959  alphas);
960  // y contribution
961  variable_step_delta(y_line_start,
962  y_line_stop,
963  y_node_start,
964  ydot,
965  y_line_defs,
966  y_ordered_nodes,
967  states,
968  rate_y,
969  alphas);
970  // z contribution
971  variable_step_delta(z_line_start,
972  z_line_stop,
973  z_node_start,
974  ydot,
975  z_line_defs,
976  z_ordered_nodes,
977  states,
978  rate_z,
979  alphas);
980  } else {
981  rate_x = 1.0 / (dx * dx);
982  rate_y = 1.0 / (dy * dy);
983  rate_z = 1.0 / (dz * dz);
984  // x contribution
985  variable_step_delta(x_line_start,
986  x_line_stop,
987  x_node_start,
988  ydot,
989  x_line_defs,
990  x_ordered_nodes,
991  states,
992  rate_x,
993  g->ics_adi_dir_x->dcgrid,
994  alphas);
995  // y contribution
996  variable_step_delta(y_line_start,
997  y_line_stop,
998  y_node_start,
999  ydot,
1000  y_line_defs,
1001  y_ordered_nodes,
1002  states,
1003  rate_y,
1004  g->ics_adi_dir_y->dcgrid,
1005  alphas);
1006  // z contribution
1007  variable_step_delta(z_line_start,
1008  z_line_stop,
1009  z_node_start,
1010  ydot,
1011  z_line_defs,
1012  z_ordered_nodes,
1013  states,
1014  rate_z,
1015  g->ics_adi_dir_z->dcgrid,
1016  alphas);
1017  }
1018 }
1019 
1020 // Inhomogeneous diffusion coefficient
1021 static void variable_step_x(ICS_Grid_node* g,
1022  int line_start,
1023  int line_stop,
1024  int node_start,
1025  double* states,
1026  double* CVodeRHS,
1027  double* RHS,
1028  double* scratchpad,
1029  double* alphas,
1030  double* dcs,
1031  double dt) {
1032  long* x_lines = g->ics_adi_dir_x->ordered_line_defs;
1033  long* ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
1034  double* l_diag = g->ics_tasks->l_diag;
1035  double* diag = g->ics_tasks->diag;
1036  double* u_diag = g->ics_tasks->u_diag;
1037  double* delta_x = g->ics_adi_dir_x->deltas;
1038  double* delta_y = g->ics_adi_dir_y->deltas;
1039  double* delta_z = g->ics_adi_dir_z->deltas;
1040 
1041  long next_index = -1;
1042  long prev_index = -1;
1043  double next;
1044  double prev;
1045  double dx = g->ics_adi_dir_x->d;
1046  double dy = g->ics_adi_dir_y->d;
1047  double dz = g->ics_adi_dir_z->d;
1048 
1049  long current_index;
1050  long ordered_index = node_start;
1051  for (int i = line_start; i < line_stop - 1; i += 2) {
1052  long N = x_lines[i + 1];
1053  long ordered_index_start = ordered_index;
1054  for (int j = 0; j < N; j++) {
1055  current_index = ordered_nodes[ordered_index];
1056  RHS[j] = CVodeRHS[current_index] -
1057  dt *
1058  (delta_x[current_index] / SQ(dx) + delta_y[current_index] / SQ(dy) +
1059  delta_z[current_index] / SQ(dz)) /
1060  alphas[current_index];
1061  ordered_index++;
1062  }
1063 
1064  ordered_index = ordered_index_start;
1065  current_index = ordered_nodes[ordered_index];
1066  ordered_index++;
1067  next_index = ordered_nodes[ordered_index];
1068  next = alphas[next_index] * dcs[next_index] / (alphas[next_index] + alphas[current_index]);
1069  diag[0] = 1.0 + dt * next / SQ(dx);
1070  u_diag[0] = -dt * next / SQ(dx);
1071  ordered_index++;
1072  for (int c = 1; c < N - 1; c++) {
1073  prev_index = current_index;
1074  current_index = next_index;
1075  next_index = ordered_nodes[ordered_index];
1076  prev = alphas[prev_index] * dcs[current_index] /
1077  (alphas[prev_index] + alphas[current_index]);
1078  next = alphas[next_index] * dcs[next_index] /
1079  (alphas[next_index] + alphas[current_index]);
1080  l_diag[c - 1] = -dt * prev / SQ(dx);
1081  diag[c] = 1. + dt * (prev + next) / SQ(dx);
1082  u_diag[c] = -dt * next / SQ(dx);
1083  ordered_index++;
1084  }
1085  prev = alphas[current_index] * dcs[current_index] /
1086  (alphas[current_index] + alphas[next_index]);
1087  diag[N - 1] = 1.0 + dt * prev / SQ(dx);
1088  l_diag[N - 2] = -dt * prev / SQ(dx);
1089 
1090  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1091 
1092  ordered_index = ordered_index_start;
1093  for (int k = 0; k < N; k++) {
1094  current_index = ordered_nodes[ordered_index];
1095  states[current_index] = RHS[k];
1096  ordered_index++;
1097  }
1098  }
1099 }
1100 static void variable_step_y(ICS_Grid_node* g,
1101  int line_start,
1102  int line_stop,
1103  int node_start,
1104  double* states,
1105  double* RHS,
1106  double* scratchpad,
1107  double* alphas,
1108  double* dcs,
1109  double dt) {
1110  double* delta = g->ics_adi_dir_y->deltas;
1111  long* lines = g->ics_adi_dir_y->ordered_line_defs;
1112 
1113  double* l_diag = g->ics_tasks->l_diag;
1114  double* diag = g->ics_tasks->diag;
1115  double* u_diag = g->ics_tasks->u_diag;
1116 
1117  long next_index = -1;
1118  long prev_index = -1;
1119  double next;
1120  double prev;
1121  double dy = g->ics_adi_dir_y->d;
1122 
1123  long current_index;
1124  long* ordered_y_nodes = g->ics_adi_dir_y->ordered_nodes;
1125  long ordered_index = node_start;
1126 
1127  for (int i = line_start; i < line_stop - 1; i += 2) {
1128  long N = lines[i + 1];
1129  long ordered_index_start = ordered_index;
1130 
1131  for (int j = 0; j < N; j++) {
1132  current_index = ordered_y_nodes[ordered_index];
1133  RHS[j] = states[current_index] -
1134  dt * delta[current_index] / (alphas[current_index] * SQ(dy));
1135  ordered_index++;
1136  }
1137 
1138  ordered_index = ordered_index_start;
1139  current_index = ordered_y_nodes[ordered_index];
1140  ordered_index++;
1141  next_index = ordered_y_nodes[ordered_index];
1142  next = alphas[next_index] * dcs[next_index] / (alphas[next_index] + alphas[current_index]);
1143  diag[0] = 1.0 + dt * next / SQ(dy);
1144  u_diag[0] = -dt * next / SQ(dy);
1145  ordered_index++;
1146  for (int c = 1; c < N - 1; c++) {
1147  prev_index = current_index;
1148  current_index = next_index;
1149  next_index = ordered_y_nodes[ordered_index];
1150  prev = alphas[prev_index] * dcs[current_index] /
1151  (alphas[prev_index] + alphas[current_index]);
1152  next = alphas[next_index] * dcs[next_index] /
1153  (alphas[next_index] + alphas[current_index]);
1154  l_diag[c - 1] = -dt * prev / SQ(dy);
1155  diag[c] = 1. + dt * (prev + next) / SQ(dy);
1156  u_diag[c] = -dt * next / SQ(dy);
1157  ordered_index++;
1158  }
1159  prev = alphas[current_index] * dcs[current_index] /
1160  (alphas[current_index] + alphas[next_index]);
1161  diag[N - 1] = 1.0 + dt * prev / SQ(dy);
1162  l_diag[N - 2] = -dt * prev / SQ(dy);
1163 
1164  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1165 
1166  ordered_index = ordered_index_start;
1167  for (int k = 0; k < N; k++) {
1168  current_index = ordered_y_nodes[ordered_index];
1169  states[current_index] = RHS[k];
1170  ordered_index++;
1171  }
1172  }
1173 }
1174 static void variable_step_z(ICS_Grid_node* g,
1175  int line_start,
1176  int line_stop,
1177  int node_start,
1178  double* states,
1179  double* RHS,
1180  double* scratchpad,
1181  double* alphas,
1182  double* dcs,
1183  double dt) {
1184  double* delta = g->ics_adi_dir_z->deltas;
1185  long* lines = g->ics_adi_dir_z->ordered_line_defs;
1186  long* ordered_z_nodes = g->ics_adi_dir_z->ordered_nodes;
1187 
1188  double* l_diag = g->ics_tasks->l_diag;
1189  double* diag = g->ics_tasks->diag;
1190  double* u_diag = g->ics_tasks->u_diag;
1191 
1192  long next_index = -1;
1193  long prev_index = -1;
1194  double next;
1195  double prev;
1196  double dz = g->ics_adi_dir_z->d;
1197 
1198  long current_index;
1199  long ordered_index = node_start;
1200  for (int i = line_start; i < line_stop - 1; i += 2) {
1201  long N = lines[i + 1];
1202  long ordered_index_start = ordered_index;
1203 
1204  for (int j = 0; j < N; j++) {
1205  current_index = ordered_z_nodes[ordered_index];
1206  RHS[j] = states[current_index] -
1207  dt * delta[current_index] / (alphas[current_index] * SQ(dz));
1208  ordered_index++;
1209  }
1210 
1211  ordered_index = ordered_index_start;
1212  current_index = ordered_z_nodes[ordered_index];
1213  ordered_index++;
1214  next_index = ordered_z_nodes[ordered_index];
1215  next = alphas[next_index] * dcs[next_index] / (alphas[next_index] + alphas[current_index]);
1216  diag[0] = 1.0 + dt * next / SQ(dz);
1217  u_diag[0] = -dt * next / SQ(dz);
1218  ordered_index++;
1219  for (int c = 1; c < N - 1; c++) {
1220  prev_index = current_index;
1221  current_index = next_index;
1222  next_index = ordered_z_nodes[ordered_index];
1223  prev = alphas[prev_index] * dcs[current_index] /
1224  (alphas[prev_index] + alphas[current_index]);
1225  next = alphas[next_index] * dcs[next_index] /
1226  (alphas[next_index] + alphas[current_index]);
1227  l_diag[c - 1] = -dt * prev / SQ(dz);
1228  diag[c] = 1. + dt * (prev + next) / SQ(dz);
1229  u_diag[c] = -dt * next / SQ(dz);
1230  ordered_index++;
1231  }
1232  prev = alphas[current_index] * dcs[current_index] /
1233  (alphas[current_index] + alphas[next_index]);
1234  diag[N - 1] = 1.0 + dt * prev / SQ(dz);
1235  l_diag[N - 2] = -dt * prev / SQ(dz);
1236  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1237 
1238  ordered_index = ordered_index_start;
1239  for (int k = 0; k < N; k++) {
1240  current_index = ordered_z_nodes[ordered_index];
1241  states[current_index] = RHS[k];
1242  ordered_index++;
1243  }
1244  }
1245 }
1246 
1247 // Homogeneous diffusion coefficient
1248 static void variable_step_x(ICS_Grid_node* g,
1249  int line_start,
1250  int line_stop,
1251  int node_start,
1252  double* states,
1253  double* CVodeRHS,
1254  double* RHS,
1255  double* scratchpad,
1256  double* alphas,
1257  double dt) {
1258  long* x_lines = g->ics_adi_dir_x->ordered_line_defs;
1259  long* ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
1260  double* l_diag = g->ics_tasks->l_diag;
1261  double* diag = g->ics_tasks->diag;
1262  double* u_diag = g->ics_tasks->u_diag;
1263  double* delta_x = g->ics_adi_dir_x->deltas;
1264  double* delta_y = g->ics_adi_dir_y->deltas;
1265  double* delta_z = g->ics_adi_dir_z->deltas;
1266  long next_index = -1;
1267  long prev_index = -1;
1268  double next;
1269  double prev;
1270  double dx = g->ics_adi_dir_x->d;
1271  double dy = g->ics_adi_dir_y->d;
1272  double dz = g->ics_adi_dir_z->d;
1273  double dc = g->ics_adi_dir_x->dc;
1274 
1275  double r = dt * dc / SQ(dx);
1276 
1277  long current_index;
1278  long ordered_index = node_start;
1279  for (int i = line_start; i < line_stop - 1; i += 2) {
1280  long N = x_lines[i + 1];
1281  long ordered_index_start = ordered_index;
1282  for (int j = 0; j < N; j++) {
1283  current_index = ordered_nodes[ordered_index];
1284  RHS[j] = CVodeRHS[current_index] -
1285  dt *
1286  (delta_x[current_index] / SQ(dx) + delta_y[current_index] / SQ(dy) +
1287  delta_z[current_index] / SQ(dz)) /
1288  alphas[current_index];
1289  ordered_index++;
1290  }
1291 
1292  ordered_index = ordered_index_start;
1293  current_index = ordered_nodes[ordered_index];
1294  ordered_index++;
1295  next_index = ordered_nodes[ordered_index];
1296  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1297  diag[0] = 1.0 + next;
1298  u_diag[0] = -next;
1299  ordered_index++;
1300  for (int c = 1; c < N - 1; c++) {
1301  prev_index = current_index;
1302  current_index = next_index;
1303  next_index = ordered_nodes[ordered_index];
1304  prev = alphas[prev_index] * r / (alphas[prev_index] + alphas[current_index]);
1305  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1306  l_diag[c - 1] = -prev;
1307  diag[c] = 1. + prev + next;
1308  u_diag[c] = -next;
1309  ordered_index++;
1310  }
1311  prev = alphas[current_index] * r / (alphas[current_index] + alphas[next_index]);
1312  diag[N - 1] = 1.0 + prev;
1313  l_diag[N - 2] = -prev;
1314 
1315  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1316 
1317  ordered_index = ordered_index_start;
1318  for (int k = 0; k < N; k++) {
1319  current_index = ordered_nodes[ordered_index];
1320  states[current_index] = RHS[k];
1321  ordered_index++;
1322  }
1323  }
1324 }
1325 static void variable_step_y(ICS_Grid_node* g,
1326  int line_start,
1327  int line_stop,
1328  int node_start,
1329  double* states,
1330  double* RHS,
1331  double* scratchpad,
1332  double* alphas,
1333  double dt) {
1334  double* delta = g->ics_adi_dir_y->deltas;
1335  long* lines = g->ics_adi_dir_y->ordered_line_defs;
1336 
1337  double* l_diag = g->ics_tasks->l_diag;
1338  double* diag = g->ics_tasks->diag;
1339  double* u_diag = g->ics_tasks->u_diag;
1340 
1341  long next_index = -1;
1342  long prev_index = -1;
1343  double next;
1344  double prev;
1345  double dy = g->ics_adi_dir_y->d;
1346  double dc = g->ics_adi_dir_y->dc;
1347 
1348  double r = dt * dc / SQ(dy);
1349 
1350  long current_index;
1351  long* ordered_y_nodes = g->ics_adi_dir_y->ordered_nodes;
1352  long ordered_index = node_start;
1353 
1354  for (int i = line_start; i < line_stop - 1; i += 2) {
1355  long N = lines[i + 1];
1356  long ordered_index_start = ordered_index;
1357 
1358  for (int j = 0; j < N; j++) {
1359  current_index = ordered_y_nodes[ordered_index];
1360  RHS[j] = states[current_index] -
1361  dt * delta[current_index] / (alphas[current_index] * SQ(dy));
1362  ordered_index++;
1363  }
1364 
1365  ordered_index = ordered_index_start;
1366  current_index = ordered_y_nodes[ordered_index];
1367  ordered_index++;
1368  next_index = ordered_y_nodes[ordered_index];
1369  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1370  diag[0] = 1.0 + next;
1371  u_diag[0] = -next;
1372  ordered_index++;
1373  for (int c = 1; c < N - 1; c++) {
1374  prev_index = current_index;
1375  current_index = next_index;
1376  next_index = ordered_y_nodes[ordered_index];
1377  prev = alphas[prev_index] * r / (alphas[prev_index] + alphas[current_index]);
1378  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1379  l_diag[c - 1] = -prev;
1380  diag[c] = 1. + prev + next;
1381  u_diag[c] = -next;
1382  ordered_index++;
1383  }
1384  prev = alphas[current_index] * r / (alphas[current_index] + alphas[next_index]);
1385  diag[N - 1] = 1.0 + prev;
1386  l_diag[N - 2] = -prev;
1387 
1388  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1389 
1390  ordered_index = ordered_index_start;
1391  for (int k = 0; k < N; k++) {
1392  current_index = ordered_y_nodes[ordered_index];
1393  states[current_index] = RHS[k];
1394  ordered_index++;
1395  }
1396  }
1397 }
1398 static void variable_step_z(ICS_Grid_node* g,
1399  int line_start,
1400  int line_stop,
1401  int node_start,
1402  double* states,
1403  double* RHS,
1404  double* scratchpad,
1405  double* alphas,
1406  double dt) {
1407  double* delta = g->ics_adi_dir_z->deltas;
1408  long* lines = g->ics_adi_dir_z->ordered_line_defs;
1409  long* ordered_z_nodes = g->ics_adi_dir_z->ordered_nodes;
1410 
1411  double* l_diag = g->ics_tasks->l_diag;
1412  double* diag = g->ics_tasks->diag;
1413  double* u_diag = g->ics_tasks->u_diag;
1414 
1415  long next_index = -1;
1416  long prev_index = -1;
1417  double next;
1418  double prev;
1419  double dz = g->ics_adi_dir_z->d;
1420  double dc = g->ics_adi_dir_z->dc;
1421  double r = dt * dc / SQ(dz);
1422 
1423 
1424  long current_index;
1425  long ordered_index = node_start;
1426  for (int i = line_start; i < line_stop - 1; i += 2) {
1427  long N = lines[i + 1];
1428  long ordered_index_start = ordered_index;
1429 
1430  for (int j = 0; j < N; j++) {
1431  current_index = ordered_z_nodes[ordered_index];
1432  RHS[j] = states[current_index] -
1433  dt * delta[current_index] / (alphas[current_index] * SQ(dz));
1434  ordered_index++;
1435  }
1436 
1437  ordered_index = ordered_index_start;
1438  current_index = ordered_z_nodes[ordered_index];
1439  ordered_index++;
1440  next_index = ordered_z_nodes[ordered_index];
1441  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1442  diag[0] = 1.0 + next;
1443  u_diag[0] = -next;
1444  ordered_index++;
1445  for (int c = 1; c < N - 1; c++) {
1446  prev_index = current_index;
1447  current_index = next_index;
1448  next_index = ordered_z_nodes[ordered_index];
1449  prev = alphas[prev_index] * r / (alphas[prev_index] + alphas[current_index]);
1450  next = alphas[next_index] * r / (alphas[next_index] + alphas[current_index]);
1451  l_diag[c - 1] = -prev;
1452  diag[c] = 1. + prev + next;
1453  u_diag[c] = -next;
1454  ordered_index++;
1455  }
1456  prev = alphas[current_index] * r / (alphas[current_index] + alphas[next_index]);
1457  diag[N - 1] = 1.0 + prev;
1458  l_diag[N - 2] = -prev;
1459  solve_dd_tridiag(N, l_diag, diag, u_diag, RHS, scratchpad);
1460 
1461  ordered_index = ordered_index_start;
1462  for (int k = 0; k < N; k++) {
1463  current_index = ordered_z_nodes[ordered_index];
1464  states[current_index] = RHS[k];
1465  ordered_index++;
1466  }
1467  }
1468 }
1469 
1470 void ics_ode_solve_helper(ICS_Grid_node* g, double dt, double* CVodeRHS) {
1471  int num_states = g->_num_nodes;
1472 
1473  long x_line_start = g->ics_adi_dir_x->line_start_stop_indices[0];
1474  long x_line_stop = g->ics_adi_dir_x->line_start_stop_indices[NUM_THREADS * 2 - 1];
1475  long x_node_start = g->ics_adi_dir_x->ordered_start_stop_indices[0];
1476  long* x_line_defs = g->ics_adi_dir_x->ordered_line_defs;
1477  long* x_ordered_nodes = g->ics_adi_dir_x->ordered_nodes;
1478  double* delta_x = g->ics_adi_dir_x->deltas;
1479 
1480  long y_line_start = g->ics_adi_dir_y->line_start_stop_indices[0];
1481  long y_line_stop = g->ics_adi_dir_y->line_start_stop_indices[NUM_THREADS * 2 - 1];
1482  long y_node_start = g->ics_adi_dir_y->ordered_start_stop_indices[0];
1483  long* y_line_defs = g->ics_adi_dir_y->ordered_line_defs;
1484  long* y_ordered_nodes = g->ics_adi_dir_y->ordered_nodes;
1485  double* delta_y = g->ics_adi_dir_y->deltas;
1486 
1487  long z_line_start = g->ics_adi_dir_z->line_start_stop_indices[0];
1488  long z_line_stop = g->ics_adi_dir_z->line_start_stop_indices[NUM_THREADS * 2 - 1];
1489  long z_node_start = g->ics_adi_dir_z->ordered_start_stop_indices[0];
1490  long* z_line_defs = g->ics_adi_dir_z->ordered_line_defs;
1491  long* z_ordered_nodes = g->ics_adi_dir_z->ordered_nodes;
1492  double* delta_z = g->ics_adi_dir_z->deltas;
1493 
1494  double* Grid_RHS = g->ics_tasks->RHS;
1495  double* Grid_ScratchPad = g->ics_tasks->scratchpad;
1496 
1497  double* CVode_states_copy = (double*) calloc(num_states, sizeof(double));
1498  memcpy(CVode_states_copy, CVodeRHS, sizeof(double) * num_states);
1499 
1500  double* alphas = g->_ics_alphas;
1501 
1502  if (g->ics_adi_dir_x->dcgrid == NULL) {
1503  // find deltas for x, y and z directions
1504  ics_find_deltas(x_line_start,
1505  x_line_stop,
1506  x_node_start,
1507  delta_x,
1508  x_line_defs,
1509  x_ordered_nodes,
1510  CVode_states_copy,
1511  g->ics_adi_dir_x->dc,
1512  alphas);
1513  ics_find_deltas(y_line_start,
1514  y_line_stop,
1515  y_node_start,
1516  delta_y,
1517  y_line_defs,
1518  y_ordered_nodes,
1519  CVode_states_copy,
1520  g->ics_adi_dir_y->dc,
1521  alphas);
1522  ics_find_deltas(z_line_start,
1523  z_line_stop,
1524  z_node_start,
1525  delta_z,
1526  z_line_defs,
1527  z_ordered_nodes,
1528  CVode_states_copy,
1529  g->ics_adi_dir_z->dc,
1530  alphas);
1531 
1532  variable_step_x(g,
1533  x_line_start,
1534  x_line_stop,
1535  x_node_start,
1536  CVode_states_copy,
1537  CVodeRHS,
1538  Grid_RHS,
1539  Grid_ScratchPad,
1540  alphas,
1541  dt);
1542  variable_step_y(g,
1543  y_line_start,
1544  y_line_stop,
1545  y_node_start,
1546  CVode_states_copy,
1547  Grid_RHS,
1548  Grid_ScratchPad,
1549  alphas,
1550  dt);
1551  variable_step_z(g,
1552  z_line_start,
1553  z_line_stop,
1554  z_node_start,
1555  CVode_states_copy,
1556  Grid_RHS,
1557  Grid_ScratchPad,
1558  alphas,
1559  dt);
1560  } else {
1561  // find deltas for x, y and z directions
1562  ics_find_deltas(x_line_start,
1563  x_line_stop,
1564  x_node_start,
1565  delta_x,
1566  x_line_defs,
1567  x_ordered_nodes,
1568  CVode_states_copy,
1569  g->ics_adi_dir_x->dcgrid,
1570  alphas);
1571  ics_find_deltas(y_line_start,
1572  y_line_stop,
1573  y_node_start,
1574  delta_y,
1575  y_line_defs,
1576  y_ordered_nodes,
1577  CVode_states_copy,
1578  g->ics_adi_dir_y->dcgrid,
1579  alphas);
1580  ics_find_deltas(z_line_start,
1581  z_line_stop,
1582  z_node_start,
1583  delta_z,
1584  z_line_defs,
1585  z_ordered_nodes,
1586  CVode_states_copy,
1587  g->ics_adi_dir_z->dcgrid,
1588  alphas);
1589 
1590  variable_step_x(g,
1591  x_line_start,
1592  x_line_stop,
1593  x_node_start,
1594  CVode_states_copy,
1595  CVodeRHS,
1596  Grid_RHS,
1597  Grid_ScratchPad,
1598  alphas,
1599  g->ics_adi_dir_x->dcgrid,
1600  dt);
1601  variable_step_y(g,
1602  y_line_start,
1603  y_line_stop,
1604  y_node_start,
1605  CVode_states_copy,
1606  Grid_RHS,
1607  Grid_ScratchPad,
1608  alphas,
1609  g->ics_adi_dir_y->dcgrid,
1610  dt);
1611  variable_step_z(g,
1612  z_line_start,
1613  z_line_stop,
1614  z_node_start,
1615  CVode_states_copy,
1616  Grid_RHS,
1617  Grid_ScratchPad,
1618  alphas,
1619  g->ics_adi_dir_z->dcgrid,
1620  dt);
1621  }
1622 
1623  memcpy(CVodeRHS, CVode_states_copy, sizeof(double) * num_states);
1624  free(CVode_states_copy);
1625 }
1626 
1627 void _ics_hybrid_helper(ICS_Grid_node* g) {
1628  double dt = *dt_ptr;
1629  long num_1d_indices = g->hybrid_data->num_1d_indices;
1630  long* indices1d = g->hybrid_data->indices1d;
1631  long* num_3d_indices_per_1d_seg = g->hybrid_data->num_3d_indices_per_1d_seg;
1632  long* indices3d = g->hybrid_data->indices3d;
1633  double* rates = g->hybrid_data->rates;
1634  double* volumes1d = g->hybrid_data->volumes1d;
1635  double* volumes3d = g->hybrid_data->volumes3d;
1636 
1637  double vol_1d, vol_3d, rate, conc_1d;
1638  int my_3d_index, my_1d_index;
1639  int vol_3d_index;
1640 
1641  int total_num_3d_indices_per_1d_seg = 0;
1642  for (int i = 0; i < num_1d_indices; i++) {
1643  total_num_3d_indices_per_1d_seg += num_3d_indices_per_1d_seg[i];
1644  }
1645 
1646  // store the state values in the order that we need them
1647  double* old_g_states = (double*) malloc(sizeof(double) * total_num_3d_indices_per_1d_seg);
1648 
1649  vol_3d_index = 0;
1650  for (int i = 0; i < num_1d_indices; i++) {
1651  for (int j = 0; j < num_3d_indices_per_1d_seg[i]; j++, vol_3d_index++) {
1652  old_g_states[vol_3d_index] = g->states[indices3d[vol_3d_index]];
1653  }
1654  }
1655 
1656  vol_3d_index = 0;
1657  for (int i = 0; i < num_1d_indices; i++) {
1658  vol_1d = volumes1d[i];
1659  my_1d_index = indices1d[i];
1660  conc_1d = states[my_1d_index];
1661 
1662 
1663  // printf("for 1d index %d: volume 1d = %g and conc1d = %g\n", my_1d_index, vol_1d,
1664  // conc_1d);
1665  for (int j = 0; j < num_3d_indices_per_1d_seg[i]; j++, vol_3d_index++) {
1666  vol_3d = volumes3d[vol_3d_index];
1667  // rate is rate of change of 3d concentration
1668  my_3d_index = indices3d[vol_3d_index];
1669  rate = (rates[vol_3d_index]) * (old_g_states[vol_3d_index] - conc_1d);
1670  // forward euler coupling
1671  g->states[my_3d_index] -= dt * rate;
1672  states[my_1d_index] += dt * rate * vol_3d / vol_1d;
1673  // printf("for 3d index %d: volume 3d = %g and states3d = %g and rate3d = %g\n",
1674  // my_3d_index, vol_3d, g->states[my_3d_index], rates[vol_3d_index]);
1675  }
1676  }
1677 
1678  free(old_g_states);
1679 }
1680 
1681 void _ics_variable_hybrid_helper(ICS_Grid_node* g,
1682  const double* cvode_states_3d,
1683  double* const ydot_3d,
1684  const double* cvode_states_1d,
1685  double* const ydot_1d) {
1686  long num_1d_indices = g->hybrid_data->num_1d_indices;
1687  long* indices1d = g->hybrid_data->indices1d;
1688  long* num_3d_indices_per_1d_seg = g->hybrid_data->num_3d_indices_per_1d_seg;
1689  long* indices3d = g->hybrid_data->indices3d;
1690  double* rates = g->hybrid_data->rates;
1691  double* volumes1d = g->hybrid_data->volumes1d;
1692  double* volumes3d = g->hybrid_data->volumes3d;
1693 
1694  double vol_1d, vol_3d, rate, conc_1d;
1695  int my_3d_index, my_1d_index;
1696  int vol_3d_index = 0;
1697 
1698  for (int i = 0; i < num_1d_indices; i++) {
1699  vol_1d = volumes1d[i];
1700  my_1d_index = indices1d[i];
1701  conc_1d = cvode_states_1d[my_1d_index];
1702  for (int j = 0; j < num_3d_indices_per_1d_seg[i]; j++, vol_3d_index++) {
1703  vol_3d = volumes3d[vol_3d_index];
1704  // rate is rate of change of 3d concentration
1705  my_3d_index = indices3d[vol_3d_index];
1706  rate = (rates[vol_3d_index]) * (cvode_states_3d[my_3d_index] - conc_1d);
1707  // forward euler coupling
1708  ydot_3d[my_3d_index] -= rate;
1709  ydot_1d[my_1d_index] += rate * vol_3d / vol_1d;
1710  }
1711  }
1712 }
Grid_node
Definition: grids.h:120
Grid_node::hybrid
bool hybrid
Definition: grids.h:139
Grid_node::next
Grid_node * next
Definition: grids.h:122
Grid_node::hybrid_data
Hybrid_data * hybrid_data
Definition: grids.h:141
ICS_Grid_node
Definition: grids.h:289
ICS_Grid_node::run_threaded_ics_dg_adi
void run_threaded_ics_dg_adi(struct ICSAdiDirection *)
Definition: rxd_intracellular.cpp:774
ICS_Grid_node::ics_tasks
struct ICSAdiGridData * ics_tasks
Definition: grids.h:320
dt
double dt
Definition: netcvode.cpp:76
dc
#define dc
c
#define c
diag
#define diag(s)
Definition: fmenu.cpp:192
dt_ptr
double * dt_ptr
Definition: grids.cpp:14
Parallel_grids
Grid_node * Parallel_grids[100]
Definition: grids.cpp:16
SQ
#define SQ(x)
Definition: grids.h:24
start
void start()
Definition: hel2mos.cpp:204
stop
void stop()
Definition: hel2mos.cpp:212
void
void
id
#define id
Definition: md1redef.h:33
i
#define i
Definition: md1redef.h:12
prev
Item * prev(Item *item)
Definition: list.cpp:93
next
Item * next(Item *item)
Definition: list.cpp:88
RHS
#define RHS(i)
Definition: multisplit.cpp:66
j
size_t j
Definition: nrngsl_real_radix2.cpp:50
k
static philox4x32_key_t k
Definition: nrnran123.cpp:11
nrnwrap_Python.h
g
#define g
Definition: passive0.cpp:21
data
#define data
Definition: rbtqueue.cpp:49
num_states
unsigned int num_states
Definition: rxd.cpp:63
TaskQueue_add_task
void TaskQueue_add_task(TaskQueue *q, void *(*task)(void *), void *args, void *result)
Definition: rxd.cpp:1087
TaskQueue_sync
void TaskQueue_sync(TaskQueue *q)
Definition: rxd.cpp:1187
variable_step_delta
static void variable_step_delta(long start, long stop, long node_start, double *ydot, long *line_defs, long *ordered_nodes, double const *const states, double r, double *alphas)
Definition: rxd_intracellular.cpp:804
do_ics_deltas
static void * do_ics_deltas(void *dataptr)
Definition: rxd_intracellular.cpp:236
ics_dg_adi_x_inhom
void ics_dg_adi_x_inhom(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:296
_ics_rhs_variable_step_helper
void _ics_rhs_variable_step_helper(ICS_Grid_node *g, double const *const states, double *ydot)
Definition: rxd_intracellular.cpp:921
run_threaded_deltas
void run_threaded_deltas(ICS_Grid_node *g, ICSAdiDirection *ics_adi_dir)
Definition: rxd_intracellular.cpp:274
ics_dg_adi_z
void ics_dg_adi_z(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:671
variable_step_z
static void variable_step_z(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double *states, double *RHS, double *scratchpad, double *alphas, double *dcs, double dt)
Definition: rxd_intracellular.cpp:1174
variable_step_y
static void variable_step_y(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double *states, double *RHS, double *scratchpad, double *alphas, double *dcs, double dt)
Definition: rxd_intracellular.cpp:1100
solve_dd_tridiag
static int solve_dd_tridiag(int N, const double *l_diag, const double *diag, const double *u_diag, double *b, double *c)
Definition: rxd_intracellular.cpp:93
AllTasks
TaskQueue * AllTasks
Definition: rxd.cpp:28
flux
static double flux(const double alphai, const double alphaj, const double statei, const double statej)
Definition: rxd_intracellular.cpp:797
_ics_variable_hybrid_helper
void _ics_variable_hybrid_helper(ICS_Grid_node *g, const double *cvode_states_3d, double *const ydot_3d, const double *cvode_states_1d, double *const ydot_1d)
Definition: rxd_intracellular.cpp:1681
states
double * states
Definition: rxd.cpp:62
ICS_PREFETCH
const int ICS_PREFETCH
Definition: rxd_intracellular.cpp:15
ics_ode_solve_helper
void ics_ode_solve_helper(ICS_Grid_node *g, double dt, double *CVodeRHS)
Definition: rxd_intracellular.cpp:1470
ics_find_deltas
void ics_find_deltas(long start, long stop, long node_start, double *delta, long *line_defs, long *ordered_nodes, double *states, double dc, double *alphas)
Definition: rxd_intracellular.cpp:119
ics_dg_adi_y
void ics_dg_adi_y(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:601
variable_step_x
static void variable_step_x(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double *states, double *CVodeRHS, double *RHS, double *scratchpad, double *alphas, double *dcs, double dt)
Definition: rxd_intracellular.cpp:1021
_ics_hybrid_helper
void _ics_hybrid_helper(ICS_Grid_node *g)
Definition: rxd_intracellular.cpp:1627
ics_dg_adi_x
void ics_dg_adi_x(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:523
NUM_THREADS
int NUM_THREADS
Definition: rxd.cpp:26
set_hybrid_data
void set_hybrid_data(int64_t *num_1d_indices_per_grid, int64_t *num_3d_indices_per_grid, int64_t *hybrid_indices1d, int64_t *hybrid_indices3d, int64_t *num_3d_indices_per_1d_seg, int64_t *hybrid_grid_ids, double *rates, double *volumes1d, double *volumes3d, double *dxs)
Definition: rxd_intracellular.cpp:20
do_ics_dg_adi
void * do_ics_dg_adi(void *dataptr)
Definition: rxd_intracellular.cpp:740
ics_dg_adi_z_inhom
void ics_dg_adi_z_inhom(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:450
ics_dg_adi_y_inhom
void ics_dg_adi_y_inhom(ICS_Grid_node *g, int line_start, int line_stop, int node_start, double, double *states, double *RHS, double *scratchpad, double *u_diag, double *diag, double *l_diag)
Definition: rxd_intracellular.cpp:377
current_state
static double current_state(void *v)
Definition: singlech.cpp:345
NULL
#define NULL
Definition: sptree.h:16
Hybrid_data::num_3d_indices_per_1d_seg
long * num_3d_indices_per_1d_seg
Definition: grids.h:64
Hybrid_data::indices3d
long * indices3d
Definition: grids.h:65
Hybrid_data::volumes1d
double * volumes1d
Definition: grids.h:67
Hybrid_data::volumes3d
double * volumes3d
Definition: grids.h:68
Hybrid_data::rates
double * rates
Definition: grids.h:66
Hybrid_data::num_1d_indices
long num_1d_indices
Definition: grids.h:62
Hybrid_data::indices1d
long * indices1d
Definition: grids.h:63
ICSAdiDirection::dcgrid
double * dcgrid
Definition: grids.h:385
ICSAdiDirection::ordered_start_stop_indices
long * ordered_start_stop_indices
Definition: grids.h:382
ICSAdiDirection::dc
double dc
Definition: grids.h:384
ICSAdiDirection::d
double d
Definition: grids.h:386
ICSAdiDirection::ordered_nodes
long * ordered_nodes
Definition: grids.h:381
ICSAdiDirection::ics_dg_adi_dir
void(* ics_dg_adi_dir)(ICS_Grid_node *g, int, int, int, double, double *, double *, double *, double *, double *, double *)
Definition: grids.h:366
ICSAdiDirection::line_start_stop_indices
long * line_start_stop_indices
Definition: grids.h:383
ICSAdiDirection::ordered_line_defs
long * ordered_line_defs
Definition: grids.h:380
ICSAdiDirection::deltas
double * deltas
Definition: grids.h:379
ICSAdiGridData::line_stop
int line_stop
Definition: grids.h:392
ICSAdiGridData::ics_adi_dir
ICSAdiDirection * ics_adi_dir
Definition: grids.h:397
ICSAdiGridData::ordered_start
int ordered_start
Definition: grids.h:394
ICSAdiGridData::line_start
int line_start
Definition: grids.h:392
TaskQueue
Definition: rxd.h:109