// The Multigrid class. Used at the coarsest level. All methods are // sequential. 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 Multigrid { static int levels; static double [1d][3d] rhos; static double [1d][3d] errors; static double [1d][3d] residuals; static double [1d][3d] factors; static RectDomain<3> [1d][1d][1d] boundaries; static Point<3>[1d][1d] redblack; static void init() { redblack = new Point<3>[[0 : 1]][1d]; redblack[0] = new Point<3>[[0 : 3]]; redblack[1] = new Point<3>[[0 : 3]]; redblack[0][0] = [0, 0, 0]; redblack[0][1] = [1, 1, 0]; redblack[0][2] = [0, 1, 1]; redblack[0][3] = [1, 0, 1]; redblack[1][0] = [1, 0, 0]; redblack[1][1] = [0, 1, 0]; redblack[1][2] = [0, 0, 1]; redblack[1][3] = [1, 1, 1]; levels = AMR.logProblemSize + 1; RectDomain<1> dlev = [0 : levels - 1]; rhos = new double[dlev][3d]; errors = new double[dlev][3d]; residuals = new double[dlev][3d]; factors = new double[dlev][3d]; boundaries = new RectDomain<3>[dlev][1d][1d]; int n = AMR.problemSize; foreach (i in dlev) { RectDomain<3> box = [[0, 0, 0] : [n - 1, n - 1, n - 1]]; rhos[i] = new double[box]; errors[i] = new double[[[-1, -1, -1] : [n, n, n]]]; residuals[i] = new double[box]; factors[i] = new double[box]; RectDomain<3> sbox = [[1, 1, 1] : [n - 2, n - 2, n - 2]]; double [3d] fff = factors[i]; foreach (p in sbox) { int x = 1; fff[p] = x / 6.0; } foreach (p in box - sbox) { int x = 1; fff[p] = x / 6.0; } boundaries[i] = new RectDomain<3>[[1 : 3]][1d]; foreach (dir in [1 : 3]) { boundaries[i][dir] = new RectDomain<3>[[-1 : 1 : 2]]; boundaries[i][dir][-1] = box.border(1, -dir[1], 1); boundaries[i][dir][1] = box.border(1, dir[1], 1); } n /= 2; } } static void relax(double [3d] field, double [3d] rhs, double [3d] residual, double [3d] factor, double h, int level) { // This is backward from AMR. Coarser levels have higher // numbers. int j = field.domain().max()[1]; String phase = "Multigrid " + j + "x" + j + "x" + j; AMR.log.begin(phase); if (AMR.verbose) System.out.println("Multigrid.relax"); if (level < levels - 1) { int coarse = level + 1; double [3d] coarseField = errors[coarse]; double [3d] coarseRho = rhos[coarse]; RectDomain<3> coarseInterior = coarseRho.domain(); // Relax on this level. AMR.gsrbRelax(field, rhs, factor, Multigrid.boundaries[level], h); // Initialize solution at next level. coarseField.set(0); // Compute fine residual. First, set dirichlet boundary // conditions like in Patch.physicalBC(). AMR.setBC(field, Multigrid.boundaries[level]); foreach (p in residual.domain()) { // 9 of 9 residual[p] = rhs[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]) / (h * h); } // Average fine residual onto coarse rhs. foreach (p in coarseInterior) { // coarseRho[p]= 0; // foreach (q in AMR.refinedCell) { // unroll // coarseRho[p] += residual[p * AMR.nRefineP + q]; // coarseRho[p] += residual[p * [2, 2, 2] + q]; // } // coarseRho[p] /= AMR.refinedCell.size(); coarseRho[p] = ( residual[p * [2, 2, 2] + [0, 0, 0]] + residual[p * [2, 2, 2] + [0, 0, 1]] + residual[p * [2, 2, 2] + [0, 1, 0]] + residual[p * [2, 2, 2] + [0, 1, 1]] + residual[p * [2, 2, 2] + [1, 0, 0]] + residual[p * [2, 2, 2] + [1, 0, 1]] + residual[p * [2, 2, 2] + [1, 1, 0]] + residual[p * [2, 2, 2] + [1, 1, 1]] ) / 8; } // Recursive call. relax(coarseField, coarseRho, residuals[coarse],factors[coarse], h * 2, coarse); // Correct fine solution by interpolating coarse // correction. // foreach (p in coarseInterior) { // foreach (q in AMR.refinedCell) { // unroll // field[p * AMR.nRefineP + q] += coarseField[p]; // } // } foreach (p in coarseInterior) { field[p * [2, 2, 2] + [0, 0, 0]] += coarseField[p]; field[p * [2, 2, 2] + [0, 0, 1]] += coarseField[p]; field[p * [2, 2, 2] + [0, 1, 0]] += coarseField[p]; field[p * [2, 2, 2] + [0, 1, 1]] += coarseField[p]; field[p * [2, 2, 2] + [1, 0, 0]] += coarseField[p]; field[p * [2, 2, 2] + [1, 0, 1]] += coarseField[p]; field[p * [2, 2, 2] + [1, 1, 0]] += coarseField[p]; field[p * [2, 2, 2] + [1, 1, 1]] += coarseField[p]; } // Relax again. AMR.gsrbRelax(field, rhs, factor, Multigrid.boundaries[level], h); } else { relax_1x1(field, rhs, h); } AMR.log.end(phase); } static void relax_1x1(double [3d] field, double [3d] rhs, double h) { field[[0, 0, 0]] += 0.5 * Math.pow(h, 2) * rhs[[0, 0, 0]]; } }