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