import java.io.*;

public class Sod2dm {

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

  protected static double shockradius = 0.25;

  protected static double gamma = 1.4;
  protected static double smallthresh = 1e-6;

  protected static double rho0 = 0.125, p0 = 0.1;
  protected static double rho1 = 1.0, p1 = 1.0;

  protected static Quantities2 q0, q1;

  protected static Vector4 [1d] q0flux, q1flux;

  single public static void main(String args[]) {

    int single numProcs = Ti.numProcs();
    RectDomain<1> single Dprocs = [0 : numProcs-1];
    int myProc = Ti.thisProc();

    boolean myprintOut;
    String myinputfile, myoutputfile;
    int mys;
    double mytotaltime;

    /*
      Get command-line arguments:  inputfile, s, totaltime, outputfile
    */

    try {
      int myarglen = args.length;
      if (myarglen < 3) exitinfo();
      myinputfile = args[0];
      mys = Integer.parseInt(args[1]);
      mytotaltime = new Double(args[2]).doubleValue();
      myprintOut = (myarglen >= 4);
      if (myprintOut) myoutputfile = args[3];
    } catch (Exception e) {
      exitinfo();
    }

    int single s = broadcast mys from 0;
    double single totaltime = broadcast mytotaltime from 0;
    boolean single printOut = broadcast myprintOut from 0;

    /*
      Set up Amr2 with global parameters in myinputfile.
    */

    DataInputStream in =
      Amr2.SetupInput(myinputfile, new Sod2dmBoundaryFlux());

    /*
      Initialize processes.
    */

    System.out.println("initializing process " + myProc); 

    // The "single" modifier makes it look as if this sets all the
    // processes to be the same.
    Amr2Process single A2p = new Amr2Process();
    Amr2.processes = (Amr2Process [1d] single)
      new Amr2Process[[0 : numProcs - 1]];
    Amr2.processes.exchange(A2p);

    /*
      Set up levels in this process.
    */

    System.out.println("getting domains for " + myProc);
    RectDomain<2> [1d][1d] patchDomains = A2p.GetDomains(in);

    System.out.println("done with getting domains for " + myProc); 
    // For each level (and each process), construct a list of patches.
    // Or for each level, construct two lists of points (min and max).


    /*
      Find cell dimensions and step sizes.
    */

    // dimensions at level 0.
    Point<2> extent = Amr2.domain.max() - Amr2.domain.min() + [1, 1];

    double [1d][1d] ends = new double[Ddims][Ddirs];
    ends[1][-1] = 0.0;
    ends[1][+1] = 1.0;
    ends[2][-1] = 0.0;
    ends[2][+1] = 1.0;

    double [1d] ds = new double[Ddims];

    foreach (pdim in Ddims) ds[pdim] =
      (ends[pdim][+1] - ends[pdim][-1]) / (double) extent[pdim[1]];

    /*
      Set initial primitive quantities.
    */

    Quantities2.gammaSet(gamma);
    Quantities2.setThresholds(smallthresh, smallthresh);

    // recall order rho, ux, uy, p

    q0 = new Quantities2(rho0, 0.0, 0.0, p0);
    q1 = new Quantities2(rho1, 0.0, 0.0, p1);

    Quantities2 [2d] q;
    Vector4 [2d] U;

    /*
      These are used for finding boundary fluxes while solving.
    */

    q0flux = new Vector4[Ddims];
    q1flux = new Vector4[Ddims];
    foreach (pd in Ddims) {
      q0flux[pd] = q0.solidBndryFlux(pd[1]);
      q1flux[pd] = q1.getFlux(pd[1]);
    }

    /*
      Now set up levels and patches.
    */

    // double [1d] dslevel = new double[Ddims];
    // dslevel.copy(ds);
    RectDomain<2> D = Amr2.domain;
    Point<2> levExtent = extent;

    for (int single height = 0; height <= Amr2.finestLevel; height++) {

      D = [Point<2>.all(0) : levExtent - Point<2>.all(1)];
      RectDomain<2> [1d] patchesHeight = patchDomains[height];
      A2p.levels[height] = new Level(patchesHeight.domain(),
				     D,       // whole domain
				     height,  // height
				     ds);     // ds / 2^lev
      Level myLevel = A2p.levels[height];
      System.out.println("Process " + myProc + " height " + height +
			 " has domain " + Format.toString(D));

      /*
	Set up patches in each level.
      */

      foreach (indpatch in patchesHeight.domain()) {
	Patch2 patch = new Patch2();

	patch.level = myLevel;
	patch.domain = patchesHeight[indpatch];
	// The line below is always true.
	patch.coarseDomain = patch.domain / Amr2.nRefineP;
        patch.coarseCoreDomain = patch.coarseDomain;

	// There are no coarser patches, because there is only one level.
	patch.coarserPatchesOld = Patch2.NoPatches;
	patch.coarserPatchesNew = Patch2.NoPatches;

	// Avoid setting patch.adjacentPatches until other patch fields are
	// set for all patches at this level.

	patch.isPhysicalSet();
	System.out.println("Process " + myProc + " height " + height +
                           " patch " + indpatch[1] +
			   " domain " + Format.toString(patch.domain) + " " +
			   patch.isPhysical[1][-1] +
			   patch.isPhysical[1][+1] +
			   patch.isPhysical[2][-1] +
			   patch.isPhysical[2][+1]);

	patch.q = qInitial(q0, q1, patch.domain, ds, 10);
	patch.U = Quantities2.QtoU(patch.q);
	
	// Make the relevant level know about this patch.
	myLevel.patches[indpatch] = patch;
      }


      // System.out.println(myProc + " finding adjacent patches");
      Ti.barrier();
      myLevel.FindAdjacent();
      myLevel.FindCoarserOld();
      Ti.barrier();

      // Careful!  What ds will you end up with?  You'll be using ds later?
      foreach (pdim in Ddims) ds[pdim] /= Amr2.dRefine;

      levExtent = levExtent * Amr2.nRefine;
    }

    /*
      Process the patches.
    */

    // start logging
    Logger.begin();
    // log = new ExecutionLog(); barrierCount = 0;

    long startTime = System.currentTimeMillis();

    //    System.out.println(myProc + " calling Amr2Process.Process");
    A2p.Process(s, totaltime);
    //    System.out.println(myProc + " called Amr2Process.Process");

    //	System.out.println(myProc + " ending log");
    // ExecutionLog.end(log);
    Logger.end();

    //	System.out.println(myProc + " checking time");
    long endTime = System.currentTimeMillis();
    //	System.out.println(myProc + " barrier");
    Ti.barrier();

    /*
      Output
    */

    System.out.println("time = " + (endTime - startTime) / 1000 + " sec");

    if (myProc == 0 && printOut) try {

      Format LongE = new Format("%26.15e");
      // Open output file.
      System.out.println("Opening file " + myoutputfile + " for output");
      OutputStream fout = new FileOutputStream(myoutputfile);
      DataOutputStream out = new DataOutputStream(fout);

      System.out.println("Writing output");
      out.writeBytes(" " + extent[1] + " " + extent[2] + '\n');
      out.writeBytes(" " + Amr2.finestLevel + '\n');

      System.out.println("Going through all processes...");
      String outputline;
      Quantities2 qp;
      for (int height = 0; height <= Amr2.finestLevel; height++) {
	foreach (indproc in Amr2.processes.domain()) {
	  Amr2Process A2proc = Amr2.processes[indproc];
	  // All processes contain the same number of Levels,
	  // even though there may be no patches at a particular Level.
	  Level lev = A2proc.levels[height];
	  foreach (indpatch in lev.indices) {
            Patch2 patch = lev.patches[indpatch];
	    RectDomain<2> dom = patch.domain;
	    System.out.println("Writing patch with height " + height +
			       " in process " + indproc[1] +
			       " with domain " +
			       Format.toString(patch.domain));
	    outputline = " " + height + "     " +
	      dom.min()[1] + " " +
	      dom.min()[2] + "   " +
	      dom.max()[1] + " " +
	      dom.max()[2] + '\n';
	    out.writeBytes(outputline);
	    foreach (p in patch.domain) {
	      qp = patch.q[p];
	      outputline = LongE.form(qp.density()) +
		LongE.form(qp.velocity(1)) +
		LongE.form(qp.velocity(2)) +
		LongE.form(qp.pressure()) + '\n';
	      out.writeBytes(outputline);
	    }
	  }
	}
      }
      out.close();
    } catch (IOException x) {
      System.out.println("Error in writing " + myoutputfile);
      System.exit(0);
    }

  }


