import java.io.*;

// Per-process data structure.

public class Amr2Process {

  public static int dimensions = 2;
  private static final RectDomain<1> Ddims = [1:dimensions], Ddirs = [-1:1:2];

  /*
    non-constant field
  */

  Level single [1d] levels;
  Level single [1d] oldLevels; // used in regridding

  public Amr2Process(Region r) {
    levels = new (r) Level single[0:Amr2.finestLevelMax];
    oldLevels = new (r) Level single[0:Amr2.finestLevelMax];  // used in regridding
  }

  single public void GetPatches(InitialPatch2 patchInit) {
    // Get initial patches.
 
     /*
      Set up level 0.
      Process 0 has one patch at this level.
      Other processes have no patches at level 0.
     */
 
    int single height = 0;
    Amr2.finestLevel = height;
    
    RectDomain<2> D = Amr2.domain;
    Point<2> extent = D.max() - D.min() + Point<2>.all(1);
    // If D == [[0, 0] : [127, 31]] then extent == [128, 32].
    double [1d] ds = patchInit.dsGet(extent);
    SharedRegion single level0Region = new SharedRegion();
    levels[0] = new (level0Region) Level(level0Region, [0:0], D, height, ds);

    // Split domain in first dimension.
    Point<2> fracextent = extent / [Ti.numProcs(), 1];  // e.g. [32, 32]
    Point<2> shifter = [fracextent[1], 0];  // e.g. [32, 0]
    RectDomain<2> myD = [Point<2>.all(0) : fracextent - Point<2>.all(1)] +
      shifter * Ti.thisProc();
    levels[0].patches[0] = new (level0Region)
      Patch2(level0Region, levels[0], myD, patchInit);
	System.out.println("Level 0 patch " + levels[0].patches[0].domain);

    Level myLevel = levels[0];

    Ti.barrier();
    myLevel.FindAdjacent();
    myLevel.FindCoarserOld();
    Ti.barrier();
 
    while (height < Amr2.finestLevelMax) {
 
      /*
	Tag cells in all patches at this height.
      */

      // System.out.println("tagging cells at height " + height);

      boolean someTagged = false;
      foreach (indpatch in myLevel.indices) {
	Patch2 patch = myLevel.patches[indpatch];
	patch.TagPoints();
	if (!patch.tagged.isNull()) someTagged = true;
	if (someTagged)
	  System.out.println("Patch " + patch.domain + " has " + patch.tagged.size() + " tagged cells");
      }
 
      boolean single anyTagged = Reduce.or(someTagged);
 
       /*
	If any cells are tagged, then refine:  set up new level height + 1.
       */
 
      if (anyTagged) {

	//System.out.println("some are tagged.  calling meshrefine.refine");

	// Refine this level to get a new level height+1.

	MeshRefine.Refine(height, height);

	height++;
	// Amr2.finestLevel = height;  <- MeshRefine has already set this

	//System.out.println("done with meshrefine.refine.  setting new Level");
	// Perhaps MeshRefine could have a flag saying whether to use
	// initial data or to interpolate.
	myLevel = levels[height];

	if (Ti.thisProc() == 0)
	  System.out.println("********* Patches in domain " + myLevel.domain +
			     " at height " + height);
	Ti.barrier();

	foreach (indpatch in myLevel.indices) {
	  Patch2 patch = myLevel.patches[indpatch];
	  patch.fillData(level0Region, patchInit);
 
	  System.out.println("Process " + Ti.thisProc() + " has " + patch.domain.toString());
	}

	Ti.barrier();
	myLevel.FindAdjacent();
	myLevel.FindCoarserOld();
	Ti.barrier();

      } else {
	// No need to refine further.
	System.out.println("no tagged cells");
	break;  //  while (height < Amr2.finestLevelMax)
      }
    }
 
    // System.out.println("found all initial patches");
    /*
    // Fill in levels that have NO patches.  Is this necessary?
    // Probably not, if you have Amr2.finestLevel anyway.
    // Solution:  Define new Levels as you need them.
    for (height = Amr2.finestLevel; height <= Amr2.finestLevelMax; height++) {
      extent = extent * Amr2.nRefine;
      D = [Point<2>.all(0) : extent - Point<2>.all(1)];
      ds = patchInit.dsGet(extent);
      SharedRegion r = new SharedRegion();
      levels[height] = new (r) Level(r, [0:-1], D, height, ds);
    }
    */
      
   }

  single public double Process(int single s, double single totaltime,
			       double single starttime,
			       boolean single justregrid) {
     // main routine
     // s = max number of steps
     // totaltime = max time to run
    // regridintvl = number of steps between regriddings (I use 2)
    // plotintvl = number of steps between plot outputs
    // returns time
 
    double single time = starttime;
     int single dimfirst = 1;
    Point<2> domainExtent = Amr2.domain.max() - Amr2.domain.min() + Point<2>.all(1);
    String outputfilePrefix = "out";
    if (starttime > 0.0) outputfilePrefix = outputfilePrefix + 'r';
    outputfilePrefix = outputfilePrefix + '.' + 
      domainExtent[1] + '.' +
      domainExtent[2] + '.' + 
      Amr2.finestLevel + '.' +
      Ti.numProcs() + '.';

    for (int single step = 1;
	 step <= s || time < totaltime;
	 step++) {

      if (justregrid || step % Amr2.regridintvl == 0) {  // Regrid.
	int single finestLevelOld = Amr2.finestLevel;
	int single finestLevelNew = finestLevelOld +
	  ((finestLevelOld >= Amr2.finestLevelMax) ? 0 : 1);
	RectDomain<1> single newlevsh = [1 : finestLevelNew];

	// Tag points in patches at levels below finestLevelNew.
	foreach (hh in newlevsh - [1]) {
	  Level lev = levels[hh];
	  foreach (indpatch in lev.indices) lev.patches[indpatch].TagPoints();
	}

	// Save the levels that will be modified.
	// Level [1d] oldLevels = new Level[newlevsh];
	foreach (hh in newlevsh) oldLevels[hh] = levels[hh];

	// Get new grids.  This will change levels[newlevsh].
	MeshRefine.Refine(0, finestLevelNew - 1);

	// Interpolate data from old hierarchy onto new hierarchy.
	for (int single h = 1; h <= finestLevelNew; h++) {
	  Level myLevel = levels[h];
	  Ti.barrier();
	  myLevel.FindAdjacent();
	  myLevel.FindCoarserOld();
	  Ti.barrier();

	  // Now interpolate data.
	  myLevel.FillRegrid();
	}

	Ti.barrier();

	// Now delete the old hierarchy
	foreach (hh in newlevsh) {
	  Level single old = oldLevels[hh];
	  oldLevels[hh] = null;
	  Util.deleteRegion(old.getRegion());
	}

	if (justregrid) return time;
      }

      // Figure out timestep to use for level 0.
      double mydt = ComputeDt();

      // Now do a reduction over all processes to find smallest dt
      double single alldt = Reduce.min(mydt);
      double single dt = alldt;

      if (totaltime > 0.0) {
	// dt = Math.min(dt, totaltime - time);
	// I really need to annotate that standard library...
	if (totaltime - time < dt) dt = totaltime - time;
      }

      TimeStep(0, time, dt, 0, dimfirst);

      time += dt;
      dimfirst = 3 - dimfirst;
      if (Ti.thisProc() == 0) {
        System.out.println("Step " + step + ":  " +
                           "timestep " + dt + ",   total " + time);
	if (step % Amr2.plotintvl == 0)
	  WriteHierarchy.write(outputfilePrefix + step, time);
      }
    }
    return time;
  }
 
 

