/*******************************************************************
 * diffusion1d_spec.c: sequential 1d-diffusion solver with constant
 * boundaries.
 *
 * To execute: cc diffusion1d_spec.c && ./a.out
 *
 * To verify: civl verify diffusion1d_spec.c
 * Author: Ziqing Luo
 ********************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>

#ifdef _CIVL
#include <civlc.cvh>
$input int NXB = 5;          // upper bound on nx
$input int nx;               // global number of points excl. boundary
$assume(1<=nx && nx<=NXB);
$input double U_INIT[nx+2];  // initial values for temperature incl. boundary
$input double k;             // the constant D*dt/(dx*dx)
$assume(k>0 && k<.5);
$input int NSTEPS_BOUND=5;   // upper bound on nsteps
$input int nsteps;           // number of time steps
$assume(1<=nsteps && nsteps<=NSTEPS_BOUND);
$input int wstep;            // write frame every this many time steps
$assume(1<=wstep && wstep<=nsteps);
$output double output[nx+2]; // output values, write only
#else
int nx;
double k;
int nsteps;
int wstep;
#endif

/* Global variables */
double *u;     /* temperature function, local */
double *u_new; /* second copy of temperature function, local */
double *buf;   /* temp. buffer used on proc 0 only */
int print_pos; /* number of cells printed on current line */
int time=0;    /* current time step */

void initialize() {
  u = (double *)malloc((nx+2) * sizeof(double));
  u_new = (double *)malloc((nx+2) * sizeof(double));
#ifdef _CIVL
  for (int i=0; i<nx+2; i++)
    u[i] = U_INIT[i];
  u_new[0] = u[0];
  u_new[nx+1] = u[nx+1];
#else
  nx = 10;
  k = 0.2;
  nsteps = 20;
  wstep = 2;
  u[nx+1] = u[0] = 0.0;
  u_new[nx+1] = u_new[0] = 0.0;
  for (int i=1; i<=nx; i++)
    u[i]=100.0;
#endif
  printf("nx=%d, k=%lf, nsteps=%d, wstep=%d\n",
	 nx, k, nsteps, wstep);
}

void write_frame(int time) {
  printf("========== Time %d ==========\n", time);
  for (int i = 0; i < nx + 2; i++) {
    printf("%lf ", u[i]);
#ifdef _CIVL
    output[i] = u[i];
#endif
  } 
  printf("\n");
}

/* Updates u_new using u, then swaps u and u_new. */
void update() {
  for(int i = 1; i <= nx; i++) 
    u_new[i] = u[i] + k*(u[i+1] + u[i-1] - 2*u[i]);
  
  double * tmp = u_new; u_new=u; u=tmp;
}

int main() {
  initialize();
  write_frame(0);
  for (int time = 1; time < nsteps; time++){
    update();
    if(time % wstep == 0)
      write_frame(time);
  }
#ifdef _CIVL
  memcpy(output, u, sizeof(double) * (nx+2));
#endif
  free(u);
  free(u_new);
  return 0;
}