  private final static void exitinfo() {
    if (Ti.thisProc() == 0) {
      System.out.println("Usage:  ./Sod2dm [inputfile] [s] [t] [[outputfile]]");
      System.out.println("        [s] is number of time steps, or 0 if using [t]")
;
      System.out.println("        [t] is total time, or 0 if using [s]");
      System.exit(0);
    }
  }


  private final static Quantities2 [2d]
    qInitial(Quantities2 q0, Quantities2 q1,
	     RectDomain<2> Dlocal,
	     double [1d] ds,
	     int nsub) {

    double radius2 = shockradius * shockradius;
    double dnsub = (double) nsub;
    Quantities2 [2d] q = new Quantities2[Dlocal];
    double rn2 = 1.0 / (dnsub * dnsub);
    double [1d] subds = new double[Ddims];
    foreach (pdim in Ddims) subds[pdim] = ds[pdim] / dnsub;
    foreach (p in Dlocal) {
      double left = (double) p[1] * ds[1];
      double right = (double) (p[1] + 1) * ds[1];
      double bottom = (double) p[2] * ds[2];
      double top = (double) (p[2] + 1) * ds[2];
      if (left * left + bottom * bottom >= radius2) {
	q[p] = q0;
      } else if (right * right + top * top <= radius2) {
	q[p] = q1;
      } else {
	int iwts = 0;
	for (int jsub = 1; jsub <= nsub; jsub++) {
	  for (int ksub = 1; ksub <= nsub; ksub++) {
	    double xsub = left + ((double) jsub - 0.5) * subds[1];
	    double ysub = bottom + ((double) ksub - 0.5) * subds[2];
	    if (xsub * xsub + ysub * ysub <= radius2) iwts++;
	  }
	}
	double wt = (double) iwts * rn2;
	q[p] = q1 * wt + q0 * (1.0 - wt);
      }
    }
    return q;
  }

}