PLUTO
fargo.h File Reference

FARGO-MHD module header file. More...

Go to the source code of this file.

Macros

#define FARGO_ORDER   3
 Set the order of interpolation during the linear transport step. More...
 
#define FARGO_NSTEP_AVERAGE   10 /* Default is 10 */
 
#define FARGO_AVERAGE_VELOCITY   YES /* Default is YES */
 
#define SDIR   KDIR
 
#define SBEG   KBEG
 
#define SEND   KEND
 
#define NS   NX3
 
#define NS_TOT   NX3_TOT
 
#define SDOM_LOOP(s)   for ((s) = SBEG; (s) <= SEND; (s)++)
 
#define FARGO_ARRAY_INDEX(A, s, k, j, i)   A[s][j][i]
 

Functions

void FARGO_AddVelocity (const Data *, Grid *)
 
void FARGO_CHECK (Data_Arr V, Data_Arr U)
 
void FARGO_ComputeVelocity (const Data *, Grid *)
 
double ** FARGO_GetVelocity (void)
 
int FARGO_HasTotalVelocity ()
 
void FARGO_SubtractVelocity (const Data *, Grid *)
 
void FARGO_ShiftSolution (Data_Arr, Data_Arr, Grid *)
 
double FARGO_SetVelocity (double, double)
 
void FARGO_Source (Data_Arr, double, Grid *)
 

Detailed Description

FARGO-MHD module header file.

Contains contains basic definitions and declarations used by the FARGO-MHD module.

Authors
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
Oct 6, 2014
Todo:

Definition in file fargo.h.

Macro Definition Documentation

#define FARGO_ARRAY_INDEX (   A,
  s,
  k,
  j,
  i 
)    A[s][j][i]

Definition at line 69 of file fargo.h.

#define FARGO_AVERAGE_VELOCITY   YES /* Default is YES */

Average background velocity is computed by averaging the orbital velocity (YES) or by prescribing the velocity analytically exactly in the user-supplied function FARGO_SetVelocity

Definition at line 33 of file fargo.h.

#define FARGO_NSTEP_AVERAGE   10 /* Default is 10 */

Set how often (in number of steps) the total azimuthal velocity should be averaged.

Definition at line 25 of file fargo.h.

#define FARGO_ORDER   3

Set the order of interpolation during the linear transport step.

Either 2 or 3. Default is 3.

Definition at line 17 of file fargo.h.

#define NS   NX3

Definition at line 58 of file fargo.h.

#define NS_TOT   NX3_TOT

Definition at line 59 of file fargo.h.

#define SBEG   KBEG

Definition at line 56 of file fargo.h.

#define SDIR   KDIR

Definition at line 55 of file fargo.h.

#define SDOM_LOOP (   s)    for ((s) = SBEG; (s) <= SEND; (s)++)

Definition at line 66 of file fargo.h.

#define SEND   KEND

Definition at line 57 of file fargo.h.

Function Documentation

void FARGO_AddVelocity ( const Data d,
Grid grid 
)

Add the mean backgroun contribution to the residual velocity in order to obtain the total velocity.

Parameters
[in,out]dpointer to PLUTO Data structure
[in]gridpointer to array of Grid structures
Returns
This function has no return value.

Definition at line 226 of file fargo_velocity.c.

236 {
237  int i,j,k;
238  double ***vphi;
239 
240  #if GEOMETRY == CARTESIAN
241  vphi = d->Vc[VX2];
242  #else
243  vphi = d->Vc[iVPHI];
244  #endif
245 
246  TOT_LOOP(k,j,i){
247  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
248  vphi[k][j][i] += wA[k][i];
249  #elif GEOMETRY == SPHERICAL
250  vphi[k][j][i] += wA[j][i];
251  #else
252  print1 ("! FARGO not supported in this geometry\n");
253  QUIT_PLUTO(1);
254  #endif
255  }
256 
258 }
#define VX2
Definition: mod_defs.h:29
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
#define YES
Definition: pluto.h:25
static double ** wA
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
#define iVPHI
Definition: mod_defs.h:67
static int vphi_is_total_velocity
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
int j
Definition: analysis.c:2
int k
Definition: analysis.c:2
int i
Definition: analysis.c:2
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

