// The Patch class. An instance of this class describes an // atomic (rectangular) section of the fields at some level. // It contains structures that specify relationships with // the grids in other patches. // // We use patch-level parallelism. Every processor owns 0 // or more patches. /****************************************************************************** CHANGE LOG. 25 Aug 98: Modified the gradient calculation so that it tries to get phi up to date first, as is done in computePhiLaplacian(). Also, removed the use of the transpose() grid method in the gradient calculation because I don't trust it. ******************************************************************************/ class Patch { // The patch index identifies the box in each field // corresponding to this patch. int index; // The level of this patch. int level; // DOMAIN is the domain of this patch. GDOMAIN is // DOMAIN + ghost cells. COARSEDOMAIN is the coarsened // DOMAIN. The boxes containing the field values for a // patch have domain DOMAIN or GDOMAIN. For instance, // PHI has boundary (ghost) cells, but RHO doesn't. RectDomain<3> domain; RectDomain<3> gdomain; RectDomain<3> coarseDomain; // Pointers to field grids. The entire Field F is also // referenced by AMR.F. double [3d] phi; double [3d] savedPhi; double [3d] residual; double [3d] rho; double [3d] error; double [3d] laplacian; double [3d] relaxFactor; int [3d] tagged; // The flux registers for this patch are defined on the // next coarser level. double [1d] [1d] [3d] fluxRegisters; int [1d] [1d] signatures; double delta; double rhoNorm; // The boundaries of DOMAIN. They are cell-centered and just // outside DOMAIN. They are indexed by [DIRECTION][SIDE]. // ISPHYSICAL is true when the boundary on that direction and side // is a physical boundary. The problem domain is a box, so a // boundary may not be only partially physical. RectDomain<3> [1d][1d] boundaries; boolean [1d][1d] isPhysical; // Red-black domains. REDBLACK[0] = RED, REDBLACK[1] = BLACK, and // RED + BLACK = DOMAIN. // Domain<3> [1d] redblack; // Interrelationships. // Patches partially covered by this one. RelatedPatch [1d] coarserPatches; // Patches which partially cover this one. RelatedPatch [1d] finerPatches; // Boundaries with coarser cells, indexed by [DIRECTION][SIDE]. CoarseFineBoundary [1d][1d] coarseBoundaries; // Boundary patches at the same level of this patch. RelatedPatch [1d] adjacentPatches; // Boundaries with finer level. Need these to access the flux // corrections during computation of the laplacian. RelatedPatch [1d] fineBoundaryPatches; static final boolean doReflux = true; static final boolean dontReflux = false; void coarseToFine(Field F, boolean reflux) { if ((AMR.debug & 1) != 0) System.out.println("coarseToFine(reflux=" + reflux + "):"); double [3d] fine = F.grids[index]; for (int dir = 1; dir <= 3; dir++) { for (int side = -1; side <= 1; side += 2) { if ((AMR.debug & 1) != 0) System.out.println("coarseToFine dir=" + dir + " side=" + side); CoarseFineBoundary cb = coarseBoundaries[dir][side]; if (cb == null) continue; // Useful numbers and vectors. // Normal directions. Point<3> normDir = Point<3>.direction(dir); Point<3> orient = normDir * side; Point<3> offset = (normDir - orient) / 2; // Tangential directions. int k1 = dir % 3 + 1; int k2 = (dir + 1) % 3 + 1; Point<3> t1 = Point<3>.direction(k1); Point<3> t2 = Point<3>.direction(k2); Point<3> tt = t1 + t2; updateCoarseBoundary(cb, F); // Interpolate the ghost cell values. double [3d] cv = cb.coarseValues; double [3d] fr = fluxRegisters[dir][side]; foreach (p in cb.coarseBoundary) { // 19 of 19 if (reflux) { // Initialize flux register. Sign convention // and centering is done in a way that the // correction is done by adding the contents // of the flux register to the value of the // operator at the coarse cell centered in the // same location: foreach (dir, side) // foreach (p in coarseBoundary) // laplacian[kCoarse][p] += // fluxRegister[dir][side][p]. fr[p] = (cv[p] - cv[p - orient]) / (AMR.nRefine * delta); } double x1 = (cv[p + t1] - cv[p - t1]) / 2; double x2 = (cv[p + t2] - cv[p - t2]) / 2; double x1x1 = (cv[p + t1] + cv[p - t1] - 2 * cv[p]); double x2x2 = (cv[p + t2] + cv[p - t2] - 2 * cv[p]); double x1x2 = (cv[p + t1 + t2] + cv[p - t1 - t2] - cv[p - t1 + t2] - cv[p + t1 - t2]) / 4; // Compute interpolant. Point<3> pf0 = (p + offset) * AMR.nRefine - offset; foreach (q in [[0, 0, 0] : tt * AMR.nRefine + normDir - [1, 1, 1]]) { double xx1 = (q[k1] + .5) / AMR.nRefine - .5; double xx2 = (q[k2] + .5) / AMR.nRefine - .5; fine[pf0 + q] = cv[p] + xx1 * x1 + xx2 * x2 + 0.5 * xx1 * xx1 * x1x1 + 0.5 * xx2 * xx2 * x2x2 + xx1 * xx2 * x1x2; fine[pf0 + q] = AMR.c1 * fine[pf0 + q] + AMR.c2 * fine[pf0 + q - orient] + AMR.c3 * fine[pf0 + q - [2, 2, 2] * orient]; if (reflux) { // Increment values in flux register with // properly normalized fine field fluxes. fr[p] -= (fine[pf0 + q] - fine[pf0 + q - orient]) / (delta * Math.pow(AMR.nRefine, 2)); } } } } } if ((AMR.debug & 1) != 0) System.out.println("coarseToFine(reflux=" + reflux + ") done"); } void updateCoarseBoundary(CoarseFineBoundary cb, Field F) { // Update with coarse values where available. cb.coarseValues.set(0); if (AMR.verbose) { System.out.print("collecting coarse boundary "); Util.printlnRD3(cb.coarseValues.domain()); } foreach (i in cb.coarsePatches.domain()) { Patch cp = cb.coarsePatches[i]; if (AMR.verbose) { System.out.print("copying from coarse field "); Util.printlnRD3(cp.domain); } cb.coarseValues.copy(F.grids[cp.index]); } // Collect fine values. if (AMR.verbose) { System.out.print("collecting fine boundary "); Util.printlnRD3(cb.fineValues.domain()); } foreach (i in cb.finePatches.domain()) { Patch fp = cb.finePatches[i]; if (AMR.verbose) { System.out.print("copying from fine field "); Util.printlnRD3(fp.domain); } cb.fineValues.copy(F.grids[fp.index].restrict(fp.domain)); } // Average and correct where covered by fine patches. foreach (p in cb.coveringFine / AMR.nRefineP) { cb.coarseValues[p] = 0; } foreach (p in cb.coveringFine) { cb.coarseValues[p / AMR.nRefine] += cb.fineValues[p] / Math.pow(AMR.nRefine, 3); } if (AMR.verbose) { System.out.println("updating by laplacian"); Util.printlnD3(cb.insideLap); } foreach (p in cb.insideLap) { double L = -6 * cb.fineValues[p]; L += cb.fineValues[p + Point<3>.direction(1)] + cb.fineValues[p - Point<3>.direction(1)] + cb.fineValues[p + Point<3>.direction(2)] + cb.fineValues[p - Point<3>.direction(2)] + cb.fineValues[p + Point<3>.direction(3)] + cb.fineValues[p - Point<3>.direction(3)]; Point<3> q = p / AMR.nRefine; cb.coarseValues[q] += L / 2 / cb.numLaps[q]; } } // Before calling this, make sure that the flux // registers at L+1 are up to date, by calling // coarseToFine(phi, doReflux) over all L+1 grids. void computePhiLaplacian(boolean reflux) { if ((AMR.debug & 1) != 0) System.out.println("computePhiLaplacian(reflux=" + reflux + ")"); // Set ghost cells on boundaries with coarser grids, // unless this is the coarsest level. if (! coarserPatches.domain().isEmpty()) { coarseToFine(AMR.phi, dontReflux); } physicalBC(AMR.phi); //// if (Ti.thisProc() == 0) AMRCanvas.dump(phi); updateGhost(AMR.phi); //// if (Ti.thisProc() == 0) AMRCanvas.dump(phi); computeLaplacian(AMR.phi); if (reflux) { foreach (i in fineBoundaryPatches.domain()) { RelatedPatch rp = fineBoundaryPatches[i]; for (int dir = 1; dir <= 3; dir++) { for (int side = -1; side <= 1; side += 2) { double [3d] fr = rp.patch.fluxRegisters[dir][side]; // This domain intersection is one // of the few things I don't bother // to precompute. It can be null. foreach (p in laplacian.domain() * fr.domain()) { laplacian[p] += fr[p]; } } } } } if ((AMR.debug & 1) != 0) System.out.println("computePhiLaplacian(reflux=" + reflux + ") done"); } // Set ghost cells to impose physical BCs. void physicalBC(Field F) { double [3d] field = F.grids[index]; for (int dir = 1; dir <= 3; dir++) { for (int side = -1; side <= 1; side += 2) { if (isPhysical[dir][side]) { Point<3> p1 = Point<3>.direction(dir, side); foreach (p in boundaries[dir][side]) { // field[p] = field[p - p1 * 2] / 3 - // field[p - p1] * 2; field[p] = - field[p - p1]; } } } } } static double [3d] XGradient(double [3d] a) { return Gradient(a, [1, 0, 0]); } static double [3d] YGradient(double [3d] a) { return Gradient(a, [0, 1, 0]); } static double [3d] ZGradient(double [3d] a) { return Gradient(a, [0, 0, 1]); } /* Compute the a component of the gradient of a. */ static double [3d] Gradient(double [3d] a, Point<3> d) { RectDomain<3> domain = [a.domain().min() + [1, 1, 1] : a.domain().max() - [1, 1, 1]]; double [3d] result = new double [domain]; foreach (p in domain) { result[p] = a[p + d] - a[p - d]; } return result; } double [1d] [3d] gradPhi() { double [1d] [3d] r = new double [[0 : 2]] [3d]; // Set ghost cells on boundaries with coarser grids, // unless this is the coarsest level. if (! coarserPatches.domain().isEmpty()) { coarseToFine(AMR.phi, dontReflux); } physicalBC(AMR.phi); updateGhost(AMR.phi); r[0] = XGradient(phi); r[1] = YGradient(phi); r[2] = ZGradient(phi); return r; } void computeLaplacian(Field F) { double [3d] field = F.grids[index]; foreach (p in laplacian.domain()) { laplacian[p] = ( field[p + [ 0, 0, 1]] + field[p + [ 0, 0,-1]] + field[p + [ 0, 1, 0]] + field[p + [ 0,-1, 0]] + field[p + [ 1, 0, 0]] + field[p + [-1, 0, 0]] - 6 * field[p]) / Math.pow(delta, 2); } } // Update field (phi or error) ghost cells from adjacent grids. void updateGhost(Field F) { double [3d] field = F.grids[index]; foreach (p in adjacentPatches.domain()) { RelatedPatch rp = adjacentPatches[p]; double [3d] fromGrid = F.grids[rp.patch.index]; if ((AMR.debug & 1) != 0 && Ti.thisProc() == 0) { System.out.print("updateGhost(field #"); System.out.print(index); System.out.print(") this patch: "); Util.printlnRD3(domain); System.out.print("Field domain: "); Util.printlnRD3(field.domain()); System.out.print(rp.patch.index); System.out.print("Related patch: "); Util.printlnRD3(rp.patch.domain); System.out.print("Related field domain: "); Util.printlnRD3(fromGrid.domain()); } // if (Ti.thisProc() == 0) AMRCanvas.dump(field); field.copy(fromGrid.restrict(rp.patch.domain)); // if (Ti.thisProc() == 0) AMRCanvas.dump(field); } } void setResidualErrorLapDifference(Field F) { double [3d] dest = F.grids[index]; updateGhost(AMR.error); physicalBC(AMR.error); coarseToFine(AMR.error, dontReflux); computeLaplacian(AMR.error); foreach (p in dest.domain()) { dest[p] = residual[p] - laplacian[p]; } } // Accrete X by 1 on both sides of all directions except DIR. static Domain<3> tAccrete(Domain<3> x, int dir) { Domain<3> r = x; foreach (d in [1:3]) { if (d[1] != dir) { r = r + (x + Point<3>.direction(d[1])) + (x - Point<3>.direction(d[1])); } } return r; } // Compute the difference between the laplacian of the coarsened // solution, and the coarsened right-hand side. Compare with norm // of fine rho, and tag points where too large. boolean estimateError() { double tolerance = rhoNorm * 0.01; boolean anyTagged = false; foreach (p in domain / AMR.nRefineP) { double error = 0; foreach (q in AMR.refinedCell) /*unroll*/ { double coarseLaplacian = phi[p * AMR.nRefineP + q + [ 0, 0, 1]] + phi[p * AMR.nRefineP + q + [ 0, 1, 0]] + phi[p * AMR.nRefineP + q + [ 1, 0, 0]] + phi[p * AMR.nRefineP + q + [ 0, 0, -1]] + phi[p * AMR.nRefineP + q + [ 0, -1, 0]] + phi[p * AMR.nRefineP + q + [-1, 0, 0]] - phi[p * AMR.nRefineP + q] * 6; error += coarseLaplacian - rho[p * AMR.nRefineP + q]; } error *= 1 / (Math.pow(delta, 2) * Math.pow(AMR.nRefine, 4)); tagged[p] = Math.abs(error) > tolerance ? 1 : 0; anyTagged = anyTagged || (tagged[p] == 1); } return anyTagged; } // Go through all other patches at level L (this level), // above, and below, and compute relationships. void computeRelationships(int l) { BoxedList_RelatedPatch adjacentList = new BoxedList_RelatedPatch(); BoxedList_RelatedPatch coarserList = new BoxedList_RelatedPatch(); BoxedList_RelatedPatch finerList = new BoxedList_RelatedPatch(); BoxedList_RelatedPatch fineBoundaryList = new BoxedList_RelatedPatch(); foreach (p in AMR.processes.domain()) { AMRProcess ap = AMR.processes[p]; foreach (ii in ap.levels[l].patches.domain()) { Patch otherPatch = ap.levels[l].patches[ii]; if (this != otherPatch) { RectDomain<3> b = gdomain * otherPatch.domain; if (! b.isEmpty()) { RelatedPatch rp = new RelatedPatch(); rp.domain = b; rp.patch = otherPatch; adjacentList.push(rp); } } } } adjacentPatches = adjacentList.toArray(); fluxRegisters = new double [[1:3]][[-1:1:2]][3d]; coarseBoundaries = new CoarseFineBoundary [[1:3]][[-1:1:2]]; foreach (p in AMR.processes.domain()) { AMRProcess ap = AMR.processes[p]; if (l > 0) { foreach (ii in ap.levels[l-1].patches.domain()) { Patch otherPatch = ap.levels[l-1].patches[ii]; RectDomain<3> x = otherPatch.domain * coarseDomain; if (! x.isEmpty()) { RelatedPatch rp = new RelatedPatch(); rp.patch = otherPatch; rp.domain = x; coarserList.push(rp); } } // Compute coarse-fine boundaries. Must // find all coarse patches that contribute // to the boundary computation. foreach (dir in [1:3]) { foreach (side in [-1:1:2]) { BoxedList_Patch patches = new BoxedList_Patch(); RectDomain<3> fb = boundaries[dir][side]; // Allocate flux register. RectDomain<3> cb = fb / AMR.nRefineP; if (fluxRegisters[dir][side] == null) fluxRegisters[dir][side] = new double[cb]; // See if any part of cb is a true // coarse-fine boundary by comparing // it with all adjacent fine patches. Domain<3> cfb = (Domain<3>) fb; Domain<3> acfb; Domain<3> nacfb; foreach (r in adjacentPatches.domain()) { RelatedPatch rp = adjacentPatches[r]; cfb = cfb - rp.domain; } if (! cfb.isEmpty() && ! isPhysical[dir][side]) { // Coarsen cfb. cfb = cfb / AMR.nRefineP; // Accrete it tangentially. acfb = tAccrete(cfb, dir[1]); // Take union with CFB accreted // on inward normal. nacfb = acfb + (cfb - Point<3>.direction(dir[1], side[1])); foreach (ii in ap.levels[l-1].patches.domain()) { Patch otherPatch = ap.levels[l-1].patches[ii]; if (! (nacfb * otherPatch.domain).isEmpty()) { patches.push(otherPatch); } } } if (patches.length() > 0) if (coarseBoundaries[dir][side] == null) coarseBoundaries[dir][side] = new CoarseFineBoundary(cfb, acfb, nacfb, patches, l, dir, side); else { System.err.println("Input hit internal bug:" + " Coarse boundary can't be" + " shared across procs"); System.exit(-1); } // coarseBoundaries[dir][side] = // patches.length() == 0 ? // null : // new CoarseFineBoundary(cfb, acfb, nacfb, patches, // l, dir, side); } } } if (l < AMR.nLevels - 1) { foreach (ii in ap.levels[l+1].patches.domain()) { Patch otherPatch = ap.levels[l+1].patches[ii]; RectDomain<3> x = domain * (otherPatch.domain / AMR.nRefineP); if (! x.isEmpty()) { RelatedPatch rp = new RelatedPatch(); rp.patch = otherPatch; rp.domain = x; finerList.push(rp); } RectDomain<3> y = gdomain * (otherPatch.domain / AMR.nRefineP); if (! x.isEmpty()) { RelatedPatch rp = new RelatedPatch(); rp.patch = otherPatch; rp.domain = y; fineBoundaryList.push(rp); } } } } coarserPatches = coarserList.toArray(); finerPatches = finerList.toArray(); fineBoundaryPatches = fineBoundaryList.toArray(); } // Allocate and initialize grids for a patch. DOMAIN and GDOMAIN // must be set before calling this. Also initialize global fields. void allocateGrids() { phi = new double[gdomain]; savedPhi = new double[gdomain]; error = new double[gdomain]; residual = new double[domain]; rho = new double[domain]; laplacian = new double[domain]; relaxFactor = new double[domain]; tagged = new int[domain / AMR.nRefineP]; RectDomain<3> sdomain = [domain.min() + [1, 1, 1] : domain.max() - [1, 1, 1]]; foreach (p in sdomain) { int x = 1; relaxFactor[p] = x / 6.0; } foreach (p in domain - sdomain) { int x = 1; relaxFactor[p] = x / 6.0; } } // Regridding. See tagErrors() for comments. private int [1d] /*local*/ [3d] /*local*/ scan; private /*local*/ final /*inline*/ boolean isTagged(Point<3> p) { return tagged[p] > 0; } private static final /*inline*/ int tagsInRect2(int [2d] /*local*/ f, Point<2> fmin, Point<2> p) { if (p[1] < fmin[1] || p[2] < fmin[2]) return 0; else return f[p]; } /* f is constructed to contain a 2d cumulative histogram, i.e., f[[q, r]] contains the number of tagged points in this patch with coordinate 1 < q and coordinate 2 < r. (The coords can be (y, z) or (x, z) or (x, y) depending on whether dim is 0, 1, or 2.) The last few lines are to compute the number of tagged points with q0 < coord 1 < q1 and r0 < coord2 < r1 given f. Return value: the number of tagged points in this patch in the box and in the plane dim=value. Note that the dim=value plane must intersect the box. */ private final /*inline*/ int tagsInRect(int dim, int value, RectDomain<3> box) { int [2d] /*local*/ f = scan[dim].slice(dim + 1, value); Point<2> fmin = f.domain().min(), boxmin = box.slice(dim + 1).min(), boxmax = box.slice(dim + 1).max(); if (fmin == boxmin) return f[boxmax]; // speeds up common case else { boxmin = boxmin - [1, 1]; return tagsInRect2(f, fmin, boxmax) + tagsInRect2(f, fmin, boxmin) - tagsInRect2(f, fmin, [boxmax[1], boxmin[2]]) - tagsInRect2(f, fmin, [boxmin[1], boxmax[2]]); } } /* Return value, ret, should have ret[dim][value] contain the number of tagged cells in plane dim=value of box intersected with this. */ public /*local*/ final int [1d] /*local*/ [1d] tagsInPlanes(RectDomain<3> box) { int [1d] /*local*/ [1d] ret = new int [[0 : 2]] [1d]; RectDomain<3> N = taggedRect(); box = box * N; if (box.isEmpty()) { RectDomain<1> R = [0 : 0]; ret[0] = new int[R]; ret[1] = new int[R]; ret[2] = new int[R]; } else { for (int dim = 0; dim < 3; dim++) { RectDomain<1> R = [box.min()[dim] : box.max()[dim] + 1]; ret[dim] = new int[R]; ///foreach (i in R) ret[dim][i] = tagsInRect(dim, i, box); } } return ret; } private static final /*inline*/ void min3(int[] x, Point<3> p) { if (p[1] < x[0]) x[0] = p[1]; if (p[2] < x[1]) x[1] = p[2]; if (p[3] < x[2]) x[2] = p[3]; } private static final /*inline*/ void max3(int[] x, Point<3> p) { if (p[1] > x[0]) x[0] = p[1]; if (p[2] > x[1]) x[1] = p[2]; if (p[3] > x[2]) x[2] = p[3]; } private static final void rescan(int [3d] /*local*/ a, int dim) { RectDomain<3> R = a.domain(); foreach (p2 in R.slice(dim + 1)) { int s = 0, q = R.min()[dim], qmax = R.max()[dim]; Point<3> p = (dim == 0) ? [q, p2[1], p2[2]] : ((dim == 1) ? [p2[1], q, p2[2]] : [p2[1], p2[2], q]); Point<3> pinc = Point<3>.direction(dim); for (; q <= qmax; q++, p = p + pinc) { s += a[p]; a[p] = s; } } } /* tagErrors() finds the cells with high error and sets up data structures to quickly query the spatial structure of the errors later. In particular, we need to answer questions of the form: Given a box b and a plane p, how many tagged cells are in the intersection of this patch, b, and p? */ public /*local*/ final void tagErrors() { RectDomain<3> R = error.domain(); RectDomain<3> N; // subset of R that contains tagged points int [1d] /*local*/ [3d] /*local*/ scan = new int [[0 : 2]][3d]; this.scan = scan; /* Define f_{x=i}(j,k) to be the number of tagged cells in the intersection of the x = i plane and the y <= j and z <= k halfspaces. Define f_{y=j}(i,k) and f_{z=k}(i,j) analogously. The goal of the code below is to produce the array scan, such that: scan[0][i,j,k] = f_{x=i}(j,k), scan[1][i,j,k] = f_{y=j}(i,k), and scan[2][i,j,k] = f_{z=k}(i,j). */ // X scan then Y, and compute N { int[] lo = new int[3], hi = new int[3]; scan[2] = new int[R]; boolean seenAny = false; foreach (yz in R.slice(1)) { for (int s = 0, x = R.min()[1], xmax = R.max()[1]; x <= xmax; x++) { Point<3> p = [x, yz[1], yz[2]]; if (isTagged(p)) { s++; if (seenAny) { min3(lo, p); max3(hi, p); } else { seenAny = true; lo[0] = p[1]; lo[1] = p[2]; lo[2] = p[3]; hi[0] = p[1]; hi[1] = p[2]; hi[2] = p[3]; } } scan[2][p] = s; } } N = [[lo[0], lo[1], lo[2]] : [hi[0], hi[1], hi[2]]]; scan[2] = scan[2].restrict(N); rescan(scan[2], 1); } // Y scan then Z scan[0] = new int[N]; foreach (xz in N.slice(2)) { for (int s = 0, y = R.min()[1], ymax = R.max()[1]; y <= ymax; y++) { Point<3> p = [xz[1], y, xz[2]]; if (isTagged(p)) s++; scan[0][p] = s; } } rescan(scan[0], 2); // Z scan then X scan[1] = new int[N]; foreach (xy in N.slice(3)) { for (int s = 0, z = R.min()[2], zmax = R.max()[2]; z <= zmax; z++) { Point<3> p = [xy[1], xy[2], z]; if (isTagged(p)) s++; scan[1][p] = s; } } rescan(scan[1], 0); } /* Returns a rectangle that just covers the tagged cells in the patch. */ public /*local*/ final RectDomain<3> taggedRect() { return scan[0].domain(); } }