// Per-process data structure. /****************************************************************************** CHANGE LOG. 23 Aug 98: Improved parsing. Added parsing of "astro" data format. [pike] 23 Nov 98: Added use of ExecutionLog. [pike] 2 Dec 98: Added some prefetching/postwriting in coarse-to-fine and fine-to-coarse calculations that, by their inter-level nature, do not use ghost regions. 28 Jan 99: Patches specified for processor p in the input file are mapped to processor (p mod Ti.numProcs()) so that we can still meaningfully run if the number of processors we have is less than the number expected by the input file. 1 Feb 99: Improved calculation of the norm of the residual. 11 Mar 99: Fixed bug in prefetching code in relaxCoarsetoFine(). ******************************************************************************/ import java.io.FileOutputStream; import java.io.FileInputStream; import java.io.PrintStream; import java.io.IOException; import java.io.StreamTokenizer; // Temporary. import ti.domains.*; // Temporary end. class AMRProcess { Level [1d] levels; double locResidualNorm; int locResidualNormCount; boolean anyTagged; // Iterate on a fixed grid hierarchy. // // Sometimes UP and DOWN stand for, respectively, FINE and COARSE. // Finest is top, coarsest is bottom. It's better to use FINE and // COARSE explicitly. // // 0 is the coarsest level, 1 is finer, etc. single boolean single mainLoop(int single maxIterations, boolean single checkConvergence) { int single iterations; for (iterations = 0; iterations < maxIterations; iterations++) { AMR.log.append("mainLoop: iteration " + iterations); AMR.log.begin("finest"); relaxFinest(); AMR.log.end("finest"); AMR.log.begin("fine to coarse"); relaxFineToCoarse(); AMR.log.end("fine to coarse"); AMR.log.begin("coarsest"); relaxCoarsest(); AMR.log.end("coarsest"); AMR.log.begin("coarse to fine"); relaxCoarseToFine(); AMR.log.end("coarse to fine"); AMR.log.begin("check for convergence"); // Check for convergence. if (checkConvergence) { boolean single converged = true; for (int single l = 0; l < AMR.nLevels; l++) { locResidualNorm = 0; foreach (i in levels[l].patches.domain()) { double [3d] local r = (double [3d] local) levels[l].patches[i].residual; locResidualNormCount += r.domain().size(); locResidualNorm += AMR.sum_of_squares(r, levels[l].delta); } boolean cv = true; AMR.barrier(); if (Ti.thisProc() == 0) { double norm = 0; int normCount = 0; foreach (i in AMR.processes.domain()) { norm += AMR.processes[i].locResidualNorm; normCount += AMR.processes[i].locResidualNormCount; } norm = Math.sqrt(norm / normCount); System.out.print("level = "); System.out.print(l); System.out.print(", residual = "); System.out.println(norm); if (norm > AMR.epsilon * AMR.rhoNorms[l]) cv = false; } converged &= broadcast cv from 0; } AMR.log.end("check for convergence"); // If requested, display RESIDUAL and PHI. if (AMR.dumping && Ti.thisProc() == 0) Dump.dump("iteration " + iterations + " "); if (AMR.interactive && Ti.thisProc() == 0) { AMRCanvas canvas = new AMRCanvas(); foreach (p in AMR.processes.domain()) { for (int l = 0; l < AMR.nLevels; l++) { Level level = AMR.processes[p].levels[l]; foreach (i in level.patches.domain()) { canvas.add(level.patches[i].residual, l); } } } canvas.interact(); canvas = new AMRCanvas(); foreach (p in AMR.processes.domain()) { for (int l = 0; l < AMR.nLevels; l++) { Level level = AMR.processes[p].levels[l]; foreach (i in level.patches.domain()) { canvas.add(level.patches[i].phi. restrict(level.patches[i].domain), l); } } } canvas.interact(); } AMR.barrier(); if (converged) break; } } // If requested, dump gradient of PHI to stdout. if (AMR.showGradPhi) { foreach (p in AMR.processes.domain()) { for (int l = 0; l < AMR.nLevels; l++) { Level level = AMR.processes[p].levels[l]; foreach (i in level.patches.domain()) { double [1d] [3d] grad = level.patches[i].gradPhi(); AMRCanvas.dump(grad[0], l, "grad phi X", System.out); AMRCanvas.dump(grad[1], l, "grad phi Y", System.out); AMRCanvas.dump(grad[2], l, "grad phi Z", System.out); AMRCanvas.dump(level.patches[i].phi, l, "phi", System.err); } } } } return iterations < maxIterations; } void relaxFinest() { Level l = levels[AMR.finest]; foreach (i in l.patches.domain()) { Patch patch = l.patches[i]; patch.computePhiLaplacian(Patch.dontReflux); foreach (p in patch.residual.domain()) { patch.residual[p] = patch.rho[p] - patch.laplacian[p]; } patch.error.set(0); } } single void relaxFineToCoarse() { for (int single l = AMR.finest; l > 0; l--) { // Clear correction at l-1. if (AMR.verbose) System.out.println("// Clear correction at l-1."); foreach (i in levels[l-1].patches.domain()) { Patch patch = levels[l-1].patches[i]; patch.error.set(0); } AMR.barrier(); foreach (i in levels[l].patches.domain()) { Patch patch = levels[l].patches[i]; // Relax at this level. if (AMR.verbose) System.out.println("// Relax at this level."); AMR.setBC(patch.error, patch.boundaries); AMR.redRelax(patch.error, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); } AMR.barrier(); foreach (i in levels[l].patches.domain()) { Patch patch = levels[l].patches[i]; AMR.setBC(patch.error, patch.boundaries); // Must update ghost cells AFTER setting BCs. patch.updateGhost(AMR.error); } AMR.barrier(); foreach (i in levels[l].patches.domain()) { Patch patch = levels[l].patches[i]; AMR.blackRelax(patch.error, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); // Save solution. if (AMR.verbose) System.out.println("// Save solution."); patch.savedPhi.copy(patch.phi); // Apply correction. if (AMR.verbose) System.out.println("// Apply correction."); foreach (p in patch.error.domain()) { patch.phi[p] += patch.error[p]; } // Update flux registers for forthcoming computation // of L(phi) one level down. patch.coarseToFine(AMR.phi, Patch.doReflux); } AMR.barrier(); foreach (i in levels[l-1].patches.domain()) { Patch patch = levels[l-1].patches[i]; patch.computePhiLaplacian(Patch.doReflux); // Set the residual at l-1, over all the // domain, then later set again the part // under the fine region. (The fine region // is small, so we don't bother to set only // the exposed part.) foreach (p in patch.domain) { patch.residual[p] = patch.rho[p] - patch.laplacian[p]; } } AMR.barrier(); { PrivateRegion local z; foreach (i in levels[l].patches.domain()) { Patch local patch = (Patch local) levels[l].patches[i]; // Compute residual - LNF(error) and temporarily store // in laplacian. patch.setResidualErrorLapDifference(AMR.laplacian); // Set the residual at l-1 under the fine region, by // simple average. foreach (j in patch.coarserPatches.domain()) { Patch rp = patch.coarserPatches[j].patch; RectDomain<3> rd = patch.coarserPatches[j].domain; boolean rpIsLocal = rp.residual.isLocal(); { if (!rpIsLocal) z = new PrivateRegion(); double [3d] local dest = rpIsLocal ? (double [3d] local) rp.residual : new (z) double [rd]; foreach (p in rd) { // rp.residual[p] = 0; // foreach (q in AMR.refinedCell) /*unroll*/ { // rp.residual[p] += // patch.laplacian[p * AMR.refineVector + q]; // } // rp.residual[p] /= 8; dest[p] = ( patch.laplacian[p * [2, 2, 2] + [0, 0, 0]] + patch.laplacian[p * [2, 2, 2] + [0, 0, 1]] + patch.laplacian[p * [2, 2, 2] + [0, 1, 0]] + patch.laplacian[p * [2, 2, 2] + [0, 1, 1]] + patch.laplacian[p * [2, 2, 2] + [1, 0, 0]] + patch.laplacian[p * [2, 2, 2] + [1, 0, 1]] + patch.laplacian[p * [2, 2, 2] + [1, 1, 0]] + patch.laplacian[p * [2, 2, 2] + [1, 1, 1]] ) / 8; } if (!rpIsLocal) rp.residual.copy(dest); } if (!rpIsLocal) try { z.delete(); } catch (RegionInUse x) { } } } AMR.barrier(); } } } single void relaxCoarsest() { if (AMR.verbose) System.out.println("relaxCoarsest"); foreach (i in levels[0].patches.domain()) { Patch patch = levels[0].patches[i]; // Clear correction. patch.error.set(0); // Relax/solve at level 0. Multigrid.relax(patch.error, patch.residual, Multigrid.residuals[0], patch.relaxFactor, patch.delta, 0); // Apply correction. foreach (p in patch.phi.domain()) { patch.phi[p] += patch.error[p]; } } AMR.barrier("after relaxCoarsest()"); } single void relaxCoarseToFine() { if (AMR.verbose) System.out.println("relaxCoarseToFine"); PrivateRegion local z; for (int single l = 1; l < AMR.nLevels; l++) { foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; double [3d] local e = (double [3d] local) patch.error; // Correct the correction by interpolation from // coarser levels. foreach (i in patch.coarserPatches.domain()) { RelatedPatch rp = patch.coarserPatches[i]; RectDomain<3> iterationDomain = rp.patch.domain * patch.coarseDomain; boolean rpIsLocal = rp.patch.error.isLocal(); if (!rpIsLocal) z = new PrivateRegion(); { double [3d] local src = rpIsLocal ? (double [3d] local) rp.patch.error : new (z) double [iterationDomain]; if (!rpIsLocal) src.copy(rp.patch.error); foreach (p in iterationDomain) { // foreach (q in AMR.refinedCell) /*unroll*/ { // patch.error[p * AMR.nRefineP + q] += // rp.patch.error[p]; // } double x = src[p]; e[p * 2 + [0, 0, 0]] += x; e[p * 2 + [0, 0, 1]] += x; e[p * 2 + [0, 1, 0]] += x; e[p * 2 + [0, 1, 1]] += x; e[p * 2 + [1, 0, 0]] += x; e[p * 2 + [1, 0, 1]] += x; e[p * 2 + [1, 1, 0]] += x; e[p * 2 + [1, 1, 1]] += x; } } if (!rpIsLocal) try { z.delete(); } catch (RegionInUse x) { System.err.println("region delete failed"); System.exit(-1); } } } AMR.barrier("after correcting the correction"); foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; // Decrement the residual by the laplacian of the // error. patch.setResidualErrorLapDifference(AMR.residual); } AMR.barrier("after decrementing the residual by the laplacian of the error"); foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; AMR.setBC(patch.phi, patch.boundaries); patch.updateGhost(AMR.phi); } AMR.barrier("after setting BCs and updating ghost regions"); foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; // Compute correction to apply to error. Use phi for // temporary space. patch.phi.set(0); AMR.redRelax(patch.phi, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); } AMR.barrier("after redrelax"); foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; AMR.setBC(patch.phi, patch.boundaries); patch.updateGhost(AMR.phi); } AMR.barrier(); foreach (j in levels[l].patches.domain()) { Patch patch = levels[l].patches[j]; AMR.blackRelax(patch.phi, patch.residual, patch.relaxFactor, patch.boundaries, patch.delta); foreach (p in patch.error.domain()) { // Correct correction. patch.error[p] += patch.phi[p]; // Compute new value of solution. patch.phi[p] = patch.savedPhi[p] + patch.error[p]; } } AMR.barrier("after black relax; done with relaxCoarseToFine"); } if (AMR.verbose) System.out.println("relaxCoarseToFine done"); } // Create new grids for hierarchy. Return TRUE for // error conditions. // // Construct new grid, from finest to coarsest level. A // new level may be added here, or several levels may // disappear. The new level is the one above level // HIGHEST-BAD-LEVEL. The regridding goes from that // extra level down to 1 (don't regrid level 0, that's // the problem domain). single boolean regridAsNeeded() { if (AMR.verbose) System.out.println("regridAsNeeded()"); int single highestBadLevel; int single newNLevels; int single newMaxLevel; for (int single l = AMR.nLevels - 1; l >= 0; l--) { // Compute tagged cells on each patch. Level level = levels[l]; anyTagged = false; foreach (i in level.patches.domain()) { anyTagged = anyTagged || level.patches[i].estimateError(); } AMR.barrier(); int hbl = -1; if (Ti.thisProc() == 0) { foreach (i in AMR.processes.domain()) { if (AMR.processes[i].anyTagged) { hbl = l; break; } } } highestBadLevel = broadcast hbl from 0; if (highestBadLevel > -1) break; } if (highestBadLevel == -1) return false; newMaxLevel = highestBadLevel + 1; newNLevels = newMaxLevel + 1; if (newNLevels > AMR.maxLevels) return true; // BUG: We may need to extend to levels[] array here. // When allocated initially, it only has enough levels for // the input hierarchy. for (int single l = newMaxLevel; l > 0; l--) levels[l].boxes = Refiner.refinement(levels[l-1], l > 1 ? levels[l-2] : null); allocatePatches(newNLevels); return false; } single void allocatePatches(int single newNLevels) { // Note: free the old patches here as needed. AMR.nLevels = newNLevels; AMR.finest = newNLevels - 1; int index = 0; for (int l = 0; l < newNLevels; l++) { levels[l].patches = new Patch[levels[l].boxes.domain()]; foreach (i in levels[l].patches.domain()) { Patch patch = new Patch(); levels[l].patches[i] = patch; patch.level = l; patch.domain = levels[l].boxes[i]; patch.coarseDomain = patch.domain / AMR.nRefineP; patch.gdomain = [patch.domain.min() - [1, 1, 1] : patch.domain.max() + [1, 1, 1]]; patch.boundaries = new RectDomain<3>[[1:3]][1d]; patch.boundaries[1] = new RectDomain<3>[[-1:1:2]]; patch.boundaries[2] = new RectDomain<3>[[-1:1:2]]; patch.boundaries[3] = new RectDomain<3>[[-1:1:2]]; patch.isPhysical = new boolean[[1:3]][1d]; patch.isPhysical[1] = new boolean[[-1:1:2]]; patch.isPhysical[2] = new boolean[[-1:1:2]]; patch.isPhysical[3] = new boolean[[-1:1:2]]; patch.delta = levels[l].delta; for (int dir = 1; dir <= 3; dir++) { for (int side = -1; side <= 1; side +=2) { RectDomain<3> b = (RectDomain<3>) ((patch.domain + Point<3>.direction(dir, side)) - patch.domain); patch.boundaries[dir][side] = b; RectDomain<3> physical = AMR.problemBoundaries[dir][side]; patch.isPhysical[dir][side] = ! ((b / Point<3>.all((int) Math.pow(AMR.nRefine, l))) * physical).isEmpty(); } } levels[l].patches[i] = patch; } } // On each level, Patch.computeRelationships() will need to // look at all other processes' patches on that level. This // barrier makes sure that all other processes have actually // finished building levels[*].patches[*] before continuing. AMR.barrier(); for (int l = 0; l < newNLevels; l++) { foreach (i in levels[l].patches.domain()) { Patch patch = levels[l].patches[i]; if (AMR.verbose) System.out.println("Allocating grids for patch " + Util.stringifyRD3(patch.domain) + " at level " + l + "."); // Allocate various grids. patch.allocateGrids(); // Compute relationships. patch.computeRelationships(l); } } // Global operations: number patches sequentially // and assign them to AMR.patches on process 0. // Also construct all fields. AMR.barrier(); if (Ti.thisProc() == 0) { index = 0; foreach (j in AMR.processes.domain()) { AMRProcess p = AMR.processes[j]; for (int l = 0; l < newNLevels; l++) { foreach (i in p.levels[l].patches.domain()) { Patch patch = p.levels[l].patches[i]; patch.index = index++; } } } AMR.patches = new Patch[[0 : index - 1]]; AMR.phi = new Field(AMR.patches); AMR.residual = new Field(AMR.patches); AMR.rho = new Field(AMR.patches); AMR.error = new Field(AMR.patches); AMR.laplacian = new Field(AMR.patches); index = 0; foreach (j in AMR.processes.domain()) { AMRProcess p = AMR.processes[j]; for (int l = 0; l < newNLevels; l++) { foreach (i in p.levels[l].patches.domain()) { Patch patch = p.levels[l].patches[i]; AMR.patches[index] = patch; AMR.phi.grids[index] = patch.phi; AMR.residual.grids[index] = patch.residual; AMR.rho.grids[index] = patch.rho; AMR.error.grids[index] = patch.error; AMR.laplacian.grids[index] = patch.laplacian; index++; } } } } AMR.barrier(); // Replicate the AMR global structures on all // processes. RectDomain<1> patchesD = broadcast AMR.patches.domain() from 0; Patch [1d] patches = broadcast AMR.patches from 0; if (Ti.thisProc() != 0) { AMR.patches = new Patch [patchesD]; AMR.patches.copy(patches); } AMR.phi = (broadcast AMR.phi from 0).replicate(); AMR.residual = (broadcast AMR.residual from 0).replicate(); AMR.rho = (broadcast AMR.rho from 0).replicate(); AMR.error = (broadcast AMR.error from 0).replicate(); AMR.laplacian = (broadcast AMR.laplacian from 0).replicate(); } static void inputError(String s) throws IOException { System.err.println(s + "\nAborting AMR."); throw new IOException(); } int inputProblemSize; RectDomain<3> [1d][1d][1d] inputBoxes; // [level][proc][i] StreamTokenizer sttok; // Read hierarchy. single void inputHierarchy(String filename) throws IOException { int l1, l2, l3, u1, u2, u3; int token; boolean p0 = Ti.thisProc() == 0; if (p0) { sttok = new StreamTokenizer(new FileInputStream(filename)); sttok.slashSlashComments(true); sttok.eolIsSignificant(false); token = sttok.nextToken(); // problemSize token = sttok.nextToken(); // = token = sttok.nextToken(); //_problemsize inputProblemSize = (int) sttok.nval; token = sttok.nextToken(); // nLevels token = sttok.nextToken(); // = token = sttok.nextToken(); //_numberoflevels int nlevels = (int) sttok.nval; token = sttok.nextToken(); // nProcs token = sttok.nextToken(); // = token = sttok.nextToken(); //_numberofprocessors int nprocs = (int) sttok.nval; if (nprocs != Ti.numProcs()) System.err.println("Warning: Running with " + Ti.numProcs() + ", but input specifies " + nprocs + "."); // Allocate levels. BoxedList_RectDomain_3 [2d] local inputBoxesList = new BoxedList_RectDomain_3[[0 : nlevels - 1, 0 : nprocs - 1]]; // Read box hierarchy. sttok.nextToken(); while (sttok.sval.equals("level")) { token = sttok.nextToken(); //_levelnumber int l = (int) sttok.nval; token = sttok.nextToken(); // processor token = sttok.nextToken(); //_processornumber int p = ((int) sttok.nval) % Ti.numProcs(); BoxedList_RectDomain_3 boxes = new BoxedList_RectDomain_3(); token = sttok.nextToken(); while (token == StreamTokenizer.TT_NUMBER) { l1 = (int) sttok.nval; token = sttok.nextToken(); l2 = (int) sttok.nval; token = sttok.nextToken(); l3 = (int) sttok.nval; token = sttok.nextToken(); u1 = (int) sttok.nval; token = sttok.nextToken(); u2 = (int) sttok.nval; token = sttok.nextToken(); u3 = (int) sttok.nval; token = sttok.nextToken(); boxes.push([[l1, l2, l3] : [u1, u2, u3]]); if (true /* AMR.verbose */) { System.out.print("Input box for p"); System.out.print(p); System.out.print(" at level "); System.out.print(l); System.out.print(": "); Util.printlnRD3([[l1, l2, l3] : [u1, u2, u3]]); } } inputBoxesList[l, p] = boxes.concat(inputBoxesList[l, p]); } inputBoxes = new RectDomain<3> [[0 : nlevels - 1]][[0 : nprocs - 1]][1d]; foreach (l in inputBoxes.domain()) { foreach (p in inputBoxes[l].domain()) { inputBoxes[l][p] = (inputBoxesList[l[1], p[1]] == null) ? new RectDomain<3>[[0 : -1]] : inputBoxesList[l[1], p[1]].toArray(); if (AMR.verbose) System.out.println( "inputBoxes[l=" + l[1] + "][p=" + p[1] + "] = " + Util.stringifyBoxedListRD3(inputBoxesList [l[1], p[1]])); } } } } // Read a rectangle of data in astro format void readBox(StreamTokenizer sttok, int level, RectDomain<3> box) throws IOException { System.out.println("Reading box at level " + level + ": " + box); int [3d] local count = new int[box]; double [3d] xgrad, ygrad, zgrad, b = new double[box]; if (AMR.inputData[level] == null) AMR.inputData[level] = new List_double_3d_ (b, null); else AMR.inputData[level] = AMR.inputData[level].cons(b); if (AMR.astroInput && AMR.showGradPhi) { xgrad = new double[box]; ygrad = new double[box]; zgrad = new double[box]; } foreach (q in box) { int l1, l2, l3; Point<3> p; Util.skipUpTo(sttok, "<"); l1 = Util.getInt(sttok); Util.expect(sttok, ","); l2 = Util.getInt(sttok); Util.expect(sttok, ","); l3 = Util.getInt(sttok); Util.expect(sttok, ">"); p = [l1, l2, l3]; b[p] = Util.getDouble(sttok); // System.out.println("Got " + b[p] + " at " + Util.stringifyP3(p)); count[p]++; if (AMR.astroInput && AMR.showGradPhi) { // Skip 4 values, then the gradient is the next three Util.getDouble(sttok); Util.getDouble(sttok); Util.getDouble(sttok); Util.getDouble(sttok); xgrad[p] = Util.getDouble(sttok); ygrad[p] = Util.getDouble(sttok); zgrad[p] = Util.getDouble(sttok); } } foreach (p in box) { if (count[p] != 1) { inputError("Error reading box " + Util.stringifyRD3(box) + " at level " + level + ":\n value at " + Util.stringifyP3(p) + " read " + count[p] + " times."); } } if (AMR.astroInput && AMR.showGradPhi) { AMRCanvas.dump(xgrad, level, "input X gradient", System.out); AMRCanvas.dump(ygrad, level, "input Y gradient", System.out); AMRCanvas.dump(zgrad, level, "input Z gradient", System.out); } } // Read hierarchy. single void inputHierarchyAstro(String filename) throws IOException { int l1, l2, l3, u1, u2, u3; boolean p0 = Ti.thisProc() == 0; if (p0) { sttok = new StreamTokenizer(new FileInputStream(filename)); // sttok.slashSlashComments(true); sttok.eolIsSignificant(false); sttok.wordChars('_', '_'); // Get problem size try { Util.skipUpTo(sttok, "n_cell"); Util.skipUpTo(sttok, "n_cell"); Util.expect(sttok, "="); inputProblemSize = Util.getInt(sttok); } catch (Throwable x) { inputError("Error reading problem size."); } // refinement ratio should be 2 try { Util.skipUpTo(sttok, "ref_ratio"); Util.expect(sttok, "="); } catch (Throwable x) { inputError("Unable to find ref_ratio in input."); } try { Util.expectInt(sttok, 2); } catch (Throwable x) { inputError("Refinement ratio must be 2."); } // Get level numbers BoxedList_int levelNumbersList = new BoxedList_int (); int maxLevel = 0; do { Util.skipUpTo(sttok, "lev"); Util.expect(sttok, "="); if (AMR.verbose) System.out.println("Got a grid at level " + (int)sttok.nval); levelNumbersList.push(Util.getInt(sttok)); if (levelNumbersList.first() > maxLevel) maxLevel = levelNumbersList.first(); Util.skipUpTo(sttok, ">"); Util.skipUpTo(sttok, ">"); Util.skipUpTo(sttok, ">"); } while (sttok.ttype == StreamTokenizer.TT_WORD && sttok.sval.equals("Constructing")); // assume level 0 is lowest int nlevels = maxLevel + 1; int k = levelNumbersList.length(); int [1d] levelNumbers = levelNumbersList.toArray(); int readingLevel = -999; BoxedList_RectDomain_3 boxes; // Allocate levels. int nprocs = Ti.numProcs(); inputBoxes = new RectDomain<3>[[0 : nlevels - 1]][1d][1d]; AMR.inputData = new List_double_3d_ [[0 : nlevels - 1]]; foreach (i in inputBoxes.domain()) { inputBoxes[i] = new RectDomain<3>[[0 : nprocs - 1]][1d]; foreach (j in inputBoxes[i].domain()) { inputBoxes[i][j] = new RectDomain<3>[[0 : -1]]; } } // Read box hierarchy and data. Fill in inputBoxes[][]. // Note that for now we put all boxes on proc 0. while (--k >= 0) { int level = levelNumbers[k]; if (readingLevel != level) { // We're starting to read the boxes for a new level if (readingLevel >= 0) { inputBoxes[readingLevel][0] = boxes.toArray(); if (AMR.verbose) System.out.println("inputBoxes[l=" + readingLevel + "][p=0] = " + Util.stringifyBoxedListRD3(boxes)); } boxes = new BoxedList_RectDomain_3 (); readingLevel = level; } Util.skipUpTo(sttok, "<"); l1 = Util.getInt(sttok); Util.expect(sttok, ","); l2 = Util.getInt(sttok); Util.expect(sttok, ","); l3 = Util.getInt(sttok); Util.expect(sttok, ">"); Util.expect(sttok, "<"); u1 = Util.getInt(sttok); Util.expect(sttok, ","); u2 = Util.getInt(sttok); Util.expect(sttok, ","); u3 = Util.getInt(sttok); Util.expect(sttok, ">"); Util.expect(sttok, "<"); boxes.push([[l1, l2, l3] : [u1, u2, u3]]); readBox(sttok, level, [[l1, l2, l3] : [u1, u2, u3]]); } if (readingLevel >= 0 && !boxes.isEmpty()) { inputBoxes[readingLevel][0] = boxes.toArray(); if (AMR.verbose) System.out.println("inputBoxes[l=" + readingLevel + "][p=0] = " + Util.stringifyBoxedListRD3(boxes)); } } } single void setup(String filename) throws IOException { boolean _die = false; // Input box hierarchy. try { if (AMR.astroInput) inputHierarchyAstro(filename); else inputHierarchy(filename); } catch (Throwable x) { System.err.println("I/O exception during hierarchy input."); _die = true; } if (!(broadcast _die from 0)) { // Set problem size. AMR.problemSize = broadcast inputProblemSize from 0; int single ps = AMR.problemSize; for (AMR.logProblemSize = -1; ps != 0; ps = ps >> 1, AMR.logProblemSize++); if (AMR.verbose) System.out.println("Problem size = " + AMR.problemSize); RectDomain<3>[1d][1d][1d] boxes = broadcast inputBoxes from 0; AMR.nLevels = (int single) boxes.domain().size(); levels = new Level[[0 : AMR.nLevels - 1]]; double d = 1.0 / AMR.problemSize; foreach (l in boxes.domain()) { levels[l] = new Level(); levels[l].boxes = boxes[l][Ti.thisProc()]; levels[l].delta = d; d /= 2; } // Build boundaries. int n = AMR.problemSize; AMR.problemDomain = [[0, 0, 0] : [n - 1, n - 1, n - 1]]; AMR.problemBoundaries = new RectDomain<3>[[0:3]][1d]; foreach (dir in [1:3]) { AMR.problemBoundaries[dir] = new RectDomain<3>[[-1:1:2]]; AMR.problemBoundaries[dir][-1] = AMR.problemDomain.border(1, -dir[1], 1); AMR.problemBoundaries[dir][1] = AMR.problemDomain.border(1, dir[1], 1); } allocatePatches(AMR.nLevels); if (AMR.astroInput) { List_double_3d_ b; foreach (ln in AMR.inputData.domain()) // ln is the level number for (b = AMR.inputData[ln]; b != null; b = b.rest()) { Level l = levels[ln]; foreach (j in l.patches.domain()) { // j is the patch number within patches at level n System.out.println("Copying from " + Util.stringifyRD3(b.first().domain()) + " to " + Util.stringifyRD3(l.patches[j].rho.domain()) + " at level " + ln[1]); l.patches[j].rho.copy(b.first()); } } } else { BoxedList_Point_3 charges; // Read charges and initialize rho. try { charges = inputRho(); } catch (Throwable x) { System.err.println("I/O exception during rho input."); _die = true; } int single i = broadcast charges.length() from 0; while (i-- > 0) setRho(broadcast charges.pop() from 0); } } if (broadcast _die from 0) throw new IOException(); } single BoxedList_Point_3 inputRho() throws IOException { int token; int l1, l2, l3, u1, u2, u3; BoxedList_Point_3 ret; boolean p0 = Ti.thisProc() == 0; if (p0) { ret = new BoxedList_Point_3(); while ((token = sttok.nextToken()) != StreamTokenizer.TT_EOF) { Point<3> p; if (Ti.thisProc() == 0) { l1 = (int) sttok.nval; token = sttok.nextToken(); l2 = (int) sttok.nval; token = sttok.nextToken(); l3 = (int) sttok.nval; p = [l1, l2, l3]; ret.push(p); } } } return ret; } // Set a point charge in rho for all relevant patches. void setRho(Point<3> p) { for (int i = 0; i < AMR.nLevels; i++, p = p * 2) { Level l = levels[i]; foreach (j in l.patches.domain()) { Patch patch = l.patches[j]; if (patch.domain.contains(p)) { patch.rho[p] = 1.23456789; } } } } // Testing methods. single void test() { // Initialize phi at all levels with quadratic function. for (int i = 0, n = AMR.problemSize; i < AMR.nLevels; i++, n *= 2) { Level l = levels[i]; foreach (j in l.patches.domain()) { Patch patch = l.patches[j]; foreach (p in patch.domain) { double x = p[1] / (double) n + 1. / (2 * n); double y = p[2] / (double) n + 1. / (2 * n); double z = p[3] / (double) n + 1. / (2 * n); patch.phi[p] = y * y; // patch.phi[p] = x * x; } } } AMR.barrier(); // Compute laplacian of phi at highest level. Level l = levels[AMR.finest]; foreach (j in l.patches.domain()) { Patch patch = l.patches[j]; // patch.coarseToFine(AMR.phi,patch.dontReflux); patch.computePhiLaplacian(patch.doReflux); } AMR.barrier(); l = levels[0]; foreach (j in l.patches.domain()) { Patch patch = l.patches[j]; // patch.coarseToFine(AMR.phi,patch.dontReflux); patch.computePhiLaplacian(patch.doReflux); } AMR.barrier(); if (Ti.thisProc() == 0) { AMRCanvas canvas = new AMRCanvas(); foreach (p in AMR.processes.domain()) { for (int ll = 0; ll < AMR.nLevels; ll++) { Level level = AMR.processes[p].levels[ll]; foreach (i in level.patches.domain()) { canvas.add(level.patches[i].laplacian // .restrict(level.patches[i].domain) , ll); } } } canvas.interact(); canvas = new AMRCanvas(); foreach (p in AMR.processes.domain()) { for (int ll = 0; ll < AMR.nLevels; ll++) { Level level = AMR.processes[p].levels[ll]; foreach (i in level.patches.domain()) { canvas.add(level.patches[i].phi // .restrict(level.patches[i].domain) , ll); } } } canvas.interact(); } AMR.barrier(); } }