PLUTO
Macros | Functions
fargo.c File Reference

Solves the linear transport step of the FARGO-MHD algorithm. More...

#include "pluto.h"
Include dependency graph for fargo.c:

Go to the source code of this file.

Macros

#define FARGO_MOD(s)   ((s) < SBEG ? (s)+NS:((s) > SEND ? (s)-NS:(s)))
 
#define MAX_BUF_SIZE   64
 
#define s   k
 

Functions

static void FARGO_Flux (double *, double *, double)
 
static void FARGO_ParallelExchange (Data_Arr U, Data_Arr Us, double ***rbl[], RBox *, double ***rbu[], RBox *, Grid *grid)
 
void FARGO_ShiftSolution (Data_Arr U, Data_Arr Us, Grid *grid)
 

Detailed Description

Solves the linear transport step of the FARGO-MHD algorithm.

Implements a uniform shift of the conservative fluid variables in the direction of orbital motion which is assumed to be $ x_2=y,\phi$ (for Cartesian or polar coordinates) or $ x_3 = \phi$ (for spherical coordinates). The shift is converted into an integer part "m" and a fractional remainder "eps".
Parallelization is performed by adding extra data layers above and below the current data set so as to enlarge the required computational stencil in the orbital direction. The additional buffers are received from the processors lying, respectively, below and above and their size changes dynamically from one step to the next.

Reference

  • "A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code"
    Mignone et al, A&A (2012) 545, A152 (–> Section 2.4)
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
G. Muscianisi (g.mus.nosp@m.cian.nosp@m.isi@c.nosp@m.inec.nosp@m.a.it)
Date
Aug 23, 2015

Definition in file fargo.c.

Macro Definition Documentation

#define FARGO_MOD (   s)    ((s) < SBEG ? (s)+NS:((s) > SEND ? (s)-NS:(s)))

Definition at line 32 of file fargo.c.

#define MAX_BUF_SIZE   64

Definition at line 33 of file fargo.c.

#define s   k

Function Documentation

void FARGO_Flux ( double *  q,
double *  flx,
double  eps 
)
static

Compute the transport flux corresponding to a fractional increment eps.

Parameters
[in]q1D array of a conserved quantity
[out]flxthe conservative flux
[in]epsthe fractional increment

The following MAPLE script may be used to check the correctness of implementation

1 restart;
2 P0 := 3/2*q[i] - (q[p]+q[m])/4;
3 P1 := q[p] - q[m];
4 P2 := 3*(q[p] - 2*q[i] + q[m]);
5 
6 FS[p] := P0 + P1/2*(1 - epsilon) + P2*(1/4 - epsilon/2 + epsilon^2/3);
7 FS[m] := P0 - P1/2*(1 - epsilon) + P2*(1/4 - epsilon/2 + epsilon^2/3);
8 
9 dq := q[p] - q[m];
10 d2q := q[p] - 2*q[i] + q[m];
11 FF[p] := q[p] - epsilon/2*(dq + d2q*(3 - 2*epsilon));
12 FF[m] := q[m] + epsilon/2*(dq - d2q*(3 - 2*epsilon));

Definition at line 770 of file fargo.c.

