diff options
| -rw-r--r-- | src/CartGrid3D.c | 562 |
1 files changed, 307 insertions, 255 deletions
diff --git a/src/CartGrid3D.c b/src/CartGrid3D.c index 384e41a..5a5a2d0 100644 --- a/src/CartGrid3D.c +++ b/src/CartGrid3D.c @@ -41,272 +41,292 @@ void CartGrid3D(CCTK_ARGUMENTS) CCTK_REAL lowerx,upperx,lowery,uppery,lowerz,upperz; int lowerxi,upperxi,loweryi,upperyi,lowerzi,upperzi; char infoline[120]; + int ierr; int domainsym[6]; int cntstag[3]; - - /* Some compilers warn about variables which might be uninitialized - because they are assigned in some if/else statements only */ - x_origin = y_origin = z_origin = 0.0; - lowerxi = loweryi = lowerzi = 0; - upperxi = cctk_gsh[0]-1; - upperyi = cctk_gsh[1]-1; - upperzi = cctk_gsh[2]-1; - this_dx = this_dy = this_dz = 0.0; - - /* Avoid origin? Default is yes */ - cntstag[0] = no_origin && no_originx && avoid_origin && avoid_originx; - cntstag[1] = no_origin && no_originy && avoid_origin && avoid_originy; - cntstag[2] = no_origin && no_originz && avoid_origin && avoid_originz; - - /* Determine symmetries of domain */ - DecodeSymParameters3D(domainsym); - - /* Calculate physical indices, using symmetries and periodicity */ - if (domainsym[0]) - { - lowerxi = cctk_nghostzones[0]; - } - else + if (cctk_levfac[0] == 1 && cctk_levfac[1] == 1 && cctk_levfac[2] == 1) { - lowerxi = 0; - } - if (domainsym[2]) - { - loweryi = cctk_nghostzones[1]; - } - else - { - loweryi = 0; - } - if (domainsym[4]) - { - lowerzi = cctk_nghostzones[2]; - } - else - { - lowerzi = 0; - } - if (periodic) - { - if (periodic_x) + /* Calculate the coordinate ranges only for the coarsest level */ + + /* Some compilers warn about variables which might be uninitialized + because they are assigned in some if/else statements only */ + x_origin = y_origin = z_origin = 0.0; + lowerxi = loweryi = lowerzi = 0; + upperxi = cctk_gsh[0]-1; + upperyi = cctk_gsh[1]-1; + upperzi = cctk_gsh[2]-1; + this_dx = this_dy = this_dz = 0.0; + + /* Avoid origin? Default is yes */ + cntstag[0] = no_origin && no_originx && avoid_origin && avoid_originx; + cntstag[1] = no_origin && no_originy && avoid_origin && avoid_originy; + cntstag[2] = no_origin && no_originz && avoid_origin && avoid_originz; + + /* Determine symmetries of domain */ + DecodeSymParameters3D(domainsym); + + /* Calculate physical indices, using symmetries and periodicity */ + if (domainsym[0]) { - lowerxi = cctk_nghostzones[0]-1; - upperxi = cctk_gsh[0]-1-cctk_nghostzones[0]; - } - if (periodic_y) + lowerxi = cctk_nghostzones[0]; + } + else { - loweryi = cctk_nghostzones[1]-1; - upperyi = cctk_gsh[1]-1-cctk_nghostzones[1]; + lowerxi = 0; } - if (periodic_z) + if (domainsym[2]) { - lowerzi = cctk_nghostzones[2]-1; - upperzi = cctk_gsh[2]-1-cctk_nghostzones[2]; - } - } - - dconv = pow(2, cctk_convlevel); - iconv = (int)dconv; - - /**************************************************************** - * - * BYRANGE - * - * User gives: minimum and maximum values of coordinates and - * the number of gridpoints on the coarse grid - * - ***************************************************************/ - /************************************************************** - * - * BOX (-0.5 to 0.5) - * - * User gives: number of gridpoints on the coarse grid - * - **************************************************************/ - - if (CCTK_Equals(type,"byrange") || CCTK_Equals(type,"box")) - { - - if (CCTK_Equals(type,"box")) + loweryi = cctk_nghostzones[1]; + } + else { - - /* Coordinates are all -0.5 to 0.5 */ - xmax1 = 0.5; - ymax1 = 0.5; - zmax1 = 0.5; - - xmin1 = -0.5; - ymin1 = -0.5; - zmin1 = -0.5; - + loweryi = 0; } + if (domainsym[4]) + { + lowerzi = cctk_nghostzones[2]; + } else { + lowerzi = 0; + } - if (xyzmax != -424242) + if (periodic) + { + if (periodic_x) { - xmax1 = xyzmax; - ymax1 = xyzmax; - zmax1 = xyzmax; + lowerxi = cctk_nghostzones[0]; + upperxi = cctk_gsh[0]-1-cctk_nghostzones[0]; } - else + if (periodic_y) { - xmax1 = xmax; - ymax1 = ymax; - zmax1 = zmax; + loweryi = cctk_nghostzones[1]; + upperyi = cctk_gsh[1]-1-cctk_nghostzones[1]; } - - if (xyzmin != -424242) + if (periodic_z) { - xmin1 = xyzmin; - ymin1 = xyzmin; - zmin1 = xyzmin; + lowerzi = cctk_nghostzones[2]; + upperzi = cctk_gsh[2]-1-cctk_nghostzones[2]; } - else + } + + dconv = pow(2, cctk_convlevel); + iconv = (int)dconv; + + /**************************************************************** + * + * BYRANGE + * + * User gives: minimum and maximum values of coordinates and + * the number of gridpoints on the coarse grid + * + ***************************************************************/ + /************************************************************** + * + * BOX (-0.5 to 0.5) + * + * User gives: number of gridpoints on the coarse grid + * + **************************************************************/ + + if (CCTK_Equals(type,"byrange") || CCTK_Equals(type,"box")) + { + + if (CCTK_Equals(type,"box")) { - xmin1 = xmin; - ymin1 = ymin; - zmin1 = zmin; + + /* Coordinates are all -0.5 to 0.5 */ + xmax1 = 0.5; + ymax1 = 0.5; + zmax1 = 0.5; + + xmin1 = -0.5; + ymin1 = -0.5; + zmin1 = -0.5; + } + else + { + + if (xyzmax != -424242) + { + xmax1 = xyzmax; + ymax1 = xyzmax; + zmax1 = xyzmax; + } + else + { + xmax1 = xmax; + ymax1 = ymax; + zmax1 = zmax; + } + + if (xyzmin != -424242) + { + xmin1 = xyzmin; + ymin1 = xyzmin; + zmin1 = xyzmin; + } + else + { + xmin1 = xmin; + ymin1 = ymin; + zmin1 = zmin; + } - } - + } - - /* Grid spacing on coarsest grid */ - /* TODO: Put the coordinates into arrays, and loop over the - dimensions. That gets ride of all the triplicated code. */ - if (domainsym[0]) - { - if (cntstag[0]) + /* Grid spacing on coarsest grid */ + /* TODO: Put the coordinates into arrays, and loop over the + dimensions. That gets ride of all the triplicated code. */ + if (domainsym[0]) { - *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 0.5); - x_origin = - (cctk_nghostzones[0] - 0.5) * *coarse_dx * iconv; + if (cntstag[0]) + { + *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 0.5); + x_origin = - (cctk_nghostzones[0] - 0.5) * *coarse_dx * iconv; + } + else + { + *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 1); + x_origin = - cctk_nghostzones[0] * *coarse_dx * iconv; + } } else { - *coarse_dx = xmax1 / (cctk_gsh[0] - cctk_nghostzones[0] - 1); - x_origin = - cctk_nghostzones[0] * *coarse_dx * iconv; + *coarse_dx = (xmax1 - xmin1) / max(cctk_gsh[0] - 1, 1); + x_origin = xmin1; } - } - else - { - *coarse_dx = (xmax1 - xmin1) / max(cctk_gsh[0] - 1, 1); - x_origin = xmin1; - } - if (domainsym[2]) - { - if (cntstag[1]) + if (domainsym[2]) { - *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 0.5); - y_origin = - (cctk_nghostzones[1] - 0.5) * *coarse_dy * iconv; + if (cntstag[1]) + { + *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 0.5); + y_origin = - (cctk_nghostzones[1] - 0.5) * *coarse_dy * iconv; + } + else + { + *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 1); + y_origin = - cctk_nghostzones[1] * *coarse_dy * iconv; + } } else { - *coarse_dy = ymax1 / (cctk_gsh[1] - cctk_nghostzones[1] - 1); - y_origin = - cctk_nghostzones[1] * *coarse_dy * iconv; + *coarse_dy = (ymax1 - ymin1) / max(cctk_gsh[1] - 1, 1); + y_origin = ymin1; } - } - else - { - *coarse_dy = (ymax1 - ymin1) / max(cctk_gsh[1] - 1, 1); - y_origin = ymin1; - } - if (domainsym[4]) - { - if (cntstag[2]) + if (domainsym[4]) { - *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 0.5); - z_origin = - (cctk_nghostzones[2] - 0.5) * *coarse_dz * iconv; + if (cntstag[2]) + { + *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 0.5); + z_origin = - (cctk_nghostzones[2] - 0.5) * *coarse_dz * iconv; + } + else + { + *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 1); + z_origin = - cctk_nghostzones[2] * *coarse_dz * iconv; + } } else { - *coarse_dz = zmax1 / (cctk_gsh[2] - cctk_nghostzones[2] - 1); - z_origin = - cctk_nghostzones[2] * *coarse_dz * iconv; + *coarse_dz = (zmax1 - zmin1) / max(cctk_gsh[2] - 1, 1); + z_origin = zmin1; } - } - else - { - *coarse_dz = (zmax1 - zmin1) / max(cctk_gsh[2] - 1, 1); - z_origin = zmin1; + + /* dx,dy,dz on the grid we are on */ + this_dx = *coarse_dx * iconv; + this_dy = *coarse_dy * iconv; + this_dz = *coarse_dz * iconv; + } - - /* dx,dy,dz on the grid we are on */ - this_dx = *coarse_dx * iconv; - this_dy = *coarse_dy * iconv; - this_dz = *coarse_dz * iconv; - } - - /************************************************************** - * BYSPACING - * - * User gives: grid spacing on the coarsest GH and - * the number of gridpoints on the coarsest GH - * - **************************************************************/ - - else if (CCTK_Equals(type,"byspacing")) - { + /************************************************************** + * BYSPACING + * + * User gives: grid spacing on the coarsest GH and + * the number of gridpoints on the coarsest GH + * + **************************************************************/ - /* Dx, Dy, Dx on the coarsest grid */ - - if (dxyz > 0) - { - *coarse_dx = dxyz; - *coarse_dy = dxyz; - *coarse_dz = dxyz; - } - else + else if (CCTK_Equals(type,"byspacing")) { - *coarse_dx = dx; - *coarse_dy = dy; - *coarse_dz = dz; - } + + /* Dx, Dy, Dx on the coarsest grid */ + + if (dxyz > 0) + { + *coarse_dx = dxyz; + *coarse_dy = dxyz; + *coarse_dz = dxyz; + } + else + { + *coarse_dx = dx; + *coarse_dy = dy; + *coarse_dz = dz; + } - /* dx, dy, dz on the grid we are on */ - this_dx = *coarse_dx * iconv; - this_dy = *coarse_dy * iconv; - this_dz = *coarse_dz * iconv; + /* dx, dy, dz on the grid we are on */ + this_dx = *coarse_dx * iconv; + this_dy = *coarse_dy * iconv; + this_dz = *coarse_dz * iconv; + /* Set minimum values of coordinates */ + if (domainsym[0]) + { + x_origin = (- cctk_nghostzones[0] + cntstag[0] * 0.5) * this_dx; + } + else + { + x_origin = - 0.5 * (cctk_gsh[0]-1 - cntstag[0] * (cctk_gsh[0])%2) + * this_dx; + } + + if (domainsym[2]) + { + y_origin = (- cctk_nghostzones[1] + cntstag[1] * 0.5) * this_dy; + } + else + { + y_origin = - 0.5 * (cctk_gsh[1]-1 - cntstag[1] * (cctk_gsh[1])%2) + * this_dy; + } + + if (domainsym[4]) + { + z_origin = (- cctk_nghostzones[2] + cntstag[2] * 0.5) * this_dz; + } + else + { + z_origin = - 0.5 * (cctk_gsh[2]-1 - cntstag[2] * (cctk_gsh[2])%2) + * this_dz; + } + } + } + else + { + /* if (not coarsest refinement level) */ + /* Use the already calculated coordinate ranges for all but the + coarsest levels */ - /* Set minimum values of coordinates */ - if (domainsym[0]) - { - x_origin = (- cctk_nghostzones[0] + cntstag[0] * 0.5) * this_dx; - } - else - { - x_origin = - 0.5 * (cctk_gsh[0]-1 - cntstag[0] * (cctk_gsh[0])%2) - * this_dx; - } + this_dx = cctk_delta_space[0] / cctk_levfac[0]; + this_dy = cctk_delta_space[1] / cctk_levfac[1]; + this_dz = cctk_delta_space[2] / cctk_levfac[2]; - if (domainsym[2]) - { - y_origin = (- cctk_nghostzones[1] + cntstag[1] * 0.5) * this_dy; - } - else - { - y_origin = - 0.5 * (cctk_gsh[1]-1 - cntstag[1] * (cctk_gsh[1])%2) - * this_dy; - } + x_origin = cctk_origin_space[0]; + y_origin = cctk_origin_space[1]; + z_origin = cctk_origin_space[2]; - if (domainsym[4]) - { - z_origin = (- cctk_nghostzones[2] + cntstag[2] * 0.5) * this_dz; - } - else - { - z_origin = - 0.5 * (cctk_gsh[2]-1 - cntstag[2] * (cctk_gsh[2])%2) - * this_dz; - } - } + } /* if (not coarsest refinement level) */ + + /* Set spatial coordinates */ @@ -326,51 +346,83 @@ void CartGrid3D(CCTK_ARGUMENTS) } } - cctk_delta_space[0] = this_dx; - cctk_delta_space[1] = this_dy; - cctk_delta_space[2] = this_dz; + if (cctk_levfac[0] == 1 && cctk_levfac[1] == 1 && cctk_levfac[2] == 1) + { + + /* Register the coordinate ranges only for the coarsest level */ + + cctk_delta_space[0] = this_dx; + cctk_delta_space[1] = this_dy; + cctk_delta_space[2] = this_dz; + + cctk_origin_space[0] = x_origin; + cctk_origin_space[1] = y_origin; + cctk_origin_space[2] = z_origin; + + lowerx = x_origin; + upperx = x_origin+this_dx*(cctk_gsh[0]-1); + ierr = CCTK_CoordRegisterRange(cctkGH,lowerx,upperx,-1,"x","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register x-coordinate computational range"); + } + + ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerxi,upperxi,-1,"x","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register x-coordinate physical range"); + } + + lowery = y_origin; + uppery = y_origin+this_dy*(cctk_gsh[1]-1); + ierr = CCTK_CoordRegisterRange(cctkGH,lowery,uppery,-1,"y","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register y-coordinate computational range"); + } + ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,loweryi,upperyi,-1,"y","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register y-coordinate physical range"); + } + + lowerz = z_origin; + upperz = z_origin+this_dz*(cctk_gsh[2]-1); + ierr = CCTK_CoordRegisterRange(cctkGH,lowerz,upperz,-1,"z","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register z-coordinate computational range"); + } + ierr = CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerzi,upperzi,-1,"z","cart3d"); + if (ierr < 0) + { + CCTK_WARN(0,"Failed to register z-coordinate physical range"); + } - cctk_origin_space[0] = x_origin; - cctk_origin_space[1] = y_origin; - cctk_origin_space[2] = z_origin; + CCTK_INFO("Grid Spacings:"); + sprintf(infoline," %s%12.7e %s%12.7e %s%12.7e ", + "dx=>",cctk_delta_space[0], + "dy=>",cctk_delta_space[1], + "dz=>",cctk_delta_space[2]); + CCTK_INFO(infoline); - lowerx = x_origin; - upperx = x_origin+this_dx*(cctk_gsh[0]-1); - CCTK_CoordRegisterRange(cctkGH,lowerx,upperx,-1,"x","cart3d"); - CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerxi,upperxi,-1,"x","cart3d"); - - lowery = y_origin; - uppery = y_origin+this_dy*(cctk_gsh[1]-1); - CCTK_CoordRegisterRange(cctkGH,lowery,uppery,-1,"y","cart3d"); - CCTK_CoordRegisterRangePhysIndex(cctkGH,loweryi,upperyi,-1,"y","cart3d"); - - lowerz = z_origin; - upperz = z_origin+this_dz*(cctk_gsh[2]-1); - CCTK_CoordRegisterRange(cctkGH,lowerz,upperz,-1,"z","cart3d"); - CCTK_CoordRegisterRangePhysIndex(cctkGH,lowerzi,upperzi,-1,"z","cart3d"); + CCTK_INFO("Computational Coordinates:"); + sprintf(infoline," %s[%6.3f,%6.3f] %s[%6.3f,%6.3f] %s[%6.3f,%6.3f] ", + "x=>",lowerx,upperx, + "y=>",lowery,uppery, + "z=>",lowerz,upperz); + CCTK_INFO(infoline); + + CCTK_INFO("Indices of Physical Coordinates:"); + sprintf(infoline," %s[%d,%d] %s[%d,%d] %s[%d,%d] ", + "x=>",lowerxi,upperxi, + "y=>",loweryi,upperyi, + "z=>",lowerzi,upperzi); + CCTK_INFO(infoline); + + } /* if (coarsest refinement level) */ - CCTK_INFO("Grid Spacings:"); - sprintf(infoline," %s%12.7e %s%12.7e %s%12.7e ", - "dx=>",cctk_delta_space[0], - "dy=>",cctk_delta_space[1], - "dz=>",cctk_delta_space[2]); - CCTK_INFO(infoline); - CCTK_INFO("Computational Coordinates:"); - sprintf(infoline," %s[%6.3f,%6.3f] %s[%6.3f,%6.3f] %s[%6.3f,%6.3f] ", - "x=>",lowerx,upperx, - "y=>",lowery,uppery, - "z=>",lowerz,upperz); - CCTK_INFO(infoline); - - CCTK_INFO("Indices of Physical Coordinates:"); - sprintf(infoline," %s[%d,%d] %s[%d,%d] %s[%d,%d] ", - "x=>",lowerxi,upperxi, - "y=>",loweryi,upperyi, - "z=>",lowerzi,upperzi); - CCTK_INFO(infoline); - - #ifdef CCTK_DEBUG printf("\n"); printf("CartGrid3D\n"); |