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

#ifdef _CIVL

#include <civlc.cvh>
$input long NXB = 3;               // nx upper bound
$input long nx;                    // global number of columns in matrix
$assume(1 <= nx && nx <= NXB);
$input long NYB = 3;               // ny upper bound
$input long ny;                    // global number of rows of matrix
$assume(1 <= ny && ny <= NYB);
$input double u_init[ny+2][nx+2];  // initial value of temperatures, including boundaries
$input double k;                   // constant coefficient  
$assume(k > 0.0 && k < 0.5);
$input int NSTEPSB = 3;            // upper bound for nsteps
$input int nsteps;                 // number of steps
$assume(1<=nsteps && nsteps<=NSTEPSB);
$input int wstep = 1;              // write frame every this many time steps
$output double output[ny+2][nx+2];

#else

long nx, ny;
int nsteps, wstep;
double constTemp;                  // value of constant boundaries for test
double initTemp;                   // value of initial temperature for test
double k;

#endif

/* dynamically-allocated 2d array of doubles specifying state of local
 * grid in current time step. Dimensions are u_curr[ny+2][nx+2].
 */
double ** u_curr;
/* dynamically-allocated 2d array of doubles specifying state of local
 * grid in next time step. Dimensions are u_next[ny+2][nx+2].
 */
double ** u_next;

/* initialize all data values owned by each process */
void initialization() {
  u_next = (double **)malloc((ny + 2) * sizeof(double *));
  u_curr = (double **)malloc((ny + 2) * sizeof(double *));
  for (int i = 0; i < ny + 2; i++) {
    u_next[i] = (double *)malloc((nx + 2) * sizeof(double));
    u_curr[i] = (double *)malloc((nx + 2) * sizeof(double));
  }
#ifdef _CIVL

  for (int i=0; i < ny + 2; i++)
    for (int j=0; j < nx + 2; j++)
      u_next[i][j] = u_curr[i][j] = u_init[i][j];

#else

  // All boundary cells are set to the same value constTemp.
  // All cells in the interior region are set to the same value 
  // initTemp.
  for (int i=0; i < ny + 2; i++)
    for (int j=0; j < nx + 2; j++)
      if (i == 0 || j == 0 || i == ny+1 || j == nx+1)
        u_next[i][j] = u_curr[i][j] = constTemp;
      else
        u_curr[i][j] = initTemp;
  nsteps = 300;
  wstep = 5;
  nx = 15;
  ny = 15;
  constTemp = 0.0;
  initTemp = 100.0;
  k = 0.13;

#endif
}

void update() {
  for (long i = 1; i <= ny; i++)
    for (long j = 1; j <= nx; j++)
      u_next[i][j] = u_curr[i][j] +
	k*(u_curr[i+1][j] + u_curr[i-1][j] + 
	   u_curr[i][j+1] + u_curr[i][j-1] - 
	   4*u_curr[i][j]);

  double ** tmp = u_curr;
  
  u_curr = u_next;
  u_next = tmp;
}

int main(int argc, char * argv[]) {
  int i,j;

  initialization();
  for (i=0; i<nsteps; i++) 
    update();
#ifdef _CIVL
  for (i = 1; i < ny+1; i++)
    memcpy(&output[i][1], &u_curr[i][1], sizeof(double) * nx);
#endif
  for (i=0; i<ny+2; i++) {
    free(u_curr[i]);
    free(u_next[i]);
  }
  free(u_curr);
  free(u_next);
  return 0;
}