797 {
798  int s;
799  double dqc, d2q, dqp, dqm;
800  static double *dq, *dql, *qp, *qm;
801 
802  if (dq == NULL){
803  dq = ARRAY_1D(NMAX_POINT, double);
804  dql = ARRAY_1D(NMAX_POINT, double);
805  }
806 
807 /* -- compute limited slopes -- */
808 
809  dq[0] = q[1] - q[0];
810  for (s = 1; s < NS_TOT-1; s++){
811  dq[s] = q[s+1] - q[s];
812  if (dq[s]*dq[s-1] > 0.0){
813  dqc = 0.5*(dq[s] + dq[s-1]);
814  d2q = 2.0*(fabs(dq[s]) < fabs(dq[s-1]) ? dq[s]:dq[s-1]);
815  dql[s] = fabs(d2q) < fabs(dqc) ? d2q: dqc;
816  }else dql[s] = 0.0;
817  }
818 
819 /* vanleer_lim(dq,dq-1,dql,JBEG-1,JEND+1,NULL); */
820 
821  #if FARGO_ORDER == 2 /* MUSCL-Hancock, 2nd order */
822  if (eps > 0.0){
823  for (s=SBEG-1; s<=SEND; s++) flx[s] = q[s] + 0.5*dql[s]*(1.0 - eps);
824  }else { /* -- remember: eps > 0 by definition -- */
825  for (s=SBEG-1; s<=SEND; s++) flx[s] = q[s+1] - 0.5*dql[s+1]*(1.0 + eps);
826  }
827  #elif FARGO_ORDER == 3 /* PPM, 3rd order */
828  if (qp == NULL){
829  qp = ARRAY_1D(NMAX_POINT, double);
830  qm = ARRAY_1D(NMAX_POINT, double);
831  }
832  for (s = SBEG-1; s <= SEND+1; s++){
833  dqp = 0.5*dq[s] - (dql[s+1] - dql[s])/6.0;
834  dqm = -0.5*dq[s-1] + (dql[s-1] - dql[s])/6.0;
835  if (dqp*dqm > 0.0) dqp = dqm = 0.0;
836  else{
837  if (fabs(dqp) >= 2.0*fabs(dqm)) dqp = -2.0*dqm;
838  if (fabs(dqm) >= 2.0*fabs(dqp)) dqm = -2.0*dqp;
839  }
840  qp[s] = q[s] + dqp;
841  qm[s] = q[s] + dqm;
842  }
843 
844  if (eps > 0.0){
845  for (s = SBEG-1; s <= SEND; s++){
846  dqc = qp[s] - qm[s];
847  d2q = qp[s] - 2.0*q[s] + qm[s];
848  flx[s] = qp[s] - 0.5*eps*(dqc + d2q*(3.0 - 2.0*eps));
849  }
850  }else{
851  for (s = SBEG-1; s <= SEND; s++){
852  dqc = qp[s+1] - qm[s+1];
853  d2q = qp[s+1] - 2.0*q[s+1] + qm[s+1];
854  flx[s] = qm[s+1] - 0.5*eps*(dqc - d2q*(3.0 + 2.0*eps));
855  }
856  }
857  #else
858  print1 ("! FARGO_ORDER should be either 2 or 3\n");
859  QUIT_PLUTO(1);
860  #endif
861 
862 }
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
SEND
#define SEND
Definition: fargo.h:57
eps
#define eps
s
#define s
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
SBEG
#define SBEG
Definition: fargo.h:56
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
NS_TOT
#define NS_TOT
Definition: fargo.h:59

Here is the call graph for this function:

Here is the caller graph for this function:

static void FARGO_ParallelExchange ( Data_Arr  U,
Data_Arr  Us,
double ***  rbl[],
RBox *  ,
double ***  rbu[],
RBox *  ,
Grid *  grid 
)
static

Here is the caller graph for this function:

void FARGO_ShiftSolution ( Data_Arr  U,
Data_Arr  Us,
Grid *  grid 
)

Shifts conserved variables in the orbital direction.

Parameters
[in,out]Ua 3D array of conserved, zone-centered values
[in,out]Usa 3D array of staggered magnetic fields
[in]gridpointer to Grid structure;
Returns
This function has no return value.
Todo:
Optimization: avoid using too many if statement like on nproc_s > 1 or == 1

Definition at line 42 of file fargo.c.

