// Titanium 3-D AMR Poisson solver. // // Written mainly by Luigi Semenzato. Algorithm by Phil // Colella. Geoff Pike wrote the initial version of the // regridding code, improving on an algorithm by // Berger/Rigoutsos. However, regridding was never finished. // // Copyright (C) The United States Department of // Energy and the University of California. // The AMR class, representing the full hierarchy of grids and all // else. /****************************************************************************** CHANGE LOG. 23 Aug 98: Improved parsing of command line arguments. Added -a flag. Added inputData to this class. [pike] 23 Nov 98: Added use of ExecutionLog. [pike] 25 Nov 98: Improved use of ExecutionLog. Added -l flag and -D flag. [pike] 1 Feb 99: Improved calculation of the norm of the residual. Added more debugging output. [pike] 9 Mar 99: Added AMR.barrier(String) for extra debugging opportunities. ******************************************************************************/ import java.io.IOException; class AMR { // Globals. static AMRProcess [1d] single processes; static int single problemSize; static int single logProblemSize; static RectDomain<3> problemDomain; static RectDomain<3> [1d][1d] problemBoundaries; static int single maxLevels; static int single nLevels; static int single finest; // nLevels - 1 static double [1d] rhoNorms; static List_double_3d_ [1d] inputData; // indexed by level // Fields. static Field phi; static Field residual; static Field rho; static Field error; static Field laplacian; // Patches. static Patch [1d] patches; // Computation control. static boolean single interactive; static boolean single oneshot; static boolean single dumping; static boolean single astroInput; static boolean single showGradPhi; static int single mainLoopIterations; static boolean single verbose; static boolean single testing; static int single debug; static ExecutionLog local log; static int barrierCount; // Constants. static final int maxRegriddings = 0; // Increase this to regrid. static final int coarsest = 0; static final int nRefine = 2; static final Point<3> nRefineP = [2, 2, 2]; static final double epsilon = 0.01; static final Point<3> refineVector = [nRefine, nRefine, nRefine]; static final RectDomain<3> refinedCell = [[0, 0, 0] : [nRefine, nRefine, nRefine] - [1, 1, 1]]; static final double c1 = 8. / ((nRefine + 1.) * (nRefine + 3.)); static final double c2 = 2. * (nRefine - 1.) / (nRefine + 1.); static final double c3 = - (nRefine - 1.) / (nRefine + 3.); static Multigrid mg = new Multigrid(); static final String bb = "begin barrier"; static final String eb = "end barrier"; // Barrier public static sglobal void barrier() { barrierCount++; log.append(2, bb); Ti.barrier(); log.append(2, eb); } public static sglobal void barrier(String local s) { if ((debug & 32) == 32) System.out.println(s); barrier(); } // Return the value of a 3d array element, for debugging. static double index3(double [3d] a, int i, int j, int k) { return a[[i, j, k]]; } // Compute the L1 norm. static double L1norm(double [3d] local field, double h) { double n = 0; foreach (p in field.domain()) n += Math.abs(field[p] * h); return n / field.domain().size(); } // Compute the norm. static double norm(double [3d] local field, double h) { return Math.sqrt(normsquared(field, h)); } // Compute the square of the norm. static double normsquared(double [3d] local field, double h) { return sum_of_squares(field, h) / field.domain().size(); } // Compute the sum over all p of the squares of h * field[p]. static double sum_of_squares(double [3d] local field, double h) { double n = 0; double hsq = h * h; foreach (p in field.domain()) n += field[p] * field[p] * hsq; return n; } // Set boundary conditions. static void setBC(double [3d] field, RectDomain<3> [1d][1d] boundaries) { for (int dir = 1; dir <= 3; dir++) { for (int side = -1; side <= 1; side += 2) { 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 void gsrbRelax(double [3d] field, double [3d] rhs, double [3d] factor, RectDomain<3> [1d][1d] boundaries, double h) { if (AMR.verbose) System.out.println("gsrbRelax"); setBC(field, boundaries); redRelax(field, rhs, factor, boundaries, h); setBC(field, boundaries); blackRelax(field, rhs, factor, boundaries, h); } static void redRelax(double [3d] field, double [3d] rhs, double [3d] factor, RectDomain<3> [1d][1d] boundaries, double h) { pointRelax(field, rhs, factor, h*h, [ 0, 0, 0]); pointRelax(field, rhs, factor, h*h, [ 0, -1, -1]); pointRelax(field, rhs, factor, h*h, [-1, -1, 0]); pointRelax(field, rhs, factor, h*h, [-1, 0, -1]); } static void blackRelax(double [3d] field, double [3d] rhs, double [3d] factor, RectDomain<3> [1d][1d] boundaries, double h) { pointRelax(field, rhs, factor, h*h, [-1, 0, 0]); pointRelax(field, rhs, factor, h*h, [-1, -1, -1]); pointRelax(field, rhs, factor, h*h, [ 0, -1, 0]); pointRelax(field, rhs, factor, h*h, [ 0, 0, -1]); } static void pointRelax(double [3d] fieldx, double [3d] rhsx, double [3d] factorx, double hsquared, Point<3> offset) { double [3d] field = fieldx.translate(offset); double [3d] rhs = rhsx.translate(offset); double [3d] factor = factorx.translate(offset); if (AMR.dumping && ((AMR.debug & 16) == 16)) { String s = "pointRelax() h=" + Math.sqrt(hsquared) + " "; Dump.dump(fieldx, s + "field " + Util.stringify(field.domain())); Dump.dump(rhsx, s + "rhs " + Util.stringify(rhs.domain())); Dump.dump(factorx, s + "factor " + Util.stringify(factor.domain())); } foreach (p in rhsx.domain() / [2, 2, 2]) { field[p * 2] += factor[p * 2] * ( field[p * 2 + [ 0, 0, 1]] + field[p * 2 + [ 0, 0,-1]] + field[p * 2 + [ 0, 1, 0]] + field[p * 2 + [ 0,-1, 0]] + field[p * 2 + [ 1, 0, 0]] + field[p * 2 + [-1, 0, 0]] - 6 * field[p * 2] - hsquared * rhs[p * 2]); } } static void xxx(double [3d] a) { foreach (p in a.domain()) { Util.printP3(p); System.out.println(a[p]); } } static void usage () { System.err.println( "Usage:\n" + "\n" + "amr [-adgitv] [-l #] [-D #] {input file} [{number of iterations}]\n" + "\n" + "All arguments except for the input file may be omitted.\n" + "\n" + "Setting the astro flag (-a) makes the program read\n" + "its input in astro format rather than `regular' format.\n" + "See inputs/sample.input for a regular example, and\n" + "inputs/sample_astro.input for an astro example.\n" + "\n" + "With the dumping flag set (-d), the program outpus\n" + "the residual and phi after each iteration into a file for\n" + "debugging.\n" + "\n" + "The show gradient flag (-g) specifies that the gradient of phi\n" + "should be dumped to stdout after each iteration.\n" + "\n" + "With the interactivity flag set (-i), the program\n" + "pauses after each iteration and outputs\n" + "residual and phi into a file for visualization.\n" + "\n" + "Setting the testing flag (-t) makes the program\n" + "evaluate the laplacian of a built-in function\n" + "(see AMRProcess.test()) and output the\n" + "result.\n" + "\n" + "Setting the verbose flag (-v) causes a whole lot\n" + "of information to be printed during\n" + "execution.\n" + "\n" + "The -l flag sets the level of logging requested. Default is 0.\n" + "\n" + "The -D flag sets the level of debug requested. Default is 0.\n" ); } public single static void main( String[] argv ) { log = new ExecutionLog(); barrierCount = 0; int _mainLoopIterations; boolean _oneshot; if (Ti.thisProc() == 0) System.out.print("Processes: "); AMR.barrier(); System.out.println(Ti.thisProc()); AMR.barrier(); if (Ti.thisProc() == 0) System.out.println("."); AMR.barrier(); // Parsing of command line arguments and // initialization of single AMR static variables. // See usage() for explanation of the command line options. boolean _interactive = false; boolean _die = false; boolean _verbose = false; boolean _dumping = false; boolean _testing = false; boolean _astroInput = false; boolean _showGradPhi = false; int _log = 0; int _debug = 0; String filename; // Note this is only set on proc 0. if (Ti.thisProc() == 0) try { int i; String s = new String(); while (argv[i].charAt(0) == '-') { String g = argv[i].substring(1); if (g.charAt(0) == 'l') { if (g.length() != 1) throw new IOException(); _log = Util.parseNthInt(argv, i + 1); i++; } else if (g.charAt(0) == 'D') { if (g.length() != 1) throw new IOException(); _debug = Util.parseNthInt(argv, i + 1); i++; } else { s += g; } i++; } // System.out.println(s); for (int j = s.length(); --j >= 0; ) { switch (s.charAt(j)) { case 'a': _astroInput = true; break; case 'g': _showGradPhi = true; break; case 'i': _interactive = true; break; case 'd': _dumping = true; break; case 'v': _verbose = true; break; case 't': _testing = true; break; case 'l': throw new IOException(); // -l handled above default: System.err.println("Ignoring bogus option (" + s.charAt(j) + ")"); } } filename = argv[i]; // System.out.println("filename = " + filename); if (argv.length >= i + 2) { _oneshot = false; _mainLoopIterations = Integer.parseInt(argv[i + 1]); System.out.println("iters = " + _mainLoopIterations); } else _oneshot = true; } catch (Throwable x) { usage(); _die = true; } if (broadcast _die from 0) return; astroInput = broadcast _astroInput from 0; dumping = broadcast _dumping from 0; showGradPhi = broadcast _showGradPhi from 0; oneshot = broadcast _oneshot from 0; debug = broadcast _debug from 0; interactive = broadcast _interactive from 0; verbose = broadcast _verbose from 0; testing = broadcast _testing from 0; mainLoopIterations = broadcast _mainLoopIterations from 0; log.loggingLevel = broadcast _log from 0; if (Ti.thisProc() == 0 && dumping) Dump.init(); long ms = System.currentTimeMillis(); maxLevels = 20; rhoNorms = new double [[0 : maxLevels - 1]]; AMRProcess single p = new AMRProcess(); // Make every process know about every other process. processes = (AMRProcess [1d] single) new AMRProcess[[0 : Ti.numProcs() - 1]]; processes.exchange(p); // Initialize process. try { p.setup(filename); } catch (Throwable x) { _die = true; } if (!(broadcast _die from 0)) { if (Ti.thisProc() == 0) { // Sequential initialization code. // The problem domain, in real-world coordinates, starts // at (0,0) and ends at (1,1). N can change, but the // problem domain doesn't. // Initialize Multigrid. Multigrid.init(); } AMR.barrier(); int single regriddings = 0; if (testing) { p.test(); return; } if (oneshot) { p.mainLoop(1, false); ExecutionLog.end(log); return; } for (;;) { if (p.mainLoop(mainLoopIterations, true)) break; regriddings++; if (regriddings >= maxRegriddings) break; p.regridAsNeeded(); } ExecutionLog.end(log); ms = System.currentTimeMillis() - ms; System.out.println(((double) ms) / 1000.0); System.out.println("Barriers: " + barrierCount); Ti.barrier(); ExecutionLog.histDump(); } } }