/* ///////////////////////////////////////////////////////////////////// */
/*! 
  \file  
  \brief Solve 1D heat equation.

  Solve the 1D heat equation using a 1st order explicit method
  on a parallel domain.
 
  \author A. Mignone (mignone@to.infn.it)
  \date   March 12, 2020
*/
/* ///////////////////////////////////////////////////////////////////// */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NX_GLOB   64  /* Global number of interior points */
#define NGHOST     1

void Write (double *, double *, int, int);

int main(int argc, char ** argv)
{
  int    i, k, beg, end;
  int    nx_loc;               /* Local grid size */
  double t, tstop, dt, cfl = 0.5;
  double *u0;  
  double *u1;
  double xbeg =  0.0;
  double xend = +1.0;
  double xglob[NX_GLOB + 2*NGHOST];   // Global grid array
  double *xloc;
  double dx;    /* Mesh spacing */

/* --------------------------------------------------------
   1. Generate global & local grids
   -------------------------------------------------------- */

  nx_loc = NX_GLOB;    
  beg    = NGHOST;
  end    = beg + nx_loc - 1;

  dx = (xend - xbeg)/(NX_GLOB+1);
  for (i = 0; i < NX_GLOB + 2*NGHOST; i++){
    xglob[i] = xbeg + (i-beg+1)*dx;
  }
  xloc = xglob;  /* Use pointer arithmetic */

/* --------------------------------------------------------
   2. Allocate memory on local grids
   -------------------------------------------------------- */

  u0 = (double *) malloc((nx_loc + 2*NGHOST)*sizeof(double));
  u1 = (double *) malloc((nx_loc + 2*NGHOST)*sizeof(double));

/* --------------------------------------------------------
   3. Set initial condition
   -------------------------------------------------------- */

  for (i = beg; i <= end; i++){
    u0[i] = sin(M_PI*xloc[i]);
  }
  
/* --------------------------------------------------------
   4. Advance solution
   -------------------------------------------------------- */

  t     = 0.0;
  tstop = 0.1;
  dt    = cfl*dx*dx;  /* dt = C*dx^2/D with D = 1 */
  k     = 0;

  Write (xloc, u0, beg, end);
  while (t <= tstop){

    printf ("step #%d; t = %8.3e\n",k,t);

  /* -- 4a. Set physical boundary conditions -- */

    u0[beg-1] = 0.0;
    u0[end+1] = 0.0;

  /* -- 4b. Advance solution by one time step -- */

    for (i = beg; i <= end; i++){
      u1[i] = u0[i] + dt/(dx*dx)*(u0[i-1] - 2.0*u0[i] + u0[i+1]);
    }
    t += dt;
    k++;

  /* -- 4c. Copy arrays for next time level -- */

    for (i = beg; i <= end; i++) u0[i] = u1[i]; 
  }
  Write (xloc, u0, beg, end);

  return 0;
}

/* ********************************************************************* */
void Write (double *x, double *u, int beg, int end)
/*
 *********************************************************************** */
{
  int    i;
  static int n = 0;  /* File number */
  FILE *fp;
  char fname[32];

  sprintf (fname,"heat_eq%02d.dat",n);
  fp = fopen (fname,"w");  
  for (i = beg; i <= end; i++) fprintf (fp, "%12.6e  %12.6e\n", x[i], u[i]);
  fclose(fp);

  n++;
}