#include "utils.h"

// Required from solver
// fields updated by time-integration
extern scalar * evolving;

///////////////////////////////////////////////////////////////////////////
// Butcher tables, simplified by Kopal and re-organized by Williamson
// See "Low-storage runge-kutta schemes", williamson,1980
// & "Fourth-order 2N-storage rk schemes", Carpenter & Kennedy, 1994, NASA
static double ** butcherd(order){
  // butcher[0][i] = ai
  // butcher[1][i] = bi
  double ** butcher;
  butcher = malloc(2*order*sizeof(double));
  butcher[0] = malloc(order*sizeof(double));
  butcher[1] = malloc(order*sizeof(double));

  butcher[0][0] = 0; //a1
  if(order == 1)
    butcher[1][0] = 1;
  else if(order == 2){
//      a[]                         b[]  
                           butcher[1][0] = 0.5;//1
  butcher[0][1] = -0.5;    butcher[1][1] = 1;//2
  }
  else if(order == 3){
//      a[]                         b[]   
                             butcher[1][0] = 1/3.; //1
  butcher[0][1] = -5/9.;     butcher[1][1] = 15/16.; //2
  butcher[0][2] = -153/128.; butcher[1][2] = 8/15.; //3
  }
  else if(order == 5){
    // 5 steps for the 4th-order scheme
//      a[]                        
  butcher[0][1] = -1 ;     
  butcher[0][2] = -1/3.+pow(2,2/3.)/6. - 2*cbrt(2)/3.;     
  butcher[0][3] = -cbrt(2) - pow(2,2/3.) - 2;      
  butcher[0][4] = -1 + cbrt(2);      

  //   b[]
  butcher[1][0] = 2/3. + cbrt(2)/3. + pow(2,2/3.)/6.; //1
  butcher[1][1] = -pow(2,2/3.)/6. + 1/6.;    //2
  butcher[1][2] = -1/3. - 2*cbrt(2)/3. - pow(2,2/3.)/3.; //3
  butcher[1][3] = 1/3. - cbrt(2)/3. - pow(2,2/3.)/6.; //3
  butcher[1][4] = 1/3. + cbrt(2)/6. + pow(2,2/3.)/12.; //3
  }
  else assert(false);
 return butcher;
}

//////////////////////////////////////////////////////
// User-provided functions
// gradient
double (* gradient)  (double, double, double) = minmod2;

// how to compute updates
double (* update) (scalar * evolving, scalar * updates, double dtmax) = NULL;

double t = 0., dt = 0.;
static void advance_generic (scalar * output, scalar * input, scalar * updates,
			     double dt){
  trash (output);  
  foreach() {
    scalar o,i,u;
    for(o,i,u in output,input,updates)
       o[] = i[]+dt*u[];
  }
  boundary (output);
}

static void (* advance) (scalar * output, scalar * input, scalar * updates,
			 double dt) = advance_generic;
event defaults (i = 0)
{
  // limiting
  for (scalar s in all)
    s.gradient = gradient;

  // default values
  foreach()
    for (scalar s in all)
      s[] = 0.;
  boundary (all);
}


//////////////////////////////////////////////////////
// Run function
void run()
{
  t = 0.;
  init_grid (N);
  timer start = timer_start();
  int i = 0; long tnc = 0;

  // Take care of the order 4 which is in 5 steps
  int order_k = order_t < 4 ? order_t : 5;
  // williamson/carpenter tables from butcher tables
  double ** butcher = butcherd(order_k);
      
  while (events (i, t, true)) {
    
    // list of updates and intermediary fields
    scalar * updates = list_clone (evolving);
    scalar * interm =  list_clone (evolving); 
    // interm = 0
    foreach(){
      for( scalar s in interm)
	s[] = 0;
    }
    
    dt = dtnext (t, update (evolving, updates, DT));
    for( int order_ = 0; order_ < order_k ; order_ ++ ){
      /* Nnd-order time-integration */
      foreach() {
	for(scalar s in interm)
	  s[] *= butcher[0][order_];
      }
      //interm = a*qj-1
      
      update(evolving,updates,dt); // update = F(Yj-1)
      advance (interm, interm, updates, dt);  // interm = qj = a*qj-1 + dt F(Yj-1) 
     
      foreach(){
	scalar s1,s2;
	for(s1,s2 in evolving,interm)
	  s1[] += s2[]*butcher[1][order_];

      }
      // Yj = Yj-1 + bj*qj
    }
    
    delete(interm);
    free(interm);
    delete (updates);
    free (updates);
    
    long nc = 0;
    foreach (reduction(+:nc)) 
      nc++;
    tnc += nc;
    i++; t = tnext;
  }
 
  timer_print (start, i, tnc);
  free_grid();
  free(butcher);
}