void FARGO_CHECK ( Data_Arr  V,
Data_Arr  U 
)
void FARGO_ComputeVelocity ( const Data d,
Grid grid 
)

Compute the mean orbital velocity as a 2D array by averaging for each pair (x,z)/(r,z)/(r,theta) (in Cartesian/polar/spherical) along the orbital coordinate y/phi/phi.

Definition at line 29 of file fargo_velocity.c.

41 {
42  int i,j,k;
43  static int first_call = 1;
44  double *x1, *x2, *x3, th, w;
45  double ***vphi;
46 
47  x1 = grid[IDIR].x;
48  x2 = grid[JDIR].x;
49  x3 = grid[KDIR].x;
50 
51  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
52  if (wA == NULL){
53  wA = ARRAY_2D(NX3_TOT, NX1_TOT, double);
54 /*
55  wAx1 = ARRAY_2D(NX3_TOT, NX1_TOT, double);
56  wAx3 = ARRAY_2D(NX3_TOT, NX1_TOT, double);
57 */
58  }
59  #else
60  if (wA == NULL){
61  wA = ARRAY_2D(NX2_TOT, NX1_TOT, double);
62 /*
63  wAx1 = ARRAY_2D(NX2_TOT, NX1_TOT, double);
64  wAx2 = ARRAY_2D(NX2_TOT, NX1_TOT, double);
65 */
66  }
67  #endif
68 
69  #if GEOMETRY == CARTESIAN
70  vphi = d->Vc[VX2];
71  #else
72  vphi = d->Vc[iVPHI];
73  #endif
74 
75  #if FARGO_AVERAGE_VELOCITY == YES
76 
77 /* ---------------------------------------------------------------------
78  average velocity along the transport direction every 10 steps
79  --------------------------------------------------------------------- */
80 
81  if (g_stepNumber%FARGO_NSTEP_AVERAGE == 0 || first_call){
82  #ifdef PARALLEL
83  double w_recv=0.0, w_sum=0.0;
84  int count, coords[3], src, dst;
85  MPI_Comm cartcomm;
86  MPI_Status status;
87  #endif
88 
89  /* -- fill ghost zones -- */
90 
91  Boundary(d, ALL_DIR, grid);
92 
93  /* ---- get ranks of the upper and lower processsors ---- */
94 
95  #ifdef PARALLEL
96  AL_Get_cart_comm(SZ, &cartcomm);
97  for (i = 0; i < DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
98  coords[SDIR] += 1;
99  MPI_Cart_rank(cartcomm, coords, &dst);
100 
101  for (i = 0; i < DIMENSIONS; i++) coords[i] = grid[i].rank_coord;
102  coords[SDIR] -= 1;
103  MPI_Cart_rank(cartcomm, coords, &src);
104  #endif
105 
106  /* ---- compute average orbital speed ---- */
107 
108  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
109  KTOT_LOOP(k) ITOT_LOOP(i) {
110  #else
111  JTOT_LOOP(j) ITOT_LOOP(i) {
112  #endif
113  w = 0.0;
114  SDOM_LOOP(s) w += vphi[k][j][i];
115  if (grid[SDIR].nproc > 1){
116  #ifdef PARALLEL
117  w_sum = w;
118  for (count=1; count < grid[SDIR].nproc; count++ ){
119  MPI_Sendrecv(&w, 1, MPI_DOUBLE, dst, 0, &w_recv, 1,
120  MPI_DOUBLE, src, 0, cartcomm, &status);
121  w = w_recv;
122  w_sum += w;
123  }
124  w = w_sum/(double)(grid[SDIR].np_int_glob);
125  #endif
126  }else{
127  w = w/(double)NS;
128  }
129  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
130  wA[k][i] = w;
131  #elif GEOMETRY == SPHERICAL
132  wA[j][i] = w;
133  #endif
134  }
135  }
136 
137  /* -- compute interface velocity -- */
138 /*
139  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
140 
141  KTOT_LOOP(k) for (i = 0; i < NX1_TOT-1; i++){
142  wAx1[k][i] = 0.5*(wA[k][i] + wA[k][i+1]);
143  }
144  for (k = 0; k < NX3_TOT-1; k++) ITOT_LOOP(i){
145  wAx3[k][i] = 0.5*(wA[k][i] + wA[k+1][i]);
146  }
147 
148  #elif GEOMETRY == SPHERICAL
149 
150  JTOT_LOOP(j) for (i = 0; i < NX1_TOT-1; i++){
151  wAx1[j][i] = 0.5*(wA[j][i] + wA[j][i+1]);
152  }
153  for (j = 0; j < NX2_TOT-1; j++) ITOT_LOOP(i){
154  wAx2[j][i] = 0.5*(wA[j][i] + wA[j+1][i]);
155  }
156  #endif
157 */
158 
159  #else
160 
161 /* ---------------------------------------------------------------------
162  define velocity analytically
163  --------------------------------------------------------------------- */
164 
165  if (first_call){
166  #if GEOMETRY == CARTESIAN
167  KTOT_LOOP(k) ITOT_LOOP(i) wA[k][i] = FARGO_SetVelocity(x1[i], x3[k]);
168  #elif GEOMETRY == POLAR
169  KTOT_LOOP(k) ITOT_LOOP(i) wA[k][i] = x1[i]*FARGO_SetVelocity(x1[i], x3[k]);
170  #elif GEOMETRY == SPHERICAL
171  JTOT_LOOP(j) ITOT_LOOP(i) {
172  th = x2[j];
173  wA[j][i] = x1[i]*sin(th)*FARGO_SetVelocity(x1[i], x2[j]);
174  }
175  #else
176  print1 ("! FARGO not supported in this geometry\n");
177  QUIT_PLUTO(1);
178  #endif
179  }
180  #endif /* FARGO_AVERAGE_VELOCITY */
181 
182  first_call = 0;
183 }
void Boundary(const Data *d, int idim, Grid *grid)
Definition: boundary.c:36
#define VX2
Definition: mod_defs.h:29
#define SDIR
Definition: fargo.h:55
#define SDOM_LOOP(s)
Definition: fargo.h:66
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
static double ** wA
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
#define KTOT_LOOP(k)
Definition: macros.h:40
#define iVPHI
Definition: mod_defs.h:67
#define KDIR
Definition: pluto.h:195
int AL_Get_cart_comm(int, MPI_Comm *)
Definition: al_sz_get.c:117
#define FARGO_NSTEP_AVERAGE
Definition: fargo.h:25
long int g_stepNumber
Gives the current integration step number.
Definition: globals.h:97
#define IDIR
Definition: pluto.h:193
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
int j
Definition: analysis.c:2
if(divB==NULL)
Definition: analysis.c:10
int k
Definition: analysis.c:2
double * x
Definition: structs.h:80
#define ITOT_LOOP(i)
Definition: macros.h:38
#define JTOT_LOOP(j)
Definition: macros.h:39
int i
Definition: analysis.c:2
#define s
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
Definition: al_hidden.h:38
#define JDIR
Definition: pluto.h:194
#define NS
Definition: fargo.h:58
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
double FARGO_SetVelocity(double, double)
Definition: init.c:502
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
int nproc
number of processors for this grid.
Definition: structs.h:120
#define DIMENSIONS
Definition: definitions_01.h:2
#define ALL_DIR
Definition: pluto.h:196

Here is the call graph for this function:

Here is the caller graph for this function:

double** FARGO_GetVelocity ( void  )

Return a pointer to the background orbital velocity wA.

Definition at line 270 of file fargo_velocity.c.

275 {
276  return wA;
277 /*
278  int i, j, k;
279 
280  if (where == CENTER){
281 
282  return wA;
283 
284  }else if (where == X1FACE){
285 
286  return wAx1;
287 
288  }else if (where == X2FACE){
289 
290  #if GEOMETRY == SPHERICAL
291  return wAx2;
292  #else
293  print1 ("! FARGO_GetVelocity: incorrect direction\n");
294  QUIT_PLUTO(1);
295  #endif
296 
297  }else if (where == X3FACE){
298 
299  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
300  return wAx3;
301  #else
302  print1 ("! FARGO_GetVelocity: incorrect direction\n");
303  QUIT_PLUTO(1);
304  #endif
305 
306  }
307 
308 */
309 }
static double ** wA

Here is the caller graph for this function:

int FARGO_HasTotalVelocity ( )

Return 1 if the velocity array contains the total velocity. Return 0 otherwise.

Definition at line 261 of file fargo_velocity.c.

266 {
267  return vphi_is_total_velocity;
268 }
static int vphi_is_total_velocity

Here is the caller graph for this function:

double FARGO_SetVelocity ( double  x1,
double  x2 
)

Compute the shear angular velocity to be subtracted from the HD or MHD equations.

Definition at line 502 of file init.c.

508 {
509  return -SB_Q*SB_OMEGA*x1;
510 }
#define SB_OMEGA
Disk local orbital frequency .
Definition: shearingbox.h:81
#define SB_Q
The shear parameter, .
Definition: shearingbox.h:76

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 }
#define FARGO_ARRAY_INDEX(A, s, k, j, i)
Definition: fargo.h:69
#define MAX(a, b)
Definition: macros.h:101
#define KDOM_LOOP(k)
Definition: macros.h:36
int jb
Lower corner index in the x2 direction.
Definition: structs.h:349
DOM_LOOP(k, j, i)
Definition: analysis.c:22
#define SDIR
Definition: fargo.h:55
#define SDOM_LOOP(s)
Definition: fargo.h:66
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
#define MAX_BUF_SIZE
Definition: fargo.c:33
double * xr
Definition: structs.h:81
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
static double ** wA
int kb
Lower corner index in the x3 direction.
Definition: structs.h:351
void FreeArrayBox(double ***t, long nrl, long ncl, long ndl)
Definition: arrays.c:403
double g_dt
The current integration time step.
Definition: globals.h:118
int dk
Directional increment (+1 or -1) when looping over the 3rd dimension of the box.
Definition: structs.h:357
double *** ArrayBox(long int nrl, long int nrh, long int ncl, long int nch, long int ndl, long int ndh)
Definition: arrays.c:341
#define BX3s
Definition: ct.h:29
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
double * dx
Definition: structs.h:83
#define KDIR
Definition: pluto.h:195
int ib
Lower corner index in the x1 direction.
Definition: structs.h:347
#define JDOM_LOOP(j)
Definition: macros.h:35
#define MIN(a, b)
Definition: macros.h:104
int dj
Directional increment (+1 or -1) when looping over the 2nd dimension of the box.
Definition: structs.h:355
void FARGO_Source(Data_Arr, double, Grid *)
Definition: fargo_source.c:28
#define SEND
Definition: fargo.h:57
#define eps
#define IDIR
Definition: pluto.h:193
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
int j
Definition: analysis.c:2
if(divB==NULL)
Definition: analysis.c:10
int nghost
Number of ghost zones.
Definition: structs.h:104
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
int k
Definition: analysis.c:2
#define BX1s
Definition: ct.h:27
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
double * x
Definition: structs.h:80
static void FARGO_Flux(double *, double *, double)
Definition: fargo.c:770
static void FARGO_ParallelExchange(Data_Arr U, Data_Arr Us, double ***rbl[], RBox *, double ***rbu[], RBox *, Grid *grid)
#define s
#define BX3
Definition: mod_defs.h:27
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
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
int i
Definition: analysis.c:2
double ** FARGO_GetVelocity(void)
int ie
Upper corner index in the x1 direction.
Definition: structs.h:348
int ke
Upper corner index in the x3 direction.
Definition: structs.h:352
int je
Upper corner index in the x2 direction.
Definition: structs.h:350
#define BX1
Definition: mod_defs.h:25
#define FARGO_MOD(s)
Definition: fargo.c:32
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
#define SBEG
Definition: fargo.h:56
#define BX2s
Definition: ct.h:28
#define BX2
Definition: mod_defs.h:26
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
int di
Directional increment (+1 or -1) when looping over the 1st dimension of the box.
Definition: structs.h:353
double g_domEnd[3]
Upper limits of the computational domain.
Definition: globals.h:126
#define JDIR
Definition: pluto.h:194
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
#define NS
Definition: fargo.h:58
#define NVAR
Definition: pluto.h:609
Definition: structs.h:346
static Runtime q
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
#define IDOM_LOOP(i)
Definition: macros.h:34
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
int nproc
number of processors for this grid.
Definition: structs.h:120
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
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:

void FARGO_Source ( Data_Arr  UU,
double  dt,
Grid grid 
)
Parameters
[in,out]UUan array of conserved variables
[in]dtthe current time increment
[in]gridpointer to an array of Grid structures

Definition at line 28 of file fargo_source.c.

35 {
36 #if (HAVE_ENERGY) && (defined SHEARINGBOX)
37  int i,j,k;
38  double scrh, rho, mx, my, Bx, By;
39 
40  scrh = SB_Q*SB_OMEGA;
41  DOM_LOOP(k,j,i){
42  rho = UU[k][j][i][RHO];
43  mx = UU[k][j][i][MX1];
44  my = UU[k][j][i][MX2];
45  Bx = UU[k][j][i][BX1];
46  By = UU[k][j][i][BX2];
47 
48  UU[k][j][i][ENG] += - dt*scrh*Bx*(By - 0.5*dt*Bx*scrh)
49  + dt*scrh*mx*my/rho;
50  }
51 #endif
52 }
#define MX1
Definition: mod_defs.h:20
DOM_LOOP(k, j, i)
Definition: analysis.c:22
#define RHO
Definition: mod_defs.h:19
tuple scrh
Definition: configure.py:200
static double Bx
Definition: hlld.c:62
#define SB_OMEGA
Disk local orbital frequency .
Definition: shearingbox.h:81
int j
Definition: analysis.c:2
#define MX2
Definition: mod_defs.h:21
int k
Definition: analysis.c:2
#define SB_Q
The shear parameter, .
Definition: shearingbox.h:76
int i
Definition: analysis.c:2
#define BX1
Definition: mod_defs.h:25
#define BX2
Definition: mod_defs.h:26

Here is the call graph for this function:

Here is the caller graph for this function:

void FARGO_SubtractVelocity ( const Data d,
Grid grid 
)

Subtract the mean background contribution from the total velocity. In other words, compute the residual velocity.

Parameters
[in,out]dpointer to PLUTO Data structure
[in]gridpointer to array of Grid structures
Returns
This function has no return value.

Definition at line 188 of file fargo_velocity.c.

198 {
199  int i,j,k;
200  double ***vphi;
201 
202  #if GEOMETRY == CARTESIAN
203  vphi = d->Vc[VX2];
204  #else
205  vphi = d->Vc[iVPHI];
206  #endif
207 
208 /* ----------------------------------------------------------------------
209  subtract velocity
210  ---------------------------------------------------------------------- */
211 
212  TOT_LOOP(k,j,i){
213  #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
214  vphi[k][j][i] -= wA[k][i];
215  #elif GEOMETRY == SPHERICAL
216  vphi[k][j][i] -= wA[j][i];
217  #else
218  print1 ("! FARGO not supported in this geometry\n");
219  QUIT_PLUTO(1);
220  #endif
221  }
222 
224 }
#define VX2
Definition: mod_defs.h:29
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
static double ** wA
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
#define iVPHI
Definition: mod_defs.h:67
static int vphi_is_total_velocity
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
int j
Definition: analysis.c:2
int k
Definition: analysis.c:2
int i
Definition: analysis.c:2
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
#define NO
Definition: pluto.h:26

Here is the call graph for this function:

Here is the caller graph for this function: