public class Patch2 { private static final int dimensions = 2; private static final RectDomain<1> Ddims = [1:dimensions], Ddirs = [-1:1:2]; public static Patch2 [1d] NoPatches = new Patch2[[0:-1]]; // coefficient of artificial viscosity private static final double viscosity = 0.1; // minimum relative difference to tag a cell private static final double MinRelDiff = 0.1; /* non-constant fields */ public Level level; // the level to which this patch belongs public RectDomain<2> domain; // domain of q and U public RectDomain<2> coarseDomain; // coarse version of same domain public RectDomain<2> coarseCoreDomain; // does not add ghost cells // Coarser patches that partially cover the patch. // Allow for different arrays of patches at old and new times. public Patch2 [1d] coarserPatchesOld; public Patch2 [1d] coarserPatchesNew; // Adjacent patches are those at the same level (other than itself) // that contain the patch's ghost cells. public Patch2 [1d] adjacentPatches; // Flux registers // public Vector4 [1d][1d][1d][2d] fluxRegisters; public FluxRegisterSegment [1d][1d][1d] frs; // [dimension 1:2][direction -1:1:2][segment] // direction refers to side of the finer patch. // Flux register will change U[e + Amr2.face[dimension][direction]] // The problem domain is a box, so a boundary may not be only // partially physical. public boolean [1d][1d] isPhysical; // [dimension 1:2][direction -1:1:2] // = new boolean[Ddims][Ddirs] ? public Quantities2 [2d] q; // primitive quantities public Vector4 [2d] U; // conserved quantities public Domain<2> tagged; /* constructors */ public Patch2(Level lev) { level = lev; } public Patch2(Level lev, RectDomain<2> D) { level = lev; domain = D; coarseDomain = domain / Amr2.nRefineP; coarseCoreDomain = coarseDomain; isPhysicalSet(); } public Patch2(Level lev, RectDomain<2> D, InitialPatch2 patchInit) { // initializes domain of a patch. level = lev; domain = D; coarseDomain = domain / Amr2.nRefineP; coarseCoreDomain = coarseDomain; // There are no coarser patches, because there is only one level. coarserPatchesOld = NoPatches; coarserPatchesNew = NoPatches; isPhysicalSet(); /* System.out.println("Process " + Ti.thisProc() + " height " + level.height + " domain " + Format.toString(domain) + " " + isPhysical[1][-1] + isPhysical[1][+1] + isPhysical[2][-1] + isPhysical[2][+1]); */ q = patchInit.qInitial(domain, level.ds); U = Quantities2.QtoU(q); } /* methods */ public final void fillData(InitialPatch2 patchInit) { q = patchInit.qInitial(domain, level.ds); U = Quantities2.QtoU(q); // System.out.println("Filled data for patch " + domain + " at height " + // level.height); } public final void copyNoFields(Patch2 patch) { // copies everything from patch, except q and U and frs. level = patch.level; domain = patch.domain; coarseDomain = patch.coarseDomain; coarseCoreDomain = patch.coarseCoreDomain; coarserPatchesOld = new Patch2[patch.coarserPatchesOld.domain()]; coarserPatchesOld.copy(patch.coarserPatchesOld); coarserPatchesNew = new Patch2[patch.coarserPatchesNew.domain()]; coarserPatchesNew.copy(patch.coarserPatchesNew); adjacentPatches = new Patch2[patch.adjacentPatches.domain()]; adjacentPatches.copy(patch.adjacentPatches); isPhysical = new boolean[Ddims][Ddirs]; foreach (pd in Ddims) isPhysical[pd].copy(patch.isPhysical[pd]); } public final void copyAllFields(Quantities2 [2d] qNew, Vector4 [2d] UNew) { // copy all of the fields q and U q = new Quantities2[qNew.domain()]; q.copy(qNew); U = new Vector4[UNew.domain()]; U.copy(UNew); } public final void copySomeFields(Quantities2 [2d] qNew, Vector4 [2d] UNew) { // copy some parts of the fields q and U q.copy(qNew); U.copy(UNew); } public final void copy(Patch2 patch) { copyNoFields(patch); copyAllFields(patch.q, patch.U); } public void isPhysicalSet() { // Sets isPhysical. Assumes that level and domain are already set isPhysical = new boolean [Ddims][Ddirs]; RectDomain<2> wholeDomain = level.domain; foreach (pdim in Ddims) { int ipdim = pdim[1]; isPhysical[pdim][-1] = (domain.min()[ipdim] == wholeDomain.min()[ipdim]); isPhysical[pdim][+1] = (domain.max()[ipdim] == wholeDomain.max()[ipdim]); } } public boolean onBoundary(Point<2> cell, int dim, int dir) { // returns true iff cell is on the boundary side [dim][dir] return domain.border(1, dim * dir, 0).contains(cell); } public boolean onBoundary(Point<2> cell, Point<1> pdim, Point<1> pdir) { return onBoundary(cell, pdim[1], pdir[1]); } public boolean onBoundary(Point<2> cell) { // returns true iff cell is on ANY boundary side boolean retval = false; foreach (pdim in Ddims) foreach (pdir in Ddirs) retval = retval || onBoundary(cell, pdim, pdir); return retval; } public boolean onPhysicalBoundary(Point<2> cell, int dim, int dir) { // returns true iff cell is on physical boundary side [dim][dir] return (isPhysical[dim][dir] && onBoundary(cell, dim, dir)); } public boolean onPhysicalBoundary(Point<2> cell, Point<1> pdim, Point<1> pdir) { return onPhysicalBoundary(cell, pdim[1], pdir[1]); } public boolean onPhysicalBoundary(Point<2> cell) { // returns true iff cell is on ANY physical boundary side boolean retval = false; foreach (pdim in Ddims) foreach (pdir in Ddirs) retval = retval || onPhysicalBoundary(cell, pdim, pdir); return retval; } public static Patch2 findCell(Patch2 [1d] patchArray, Point<2> cell) { // Returns the patch in patchArray whose grid contains cell. // Error if not found. Patch2 patch = null; foreach (ind in patchArray.domain()) { if (patchArray[ind].domain.contains(cell)) { patch = patchArray[ind]; break; } } if (patch == null) { // Could not find cell in any of the patches. System.out.println("Patch2.findCell: cell " + cell + " not found in any of the " + patchArray.domain().size() + " patches given:"); foreach (ind in patchArray.domain()) System.out.println("-> " + patchArray[ind].domain); } return patch; } public final void LevelGodunov2(Patch2 patch, double wtOld, double wtNew, double time, double dt, int dimfirst) { // The Patch2 "patch" contains data at the old timestep. domain = patch.domain; // to be expanded with ghost cells coarseDomain = patch.coarseDomain; // ditto coarseCoreDomain = patch.coarseDomain; // do not expand // Find flux registers for patchTemp. // If this is the finest level, then set to null. // Otherwise, flux registers go where a finer patch has an // interface with patch. // Do you know anything about the finer level? // You have to take the old version. // Another idea: loop over all finer patches, find their coarser patches // and set flux registers there if necessary. // Problem with this idea: you can't assign all flux registers on // a level at once, the way we've organized this. isPhysicalSet(); getFluxRegisters(); // needs level and domain and isPhysical // System.out.println("got flux registers"); // Add two cells in each direction unless stopped by physical boundary. foreach (pdim in Ddims) foreach (pdir in Ddirs) { if (!patch.isPhysical[pdim][pdir]) { int side = pdir[1] * pdim[1]; domain = domain.accrete(2, side); coarseDomain = coarseDomain.accrete(1, side); } } // System.out.println("added ghost cells"); coarserPatchesOld = getCoarserPatches(patch.coarserPatchesOld, coarseDomain); // System.out.println("got coarserPatchesOld"); coarserPatchesNew = getCoarserPatches(patch.coarserPatchesNew, coarseDomain); // System.out.println("got coarserPatchesNew"); // adjacentPatches is not used anyway. // Leave it undefined? isPhysicalSet(); // System.out.println("got isPhysical"); // Need patchTemp.q and patchTemp.U to be new, because // you don't want to mess with patch.q and patch.U q = new Quantities2[domain]; U = new Vector4[domain]; copySomeFields(patch.q, patch.U); // System.out.println("copied q and U"); // Logger.append("LevelGodunov filling in ghosts"); // Keep track of which coarse cells contain unfilled fine cells. Domain<2> unfilledC = coarseDomain - patch.coarseDomain; // Fill in data from adjacent patches at the fine level. // The method fills in as much as possible of patchTemp // and returns the new domain of unfilled coarse cells. unfilledC = FillFine(unfilledC, patch.adjacentPatches); // System.out.println("did FillFine"); /* Now fill the still-unfilled ghost cells of the patch with interpolated data from the coarse grid. */ FillInterp(unfilledC, wtOld, wtNew); // System.out.println("did FillInterp on domain " + domain); /* Advance one timestep. */ SolveStep2(dimfirst, dt, time); // System.out.println("did SolveStep2"); /* Restrict domain again. (Remove the ghost cells.) */ domain = patch.domain; coarseDomain = coarseCoreDomain; // Reset coarser patches, eliminating the ones that cover only // ghost cells. coarserPatchesOld = new Patch2[patch.coarserPatchesOld.domain()]; coarserPatchesOld.copy(patch.coarserPatchesOld); coarserPatchesNew = new Patch2[patch.coarserPatchesNew.domain()]; coarserPatchesNew.copy(patch.coarserPatchesNew); // Re-check physical boundary conditions. isPhysicalSet(); /* Copy patch, omitting ghost cells. */ //System.out.println("resetting levelOut.patches[" + ind[1] + "]"); // copy domain, coarseDomain, coarser patches //levelOut.patches[ind].copyNoFields(patch); // Now the only field that hasn't been set is adjacentPatches, // which is found after return from this method. copyAllFields(q.restrict(patch.domain), U.restrict(patch.domain)); // System.out.println("did LevelGodunov2"); } public final static Patch2 [1d] getCoarserPatches(Patch2 [1d] coarserPatches, RectDomain<2> coarseDomain) { // Return list of coarser patches covering coarseDomain, // starting with the list in coarserPatches and then searching // through their neighbors to find anything else uncovered. // See Level.FindCoarserOld(). // To get the list of coarser patches covering coarseDomain, // make a list. // First include the input coarserPatches. // Then, if necessary, look through the adjacent patches of // the coarser patches. BoxedList_Patch2 coarseList = new BoxedList_Patch2(); // Keep track of what part of the coarse domain is not yet covered // by coarse patches in the list. // You need this because you will be seeing the same coarse // patches more than once, going through the adjacent patches. Domain<2> uncoveredC = coarseDomain; foreach (indC in coarserPatches.domain()) { Patch2 patchC = coarserPatches[indC]; coarseList.push(patchC); uncoveredC = uncoveredC - patchC.domain; } if (!uncoveredC.isNull()) { searchAdjacentC: foreach (indC in coarserPatches.domain()) { Patch2 patchC = coarserPatches[indC]; foreach (indAdjC in patchC.adjacentPatches.domain()) { Patch2 patchAdjC = patchC.adjacentPatches[indAdjC]; Domain<2> intersect = (Domain<2>) patchAdjC.domain * uncoveredC; if (!intersect.isNull()) { coarseList.push(patchAdjC); uncoveredC = uncoveredC - patchAdjC.domain; if (uncoveredC.isNull()) break searchAdjacentC; } } } } return coarseList.toArray(); } // These are used in getFluxRegisters, which is not called recursively. private static Domain<2> [1d][1d] fluxEdges = new Domain<2>[Ddims][Ddirs]; private static Domain<2> [1d][1d] commonBoundary = new Domain<2>[Ddims][Ddirs]; private static RectDomain<2> [1d][1d] myBoundary = new RectDomain<2>[Ddims][Ddirs]; private static RectDomain<2> [1d][1d] interiorEdges = new RectDomain<2>[Ddims][Ddirs]; public final void getFluxRegisters() { /* Get flux registers for this patch. A flux register exists where a cell borders a patch at the next finer level. Need: * domain of this patch; * domain of all patches at next finer level */ frs = new FluxRegisterSegment[Ddims][Ddirs][1d]; if (level.height == Amr2.finestLevel) { // Patch is at the finest level: set flux registers to null. foreach (pdim in Ddims) foreach (pdir in Ddirs) frs[pdim][pdir] = new FluxRegisterSegment[0:-1]; } else { /* Find all edges in this patch in each dimension. Flux registers occur only at intersection of these edges with edges of finer-level patches (using coarseDomain). */ foreach (pdim in Ddims) { Point<2> unit = Point<2>.direction(pdim[1]); foreach (pdir in Ddirs) { // Include the boundary edges if and only if // (not physical) and (different side) myBoundary[pdim][pdir] = domain.border(1, pdim[1]*pdir[1], 0) + Amr2.edge[pdim][pdir]; Point<2> interiorEdgesMin = domain.min() + Amr2.edge[pdim][+1]; if (!isPhysical[pdim][-1] && pdir == [+1]) interiorEdgesMin = interiorEdgesMin - unit; Point<2> interiorEdgesMax = domain.max() + Amr2.edge[pdim][-1]; if (!isPhysical[pdim][+1] && pdir == [-1]) interiorEdgesMax = interiorEdgesMax + unit; interiorEdges[pdim][pdir] = [interiorEdgesMin : interiorEdgesMax]; } } /* fluxEdges[pdim][pdir] will contain all possible edges where flux registers can occur. commonBoundary[pdim][pdir] will contain all edges that lie on side [pdim][pdir] of the patch and also lie on boundaries of finer patches. These should be removed from fluxEdges[pdim][-pdir]. */ // Initialize to the empty domains. foreach (pdim in Ddims) foreach (pdir in Ddirs) { fluxEdges[pdim][pdir] = [[0, 0] : [-1, -1]]; commonBoundary[pdim][pdir] = [[0, 0] : [-1, -1]]; } // Loop over ALL patches at the next finer level. foreach (proc in Amr2.processes.domain()) { Level lev = Amr2.processes[proc].levels[level.height+1]; foreach (indF in lev.indices) { Patch2 patchF = lev.patches[indF]; foreach (pdim in Ddims) foreach (pdir in Ddirs) { RectDomain<2> patchFboundary = patchF.coarseDomain.border(1, pdim[1] * pdir[1], 0) + Amr2.edge[pdim][pdir]; Domain<2> boundaryI = interiorEdges[pdim][pdir] * patchFboundary; fluxEdges[pdim][pdir] = fluxEdges[pdim][pdir] + boundaryI; Domain<2> boundaryC = patchFboundary * myBoundary[pdim][pdir]; commonBoundary[pdim][pdir] = commonBoundary[pdim][pdir]+boundaryC; } } } // Remove edges that are boundaries of other fine patches. foreach (pdim in Ddims) { Domain<2> fluxEdgesBothSides = fluxEdges[pdim][+1] * fluxEdges[pdim][-1]; foreach (pdir in Ddirs) fluxEdges[pdim][pdir] = fluxEdges[pdim][pdir] - fluxEdgesBothSides - commonBoundary[pdim][-pdir[1]]; } /* Set up flux register segments. */ foreach (pdim in Ddims) foreach (pdir in Ddirs) { // Split up fluxEdges[pdim][pdir] into rectangles. // Each rectangle should have width 1, because any pair of non-adjacent // fine patches are separated by at least two coarse cells. RectDomain<2> [1d] segments = Fixme2.RectDomainList(fluxEdges[pdim][pdir]); RectDomain<1> segmentIndices = segments.domain(); frs[pdim][pdir] = new FluxRegisterSegment[segmentIndices]; foreach (indSegment in segmentIndices) { RectDomain<2> segment = segments[indSegment]; frs[pdim][pdir][indSegment] = new FluxRegisterSegment(segment); /* System.out.println("Process " + Ti.thisProc() + " height " + level.height + " patch " + Format.toString(domain) + " side " + (pdim[1] * pdir[1]) + " segment " + indSegment[1] + " is " + Format.toString(segment, pdim[1])); */ } } } } public final Domain<2> FillFine(Domain<2> unfilledCin, Patch2 [1d] adjacentPatches) { // Fill the patch's cells in unfilledCin (indexed at coarse level) // with data from patches in adjacentPatches. Domain<2> unfilledC = unfilledCin; foreach (indAdj in adjacentPatches.domain()) { Patch2 patchAdj = adjacentPatches[indAdj]; Domain<2> intersect = patchAdj.domain * domain; // Recall that grids at the same level do not overlap. // Hence nothing is copied twice. if (!intersect.isNull()) { copySomeFields(patchAdj.q, patchAdj.U); // Boundaries of patches must coincide with boundaries of // coarse cells. Hence by copying from a patch, we have just // filled up all the subcells of some coarse cells. unfilledC = unfilledC - (patchAdj.coarseDomain * coarseDomain); if (unfilledC.isNull()) break; } } return unfilledC; } public void FillInterp(Domain<2> unfilledC, double wtOld, double wtNew) { /* System.out.println("FillInterp: " + unfilledC.size() + " unfilled cells in patch " + domain + " at height " + level.height); */ foreach (cellC in unfilledC) { // Find the old and new coarse grid patches that contain cellC. Patch2 patchOldC = findCell(coarserPatchesOld, cellC); Patch2 patchNewC = findCell(coarserPatchesNew, cellC); // Get conserved quantities at cellC at the current time. /* if (level.height == 2) { System.out.println("found patchOldC, patchNewC for " + cellC + " in " + domain); System.out.println("patchOldC.U" + cellC + " domain is " + patchOldC.U.domain()); } */ Vector4 UCold = patchOldC.U[cellC]; Vector4 UCnew = patchNewC.U[cellC]; Vector4 UC = UCold * wtOld + UCnew * wtNew; // Vector4 UC = patchOldC.U[cellC] * wtOld + patchNewC.U[cellC] * wtNew; // if (level.height == 2) System.out.println("found conserved quantities for " + cellC + " in " + domain); /* Find slopes in each dimension. */ Vector4 [1d] UCslopes = CoarseSlopes(patchOldC, patchNewC, cellC, UC, wtOld, wtNew); // if (level.height == 2) System.out.println("found slopes for " + cellC + " in " + domain); /* Fill in all ghost cells in cellC. Recall that by proper nesting, the boundaries of domain occur at coarse cell boundaries. Therefore cellC will not contain any cells in domain; it'll contain only the ghost cells. */ Domain<2> ghosts = Amr2.subcells(cellC) * domain; // System.out.println("Interpolating"); foreach (p in ghosts) { U[p] = UC; foreach (pd in Ddims) { int dim = pd[1]; U[p] = U[p] + UCslopes[pd] * (((double) p[dim] + 0.5) / Amr2.dRefine - ((double) cellC[dim] + 0.5)); } q[p] = Quantities2.UtoQ(U[p]); } // if (level.height == 2) System.out.println("interpolated for " + cellC + " in " + domain); } // end cellC in unfilledC } // These can be globals because CoarseSlopes doesn't call itself. private static Vector4 [1d] UCslopes = new Vector4[Ddims]; private static Vector4 [1d] UCadj = new Vector4[Ddirs]; private static boolean [1d] foundAdjC = new boolean[Ddirs]; public final static Vector4 [1d] CoarseSlopes(Patch2 patchOldC, Patch2 patchNewC, Point<2> cellC, Vector4 UC, double wtOld, double wtNew) { // Check all directions for neighbors. foreach (pdim in Ddims) { foreach (pdir in Ddirs) { // pdir in {-1, 1} Point<2> adjC = cellC + Point<2>.direction(pdir[1] * pdim[1]); foundAdjC[pdir] = (patchOldC.domain.contains(adjC) && patchNewC.domain.contains(adjC)); if (foundAdjC[pdir]) { UCadj[pdir] = patchOldC.U[adjC] * wtOld + patchNewC.U[adjC] * wtNew; } else { // Look in adjacent patches. (Hope you've found them!) findAdj: foreach (indAdjOld in patchOldC.adjacentPatches.domain()) { Patch2 patchAdjOld = patchOldC.adjacentPatches[indAdjOld]; if (patchAdjOld.domain.contains(adjC)) { foreach (indAdjNew in patchNewC.adjacentPatches.domain()) { Patch2 patchAdjNew = patchNewC.adjacentPatches[indAdjNew]; if (patchAdjNew.domain.contains(adjC)) { foundAdjC[pdir] = true; UCadj[pdir] = patchAdjOld.U[adjC] * wtOld + patchAdjNew.U[adjC] * wtNew; break findAdj; } } } } } } // end pdir if (!foundAdjC[-1] && !foundAdjC[+1]) { // This shouldn't happen. UCslopes[pdim] = new Vector4(0.0, 0.0, 0.0, 0.0); } else if (foundAdjC[-1] && !foundAdjC[+1]) { UCslopes[pdim] = UC - UCadj[-1]; } else if (!foundAdjC[-1] && foundAdjC[+1]) { UCslopes[pdim] = UCadj[+1] - UC; } else { UCslopes[pdim] = Vector4.vanLeer(UCadj[-1], UC, UCadj[+1]); } } // end pdim return UCslopes; } // Used in SolveStep2 private static RectDomain<2> [1d] coarseCoreBoundary = new RectDomain<2>[Ddirs]; private static Point<2> [1d] edge = new Point<2>[Ddirs]; private static Point<2> [1d] face = new Point<2>[Ddirs]; public final void SolveStep2(int dimfirst, // dimension of first pass double dt, double time) { // Solver using 2nd-order slopes RectDomain<2> D = domain; Logger.append("Entering SolveStep2 level " + level.height + " patch " + "[" + D.min()[1] + "," + D.min()[2] + "]:" + "[" + D.max()[1] + "," + D.max()[2] + "]"); Quantities2 [1d][2d] qhext = new Quantities2[Ddirs][D]; Quantities2 [2d] qp = new Quantities2[D]; // D should contain Dh RectDomain<2> DhoutBounds = [D.min() + Amr2.edge[1][-1] + Amr2.edge[2][-1] : D.max() + Amr2.edge[1][+1] + Amr2.edge[2][+1]]; Vector4 [2d] flux = new Vector4[DhoutBounds]; Quantities2 [2d] deltaq = new Quantities2[D]; int dims[] = {dimfirst, 3 - dimfirst}; //System.out.println("Entering SolveStep2 with domain " + // Format.toString(D) + ", time=" + time + ", dt=" + dt); for (int idim = 0; idim < dimensions; idim++) { int dim = dims[idim]; foreach (pdir in Ddirs) { face[pdir] = Amr2.face[dim][pdir]; edge[pdir] = Amr2.edge[dim][pdir]; } // System.out.println("SolveStep2 pass dim=" + dim + " on patch " + // Format.toString(D) + " at level " + level.height); /* (1) Get slopes of primitive quantities. */ // System.out.println("(1) slopes"); // Logger.append("(1) Getting slopes in dim " + dim); Quantities2.Slopes2(q, dim, deltaq); /* (2) Extrapolate primitive quantities from centers to edges. */ // Logger.append("(2) Extrapolating in dim " + dim); // System.out.println("(2) extrapolation"); double dtdx = dt / level.ds[dim]; Quantities2.extrapolate(q, deltaq, dim, dtdx, qhext); /* (3) Solve Riemann problems at cell boundaries. */ // Logger.append("(3) Solving Riemann problems in dim " + dim); // System.out.println("(3) Riemann"); RectDomain<2> Dh = [D.min() + edge[+1] : D.max() + edge[-1]]; // space in qp is reused foreach (e in Dh) { qp[e] = Quantities2.Riemann0(dim, qhext[+1][e + face[-1]], qhext[-1][e + face[+1]]); if (qp[e].norm() > 1e+6 || qp[e].density() < 0.0) { System.out.println("qp[" + Format.toString(e) + "] problem: " + qp[e].density() + " " + qp[e].velocity(1) + " " + qp[e].velocity(2) + " " + qp[e].pressure()); System.out.println("problem from left at " + Format.toString(e + face[-1]) + ": " + qhext[+1][e+face[-1]].density()+" " + qhext[+1][e+face[-1]].velocity(1)+" "+ qhext[+1][e+face[-1]].velocity(2)+" "+ qhext[+1][e+face[-1]].pressure()); System.out.println("problem from right at " + Format.toString(e + face[+1]) + ": " + qhext[-1][e+face[+1]].density()+" " + qhext[-1][e+face[+1]].velocity(1)+" "+ qhext[-1][e+face[+1]].velocity(2)+" "+ qhext[-1][e+face[+1]].pressure()); } } /* (4) Conservative differencing of fluxes. */ // System.out.println("(4) flux"); // Logger.append("(4) Flux differencing in dim " + dim); RectDomain<2> Dhout = [D.min() + edge[-1] : D.max() + edge[+1]]; // space in flux is reused foreach (e in Dh) { flux[e] = qp[e].getFlux(dim); // artificial viscosity double vchange = q[e + face[-1]].velocityN(dim) - q[e + face[+1]].velocityN(dim); if (vchange > 0.0) flux[e] = flux[e] - (U[e + face[+1]] - U[e + face[-1]]) * vchange * viscosity; } // Get flux at boundaries, using user-supplied routine. // If the boundary is not physical, then this is not correct. // System.out.println("(4.) Flux at boundaries"); // Logger.append("(4.) Flux at boundaries"); foreach (pd in Ddirs) { RectDomain<2> borderSide = D.border(1, pd[1] * dim, 0); Amr2.boundaryFluxFn.f(qhext[pd].restrict(borderSide), level, pd[1] * dim, time, flux); } foreach (p in D) { U[p] = U[p] - ((flux[p + edge[+1]] - flux[p + edge[-1]]) * dtdx); // System.out.println("UtoQ"); q[p] = Quantities2.UtoQ(U[p]); } // System.out.println("(5) Filling in flux registers for this patch"); // Logger.append("(5) Filling in flux registers for this patch"); // Fill in the flux registers for this patch. FluxRegisterSegment [1d][1d] fluxRegistersDim = frs[dim]; foreach (pdir in Ddirs) { foreach (indSegment in fluxRegistersDim[pdir].domain()) { FluxRegisterSegment fluxSeg = fluxRegistersDim[pdir][indSegment]; foreach (e in fluxSeg.domain) fluxSeg.Set(e, flux[e] * dt, "level " + level.height + " coarse, side " + (pdir[1]*dim)); } } //System.out.println("(6) Updating flux registers for coarser patches."); //Logger.append("(6) Updating flux registers for coarser patches."); // Update the flux registers at the coarser level. // The flux registers will be at the boundary of the (finer) patch. // CAUTION: Avoid coarser patches that cover ghost cells of // the patch domain but not non-ghost cells. Is this a problem? // Not the way I have it set up, where you omit the ghost cells // when finding flux registers. foreach (pdir in Ddirs) coarseCoreBoundary[pdir] = coarseCoreDomain.border(1, dim*pdir[1], 0) + edge[pdir]; foreach (indpatchC in coarserPatchesNew.domain()) { fluxRegistersDim = coarserPatchesNew[indpatchC].frs[dim]; foreach (pdir in Ddirs) { foreach (indSegment in fluxRegistersDim[pdir].domain()) { FluxRegisterSegment fluxSeg = fluxRegistersDim[pdir][indSegment]; // We do not update flux registers in the ghost region. Domain<2> fluxEdges = fluxSeg.domain * coarseCoreBoundary[pdir]; foreach (e in fluxEdges) { foreach (esub in Amr2.subedges(e, dim)) fluxSeg.Set(e, fluxSeg[e] - flux[esub] * (dt / Amr2.dRefine), "level " + (level.height-1) + " fine, side " + (pdir[1]*dim)); } } } } } // System.out.println("Leaving SolveStep2"); Logger.append("Leaving SolveStep2 level " + level.height + " patch " + "[" + D.min()[1] + "," + D.min()[2] + "]:" + "[" + D.max()[1] + "," + D.max()[2] + "]"); } public final void FindAdjacent() { // Find patches at the same level that are adjacent to this one. // "Adjacent" means that its ghost cells overlap with the other patch. // Easier: expandedDomain = domain.accrete(2) * level.domain; RectDomain<2> expandedDomain = domain; foreach (pdim in Ddims) foreach (pdir in Ddirs) // Add two cells in each direction unless stopped by physical boundary. if (!isPhysical[pdim][pdir]) expandedDomain = expandedDomain.accrete(2, pdir[1] * pdim[1]); // System.out.println("setting up adjList"); BoxedList_Patch2 adjList = new BoxedList_Patch2(); // Loop over ALL patches at this level. foreach (proc in Amr2.processes.domain()) { Level lev = Amr2.processes[proc].levels[level.height]; foreach (indOther in lev.indices) {// loop again over all patches... Patch2 patchOther = lev.patches[indOther]; // ...other than this one if (patchOther.domain != domain) { Domain<2> intersect = expandedDomain * patchOther.domain; if (!intersect.isNull()) adjList.push(patchOther); } } } adjacentPatches = adjList.toArray(); } public final void AverageDown() { // Average solution at this level to get a new coarse solution // on all coarse cells (at next coarser level) covered by this level. foreach (indC in coarserPatchesNew.domain()) { // patches at coarser level Patch2 coarserPatch = coarserPatchesNew[indC]; Domain<2> intersect = coarserPatch.domain * coarseDomain; foreach (cellC in intersect) { coarserPatch.U[cellC] = Vector4.Average(U.restrict(Amr2.subcells(cellC))); coarserPatch.q[cellC] = Quantities2.UtoQ(coarserPatch.U[cellC]); } } } public final void Reflux() { // Correct the solution (U and q) at coarse-fine boundaries // due to flux discrepancy. foreach (pdim in Ddims) { double dx = level.ds[pdim]; foreach (pdir in Ddirs) { /* System.out.println("reflux for side " + (pdim[1] * pdir[1]) + " of level " + level.height + " patch " + Format.toString(domain)); */ foreach (indSegment in frs[pdim][pdir].domain()) { // System.out.println("segment " + indSegment[1]); FluxRegisterSegment fluxSegment = frs[pdim][pdir][indSegment]; // System.out.println("domain " + // Format.toString(fluxSegment.domain)); foreach (e in fluxSegment.domain) { Point<2> f = e + Amr2.face[pdim][pdir]; /* System.out.println("Reflux face " + Format.toString(e) + " from " + U[f].elem(1) + " " + U[f].elem(2) + " " + U[f].elem(3) + " " + U[f].elem(4)); */ U[f] = U[f] + fluxSegment[e] / dx; /* System.out.println("Reflux face " + Format.toString(e) + " to " + U[f].elem(1) + " " + U[f].elem(2) + " " + U[f].elem(3) + " " + U[f].elem(4)); */ q[f] = Quantities2.UtoQ(U[f]); // q[f].thresholded(); -- this is part of UtoQ now // U[f] = q[f].QtoU(); -- should we use this?? } } } } } public final void TagPoints() { // Tag appropriate cells of domain. // Needs q to be defined. // First find points where relative difference in density // in some direction exceeds MinRelDiff. boolean [2d] bigdiff = new boolean[domain]; foreach (p in domain) { double dens0 = q[p].density(); foreach (pdim in Ddims) foreach (pdir in Ddirs) { Point<2> pnew = p + Point<2>.direction(pdim[1] * pdir[1]); if (domain.contains(pnew)) { double dens1 = q[pnew].density(); double relDiff = 2.0 * Math.abs(dens0 - dens1) / (Math.abs(dens0) + Math.abs(dens1)); if (relDiff > MinRelDiff) bigdiff[p] = true; } } } // Now tag points of large relative difference, // as well as the neighbors of these points. // You may go beyond the domain, but that's OK, because you // clip tags outside the proper nesting domain anyway. tagged = [Point<2>.all(0) : Point<2>.all(-1)]; RectDomain<2> baseSquare = [Point<2>.all(-1) : Point<2>.all(+1)]; foreach (p in domain) if (bigdiff[p]) tagged = tagged + (baseSquare + p); } public final void FillRegrid() { q = new Quantities2[domain]; U = new Vector4[domain]; // ASSUME patch boundaries coincide with coarser cell boundaries. Domain<2> unfilledC = coarseDomain; /* First try to copy from old patches. */ // Loop over ALL patches at this level. foreach (proc in Amr2.processes.domain()) { Level oldLev = Amr2.processes[proc].oldLevels[level.height]; foreach (indOld in oldLev.indices) {// loop again over all patches... Patch2 patchOld = oldLev.patches[indOld]; RectDomain<2> intersect = patchOld.coarseDomain * coarseDomain; if (!intersect.isNull()) { copySomeFields(patchOld.q, patchOld.U); } } } /* Now fill in the rest by interpolation. */ foreach (cellC in unfilledC) { // Find the old and new coarse grid patches that contain cellC. Patch2 patchOldC = findCell(coarserPatchesOld, cellC); Vector4 UC = patchOldC.U[cellC]; /* Find slopes in each dimension. */ Vector4 [1d] UCslopes = CoarseSlopes(patchOldC, patchOldC, cellC, UC, 1.0, 0.0); /* Fill in all subcells of cellC. Recall that by proper nesting, the boundaries of domain occur at coarse cell boundaries. Therefore cellC will not contain any cells in domain; it'll contain only the ghost cells. */ // System.out.println("Interpolating"); foreach (p in Amr2.subcells(cellC) * domain) { U[p] = UC; foreach (pd in Ddims) { int dim = pd[1]; U[p] = U[p] + UCslopes[pd] * (((double) p[dim] + 0.5) / Amr2.dRefine - ((double) cellC[dim] + 0.5)); } q[p] = Quantities2.UtoQ(U[p]); } // if (level.height == 2) System.out.println("interpolated for " + cellC + " in " + domain); } // end cellC in unfilledC } }