NEURON
rxd_vol.cpp
Go to the documentation of this file.
1 #include <../../nrnconf.h>
2 #include <stdio.h>
3 #include <assert.h>
4 #include "grids.h"
5 #include "rxd.h"
6 #include "nrnwrap_Python.h"
7 
8 /*Tortuous diffusion coefficients*/
9 #define DcX(x, y, z) (g->dc_x * PERM(x, y, z))
10 #define DcY(x, y, z) (g->dc_y * PERM(x, y, z))
11 #define DcZ(x, y, z) (g->dc_z * PERM(x, y, z))
12 
13 /*Flux in the x,y,z directions*/
14 // TODO: Refactor to avoid calculating indices
15 #define Fxx(x1, x2) \
16  (ALPHA(x1, y, z) * ALPHA(x2, y, z) * DcX(x1, y, z) * \
17  (g->states[IDX(x1, y, z)] - g->states[IDX(x2, y, z)]) / \
18  ((ALPHA(x1, y, z) + ALPHA(x2, y, z))))
19 #define Fxy(y1, y1d, y2) \
20  (ALPHA(x, y1, z) * ALPHA(x, y2, z) * DcY(x, y1d, z) * \
21  (g->states[IDX(x, y1, z)] - g->states[IDX(x, y2, z)]) / \
22  ((ALPHA(x, y1, z) + ALPHA(x, y2, z))))
23 #define Fxz(z1, z1d, z2) \
24  (ALPHA(x, y, z1) * ALPHA(x, y, z2) * DcZ(x, y, z1d) * \
25  (g->states[IDX(x, y, z1)] - g->states[IDX(x, y, z2)]) / \
26  ((ALPHA(x, y, z1) + ALPHA(x, y, z2))))
27 #define Fyy(y1, y2) \
28  (ALPHA(x, y1, z) * ALPHA(x, y2, z) * DcY(x, y1, z) * \
29  (g->states[IDX(x, y1, z)] - g->states[IDX(x, y2, z)]) / \
30  ((ALPHA(x, y1, z) + ALPHA(x, y2, z))))
31 #define Fzz(z1, z2) \
32  (ALPHA(x, y, z1) * ALPHA(x, y, z2) * DcZ(x, y, z1) * \
33  (g->states[IDX(x, y, z1)] - g->states[IDX(x, y, z2)]) / \
34  ((ALPHA(x, y, z1) + ALPHA(x, y, z2))))
35 
36 /*Flux for used by variable step inhomogeneous volume fraction*/
37 #define FLUX(pidx, idx) \
38  (VOLFRAC(pidx) * VOLFRAC(idx) * (states[pidx] - states[idx])) / \
39  (0.5 * (VOLFRAC(pidx) + VOLFRAC(idx)))
40 
41 extern int NUM_THREADS;
42 /* solve Ax=b where A is a diagonally dominant tridiagonal matrix
43  * N - length of the matrix
44  * l_diag - pointer to the lower diagonal (length N-1)
45  * diag - pointer to diagonal (length N)
46  * u_diag - pointer to the upper diagonal (length N-1)
47  * B - pointer to the RHS (length N)
48  * The solution (x) will be stored in B.
49  * l_diag, diag and u_diag are not changed.
50  * c - scratch pad, preallocated memory for (N-1) doubles
51  */
52 static int solve_dd_tridiag(int N,
53  const double* l_diag,
54  const double* diag,
55  const double* u_diag,
56  double* b,
57  double* c) {
58  int i;
59 
60  c[0] = u_diag[0] / diag[0];
61  b[0] = b[0] / diag[0];
62 
63  for (i = 1; i < N - 1; i++) {
64  c[i] = u_diag[i] / (diag[i] - l_diag[i - 1] * c[i - 1]);
65  b[i] = (b[i] - l_diag[i - 1] * b[i - 1]) / (diag[i] - l_diag[i - 1] * c[i - 1]);
66  }
67  b[N - 1] = (b[N - 1] - l_diag[N - 2] * b[N - 2]) / (diag[N - 1] - l_diag[N - 2] * c[N - 2]);
68 
69 
70  /*back substitution*/
71  for (i = N - 2; i >= 0; i--) {
72  b[i] = b[i] - c[i] * b[i + 1];
73  }
74  return 0;
75 }
76 
77 
78 /* dg_adi_vol_x performs the first of 3 steps in DG-ADI
79  * g - the parameters and state of the grid
80  * dt - the time step
81  * y - the index for the y plane
82  * z - the index for the z plane
83  * state - where the output of this step is stored
84  * scratch - scratchpad array of doubles, length g.size_x - 1
85  * like dg_adi_x except the grid has a variable volume fraction
86  * g.alpha and my have variable tortuosity g.lambda
87  */
88 static void ecs_dg_adi_vol_x(ECS_Grid_node* g,
89  const double dt,
90  const int y,
91  const int z,
92  double const* const state,
93  double* const RHS,
94  double* const scratch) {
95  int yp, ypd, ym, ymd, zp, zpd, zm, zmd;
96  int x;
97  double* diag;
98  double* l_diag;
99  double* u_diag;
100  double prev, next, div_y, div_z;
101 
102  /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
103  if (g->bc->type == DIRICHLET &&
104  (y == 0 || z == 0 || y == g->size_y - 1 || z == g->size_z - 1)) {
105  for (x = 0; x < g->size_x; x++)
106  RHS[x] = g->bc->value;
107  return;
108  }
109 
110  /*zero flux boundary condition*/
111  div_y = ((y == 0 || y == g->size_y - 1) ? 1.0 : 0.5) * SQ(g->dy);
112  div_z = ((z == 0 || z == g->size_z - 1) ? 1.0 : 0.5) * SQ(g->dz);
113  if (g->size_y == 1) {
114  yp = 0;
115  ypd = 0;
116  ym = 0;
117  ymd = 0;
118  } else {
119  yp = (y == g->size_y - 1) ? y - 1 : y + 1;
120  ypd = (y == g->size_y - 1) ? y : y + 1;
121  ym = (y == 0) ? y + 1 : y - 1;
122  ymd = (y == 0) ? 1 : y;
123  }
124  if (g->size_z == 1) {
125  zp = 0;
126  zpd = 0;
127  zm = 0;
128  zmd = 0;
129  } else {
130  zp = (z == g->size_z - 1) ? z - 1 : z + 1;
131  zpd = (z == g->size_z - 1) ? z : z + 1;
132  zm = (z == 0) ? z + 1 : z - 1;
133  zmd = (z == 0) ? 1 : z;
134  }
135 
136  if (g->size_x == 1) {
137  if (g->bc->type == DIRICHLET) {
138  RHS[0] = g->bc->value;
139  } else {
140  x = 0;
141  RHS[0] = 0;
142  if (g->size_y > 1)
143  RHS[0] += (Fxy(yp, ypd, y) - Fxy(y, ymd, ym)) / div_y;
144  if (g->size_z > 1)
145  RHS[0] += (Fxz(zp, zpd, z) - Fxz(z, zmd, zm)) / div_z;
146  RHS[0] *= (dt / ALPHA(0, y, z));
147  RHS[0] += state[IDX(0, y, z)] + g->states_cur[IDX(0, y, z)];
148  }
149  return;
150  }
151 
152  /*TODO: move allocation out of loop*/
153  diag = (double*) malloc(g->size_x * sizeof(double));
154  l_diag = (double*) malloc((g->size_x - 1) * sizeof(double));
155  u_diag = (double*) malloc((g->size_x - 1) * sizeof(double));
156 
157  for (x = 1; x < g->size_x - 1; x++) {
158  prev = ALPHA(x - 1, y, z) * DcX(x, y, z) / (ALPHA(x - 1, y, z) + ALPHA(x, y, z));
159  next = ALPHA(x + 1, y, z) * DcX(x + 1, y, z) / (ALPHA(x + 1, y, z) + ALPHA(x, y, z));
160 
161  l_diag[x - 1] = -dt * prev / SQ(g->dx);
162  diag[x] = 1. + dt * (prev + next) / SQ(g->dx);
163  u_diag[x] = -dt * next / SQ(g->dx);
164  }
165 
166 
167  if (g->bc->type == NEUMANN) {
168  /*zero flux boundary condition*/
169  next = ALPHA(1, y, z) * DcX(1, y, z) / (ALPHA(1, y, z) + ALPHA(0, y, z));
170  diag[0] = 1. + dt * next / SQ(g->dx);
171  u_diag[0] = -dt * next / SQ(g->dx);
172 
173  prev = ALPHA(g->size_x - 2, y, z) * DcX(g->size_x - 1, y, z) /
174  (ALPHA(g->size_x - 1, y, z) + ALPHA(g->size_x - 2, y, z));
175  diag[g->size_x - 1] = 1. + dt * prev / SQ(g->dx);
176  l_diag[g->size_x - 2] = -dt * prev / SQ(g->dx);
177 
178  x = 0;
179  RHS[x] = state[IDX(x, y, z)] +
180  (dt / ALPHA(x, y, z)) *
181  ((Fxx(x + 1, x) / SQ(g->dx)) + (Fxy(yp, ypd, y) - Fxy(y, ymd, ym)) / div_y +
182  (Fxz(zp, zpd, z) - Fxz(z, zmd, zm)) / div_z) +
183  g->states_cur[IDX(x, y, z)];
184 
185  x = g->size_x - 1;
186  RHS[x] = state[IDX(x, y, z)] +
187  (dt / ALPHA(x, y, z)) *
188  ((Fxx(x - 1, x)) / SQ(g->dx) + (Fxy(yp, ypd, y) - Fxy(y, ymd, ym)) / div_y +
189  (Fxz(zp, zpd, z) - Fxz(z, zmd, zm)) / div_z) +
190  g->states_cur[IDX(x, y, z)];
191  } else {
192  diag[0] = 1.0;
193  u_diag[0] = 0;
194  diag[g->size_x - 1] = 1.0;
195  l_diag[g->size_x - 2] = 0;
196  RHS[0] = g->bc->value;
197  RHS[g->size_x - 1] = g->bc->value;
198  }
199 
200  for (x = 1; x < g->size_x - 1; x++) {
201 #ifndef __PGI
202  __builtin_prefetch(&(state[IDX(x + PREFETCH, y, z)]), 0, 1);
203  __builtin_prefetch(&(state[IDX(x + PREFETCH, yp, z)]), 0, 0);
204  __builtin_prefetch(&(state[IDX(x + PREFETCH, ym, z)]), 0, 0);
205  __builtin_prefetch(&(state[IDX(x + PREFETCH, ypd, z)]), 0, 0);
206  __builtin_prefetch(&(state[IDX(x + PREFETCH, ymd, z)]), 0, 0);
207 #endif
208  RHS[x] = state[IDX(x, y, z)] +
209  (dt / ALPHA(x, y, z)) * ((Fxx(x + 1, x) - Fxx(x, x - 1)) / SQ(g->dx) +
210  (Fxy(yp, ypd, y) - Fxy(y, ymd, ym)) / div_y +
211  (Fxz(zp, zpd, z) - Fxz(z, zmd, zm)) / div_z) +
212  g->states_cur[IDX(x, y, z)];
213  }
214 
215  solve_dd_tridiag(g->size_x, l_diag, diag, u_diag, RHS, scratch);
216 
217  free(diag);
218  free(l_diag);
219  free(u_diag);
220 }
221 
222 /* dg_adi_vol_y performs the second of 3 steps in DG-ADI
223  * g - the parameters and state of the grid
224  * dt - the time step
225  * y - the index for the y plane
226  * z - the index for the z plane
227  * state - where the output of this step is stored
228  * scratch - scratchpad array of doubles, length g->size_y
229  * like dg_adi_y except the grid has a variable volume fraction
230  * g->alpha and my have variable tortuosity g->lambda
231  */
232 static void ecs_dg_adi_vol_y(ECS_Grid_node* g,
233  double const dt,
234  int const x,
235  int const z,
236  double const* const state,
237  double* const RHS,
238  double* const scratch) {
239  int y;
240  double* diag;
241  double* l_diag;
242  double* u_diag;
243  double prev, next;
244 
245  /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
246  if (g->bc->type == DIRICHLET &&
247  (x == 0 || z == 0 || x == g->size_x - 1 || z == g->size_z - 1)) {
248  for (y = 0; y < g->size_y; y++)
249  RHS[y] = g->bc->value;
250  return;
251  }
252 
253  if (g->size_y == 1) {
254  if (g->bc->type == DIRICHLET) {
255  RHS[0] = g->bc->value;
256  } else {
257  RHS[0] = state[x + z * g->size_x];
258  }
259  return;
260  }
261 
262  diag = (double*) malloc(g->size_y * sizeof(double));
263  l_diag = (double*) malloc((g->size_y - 1) * sizeof(double));
264  u_diag = (double*) malloc((g->size_y - 1) * sizeof(double));
265 
266  for (y = 1; y < g->size_y - 1; y++) {
267  prev = ALPHA(x, y - 1, z) * DcY(x, y, z) / (ALPHA(x, y - 1, z) + ALPHA(x, y, z));
268  next = ALPHA(x, y + 1, z) * DcY(x, y + 1, z) / (ALPHA(x, y + 1, z) + ALPHA(x, y, z));
269 
270  l_diag[y - 1] = -dt * prev / SQ(g->dy);
271  diag[y] = 1. + dt * (prev + next) / SQ(g->dy);
272  u_diag[y] = -dt * next / SQ(g->dy);
273  }
274 
275 
276  if (g->bc->type == NEUMANN) {
277  /*zero flux boundary condition*/
278  next = ALPHA(x, 1, z) * DcY(x, 1, z) / (ALPHA(x, 1, z) + ALPHA(x, 0, z));
279  diag[0] = 1. + dt * next / SQ(g->dy);
280  u_diag[0] = -dt * next / SQ(g->dy);
281 
282  prev = ALPHA(x, g->size_y - 2, z) * DcY(x, g->size_y - 1, z) /
283  (ALPHA(x, g->size_y - 1, z) + ALPHA(x, g->size_y - 2, z));
284 
285  diag[g->size_y - 1] = 1. + dt * prev / SQ(g->dy);
286  l_diag[g->size_y - 2] = -dt * prev / SQ(g->dy);
287 
288  RHS[0] = state[x + z * g->size_x] - dt * Fyy(1, 0) / (SQ(g->dy) * ALPHA(x, 0, z));
289  y = g->size_y - 1;
290  RHS[y] = state[x + (z + y * g->size_z) * g->size_x] +
291  (dt / ALPHA(x, y, z)) * Fyy(y, y - 1) / SQ(g->dy);
292  } else {
293  diag[0] = 1.0;
294  u_diag[0] = 0;
295  diag[g->size_y - 1] = 1.0;
296  l_diag[g->size_y - 2] = 0;
297  RHS[0] = g->bc->value;
298  RHS[g->size_y - 1] = g->bc->value;
299  }
300  for (y = 1; y < g->size_y - 1; y++) {
301 #ifndef __PGI
302  __builtin_prefetch(&state[x + (z + (y + PREFETCH) * g->size_z) * g->size_x], 0, 0);
303  __builtin_prefetch(&(g->states[IDX(x, y + PREFETCH, z)]), 0, 1);
304 #endif
305  RHS[y] = state[x + (z + y * g->size_z) * g->size_x] -
306  (dt / ALPHA(x, y, z)) * (Fyy(y + 1, y) - Fyy(y, y - 1)) / SQ(g->dy);
307  }
308 
309  solve_dd_tridiag(g->size_y, l_diag, diag, u_diag, RHS, scratch);
310 
311  free(diag);
312  free(l_diag);
313  free(u_diag);
314 }
315 
316 /* dg_adi_vol_z performs the third of 3 steps in DG-ADI
317  * g - the parameters and state of the grid
318  * dt - the time step
319  * y - the index for the y plane
320  * z - the index for the z plane
321  * state - where the output of this step is stored
322  * scratch - scratchpad array of doubles, length g->size_y
323  * like dg_adi_z except the grid has a variable volume fraction
324  * g->alpha and my have variable tortuosity g->lambda
325  */
326 static void ecs_dg_adi_vol_z(ECS_Grid_node* g,
327  double const dt,
328  int const x,
329  int const y,
330  double const* const state,
331  double* const RHS,
332  double* const scratch) {
333  int z;
334  double* diag;
335  double* l_diag;
336  double* u_diag;
337  double prev, next;
338 
339  /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
340  if (g->bc->type == DIRICHLET &&
341  (x == 0 || y == 0 || x == g->size_x - 1 || y == g->size_y - 1)) {
342  for (z = 0; z < g->size_z; z++)
343  RHS[z] = g->bc->value;
344  return;
345  }
346  if (g->size_z == 1) {
347  if (g->bc->type == DIRICHLET) {
348  RHS[0] = g->bc->value;
349  } else {
350  RHS[0] = state[y + g->size_y * x];
351  }
352  return;
353  }
354 
355  diag = (double*) malloc(g->size_z * sizeof(double));
356  l_diag = (double*) malloc((g->size_z - 1) * sizeof(double));
357  u_diag = (double*) malloc((g->size_z - 1) * sizeof(double));
358 
359  for (z = 1; z < g->size_z - 1; z++) {
360  prev = ALPHA(x, y, z - 1) * DcZ(x, y, z) / (ALPHA(x, y, z - 1) + ALPHA(x, y, z));
361  next = ALPHA(x, y, z + 1) * DcZ(x, y, z + 1) / (ALPHA(x, y, z + 1) + ALPHA(x, y, z));
362 
363  l_diag[z - 1] = -dt * prev / SQ(g->dz);
364  diag[z] = 1. + dt * (prev + next) / SQ(g->dz);
365  u_diag[z] = -dt * next / SQ(g->dz);
366  }
367 
368  if (g->bc->type == NEUMANN) {
369  /*zero flux boundary condition*/
370  next = ALPHA(x, y, 1) * DcZ(x, y, 1) / (ALPHA(x, y, 1) + ALPHA(x, y, 0));
371  diag[0] = 1. + dt * next / SQ(g->dz);
372  u_diag[0] = -dt * next / SQ(g->dz);
373 
374  prev = ALPHA(x, y, g->size_z - 2) * DcZ(x, y, g->size_z - 1) /
375  (ALPHA(x, y, g->size_z - 1) + ALPHA(x, y, g->size_z - 2));
376 
377  diag[g->size_z - 1] = 1. + dt * prev / SQ(g->dz);
378  l_diag[g->size_z - 2] = -dt * prev / SQ(g->dz);
379 
380  RHS[0] = state[y + g->size_y * (x * g->size_z)] -
381  dt * Fzz(1, 0) / (SQ(g->dz) * ALPHA(x, y, 0));
382  z = g->size_z - 1;
383  RHS[z] = state[y + g->size_y * (x * g->size_z + z)] +
384  (dt / ALPHA(x, y, z)) * Fzz(z, z - 1) / SQ(g->dz);
385  } else {
386  diag[0] = 1.0;
387  u_diag[0] = 0;
388  diag[g->size_z - 1] = 1.0;
389  l_diag[g->size_z - 2] = 0;
390  RHS[0] = g->bc->value;
391  RHS[g->size_z - 1] = g->bc->value;
392  }
393  for (z = 1; z < g->size_z - 1; z++) {
394  RHS[z] = state[y + g->size_y * (x * g->size_z + z)] -
395  (dt / ALPHA(x, y, z)) * (Fzz(z + 1, z) - Fzz(z, z - 1)) / SQ(g->dz);
396  }
397 
398  solve_dd_tridiag(g->size_z, l_diag, diag, u_diag, RHS, scratch);
399 
400  free(diag);
401  free(l_diag);
402  free(u_diag);
403 }
404 
405 /*DG-ADI implementation the 3 step process to diffusion species in grid g by time step *dt_ptr
406  * g - the state and parameters
407  * like dg_adi execpt the Grid_node g has variable volume fraction g.alpha and may have
408  * variable tortuosity g.lambda
409  */
410 void ecs_set_adi_vol(ECS_Grid_node* g) {
411  g->ecs_adi_dir_x->ecs_dg_adi_dir = ecs_dg_adi_vol_x;
412  g->ecs_adi_dir_y->ecs_dg_adi_dir = ecs_dg_adi_vol_y;
413  g->ecs_adi_dir_z->ecs_dg_adi_dir = ecs_dg_adi_vol_z;
414 }
415 
416 
417 /* dg_adi_tort_x performs the first of 3 steps in DG-ADI
418  * g - the parameters and state of the grid
419  * dt - the time step
420  * y - the index for the y plane
421  * z - the index for the z plane
422  * state - where the output of this step is stored
423  * scratch - scratchpad array of doubles, length g.size_x - 1
424  * like dg_adi_x except the grid has a variable tortuosity
425  * g.lambda (but it still has fixed volume fraction)
426  */
427 static void ecs_dg_adi_tort_x(ECS_Grid_node* g,
428  const double dt,
429  const int y,
430  const int z,
431  double const* const state,
432  double* const RHS,
433  double* const scratch) {
434  int yp, ypd, ym, ymd, zp, zpd, zm, zmd, div_y, div_z;
435  int x;
436  double* diag;
437  double* l_diag;
438  double* u_diag;
439 
440  if (g->bc->type == DIRICHLET &&
441  (y == 0 || z == 0 || y == g->size_y - 1 || z == g->size_z - 1)) {
442  for (x = 0; x < g->size_x; x++)
443  RHS[x] = g->bc->value;
444  return;
445  }
446 
447  /*zero flux boundary condition*/
448  div_y = (y == 0 || y == g->size_y - 1) ? 2 : 1;
449  div_z = (z == 0 || z == g->size_z - 1) ? 2 : 1;
450  if (g->size_y == 1) {
451  yp = 0;
452  ypd = 0;
453  ym = 0;
454  ymd = 0;
455  } else {
456  yp = (y == g->size_y - 1) ? y - 1 : y + 1;
457  ypd = (y == g->size_y - 1) ? y : y + 1;
458  ym = (y == 0) ? y + 1 : y - 1;
459  ymd = (y == 0) ? 1 : y;
460  }
461  if (g->size_z == 1) {
462  zp = 0;
463  zpd = 0;
464  zm = 0;
465  zmd = 0;
466  } else {
467  zp = (z == g->size_z - 1) ? z - 1 : z + 1;
468  zpd = (z == g->size_z - 1) ? z : z + 1;
469  zm = (z == 0) ? z + 1 : z - 1;
470  zmd = (z == 0) ? 1 : z;
471  }
472 
473  if (g->size_x == 1) {
474  if (g->bc->type == DIRICHLET) {
475  RHS[0] = g->bc->value;
476  } else {
477  RHS[0] = 0;
478  if (g->size_y > 1)
479  RHS[0] += (DcY(0, ypd, z) * state[IDX(0, yp, z)] -
480  (DcY(0, ypd, z) + DcY(0, ymd, z)) * state[IDX(0, y, z)] +
481  DcY(0, ymd, z) * state[IDX(0, ym, z)]) /
482  (div_y * SQ(g->dy));
483  if (g->size_z > 1)
484  RHS[0] += (DcZ(0, y, zpd) * state[IDX(0, y, zp)] -
485  (DcZ(0, y, zpd) + DcZ(0, y, zmd)) * state[IDX(0, y, z)] +
486  DcZ(0, y, zmd) * state[IDX(0, y, zm)]) /
487  (div_z * SQ(g->dz));
488  RHS[0] *= dt;
489  RHS[0] += state[IDX(0, y, z)] + g->states_cur[IDX(0, y, z)];
490  }
491  return;
492  }
493 
494  diag = (double*) malloc(g->size_x * sizeof(double));
495  l_diag = (double*) malloc((g->size_x - 1) * sizeof(double));
496  u_diag = (double*) malloc((g->size_x - 1) * sizeof(double));
497 
498  for (x = 1; x < g->size_x - 1; x++) {
499  l_diag[x - 1] = -dt * DcX(x, y, z) / (2. * SQ(g->dx));
500  diag[x] = 1. + dt * (DcX(x, y, z) + DcX(x + 1, y, z)) / (2. * SQ(g->dx));
501  u_diag[x] = -dt * DcX(x + 1, y, z) / (2. * SQ(g->dx));
502  }
503 
504 
505  if (g->bc->type == NEUMANN) {
506  diag[0] = 1. + 0.5 * dt * DcX(1, y, z) / SQ(g->dx);
507  u_diag[0] = -0.5 * dt * DcX(1, y, z) / SQ(g->dx);
508  diag[g->size_x - 1] = 1. + 0.5 * dt * DcX(g->size_x - 1, y, z) / SQ(g->dx);
509  l_diag[g->size_x - 2] = -0.5 * dt * DcX(g->size_x - 1, y, z) / SQ(g->dx);
510 
511 
512  RHS[0] = state[IDX(0, y, z)] +
513  dt * ((DcX(1, y, z) * state[IDX(1, y, z)] - DcX(1, y, z) * state[IDX(0, y, z)]) /
514  (2. * SQ(g->dx)) +
515  (DcY(0, ypd, z) * state[IDX(0, yp, z)] -
516  (DcY(0, ypd, z) + DcY(0, ymd, z)) * state[IDX(0, y, z)] +
517  DcY(0, ymd, z) * state[IDX(0, ym, z)]) /
518  (div_y * SQ(g->dy)) +
519  (DcZ(0, y, zpd) * state[IDX(0, y, zp)] -
520  (DcZ(0, y, zpd) + DcZ(0, y, zmd)) * state[IDX(0, y, z)] +
521  DcZ(0, y, zmd) * state[IDX(0, y, zm)]) /
522  (div_z * SQ(g->dz))) +
523  g->states_cur[IDX(0, y, z)];
524  x = g->size_x - 1;
525  RHS[x] =
526  state[IDX(x, y, z)] +
527  dt * ((DcX(x, y, z) * state[IDX(x - 1, y, z)] - DcX(x, y, z) * state[IDX(x, y, z)]) /
528  (2. * SQ(g->dx)) +
529  (DcY(x, ypd, z) * state[IDX(x, yp, z)] -
530  (DcY(x, ypd, z) + DcY(x, ymd, z)) * state[IDX(x, y, z)] +
531  DcY(x, ymd, z) * state[IDX(x, ym, z)]) /
532  (div_y * SQ(g->dy)) +
533  (DcZ(x, y, zpd) * state[IDX(x, y, zp)] -
534  (DcZ(x, y, zpd) + DcZ(x, y, zmd)) * state[IDX(x, y, z)] +
535  DcZ(x, y, zmd) * state[IDX(x, y, zm)]) /
536  (div_z * SQ(g->dz))) +
537  g->states_cur[IDX(x, y, z)];
538 
539  } else {
540  diag[0] = 1.0;
541  u_diag[0] = 0.0;
542  diag[g->size_x - 1] = 1.0;
543  l_diag[g->size_x - 2] = 0.0;
544  RHS[0] = g->bc->value;
545  RHS[g->size_x - 1] = g->bc->value;
546  }
547 
548  for (x = 1; x < g->size_x - 1; x++) {
549 #ifndef __PGI
550  __builtin_prefetch(&(state[IDX(x + PREFETCH, y, z)]), 0, 1);
551  __builtin_prefetch(&(state[IDX(x + PREFETCH, yp, z)]), 0, 0);
552  __builtin_prefetch(&(state[IDX(x + PREFETCH, ym, z)]), 0, 0);
553 #endif
554 
555  RHS[x] = state[IDX(x, y, z)] +
556  dt * ((DcX(x + 1, y, z) * state[IDX(x + 1, y, z)] -
557  (DcX(x + 1, y, z) + DcX(x, y, z)) * state[IDX(x, y, z)] +
558  DcX(x, y, z) * state[IDX(x - 1, y, z)]) /
559  (2. * SQ(g->dx)) +
560  (DcY(x, ypd, z) * state[IDX(x, yp, z)] -
561  (DcY(x, ypd, z) + DcY(x, ymd, z)) * state[IDX(x, y, z)] +
562  DcY(x, ymd, z) * state[IDX(x, ym, z)]) /
563  (div_y * SQ(g->dy)) +
564  (DcZ(x, y, zpd) * state[IDX(x, y, zp)] -
565  (DcZ(x, y, zpd) + DcZ(x, y, zmd)) * state[IDX(x, y, z)] +
566  DcZ(x, y, zmd) * state[IDX(x, y, zm)]) /
567  (div_z * SQ(g->dz))) +
568  g->states_cur[IDX(x, y, z)];
569  }
570 
571  solve_dd_tridiag(g->size_x, l_diag, diag, u_diag, RHS, scratch);
572 
573  free(diag);
574  free(l_diag);
575  free(u_diag);
576 }
577 
578 
579 /* dg_adi_tort_y performs the second of 3 steps in DG-ADI
580  * g - the parameters and state of the grid
581  * dt - the time step
582  * y - the index for the y plane
583  * z - the index for the z plane
584  * state - where the output of this step is stored
585  * scratch - scratchpad array of doubles, length g->size_y - 1
586  * like dg_adi_y except the grid has a variable tortuosity
587  * g->lambda (but it still has fixed volume fraction)
588  */
589 static void ecs_dg_adi_tort_y(ECS_Grid_node* g,
590  double const dt,
591  int const x,
592  int const z,
593  double const* const state,
594  double* const RHS,
595  double* const scratch) {
596  int y;
597  double* diag;
598  double* l_diag;
599  double* u_diag;
600 
601  /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
602  if (g->bc->type == DIRICHLET &&
603  (x == 0 || z == 0 || x == g->size_x - 1 || z == g->size_z - 1)) {
604  for (y = 0; y < g->size_y; y++)
605  RHS[y] = g->bc->value;
606  return;
607  }
608 
609  if (g->size_y == 1) {
610  if (g->bc->type == DIRICHLET) {
611  RHS[0] = g->bc->value;
612  } else {
613  RHS[0] = state[x + z * g->size_x];
614  }
615  return;
616  }
617 
618  diag = (double*) malloc(g->size_y * sizeof(double));
619  l_diag = (double*) malloc((g->size_y - 1) * sizeof(double));
620  u_diag = (double*) malloc((g->size_y - 1) * sizeof(double));
621 
622  for (y = 1; y < g->size_y - 1; y++) {
623  l_diag[y - 1] = -dt * DcY(x, y, z) / (2. * SQ(g->dy));
624  diag[y] = 1. + dt * (DcY(x, y, z) + DcY(x, y + 1, z)) / (2. * SQ(g->dy));
625  u_diag[y] = -dt * DcY(x, y + 1, z) / (2. * SQ(g->dy));
626  }
627 
628  if (g->bc->type == NEUMANN) {
629  /*zero flux boundary condition*/
630  diag[0] = 1. + 0.5 * dt * DcY(x, 1, z) / SQ(g->dy);
631  u_diag[0] = -0.5 * dt * DcY(x, 1, z) / SQ(g->dy);
632  diag[g->size_y - 1] = 1. + 0.5 * dt * DcY(x, g->size_y - 1, z) / SQ(g->dy);
633  l_diag[g->size_y - 2] = -0.5 * dt * DcY(x, g->size_y - 1, z) / SQ(g->dy);
634 
635 
636  RHS[0] = state[x + z * g->size_x] - dt * ((DcY(x, 1, z) * g->states[IDX(x, 1, z)] -
637  DcY(x, 1, z) * g->states[IDX(x, 0, z)]) /
638  (2. * SQ(g->dy)));
639  y = g->size_y - 1;
640  RHS[y] = state[x + (z + y * g->size_z) * g->size_x] -
641  dt *
642  (DcY(x, y, z) * g->states[IDX(x, y - 1, z)] -
643  DcY(x, y, z) * g->states[IDX(x, y, z)]) /
644  (2. * SQ(g->dy));
645  } else {
646  diag[0] = 1.0;
647  u_diag[0] = 0.0;
648  diag[g->size_y - 1] = 1.0;
649  l_diag[g->size_y - 2] = 0.0;
650  RHS[0] = g->bc->value;
651  RHS[g->size_y - 1] = g->bc->value;
652  }
653 
654 
655  for (y = 1; y < g->size_y - 1; y++) {
656 #ifndef __PGI
657  __builtin_prefetch(&state[x + (z + (y + PREFETCH) * g->size_z) * g->size_x], 0, 0);
658  __builtin_prefetch(&(g->states[IDX(x, y + PREFETCH, z)]), 0, 1);
659 #endif
660  RHS[y] = state[x + (z + y * g->size_z) * g->size_x] -
661  dt *
662  (DcY(x, y + 1, z) * g->states[IDX(x, y + 1, z)] -
663  (DcY(x, y + 1, z) + DcY(x, y, z)) * g->states[IDX(x, y, z)] +
664  DcY(x, y, z) * g->states[IDX(x, y - 1, z)]) /
665  (2. * SQ(g->dy));
666  }
667 
668  solve_dd_tridiag(g->size_y, l_diag, diag, u_diag, RHS, scratch);
669 
670  free(diag);
671  free(l_diag);
672  free(u_diag);
673 }
674 
675 /* dg_adi_tort_z performs the second of 3 steps in DG-ADI
676  * g - the parameters and state of the grid
677  * dt - the time step
678  * y - the index for the y plane
679  * z - the index for the z plane
680  * state - where the output of this step is stored
681  * scratch - scratchpad array of doubles, length g->size_z - 1
682  * like dg_adi_z except the grid has a variable tortuosity
683  * g->lambda (but it still has fixed volume fraction)
684  */
685 static void ecs_dg_adi_tort_z(ECS_Grid_node* g,
686  double const dt,
687  int const x,
688  int const y,
689  double const* const state,
690  double* const RHS,
691  double* const scratch) {
692  int z;
693  double* diag;
694  double* l_diag;
695  double* u_diag;
696 
697  /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
698  if (g->bc->type == DIRICHLET &&
699  (x == 0 || y == 0 || x == g->size_x - 1 || y == g->size_y - 1)) {
700  for (z = 0; z < g->size_z; z++)
701  RHS[z] = g->bc->value;
702  return;
703  }
704 
705  if (g->size_z == 1) {
706  if (g->bc->type == DIRICHLET) {
707  RHS[0] = g->bc->value;
708  } else {
709  RHS[0] = state[y + g->size_y * x];
710  }
711  return;
712  }
713  diag = (double*) malloc(g->size_z * sizeof(double));
714  l_diag = (double*) malloc((g->size_z - 1) * sizeof(double));
715  u_diag = (double*) malloc((g->size_z - 1) * sizeof(double));
716 
717 
718  for (z = 1; z < g->size_z - 1; z++) {
719  l_diag[z - 1] = -dt * DcZ(x, y, z) / (2. * SQ(g->dz));
720  diag[z] = 1. + dt * (DcZ(x, y, z) + DcZ(x, y, z + 1)) / (2. * SQ(g->dz));
721  u_diag[z] = -dt * DcZ(x, y, z + 1) / (2. * SQ(g->dz));
722  }
723 
724  if (g->bc->type == NEUMANN) {
725  /*zero flux boundary condition*/
726  diag[0] = 1. + 0.5 * dt * DcZ(x, y, 1) / SQ(g->dz);
727  u_diag[0] = -0.5 * dt * DcZ(x, y, 1) / SQ(g->dz);
728  diag[g->size_z - 1] = 1. + 0.5 * dt * DcZ(x, y, g->size_z - 1) / SQ(g->dz);
729  l_diag[g->size_z - 2] = -0.5 * dt * DcZ(x, y, g->size_z - 1) / SQ(g->dz);
730 
731  RHS[0] = state[y + g->size_y * (x * g->size_z)] -
732  dt * ((DcZ(x, y, 1) * g->states[IDX(x, y, 1)] -
733  DcZ(x, y, 1) * g->states[IDX(x, y, 0)]) /
734  (2. * SQ(g->dz)));
735  z = g->size_z - 1;
736  RHS[z] = state[y + g->size_y * (x * g->size_z + z)] -
737  dt *
738  (DcZ(x, y, z) * g->states[IDX(x, y, z - 1)] -
739  DcZ(x, y, z) * g->states[IDX(x, y, z)]) /
740  (2. * SQ(g->dz));
741  } else {
742  diag[0] = 1.0;
743  u_diag[0] = 0.0;
744  diag[g->size_z - 1] = 1.0;
745  l_diag[g->size_z - 2] = 0.0;
746  RHS[0] = g->bc->value;
747  RHS[g->size_z - 1] = g->bc->value;
748  }
749 
750 
751  for (z = 1; z < g->size_z - 1; z++) {
752  RHS[z] = state[y + g->size_y * (x * g->size_z + z)] -
753  dt *
754  (DcZ(x, y, z + 1) * g->states[IDX(x, y, z + 1)] -
755  (DcZ(x, y, z + 1) + DcZ(x, y, z)) * g->states[IDX(x, y, z)] +
756  DcZ(x, y, z) * g->states[IDX(x, y, z - 1)]) /
757  (2. * SQ(g->dz));
758  }
759 
760  solve_dd_tridiag(g->size_z, l_diag, diag, u_diag, RHS, scratch);
761 
762  free(diag);
763  free(l_diag);
764  free(u_diag);
765 }
766 
767 /*DG-ADI implementation the 3 step process to diffusion species in grid g by time step *dt_ptr
768  * g - the state and parameters
769  * like dg_adi except the grid node g has variable tortuosity g.lambda (but fixed volume
770  * fraction)
771  */
772 void ecs_set_adi_tort(ECS_Grid_node* g) {
773  g->ecs_adi_dir_x->ecs_dg_adi_dir = ecs_dg_adi_tort_x;
774  g->ecs_adi_dir_y->ecs_dg_adi_dir = ecs_dg_adi_tort_y;
775  g->ecs_adi_dir_z->ecs_dg_adi_dir = ecs_dg_adi_tort_z;
776 }
777 
778 
779 /*****************************************************************************
780  *
781  * Begin variable step code
782  *
783  *****************************************************************************/
784 
785 void _rhs_variable_step_helper_tort(Grid_node* g, double const* const states, double* ydot) {
786  int num_states_x = g->size_x, num_states_y = g->size_y, num_states_z = g->size_z;
787  double dx = g->dx, dy = g->dy, dz = g->dz;
788  int i, j, k, stop_i, stop_j, stop_k;
789  double rate_x = 1. / (dx * dx);
790  double rate_y = 1. / (dy * dy);
791  double rate_z = 1. / (dz * dz);
792 
793  /*indices*/
794  int index, prev_i, prev_j, prev_k, next_i, next_j, next_k;
795  int xpd, xmd, ypd, ymd, zpd, zmd;
796  double div_x, div_y, div_z;
797 
798  /* Euler advance x, y, z (all points) */
799  stop_i = num_states_x - 1;
800  stop_j = num_states_y - 1;
801  stop_k = num_states_z - 1;
802  if (g->bc->type == NEUMANN) {
803  for (i = 0,
804  index = 0,
805  prev_i = num_states_y * num_states_z,
806  next_i = num_states_y * num_states_z;
807  i < num_states_x;
808  i++) {
809  /*Zero flux boundary conditions*/
810  div_x = (i == 0 || i == stop_i) ? 2. : 1.;
811  xpd = (i == stop_i) ? i : i + 1;
812  xmd = (i == 0) ? 1 : i;
813 
814  for (j = 0, prev_j = index + num_states_z, next_j = index + num_states_z;
815  j < num_states_y;
816  j++) {
817  div_y = (j == 0 || j == stop_j) ? 2. : 1.;
818  ypd = (j == stop_j) ? j : j + 1;
819  ymd = (j == 0) ? 1 : j;
820 
821  for (k = 0, prev_k = index + 1, next_k = index + 1; k < num_states_z;
822  k++, index++, prev_i++, next_i++, prev_j++, next_j++) {
823  div_z = (k == 0 || k == stop_k) ? 2. : 1.;
824  zpd = (k == stop_k) ? k : k + 1;
825  zmd = (k == 0) ? 1 : k;
826 
827  if (stop_i > 0) {
828  ydot[index] += rate_x *
829  (DcX(xmd, j, k) * states[prev_i] -
830  (DcX(xmd, j, k) + DcX(xpd, j, k)) * states[index] +
831  DcX(xpd, j, k) * states[next_i]) /
832  div_x;
833  }
834  if (stop_j > 0) {
835  ydot[index] += rate_y *
836  (DcY(i, ymd, k) * states[prev_j] -
837  (DcY(i, ymd, k) + DcY(i, ypd, k)) * states[index] +
838  DcY(i, ypd, k) * states[next_j]) /
839  div_y;
840  }
841 
842  if (stop_k > 0) {
843  ydot[index] += rate_z *
844  (DcZ(i, j, zmd) * states[prev_k] -
845  (DcZ(i, j, zmd) + DcZ(i, j, zpd)) * states[index] +
846  DcZ(i, j, zpd) * states[next_k]) /
847  div_z;
848  }
849 
850  next_k = (k == stop_k - 1) ? index : index + 2;
851  prev_k = index;
852  }
853  prev_j = index - num_states_z;
854  next_j = (j == stop_j - 1) ? prev_j : index + num_states_z;
855  }
856  prev_i = index - num_states_y * num_states_z;
857  next_i = (i == stop_i - 1) ? prev_i : index + num_states_y * num_states_z;
858  }
859  } else {
860  for (i = 0, index = 0, prev_i = 0, next_i = num_states_y * num_states_z; i < num_states_x;
861  i++) {
862  for (j = 0, prev_j = index - num_states_z, next_j = index + num_states_z;
863  j < num_states_y;
864  j++) {
865  for (k = 0, prev_k = index - 1, next_k = index + 1; k < num_states_z;
866  k++, index++, prev_i++, next_i++, prev_j++, next_j++, next_k++, prev_k++) {
867  if (i == 0 || i == stop_i || j == 0 || j == stop_j || k == 0 || k == stop_k) {
868  // set to zero to prevent currents altering concentrations at the boundary
869  ydot[index] = 0;
870  } else {
871  ydot[index] += rate_x * (DcX(i, j, k) * states[prev_i] -
872  (DcX(i, j, k) + DcX(i + 1, j, k)) * states[index] +
873  DcX(i + 1, j, k) * states[next_i]);
874 
875  ydot[index] += rate_y * (DcY(i, j, k) * states[prev_j] -
876  (DcY(i, j, k) + DcY(i, j + 1, k)) * states[index] +
877  DcY(i, j + 1, k) * states[next_j]);
878 
879  ydot[index] += rate_z * (DcZ(i, j, k) * states[prev_k] -
880  (DcZ(i, j, k) + DcZ(i, j, k + 1)) * states[index] +
881  DcZ(i, j, k + 1) * states[next_k]);
882  }
883  }
884  prev_j = index - num_states_z;
885  next_j = index + num_states_z;
886  }
887  prev_i = index - num_states_y * num_states_z;
888  next_i = index + num_states_y * num_states_z;
889  }
890  }
891 }
892 
893 
894 void _rhs_variable_step_helper_vol(Grid_node* g, double const* const states, double* ydot) {
895  int num_states_x = g->size_x, num_states_y = g->size_y, num_states_z = g->size_z;
896  double dx = g->dx, dy = g->dy, dz = g->dz;
897  int i, j, k, stop_i, stop_j, stop_k;
898  double rate_x = g->dc_x / (dx * dx);
899  double rate_y = g->dc_y / (dy * dy);
900  double rate_z = g->dc_z / (dz * dz);
901 
902 
903  /*indices*/
904  int index, prev_i, prev_j, prev_k, next_i, next_j, next_k;
905 
906  /*zero flux boundary conditions*/
907  int xpd, xmd, ypd, ymd, zpd, zmd;
908  double div_x, div_y, div_z;
909 
910  /* Euler advance x, y, z (all points) */
911  stop_i = num_states_x - 1;
912  stop_j = num_states_y - 1;
913  stop_k = num_states_z - 1;
914  if (g->bc->type == NEUMANN) {
915  for (i = 0,
916  index = 0,
917  prev_i = num_states_y * num_states_z,
918  next_i = num_states_y * num_states_z;
919  i < num_states_x;
920  i++) {
921  /*Zero flux boundary conditions*/
922  div_x = (i == 0 || i == stop_i) ? 2. : 1.;
923  xpd = (i == stop_i) ? i : i + 1;
924  xmd = (i == 0) ? 1 : i;
925 
926  for (j = 0, prev_j = index + num_states_z, next_j = index + num_states_z;
927  j < num_states_y;
928  j++) {
929  div_y = (j == 0 || j == stop_j) ? 2. : 1.;
930  ypd = (j == stop_j) ? j : j + 1;
931  ymd = (j == 0) ? 1 : j;
932 
933  for (k = 0, prev_k = index + 1, next_k = index + 1; k < num_states_z;
934  k++, index++, prev_i++, next_i++, prev_j++, next_j++) {
935  div_z = (k == 0 || k == stop_k) ? 2. : 1.;
936  zpd = (k == stop_k) ? k : k + 1;
937  zmd = (k == 0) ? 1 : k;
938 
939  /*x-direction*/
940  if (stop_i > 0) {
941  ydot[index] += rate_x *
942  (FLUX(next_i, index) * PERM(xpd, j, k) +
943  FLUX(prev_i, index) * PERM(xmd, j, k)) /
944  (VOLFRAC(index) * div_x);
945  }
946 
947  /*y-direction*/
948  if (stop_j > 0) {
949  ydot[index] += rate_y *
950  (FLUX(next_j, index) * PERM(i, ypd, k) +
951  FLUX(prev_j, index) * PERM(i, ymd, k)) /
952  (VOLFRAC(index) * div_y);
953  }
954 
955  /*z-direction*/
956  if (stop_k > 0) {
957  ydot[index] += rate_z *
958  (FLUX(next_k, index) * PERM(i, j, zpd) +
959  FLUX(prev_k, index) * PERM(i, j, zmd)) /
960  (VOLFRAC(index) * div_z);
961  }
962 
963  next_k = (k == stop_k - 1) ? index : index + 2;
964  prev_k = index;
965  }
966  prev_j = index - num_states_z;
967  next_j = (j == stop_j - 1) ? prev_j : index + num_states_z;
968  }
969  prev_i = index - num_states_y * num_states_z;
970  next_i = (i == stop_i - 1) ? prev_i : index + num_states_y * num_states_z;
971  }
972  } else {
973  for (i = 0, index = 0, prev_i = 0, next_i = num_states_y * num_states_z; i < num_states_x;
974  i++) {
975  for (j = 0, prev_j = index - num_states_z, next_j = index + num_states_z;
976  j < num_states_y;
977  j++) {
978  for (k = 0, prev_k = index - 1, next_k = index + 1; k < num_states_z;
979  k++, index++, prev_i++, next_i++, prev_j++, next_j++, next_k++, prev_k++) {
980  if (i == 0 || i == stop_i || j == 0 || j == stop_j || k == 0 || k == stop_k) {
981  // set to zero to prevent currents altering concentrations at the boundary
982  ydot[index] = 0;
983  } else {
984  /*x-direction*/
985  ydot[index] += rate_x *
986  (FLUX(next_i, index) * PERM(i + 1, j, k) +
987  FLUX(prev_i, index) * PERM(i, j, k)) /
988  (VOLFRAC(index));
989 
990  /*y-direction*/
991  ydot[index] += rate_y *
992  (FLUX(next_j, index) * PERM(i, j + 1, k) +
993  FLUX(prev_j, index) * PERM(i, j, k)) /
994  (VOLFRAC(index));
995 
996  /*z-direction*/
997  ydot[index] += rate_z *
998  (FLUX(next_k, index) * PERM(i, j, k + 1) +
999  FLUX(prev_k, index) * PERM(i, j, k)) /
1000  (VOLFRAC(index));
1001  }
1002  }
1003  prev_j = index - num_states_z;
1004  next_j = index + num_states_z;
1005  }
1006  prev_i = index - num_states_y * num_states_z;
1007  next_i = index + num_states_y * num_states_z;
1008  }
1009  }
1010 }
ECS_Grid_node
Definition: grids.h:161
ECS_Grid_node::ecs_adi_dir_z
struct ECSAdiDirection * ecs_adi_dir_z
Definition: grids.h:184
ECS_Grid_node::ecs_adi_dir_y
struct ECSAdiDirection * ecs_adi_dir_y
Definition: grids.h:183
ECS_Grid_node::ecs_adi_dir_x
struct ECSAdiDirection * ecs_adi_dir_x
Definition: grids.h:182
Grid_node
Definition: grids.h:82
Grid_node::size_y
int size_y
Definition: grids.h:92
Grid_node::size_x
int size_x
Definition: grids.h:91
Grid_node::states_cur
double * states_cur
Definition: grids.h:90
Grid_node::states
double * states
Definition: grids.h:86
Grid_node::bc
BoundaryConditions * bc
Definition: grids.h:102
Grid_node::dy
double dy
Definition: grids.h:98
Grid_node::dc_z
double dc_z
Definition: grids.h:96
Grid_node::dc_y
double dc_y
Definition: grids.h:95
Grid_node::size_z
int size_z
Definition: grids.h:93
Grid_node::dc_x
double dc_x
Definition: grids.h:94
Grid_node::dz
double dz
Definition: grids.h:99
Grid_node::dx
double dx
Definition: grids.h:97
i
#define i
Definition: md1redef.h:19
k
static RNG::key_type k
Definition: nrnran123.cpp:9
IDX
#define IDX(x, y, z)
Definition: grids.h:14
VOLFRAC
#define VOLFRAC(idx)
Definition: grids.h:17
NEUMANN
#define NEUMANN
Definition: grids.h:28
PERM
#define PERM(x, y, z)
Definition: grids.h:19
DIRICHLET
#define DIRICHLET
Definition: grids.h:29
SQ
#define SQ(x)
Definition: grids.h:20
ALPHA
#define ALPHA(x, y, z)
Definition: grids.h:16
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
j
size_t j
Definition: nrngsl_real_radix2.cpp:50
index
short index
Definition: cabvars.h:11
nrnwrap_Python.h
states
double * states
Definition: rxd.cpp:75
PREFETCH
#define PREFETCH
Definition: rxd.h:8
ecs_dg_adi_vol_y
static void ecs_dg_adi_vol_y(ECS_Grid_node *g, double const dt, int const x, int const z, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:232
FLUX
#define FLUX(pidx, idx)
Definition: rxd_vol.cpp:37
ecs_dg_adi_vol_x
static void ecs_dg_adi_vol_x(ECS_Grid_node *g, const double dt, const int y, const int z, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:88
ecs_set_adi_tort
void ecs_set_adi_tort(ECS_Grid_node *g)
Definition: rxd_vol.cpp:772
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_vol.cpp:52
Fzz
#define Fzz(z1, z2)
Definition: rxd_vol.cpp:31
DcX
#define DcX(x, y, z)
Definition: rxd_vol.cpp:9
Fyy
#define Fyy(y1, y2)
Definition: rxd_vol.cpp:27
_rhs_variable_step_helper_tort
void _rhs_variable_step_helper_tort(Grid_node *g, double const *const states, double *ydot)
Definition: rxd_vol.cpp:785
ecs_dg_adi_tort_z
static void ecs_dg_adi_tort_z(ECS_Grid_node *g, double const dt, int const x, int const y, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:685
Fxx
#define Fxx(x1, x2)
Definition: rxd_vol.cpp:15
_rhs_variable_step_helper_vol
void _rhs_variable_step_helper_vol(Grid_node *g, double const *const states, double *ydot)
Definition: rxd_vol.cpp:894
Fxy
#define Fxy(y1, y1d, y2)
Definition: rxd_vol.cpp:19
DcZ
#define DcZ(x, y, z)
Definition: rxd_vol.cpp:11
NUM_THREADS
int NUM_THREADS
Definition: rxd.cpp:34
ecs_dg_adi_tort_x
static void ecs_dg_adi_tort_x(ECS_Grid_node *g, const double dt, const int y, const int z, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:427
Fxz
#define Fxz(z1, z1d, z2)
Definition: rxd_vol.cpp:23
ecs_dg_adi_vol_z
static void ecs_dg_adi_vol_z(ECS_Grid_node *g, double const dt, int const x, int const y, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:326
ecs_set_adi_vol
void ecs_set_adi_vol(ECS_Grid_node *g)
Definition: rxd_vol.cpp:410
ecs_dg_adi_tort_y
static void ecs_dg_adi_tort_y(ECS_Grid_node *g, double const dt, int const x, int const z, double const *const state, double *const RHS, double *const scratch)
Definition: rxd_vol.cpp:589
DcY
#define DcY(x, y, z)
Definition: rxd_vol.cpp:10
BoundaryConditions::value
double value
Definition: grids.h:79
BoundaryConditions::type
unsigned char type
Definition: grids.h:78
ECSAdiDirection::ecs_dg_adi_dir
void(* ecs_dg_adi_dir)(ECS_Grid_node *, const double, const int, const int, double const *const, double *const, double *const)
Definition: grids.h:229