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