  public double ComputeDt() {
    // Compute the best dt for this process
    // See HSCL::estTimeStep() in appsrc/PDE/HSCL/HSCL/HSCL.cpp
    // CONS2PRM there converts from conserved to primitive quantities, and
    // MAXWAVESPEED there gives max of abs(q.velocity(dim)) + c.
    // But the HSCL routine assumes ds[1]==ds[2].
    double maxrat = 0.0;
    // Use coarsest level ds, because you're finding only coarsest level dt.
    double [1d] dsl = levels[0].ds;
    foreach (l in levels.domain()) {
      Level level = levels[l];
      foreach (ind in level.indices) {
	Patch2 patch = level.patches[ind];
	foreach (p in patch.domain) {
	  Quantities2 q = patch.q[p];
	  double c = q.soundspeed();
	  for (int dim = 1; dim <= dimensions; dim++) {
	    double rat = (Math.abs(q.velocity(dim)) + c) / dsl[dim];
	    maxrat = Math.max(maxrat, rat);
	  }
	}
      }
    }
    return (Amr2.maxcfl / maxrat);
  }

  single public final void
    TimeStep(int single levi, double time, double dt, int iStep, int dimfirst) {

    // Advance this level in time using the level operator.

    Logger.barrier();
    SharedRegion single levelRegion = new SharedRegion();

    // LevelGodunov2 calls FillInterp, which needs coarserPatchesNew if levi>0.

    // Can't call this way because of singleness issues.  Hmmm.
    //    Level newLevel = myLevel.LevelGodunov2(dt,time,iStep,dimfirst);
    Level single oldLevel = levels[levi];
    Level single newLevel = Level.LevelGodunov2(levelRegion, oldLevel,
						time, dt, iStep, dimfirst);

    Logger.barrier();
    levels[levi] = newLevel;
    Logger.barrier();

    //	System.out.println("Calling FindAdjacent");
    // Logger.append("Doing FindAdjacent");
    levels[levi].FindAdjacent();
    Logger.barrier();


    // advance the next finer level to time+dt
    if (levi < Amr2.finestLevel) {

      dt /= Amr2.dRefine;

      // Don't call TimeStep until you have filled in coarserPatchesNew
      // for levels[levi + 1].
      // Run through something like FindCoarserNew?
      // Go through processes[ALL].levels[levi] to search for
      // coarserPatchesNew?
      Logger.barrier();

      levels[levi + 1].FindCoarserNew();
      // The pointers to old and new coarser patches
      // will be copied as you progress nRefine steps.
      Logger.barrier();

      // At finer levels, start with dimension 1 if Amr2.nRefine is even.
      // Then end up with the other dimension after Amr2.nRefine steps.
      int dimfirstsub = 1;
      for (int single iStepNew = 0; iStepNew < Amr2.nRefine; iStepNew++) {
	// TimeStep will modify levels[levi + 1]
	TimeStep(levi + 1, time, dt, iStepNew, dimfirstsub);
        dimfirstsub = 3 - dimfirstsub;
	time += dt;
      }

      // Correct the solution at coarse-fine boundaries
      // due to flux discrepancy.
      // System.out.println("reflux");
      levels[levi].Reflux();

      // Average fine solution to coarse solution on all cells covered
      // by finer level.
      // Need coarserPatchesNew to be able to do AverageDown.
      // AverageDown changes levels[levi] but not levels[levi + 1].
      // System.out.println("averaging down");
      Ti.barrier();
      levels[levi + 1].AverageDown();
      Ti.barrier();
      // Reset coarserPatchesOld to coarserPatchesNew
      foreach (indpatch in levels[levi+1].indices) {
	Patch2 patch = levels[levi+1].patches[indpatch];
	patch.coarserPatchesOld = patch.coarserPatchesNew;
      }
      Ti.barrier();
    }
    // The old level is now no longer needed. Delete it.
    Util.deleteRegion(oldLevel.getRegion());
  }
}