/*-----------------------------------------------------------------------*/
// initial.C
// 
//
/*-----------------------------------------------------------------------*/

#include "DAGH.h"
#include "waveamr2d.h"
#include "wavefortran2d.h"
#define VizServer 0
//
// Routine to take initial time step for wave equation given that the 
// initial time slice is defined on the t=0 step
//
// To take the first step in leapfrog we
// 1. Copy level 0 to level 1
// 2. Average levels 0 and 1 to obtain level 1/2.
// 3. Using levels 0 and 1/2 we get level 1 using leapfrog.
// 4. Compare level 1 before and after by computing a norm of the difference.
//    If the difference is smaller than a stopping criterion then return.
// 5. If the stopping criterion is not met go to step 2
//

void initial_step(INTEGER const IDENT,
                GridHierarchy &GH,
                GridFunction(DIM)<DOUBLE> &A,
                GridFunction(DIM)<DOUBLE> &B,
                GridFunction(DIM)<DOUBLE> &C,
                GridFunction(DIM)<DOUBLE> &phi,
                DOUBLE h, DOUBLE dt, const DOUBLE stp, const INTEGER maxsteps, 
                const INTEGER Level)
{ /* initial_step */
	
	INTEGER me = MY_PROC(GH);

	INTEGER t, l, idt, nsteps;
	DOUBLE  tnorm, norma, normb, normc;

	DOUBLE h2, dt2;

	l = Level;

	t = CurrentTime(GH,l);
	idt = TimeStep(GH,l);

	h2 = h * 0.5;
	dt2 = -dt * 0.5;

	INTEGER t_sten = 0, s_sten = 1;

	GridFunction(DIM)<DOUBLE> temp("temp", t_sten, s_sten, GH, COMM, HAS_SHADOW);
	GridFunction(DIM)<DOUBLE> temp1("temp1", t_sten, s_sten, GH, COMM, HAS_SHADOW);
	GridFunction(DIM)<DOUBLE> temp2("temp2", t_sten, s_sten, GH, COMM, HAS_SHADOW);
	GridFunction(DIM)<DOUBLE> temp3("temp3", t_sten, s_sten, GH, COMM, HAS_SHADOW);


//	Copy time level 0 to time level 1 and -1

	forall(A,t,l,c)
		A(t+idt,l,c,IDENT) = A(t,l,c,IDENT);
	end_forall
	forall(A,t,l,c)
		A(t-idt,l,c,IDENT) = A(t,l,c,IDENT);
	end_forall

	forall(B,t,l,c)
		B(t+idt,l,c,IDENT) = B(t,l,c,IDENT);
	end_forall
	forall(B,t,l,c)
		B(t-idt,l,c,IDENT) = B(t,l,c,IDENT);
	end_forall

	forall(C,t,l,c)
		C(t+idt,l,c,IDENT) = C(t,l,c,IDENT);
	end_forall
	forall(C,t,l,c)
	        C(t-idt,l,c,IDENT) = C(t,l,c,IDENT);
	end_forall

	forall(phi,t,l,c)
		phi(t-idt,l,c,IDENT) = phi(t,l,c,IDENT);
		phi(t+idt,l,c,IDENT) = phi(t,l,c,IDENT);
	end_forall

  if (0 && IDENT==MAIN) {
       char  fname[40];
       sprintf(fname,"Initial0.dat");
       writeslice(MAIN, fname, GH, phi, t, l, VizServer);
  }

// Repeat the averaging process until the stopping criteria is met
	tnorm = 999999.0;
	nsteps=0;
	while((tnorm > stp) && (nsteps < maxsteps)){
        
        nsteps++;

	// Carry out an average 
	  forall(A,t,l,c)
                  f_average(FBA(A(t,l,c,IDENT)),
                            FDA(A(t,l,c,IDENT)), 
                            FDA(A(t-idt,l,c,IDENT)), 
                            FDA(A(t+idt,l,c,IDENT))); 
	  end_forall

	  forall(B,t,l,c)
                  f_average(FBA(B(t,l,c,IDENT)),
                            FDA(B(t,l,c,IDENT)), 
                            FDA(B(t-idt,l,c,IDENT)), 
                            FDA(B(t+idt,l,c,IDENT))); 
	  end_forall

	  forall(C,t,l,c)
                  f_average(FBA(C(t,l,c,IDENT)),
                            FDA(C(t,l,c,IDENT)), 
                            FDA(C(t-idt,l,c,IDENT)), 
                            FDA(C(t+idt,l,c,IDENT))); 
	  end_forall


	// Store +1 time level before embarking on leapfrog

		forall(temp1,t,l,c)
			temp1(t,l,c,IDENT) = A(t+idt,l,c,IDENT);
		end_forall

		forall(temp2,t,l,c)
			temp2(t,l,c,IDENT) = B(t+idt,l,c,IDENT);
		end_forall

		forall(temp3,t,l,c)
			temp3(t,l,c,IDENT) = C(t+idt,l,c,IDENT);
		end_forall

      if(me == 0){
         cout << "MaxVal on before Initial leapfrog (t) " 
                                   << MaxVal(A, t, l, IDENT) << " "
                                   << MaxVal(B, t, l, IDENT) << " "
                                   << MaxVal(C, t, l, IDENT) << " "
                                   << MaxVal(phi, t, l, IDENT) << " "
                                   << endl << flush;
         cout << "MaxVal on before Initial leapfrog (t-idt) " 
                                   << MaxVal(A, t-idt, l, IDENT) << " "
                                   << MaxVal(B, t-idt, l, IDENT) << " "
                                   << MaxVal(C, t-idt, l, IDENT) << " "
                                   << MaxVal(phi, t-idt, l, IDENT) << " "
                                   << endl << flush;
      }

	// carry out leapfrog step
      forall(phi, t, l, c)
         f_leapfrog(FBA(A(t,l,c,IDENT)),
         FDA(A(t+idt,l,c,IDENT)), FDA(A(t-idt,l,c,IDENT)), FDA(A(t,l,c,IDENT)),
         FDA(B(t+idt,l,c,IDENT)), FDA(B(t-idt,l,c,IDENT)), FDA(B(t,l,c,IDENT)),
         FDA(C(t+idt,l,c,IDENT)), FDA(C(t-idt,l,c,IDENT)), FDA(C(t,l,c,IDENT)),
	    	FDA(phi(t+idt,l,c,IDENT)), FDA(phi(t-idt,l,c,IDENT)), 
         FDA(phi(t,l,c,IDENT)), &dt2, &h2);
      end_forall

	// Sync'ing
        Sync(A,t+idt,l,IDENT);
        Sync(B,t+idt,l,IDENT);
        Sync(C,t+idt,l,IDENT);
        Sync(phi,t+idt,l,IDENT);

      if(me == 0){
         cout << "MaxVal on after Initial leapfrog (t+idt) " 
                                   << MaxVal(A, t+idt, l, IDENT) << " "
                                   << MaxVal(B, t+idt, l, IDENT) << " "
                                   << MaxVal(C, t+idt, l, IDENT) << " "
                                   << MaxVal(phi, t+idt, l, IDENT) << " "
                                   << endl << flush;
      }

	// Calculate the change in the +1 level

		forall(temp1, t, l, c)
			temp1(t,l,c,IDENT) -= A(t+idt,l,c,IDENT); 
		end_forall
		Sync(temp1,t,l,IDENT);

		forall(temp2, t, l, c)
			temp2(t,l,c,IDENT) -= B(t+idt,l,c,IDENT); 
		end_forall
		Sync(temp2,t,l,IDENT);

		forall(temp3, t, l, c)
			temp3(t,l,c,IDENT) -= C(t+idt,l,c,IDENT); 
		end_forall
		Sync(temp3,t,l,IDENT);

		norma = Norm2(temp1,t,l,IDENT);
		normb = Norm2(temp2,t,l,IDENT);
		normc = Norm2(temp3,t,l,IDENT);

		tnorm = sqrt(norma * norma + normb*normb + normc*normc);

		if( MY_PROC(GH) == 0 ){
			cout << "Output from first step averaging procedure\n" << flush;
			cout << norma << " " << normb << " " << normc << "\n" << flush ;
		}
	}

// Copy level 0 to n	with -1 as 1
	forall(A, t, l, c)
		A(t,l,c,IDENT) = A(t-idt,l,c,IDENT);
	end_forall
	
	forall(B,t,l,c)
		B(t,l,c,IDENT) = B(t-idt,l,c,IDENT);
	end_forall
	
	forall(C,t,l,c)
		C(t,l,c,IDENT) = C(t-idt,l,c,IDENT);
	end_forall


	forall(phi,t,l,c)
		phi(t,l,c,IDENT) = phi(t-idt,l,c,IDENT);
	end_forall


//	Sync t=0 level or nth level
	Sync(A,t,l,IDENT);
	Sync(B,t,l,IDENT);
	Sync(C,t,l,IDENT);

	Sync(phi,t,l,IDENT);

	forall(A, t-idt, l, c)
		A(t-idt,l,c,IDENT) = A(t+idt,l,c,IDENT);
	end_forall
	
	forall(B,t-idt,l,c)
		B(t-idt,l,c,IDENT) = B(t+idt,l,c,IDENT);
	end_forall
	
	forall(C,t-idt,l,c)
		C(t-idt,l,c,IDENT) = C(t+idt,l,c,IDENT);
	end_forall


	forall(phi,t-idt,l,c)
		phi(t-idt,l,c,IDENT) = phi(t+idt,l,c,IDENT);
	end_forall

	
} /* end initial_step */
/*-----------------------------------------------------------------------*/
// Routine to take evolve one timestep on the SHADOW and two on the MAIN.
// This way, going into recursive integrate we have two levels with the
// proper data for doing truncation error estimation.
//
void TakeTwoSteps( GridHierarchy &GH,
              GridFunction(DIM)<DOUBLE> &A,
              GridFunction(DIM)<DOUBLE> &B,
              GridFunction(DIM)<DOUBLE> &C,
              GridFunction(DIM)<DOUBLE> &phi,
              const DOUBLE cfl , const INTEGER Level)
{
      INTEGER l = Level;
      INTEGER t = CurrentTime(GH,Level);
      INTEGER idt = TimeStep(GH,Level);

      // Evolve one step on the SHADOW
      DOUBLE dagh_h = DeltaX(GH,X_AXIS,l,SHADOW);
      DOUBLE dagh_dt = dagh_h * cfl;

      if(MY_PROC(GH) == 0 )cout << "Leapfrog on SHADOW \n" << flush;

      // Leapfrog on SHADOW

      forall(phi, t, l, c)
         f_leapfrog(FBA(A(t,l,c,SHADOW)),
           FDA(A(t+idt,l,c,SHADOW)), FDA(A(t-idt,l,c,SHADOW)), 
             FDA(A(t,l,c,SHADOW)),
           FDA(B(t+idt,l,c,SHADOW)), FDA(B(t-idt,l,c,SHADOW)), 
             FDA(B(t,l,c,SHADOW)),
           FDA(C(t+idt,l,c,SHADOW)), FDA(C(t-idt,l,c,SHADOW)), 
             FDA(C(t,l,c,SHADOW)),
           FDA(phi(t+idt,l,c,SHADOW)), FDA(phi(t-idt,l,c,SHADOW)),
           FDA(phi(t,l,c,SHADOW)), &dagh_dt, &dagh_h);
      end_forall

      Sync(A, t+idt, l, SHADOW);
      Sync(B, t+idt, l, SHADOW);
      Sync(C, t+idt, l, SHADOW);

      Sync(phi, t+idt, l, SHADOW);

      // At this point t+idt, t, t-idt are defined. We do not move the
      // stencil up until we are down with children.

      // ******* Update Boundaries on SHADOW ******
      if(MY_PROC(GH) == 0)cout << "Updating boundaries on SHADOW\n" << flush ; 

      UpdateExternalBoundary(phi,t+idt,l,SHADOW);
      UpdateExternalBoundary(A,t+idt,l,SHADOW);
      UpdateExternalBoundary(B,t+idt,l,SHADOW);
      UpdateExternalBoundary(C,t+idt,l,SHADOW);

      Sync(A, t+idt, l, SHADOW);
      Sync(B, t+idt, l, SHADOW);
      Sync(C, t+idt, l, SHADOW);

      Sync(phi, t+idt, l, SHADOW);

      // Move stencil up
      MoveStencilUp(SHADOW, A, t, l,idt);
      MoveStencilUp(SHADOW, B, t, l,idt);
      MoveStencilUp(SHADOW, C, t, l,idt);
      MoveStencilUp(SHADOW, phi, t, l,idt);

      // Evolve two steps on MAIN
       l = Level;
       t = CurrentTime(GH,Level);
       idt = TimeStep(GH,Level);

       dagh_h = DeltaX(GH,X_AXIS,l,MAIN);
       dagh_dt = dagh_h * cfl;

      for(int i = 0; i < 2 ; i++){

        if(MY_PROC(GH) == 0 )cout << "Leapfrog on MAIN " << i << "\n" << flush;

        // Leapfrog on MAIN

        t = CurrentTime(GH,Level);
        forall(phi, t, l, c)
         f_leapfrog(FBA(A(t,l,c,MAIN)),
           FDA(A(t+idt,l,c,MAIN)), FDA(A(t-idt,l,c,MAIN)), FDA(A(t,l,c,MAIN)),
           FDA(B(t+idt,l,c,MAIN)), FDA(B(t-idt,l,c,MAIN)), FDA(B(t,l,c,MAIN)),
           FDA(C(t+idt,l,c,MAIN)), FDA(C(t-idt,l,c,MAIN)), FDA(C(t,l,c,MAIN)),
           FDA(phi(t+idt,l,c,MAIN)), FDA(phi(t-idt,l,c,MAIN)),
           FDA(phi(t,l,c,MAIN)), &dagh_dt, &dagh_h);
        end_forall

        Sync(A, t+idt, l, MAIN);
        Sync(B, t+idt, l, MAIN);
        Sync(C, t+idt, l, MAIN);

        Sync(phi, t+idt, l, MAIN);
  
        // ******* Update Boundaries on MAIN ******
        if(MY_PROC(GH) == 0)cout << "Updating boundaries on MAIN\n" << flush ; 

        UpdateExternalBoundary(phi,t+idt,l,MAIN);
        UpdateExternalBoundary(A,t+idt,l,MAIN);
        UpdateExternalBoundary(B,t+idt,l,MAIN);
        UpdateExternalBoundary(C,t+idt,l,MAIN);

        Sync(A, t+idt, l, MAIN);
        Sync(B, t+idt, l, MAIN);
        Sync(C, t+idt, l, MAIN);

        Sync(phi, t+idt, l, MAIN);

        MoveStencilUp(MAIN, A, t, l,idt);
        MoveStencilUp(MAIN, B, t, l,idt);
        MoveStencilUp(MAIN, C, t, l,idt);
        MoveStencilUp(MAIN, phi, t, l,idt);

        // Increment time step on current level
//        IncrCurrentTime(GH,Level);
      }

} // End TwoStep