59 {
60  int i, j, k, nv, ngh, nvar_c;
61  int sm, m;
62  static int mmax, mmin; /* Used to determine buffer size from time to time */
63  int nbuf_lower, nbuf_upper, nproc_s = 1;
64  double *dx, *dy, *dz, dphi;
65  double *x, *xp, *y, *yp, *z, *zp, w, Lphi, dL, eps;
66  double **wA, ***Uc[NVAR];
67  static double *q00, *flux00;
68  double *q, *flux;
69  static double ***Ez, ***Ex;
70  #ifdef PARALLEL
71  double ***upper_buf[NVAR]; /* upper layer buffer */
72  double ***lower_buf[NVAR]; /* lower layer buffer */
73  RBox lbox, ubox; /* lower and upper box structures */
74  #endif
75 
76  #if GEOMETRY == CYLINDRICAL
77  print1 ("! FARGO cannot be used in this geometry\n");
78  QUIT_PLUTO(1);
79  #endif
80 
81 /* -----------------------------------------------------------------
82  1. Allocate memory for static arrays.
83  In parallel mode, one-dimensional arrays must be large enough
84  to contain data values coming from neighbour processors.
85  For this reason we augment them with an extra zone buffer
86  containing MAX_BUF_SIZE cells on both sides.
87  In this way, the array indices for q can be negative.
88 
89 
90  <..MAX_BUF_SIZE..>|---|---| ... |---|----|<..MAX_BUF_SIZE..>
91  0 1 NX2-1
92  ----------------------------------------------------------------- */
93 
94  if (q00 == NULL){
95  #if GEOMETRY == SPHERICAL && DIMENSIONS == 3
96  q00 = ARRAY_1D(NX3_TOT + 2*MAX_BUF_SIZE,double);
97  flux00 = ARRAY_1D(NX3_TOT + 2*MAX_BUF_SIZE,double);
98  #else
99  q00 = ARRAY_1D(NX2_TOT + 2*MAX_BUF_SIZE,double);
100  flux00 = ARRAY_1D(NX2_TOT + 2*MAX_BUF_SIZE,double);
101  #endif
102  #ifdef STAGGERED_MHD
103  Ez = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
104  Ex = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
105  #endif
106 
107  /* ------------------------------------------------------
108  The value of mmax and mmin is kept from call to call
109  because they are used to determine the size of the
110  send/receive buffers each time step.
111  At the very beginning, we initialize them to the
112  largest possible value.
113  ------------------------------------------------------ */
114 
115  mmax = MIN(NS, MAX_BUF_SIZE-1) - grid[SDIR].nghost - 2;
116  mmin = mmax;
117  }
118 
119  q = q00 + MAX_BUF_SIZE;
120  flux = flux00 + MAX_BUF_SIZE;
121 
122 /* ------------------------------------------------------------
123  2. Get the pointer to the orbital speed array
124  ------------------------------------------------------------ */
125 
126  wA = FARGO_GetVelocity();
127 
128  nproc_s = grid[SDIR].nproc; /* -- # of processors in shift dir -- */
129  ngh = grid[SDIR].nghost;
130  Lphi = g_domEnd[SDIR] - g_domBeg[SDIR];
131 
132  x = grid[IDIR].x; xp = grid[IDIR].xr; dx = grid[IDIR].dx;
133  y = grid[JDIR].x; yp = grid[JDIR].xr; dy = grid[JDIR].dx;
134  z = grid[KDIR].x; zp = grid[KDIR].xr; dz = grid[KDIR].dx;
135 
136  dphi = grid[SDIR].dx[SBEG];
137 
138 /* -----------------------------------------------------------
139  3. Add source term (if any)
140  ----------------------------------------------------------- */
141 
142 #if (HAVE_ENERGY) && (defined SHEARINGBOX)
143  FARGO_Source(U, g_dt, grid);
144 #endif
145 
146 /* -----------------------------------------------------------
147  4. Fargo parallelization is performed by adding extra
148  data layers above and below the current data set so
149  as to enlarge the required computational stencil in
150  the orbital direction:
151 
152  ^
153  | |recv_lower|
154  | +----------+
155  | | |
156  | | U |
157  | | |
158  | +----------+
159  | |recv_upper|
160 
161  The buffers recv_upper and recv_lower are received from
162  the processors lying, respectively, below and above and
163  their size changes dynamically from one step to the next.
164  ----------------------------------------------------------- */
165 
166  #ifdef PARALLEL
167  if (nproc_s > 1){
168  nbuf_lower = abs(mmin) + ngh + 1;
169  nbuf_upper = mmax + ngh + 1;
170 
171  /* -- check that buffer is not larger than processor size -- */
172 
173  if (nbuf_upper > NS || nbuf_lower > NS){
174  print ("! FARGO_ShiftSolution: buffer size exceeds processsor size\n");
175  QUIT_PLUTO(1);
176  }
177 
178  /* -- check that buffer is less than MAX_BUF_SIZE -- */
179 
180  if (nbuf_upper > MAX_BUF_SIZE || nbuf_lower > MAX_BUF_SIZE){
181  print ("! FARGO_ShiftSolution: buffer size exceeds maximum length\n");
182  QUIT_PLUTO(1);
183  }
184 
185  /* -- set grid indices of the lower & upper boxes -- */
186 
187  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
188 
189  lbox.ib = IBEG; lbox.jb = 0; lbox.kb = KBEG;
190  lbox.ie = IEND; lbox.je = nbuf_lower-1; lbox.ke = KEND;
191  lbox.di = lbox.dj = lbox.dk = 1;
192 
193  ubox.ib = IBEG; ubox.jb = 0; ubox.kb = KBEG;
194  ubox.ie = IEND; ubox.je = nbuf_upper-1; ubox.ke = KEND;
195  ubox.di = ubox.dj = ubox.dk = 1;
196 
197  #ifdef STAGGERED_MHD /* -- enlarge boxes to fit staggered fields --*/
198  D_EXPAND(lbox.ib--; ubox.ib--; ,
199  ; ,
200  lbox.kb--; ubox.kb--;)
201  #endif
202 
203  #elif GEOMETRY == SPHERICAL
204 
205  lbox.ib = IBEG; lbox.jb = JBEG; lbox.kb = 0;
206  lbox.ie = IEND; lbox.je = JEND; lbox.ke = nbuf_lower-1;
207  lbox.di = lbox.dj = lbox.dk = 1;
208 
209  ubox.ib = IBEG; ubox.jb = JBEG; ubox.kb = 0;
210  ubox.ie = IEND; ubox.je = JEND; ubox.ke = nbuf_upper-1;
211  ubox.di = ubox.dj = ubox.dk = 1;
212 
213  #ifdef STAGGERED_MHD /* -- enlarge boxes to fit staggered fields --*/
214  D_EXPAND(lbox.ib--; ubox.ib--; ,
215  lbox.jb--; ubox.jb--; ,
216  ; )
217  #endif
218 
219  #endif /* GEOMETRY */
220 
221  /* -- allocate memory -- */
222 
223  for (nv = 0; nv < NVAR; nv++){
224  upper_buf[nv] = ArrayBox(ubox.kb, ubox.ke,
225  ubox.jb, ubox.je,
226  ubox.ib, ubox.ie);
227  lower_buf[nv] = ArrayBox(lbox.kb, lbox.ke,
228  lbox.jb, lbox.je,
229  lbox.ib, lbox.ie);
230  }
231 
232  /* -- exchange data values between neighbour processes -- */
233 
234  FARGO_ParallelExchange(U, Us, lower_buf, &lbox, upper_buf, &ubox, grid);
235  }
236  #endif /* PARALLEL */
237 
238 /* -------------------------------------------------
239  5. Shift cell-centered quantities
240  ------------------------------------------------- */
241 
242  mmax = mmin = 0;
243  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
244  KDOM_LOOP(k) IDOM_LOOP(i) {
245  #else
246  JDOM_LOOP(j) IDOM_LOOP(i) {
247  #endif
248 
249  #if GEOMETRY == CARTESIAN
250  w = wA[k][i];
251  #elif GEOMETRY == POLAR
252  w = wA[k][i]/x[i];
253  #elif GEOMETRY == SPHERICAL
254  w = wA[j][i]/(x[i]*sin(y[j]));
255  #endif
256 
257  /* --------------------------------------------------
258  5a. Obtain shift in terms of an integer number of
259  cells (m) plus a remainder eps.
260  The integer part is obtained by rounding
261  dL/dphi to the nearest integer.
262  Examples:
263 
264  dL/dphi = 3.4 --> m = 3, eps = 0.4
265  dL/dphi = 4.7 --> m = 5, eps = -0.3
266  dL/dphi = -0.8 --> m = -1, eps = 0.2
267  -------------------------------------------------- */
268 
269  dL = w*g_dt; /* spatial shift */
270  dL = fmod(dL, Lphi);
271  m = (int)floor(dL/dphi + 0.5); /* nearest integer part */
272  eps = dL/dphi - m; /* fractional remainder */
273 
274  /* -- check if shift is compatible with estimated buffer size -- */
275 
276  if (nproc_s > 1){
277  if (w > 0.0 && (m + ngh) > nbuf_upper){
278  print1 ("! FARGO_ShiftSolution: positive shift too large\n");
279  QUIT_PLUTO(1);
280  }
281  if (w < 0.0 && (-m + ngh) > nbuf_lower){
282  print1 ("! FARGO_ShiftSolution: negative shift too large\n");
283  QUIT_PLUTO(1);
284  }
285  }
286 
287  mmax = MAX(m, mmax);
288  mmin = MIN(m, mmin);
289 
290  /* ------------------------------------
291  5b. Start main loop on variables
292  ------------------------------------ */
293 
294  for (nv = NVAR; nv--; ){
295 
296  #ifdef STAGGERED_MHD
297  D_EXPAND(if (nv == BX1) continue; ,
298  if (nv == BX2) continue; ,
299  if (nv == BX3) continue;)
300  #endif
301 
302  /* -- copy 3D data into 1D array -- */
303 
304  SDOM_LOOP(s) q[s] = U[k][j][i][nv];
305 
306  if (nproc_s > 1){ /* -- copy values from lower and upper buffers -- */
307  #ifdef PARALLEL
308  for (s = SBEG-1; s >= SBEG - nbuf_upper; s--){
309  q[s] = FARGO_ARRAY_INDEX(upper_buf[nv], SBEG-1-s, k, j, i);
310  }
311  for (s = SEND+1; s <= SEND + nbuf_lower; s++){
312  q[s] = FARGO_ARRAY_INDEX(lower_buf[nv], s-SEND-1, k, j, i);
313  }
314  FARGO_Flux (q-m, flux-m, eps);
315  #endif
316  }else{ /* -- impose periodic b.c. -- */
317  for (s = 1; s <= ngh; s++) {
318  q[SBEG - s] = q[SBEG - s + NS];
319  q[SEND + s] = q[SEND + s - NS];
320  }
321  FARGO_Flux (q, flux, eps);
322  }
323 
324  /* -- update main solution array -- */
325 
326  if (nproc_s > 1){
327  #ifdef PARALLEL
328  for (s = SBEG; s <= SEND; s++){
329  sm = s - m;
330  U[k][j][i][nv] = q[sm] - eps*(flux[sm] - flux[sm-1]);
331  }
332  #endif
333  }else{
334  for (s = SBEG; s <= SEND; s++){
335  sm = FARGO_MOD(s - m);
336  U[k][j][i][nv] = q[sm] - eps*(flux[sm] - flux[sm-1]);
337  }
338  }
339  }
340  } /* -- end main spatial loop -- */
341 
342 /* -------------------------------------------------
343  6. Shift cell-centered quantities
344  ------------------------------------------------- */
345 
346 #ifdef STAGGERED_MHD
347 {
348  int ss;
349  double *Ar, *Ath, *dr, *r;
350  double ***bx1, ***bx2, ***bx3;
351 
352 /* -- pointer shortcuts -- */
353 
354  D_EXPAND(bx1 = Us[BX1s]; ,
355  bx2 = Us[BX2s]; ,
356  bx3 = Us[BX3s];)
357  r = x;
358  dr = dx;
359  Ar = grid[IDIR].A;
360  Ath = grid[JDIR].A;
361 
362 /* ------------------------------------------------------------------
363  6a. Compute Ez at x(i+1/2), y(j+1/2), z(k) (Cartesian or Polar)
364  Compute -Eth at x(i+1/2), y(j), z(k+1/2) (Spherical)
365  ------------------------------------------------------------------ */
366 
367  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
368  KDOM_LOOP(k) for (i = IBEG-1; i <= IEND; i++){
369  #else
370  JDOM_LOOP(j) for (i = IBEG-1; i <= IEND; i++){
371  #endif
372 
373  #if GEOMETRY == CARTESIAN
374  w = 0.5*(wA[k][i] + wA[k][i+1]);
375  #elif GEOMETRY == POLAR
376  w = 0.5*(wA[k][i] + wA[k][i+1])/xp[i];
377  #elif GEOMETRY == SPHERICAL
378  w = 0.5*(wA[j][i] + wA[j][i+1])/(xp[i]*sin(y[j]));
379  #endif
380 
381  /* --------------------------------------------------
382  6b. Obtain shift in terms of an integer number of
383  cells (m) plus a remainder eps.
384  -------------------------------------------------- */
385 
386  dL = w*g_dt; /* spatial shift */
387  dL = fmod(dL, Lphi);
388  m = (int)floor(dL/dphi + 0.5); /* nearest integer part */
389  eps = dL/dphi - m; /* remainder in units of dphi */
390 
391  /* -- copy 3D array into 1D buffer -- */
392 
393  SDOM_LOOP(s) q[s] = bx1[k][j][i];
394  if (nproc_s > 1){ /* -- copy values from lower and upper buffers -- */
395  #ifdef PARALLEL
396  for (s = SBEG-1; s >= SBEG - nbuf_upper; s--){
397  q[s] = FARGO_ARRAY_INDEX(upper_buf[BX1], SBEG-1-s, k, j, i);
398  }
399  for (s = SEND+1; s <= SEND + nbuf_lower; s++){
400  q[s] = FARGO_ARRAY_INDEX(lower_buf[BX1], s-SEND-1, k, j, i);
401  }
402  FARGO_Flux (q-m, flux-m, eps);
403  #endif
404  }else{ /* -- assign periodic b.c. -- */
405  for (s = 1; s <= ngh; s++) {
406  q[SBEG - s] = q[SBEG - s + NS];
407  q[SEND + s] = q[SEND + s - NS];
408  }
409  FARGO_Flux (q, flux, eps);
410  }
411 
412  /* -- update main solution array -- */
413 
414  if (nproc_s > 1){
415  #ifdef PARALLEL
416  if (m >= 0){
417  for (s = SBEG - 1; s <= SEND; s++) {
418  sm = s - m;
419  Ez[k][j][i] = eps*flux[sm];
420  for (ss = s; ss > (s-m); ss--) Ez[k][j][i] += q[ss];
421  Ez[k][j][i] *= dphi;
422  }
423  }else{
424  for (s = SBEG-1; s <= SEND; s++) {
425  sm = s - m;
426  Ez[k][j][i] = eps*flux[sm];
427  for (ss = s+1; ss < (s-m+1); ss++) Ez[k][j][i] -= q[ss];
428  Ez[k][j][i] *= dphi;
429  }
430  }
431  #endif
432  }else{
433  if (m >= 0){
434  for (s = SBEG - 1; s <= SEND; s++) {
435  sm = FARGO_MOD(s - m);
436  Ez[k][j][i] = eps*flux[sm];
437  for (ss = s; ss > (s-m); ss--) Ez[k][j][i] += q[FARGO_MOD(ss)];
438  Ez[k][j][i] *= dphi;
439  }
440  }else {
441  for (s = SBEG-1; s <= SEND; s++) {
442  sm = FARGO_MOD(s - m);
443  Ez[k][j][i] = eps*flux[sm];
444  for (ss = s+1; ss < (s-m+1); ss++) Ez[k][j][i] -= q[FARGO_MOD(ss)];
445  Ez[k][j][i] *= dphi;
446  }
447  }
448  }
449  }
450 
451 #if DIMENSIONS == 3
452 
453 /* -----------------------------------------------------------------
454  6c. Compute Ex at x(i), y(j+1/2), z(k+1/2) (Cartesian or Polar)
455  Compute -Er at x(i), y(j+1/2), z(k+1/2) (Spherical)
456  ----------------------------------------------------------------- */
457 
458  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
459  for (k = KBEG-1; k <= KEND; k++) IDOM_LOOP(i){
460  int BS = BX3;
461  double ***bs = Us[BX3s];
462  #else
463  for (j = JBEG-1; j <= JEND; j++) IDOM_LOOP(i){
464  int BS = BX2;
465  double ***bs = Us[BX2s];
466  #endif
467 
468  #if GEOMETRY == CARTESIAN
469  w = 0.5*(wA[k][i] + wA[k+1][i]);
470  #elif GEOMETRY == POLAR
471  w = 0.5*(wA[k][i] + wA[k+1][i])/x[i];;
472  #elif GEOMETRY == SPHERICAL
473  w = 0.5*(wA[j][i] + wA[j+1][i])/(x[i]*sin(yp[j]));
474  #endif
475 
476  /* --------------------------------------------------
477  6d. Obtain shift in terms of an integer number of
478  cells (m) plus a remainder eps.
479  -------------------------------------------------- */
480 
481  dL = w*g_dt; /* spatial shift */
482  dL = fmod(dL, Lphi);
483  m = (int)floor(dL/dphi + 0.5); /* nearest integer part */
484  eps = dL/dphi - m; /* remainder in units of dphi */
485 
486  /* -- copy 3D array into 1D buffer -- */
487 
488  SDOM_LOOP(s) q[s] = -bs[k][j][i];
489  if (nproc_s > 1){ /* -- copy values from lower and upper buffers -- */
490  #ifdef PARALLEL
491  for (s = SBEG-1; s >= SBEG - nbuf_upper; s--){
492  q[s] = -FARGO_ARRAY_INDEX(upper_buf[BS], SBEG-1-s, k, j, i);
493  }
494  for (s = SEND+1; s <= SEND + nbuf_lower; s++){
495  q[s] = -FARGO_ARRAY_INDEX(lower_buf[BS], s-SEND-1, k, j, i);
496  }
497  FARGO_Flux (q-m, flux-m, eps);
498  #endif
499  }else{ /* -- assign periodic b.c. -- */
500  for (s = 1; s <= ngh; s++) {
501  q[SBEG - s] = q[SBEG - s + NS];
502  q[SEND + s] = q[SEND + s - NS];
503  }
504  FARGO_Flux (q, flux, eps);
505  }
506 
507  /* -- update main solution array -- */
508 
509  if (nproc_s > 1){
510  #ifdef PARALLEL
511  if (m >= 0){
512  for (s = SBEG-1; s <= SEND; s++) {
513  sm = s - m;
514  Ex[k][j][i] = eps*flux[sm];
515  for (ss = s; ss > (s-m); ss--) Ex[k][j][i] += q[ss];
516  Ex[k][j][i] *= dphi;
517  }
518  }else {
519  for (s = SBEG-1; j <= SEND; s++) {
520  sm = s - m;
521  Ex[k][j][i] = eps*flux[sm];
522  for (ss = s+1; ss < (s-m+1); ss++) Ex[k][j][i] -= q[ss];
523  Ex[k][j][i] *= dphi;
524  }
525  }
526  #endif
527  }else{
528  if (m >= 0){
529  for (s = SBEG-1; s <= SEND; s++) {
530  sm = FARGO_MOD(s - m);
531  Ex[k][j][i] = eps*flux[sm];
532  for (ss = s; ss > (s-m); ss--) Ex[k][j][i] += q[FARGO_MOD(ss)];
533  Ex[k][j][i] *= dphi;
534  }
535  }else {
536  for (s = SBEG-1; j <= SEND; s++) {
537  sm = FARGO_MOD(s - m);
538  Ex[k][j][i] = eps*flux[sm];
539  for (ss = s+1; ss < (s-m+1); ss++) Ex[k][j][i] -= q[FARGO_MOD(ss)];
540  Ex[k][j][i] *= dphi;
541  }
542  }
543  }
544  }
545 #endif /* DIMENSIONS == 3 */
546 
547 /* ------------------------------------------
548  6e. Update bx1 staggered magnetic field
549  ------------------------------------------ */
550 
551  KDOM_LOOP(k) JDOM_LOOP(j) for (i = IBEG-1; i <= IEND; i++){
552  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
553  bx1[k][j][i] -= (Ez[k][j][i] - Ez[k][j-1][i])/dphi;
554  #elif GEOMETRY == SPHERICAL /* in spherical coordinates Ez = -Eth */
555  bx1[k][j][i] -= (Ez[k][j][i] - Ez[k-1][j][i])/dphi;
556  #endif
557  }
558 
559 /* ------------------------------------------
560  6f. Update bx2 staggered magnetic field
561  ------------------------------------------ */
562 
563  KDOM_LOOP(k) for (j = JBEG-1; j <= JEND; j++) IDOM_LOOP(i){
564  #if GEOMETRY == CARTESIAN
565  bx2[k][j][i] += D_EXPAND( ,
566  (Ez[k][j][i] - Ez[k][j][i-1])/dx[i] ,
567  - (Ex[k][j][i] - Ex[k-1][j][i])/dz[k]);
568  #elif GEOMETRY == POLAR
569  bx2[k][j][i] += D_EXPAND( ,
570  (Ar[i]*Ez[k][j][i] - Ar[i-1]*Ez[k][j][i-1])/dr[i] ,
571  - x[i]*(Ex[k][j][i] - Ex[k-1][j][i])/dz[k]);
572 
573  #elif GEOMETRY == SPHERICAL /* in spherical coordinates Ex = -Er */
574  bx2[k][j][i] += (Ex[k][j][i] - Ex[k-1][j][i])/dphi;
575  #endif
576  }
577 
578 /* ------------------------------------------
579  6g. Update bx3 staggered magnetic field
580  ------------------------------------------ */
581 
582  #if DIMENSIONS == 3
583  for (k = KBEG-1; k <= KEND; k++) JDOM_LOOP(j) IDOM_LOOP(i){
584  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
585  bx3[k][j][i] += (Ex[k][j][i] - Ex[k][j-1][i])/dphi;
586  #elif GEOMETRY == SPHERICAL
587  bx3[k][j][i] += sin(y[j])*( Ar[i]*Ez[k][j][i]
588  - Ar[i-1]*Ez[k][j][i-1])/(r[i]*dr[i])
589  - (Ath[j]*Ex[k][j][i] - Ath[j-1]*Ex[k][j-1][i])/dy[j];
590  #endif
591  }
592  #endif
593 
594 /* -- redefine cell-centered field -- */
595 
596  DOM_LOOP(k,j,i){
597  D_EXPAND(
598  U[k][j][i][BX1] = 0.5*(Us[BX1s][k][j][i] + Us[BX1s][k][j][i-1]); ,
599  U[k][j][i][BX2] = 0.5*(Us[BX2s][k][j][i] + Us[BX2s][k][j-1][i]); ,
600  U[k][j][i][BX3] = 0.5*(Us[BX3s][k][j][i] + Us[BX3s][k-1][j][i]);
601  )
602  }
603 }
604 #endif /* STAGGERED_MHD */
605 
606  #ifdef PARALLEL
607  if (nproc_s > 1){
608  for (nv = 0; nv < NVAR; nv++){
609  FreeArrayBox(upper_buf[nv], ubox.kb, ubox.jb, ubox.ib);
610  FreeArrayBox(lower_buf[nv], lbox.kb, lbox.jb, lbox.ib);
611  }
612  }
613  #endif
614 }
FARGO_ARRAY_INDEX
#define FARGO_ARRAY_INDEX(A, s, k, j, i)
Definition: fargo.h:69
MAX
#define MAX(a, b)
Definition: macros.h:101
KDOM_LOOP
#define KDOM_LOOP(k)
Definition: macros.h:36
RBOX::jb
int jb
Lower corner index in the x2 direction.
Definition: structs.h:349
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
SDIR
#define SDIR
Definition: fargo.h:55
SDOM_LOOP
#define SDOM_LOOP(s)
Definition: fargo.h:66
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
MAX_BUF_SIZE
#define MAX_BUF_SIZE
Definition: fargo.c:33
GRID::xr
double * xr
Definition: structs.h:81
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
wA
static double ** wA
Definition: fargo_velocity.c:24
RBOX::kb
int kb
Lower corner index in the x3 direction.
Definition: structs.h:351
FreeArrayBox
void FreeArrayBox(double ***t, long nrl, long ncl, long ndl)
Definition: arrays.c:403
g_dt
double g_dt
The current integration time step.
Definition: globals.h:118
RBOX::dk
int dk
Directional increment (+1 or -1) when looping over the 3rd dimension of the box.
Definition: structs.h:357
ArrayBox
double *** ArrayBox(long int nrl, long int nrh, long int ncl, long int nch, long int ndl, long int ndh)
Definition: arrays.c:341
BX3s
#define BX3s
Definition: ct.h:29
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
RBOX::ib
int ib
Lower corner index in the x1 direction.
Definition: structs.h:347
JDOM_LOOP
#define JDOM_LOOP(j)
Definition: macros.h:35
MIN
#define MIN(a, b)
Definition: macros.h:104
RBOX::dj
int dj
Directional increment (+1 or -1) when looping over the 2nd dimension of the box.
Definition: structs.h:355
FARGO_Source
void FARGO_Source(Data_Arr, double, Grid *)
Definition: fargo_source.c:28
SEND
#define SEND
Definition: fargo.h:57
eps
#define eps
IDIR
#define IDIR
Definition: pluto.h:193
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
GRID::nghost
int nghost
Number of ghost zones.
Definition: structs.h:104
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
g_domBeg
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
GRID::x
double * x
Definition: structs.h:80
FARGO_Flux
static void FARGO_Flux(double *, double *, double)
Definition: fargo.c:770
FARGO_ParallelExchange
static void FARGO_ParallelExchange(Data_Arr U, Data_Arr Us, double ***rbl[], RBox *, double ***rbu[], RBox *, Grid *grid)
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
FARGO_GetVelocity
double ** FARGO_GetVelocity(void)
Definition: fargo_velocity.c:270
RBOX::ie
int ie
Upper corner index in the x1 direction.
Definition: structs.h:348
RBOX::ke
int ke
Upper corner index in the x3 direction.
Definition: structs.h:352
RBOX::je
int je
Upper corner index in the x2 direction.
Definition: structs.h:350
BX1
#define BX1
Definition: mod_defs.h:25
FARGO_MOD
#define FARGO_MOD(s)
Definition: fargo.c:32
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
SBEG
#define SBEG
Definition: fargo.h:56
BX2s
#define BX2s
Definition: ct.h:28
BX2
#define BX2
Definition: mod_defs.h:26
KEND
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
RBOX::di
int di
Directional increment (+1 or -1) when looping over the 1st dimension of the box.
Definition: structs.h:353
g_domEnd
double g_domEnd[3]
Upper limits of the computational domain.
Definition: globals.h:126
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
NS
#define NS
Definition: fargo.h:58
NVAR
#define NVAR
Definition: pluto.h:609
RBOX
Definition: structs.h:346
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
IDOM_LOOP
#define IDOM_LOOP(i)
Definition: macros.h:34
JEND
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
GRID::nproc
int nproc
number of processors for this grid.
Definition: structs.h:120
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

Here is the caller graph for this function: