18 static void stretch_fun (
double,
double *,
double *,
double *);
32 double *
dx, *xlft, *xrgt;
35 double xpatch_lft, xpatch_rgt;
44 for (idim = 0; idim < 3; idim++) {
50 MakeGrid (idim, rtime, xlft, xrgt, dx);
55 for (i = 1; i <= rtime->
npoint[idim]; i++) {
62 iR = rtime->
npoint[idim] + ngh - 1;
69 for (i = 0; i < ngh; i++) {
86 for (i = 0; i <= iR + ngh; i++) {
109 sprintf (fname,
"%s/grid.out", rtime->
output_dir);
111 fg = fopen(fname,
"w");
112 for (idim = 0; idim < 3; idim++) {
118 fprintf (fg,
"%d \n", iR - iL + 1);
119 for (i = iL; i <= iR; i++) {
120 fprintf (fg,
" %d %12.6e %12.6e %12.6e %12.6e\n",
130 fg = fopen(fname,
"w");
131 fprintf (fg,
"# ******************************************************\n");
133 fprintf (fg,
"# PLUTO-Chombo %s (Base) Grid File\n",
PLUTO_VERSION);
138 fprintf (fg,
"# Generated on %s",asctime(localtime(&time_now)));
142 fprintf (fg,
"# \n");
144 #if GEOMETRY == CARTESIAN
145 fprintf (fg,
"# GEOMETRY: CARTESIAN\n");
146 #elif GEOMETRY == CYLINDRICAL
147 fprintf (fg,
"# GEOMETRY: CYLINDRICAL\n");
148 #elif GEOMETRY == POLAR
149 fprintf (fg,
"# GEOMETRY: POLAR\n");
150 #elif GEOMETRY == SPHERICAL
151 fprintf (fg,
"# GEOMETRY: SPHERICAL\n");
154 fprintf (fg,
"# X%d: [% f, % f], %d point(s), %d ghosts\n", idim+1,
158 fprintf (fg,
"# ******************************************************\n");
159 for (idim = 0; idim < 3; idim++) {
165 fprintf (fg,
"%d \n", iR - iL + 1);
166 for (i = iL; i <= iR; i++) {
167 fprintf (fg,
" %d %18.12e %18.12e\n",
184 print1 (
" X%d: [% f, % f], %d point(s), %d ghosts\n", idim+1,
199 for (dir = 0; dir < 3; dir++){
227 for (idim = 0; idim < 3; idim++) {
316 int i, iseg,
n, nstart, nseg, npoint;
317 int log_inc, log_dec, next_seg_is;
318 int done_with_segment[16], all_segments_done;
320 int i_patch_lft[16], i_patch_rgt[16];
321 double x_patch_lft[16], x_patch_rgt[16];
324 double dalpha,
alpha, f, df;
327 nseg = rtime->
npatch[idim];
339 for (iseg = 2; iseg <= nseg; iseg++) {
340 i_patch_lft[iseg] = i_patch_rgt[iseg - 1] + 1;
341 i_patch_rgt[iseg] = i_patch_lft[iseg] + rtime->
patch_npoint[idim][iseg] - 1;
346 done_with_segment[0] = 0;
347 done_with_segment[nseg + 1] = 0;
349 for (iseg = 1; iseg <= nseg; iseg++) {
351 done_with_segment[iseg] = 0;
352 xL = x_patch_lft[iseg];
353 xR = x_patch_rgt[iseg];
354 iL = i_patch_lft[iseg];
355 iR = i_patch_rgt[iseg];
363 for (i = iL; i <= iR; i++) {
364 dx[
i] = (xR - xL)/(
double)(npoint);
365 xlft[
i] = xL + (double)(i - iL)*dx[
i];
366 xrgt[
i] = xlft[
i] + dx[
i];
368 done_with_segment[iseg] = 1;
376 dy = log10( (xR + fabs(xL) - xL)/fabs(xL));
377 dy /= (double) (npoint);
379 if (log_dec) dy *= -1.0;
383 for (i = iL; i <= iR; i++) {
385 dx[
i] = (xlft[
i] + fabs(xL) - xL)*(pow(10.0,dy) - 1.0);
387 dx[
i] = (xlft[
i] - fabs(xL) - xR)*(pow(10.0,dy) - 1.0);
389 xrgt[
i] = xlft[
i] + dx[
i];
390 xlft[i + 1] = xrgt[
i];
392 done_with_segment[iseg] = 1;
399 all_segments_done = 1;
400 for (iseg = 1; iseg <= nseg; iseg++) {
401 all_segments_done = all_segments_done && done_with_segment[iseg];
408 while (!all_segments_done) {
416 for (iseg = 1; iseg <= nseg; iseg++) {
418 xL = x_patch_lft[iseg];
419 xR = x_patch_rgt[iseg];
420 iL = i_patch_lft[iseg];
421 iR = i_patch_rgt[iseg];
430 if (done_with_segment[iseg + 1]) next_seg_is = iseg + 1;
431 if (done_with_segment[iseg - 1]) next_seg_is = iseg - 1;
440 next_seg_is != 0 && !done_with_segment[iseg]) {
451 nstart = (next_seg_is > iseg ? iR:iL);
461 par[0] = (double)npoint;
463 if (next_seg_is > iseg)
464 par[1] = (xR - xL)/dx[nstart + 1];
466 par[1] = (xR - xL)/dx[nstart - 1];
472 print1 (
"Too many iterations during grid (%d) generation!\n",idim);
481 if (fabs (dalpha) < 1.e-14*alpha)
486 print1 (
" ! WARNING: alpha=%12.6e > 1.05 , dimension %d\n", alpha, idim);
487 print1 (
" ! While stretching segment %d\n", iseg);
491 print1 (
" - Stretched grid on dim %d (ratio: %f)\n",idim, alpha);
493 if (next_seg_is > iseg) {
495 for (i = nstart; i >= iL; i--) {
496 dx[
i] = pow (alpha, nstart - i + 1)*dx[nstart + 1];
497 xlft[
i] = xrgt[
i] - dx[
i];
498 xrgt[i - 1] = xlft[
i];
502 for (i = nstart; i <= iR; i++) {
503 dx[
i] = pow (alpha, i - nstart + 1)*dx[nstart - 1];
504 xlft[i + 1] = xrgt[
i] = xlft[
i] + dx[
i];
510 done_with_segment[iseg] = 1;
514 all_segments_done = 1;
515 for (i = 1; i <= nseg; i++)
516 all_segments_done = all_segments_done && done_with_segment[i];
535 double xN, L_dx,
scrh;
540 scrh = x*(pow (x, xN) - 1.0)/(x - 1.0);
541 *f = log (scrh) - log (L_dx);
542 *dfdx = 1.0/x + xN*pow (x, xN - 1.0)/(pow (x, xN) - 1.0) -
void FreeArray1D(void *v)
void MakeGeometry(Grid *)
int end
Global end index for the local array.
int npoint[3]
Global number of zones in the interior domain.
void print1(const char *fmt,...)
static void InitializeGrid(Runtime *, Grid *)
char output_dir[256]
The name of the output directory (output_dir for static PLUTO, Output_dir for PLUTO-Chombo) ...
int np_int_glob
Total number of points in the global domain (boundaries excluded).
#define LOGARITHMIC_INC_GRID
double * x_glob
Cell geometrical central points.
double * dx_glob
Cell size.
int lend
Local end index for the local array.
int lbeg
Local start index for the local array.
double * xr_glob
Cell right interface.
double xf
Leftmost and rightmost point in the local domain.
int gend
Global end index for the global array.
int np_tot_glob
Total number of points in the global domain (boundaries included).
#define LOGARITHMIC_DEC_GRID
double patch_left_node[5][16]
int beg
Global start index for the local array.
static void MakeGrid(int, Runtime *, double *, double *, double *)
int gbeg
Global start index for the global array.
void SetGrid(Runtime *rtime, Grid *GXYZ)
int nghost
Number of ghost zones.
double * xl_glob
Cell left interface.
double g_domBeg[3]
Lower limits of the computational domain.
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
double * r_1
Geometrical factor 1/r.
#define ARRAY_1D(nx, type)
double * ct
Geometrical factor cot(theta).
int npatch[5]
The number of grid patches.
int np_tot
Total number of points in the local domain (boundaries included).
static void stretch_fun(double, double *, double *, double *)
double g_domEnd[3]
Upper limits of the computational domain.
void FreeGrid(Grid *grid)
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
#define QUIT_PLUTO(e_code)
int nproc
number of processors for this grid.
double * A
Right interface area, A[i] = .
int np_int
Total number of points in the local domain (boundaries excluded).