/****************************************************************************
 *  * FILE: ser_heat2D.c
 *   * DESCRIPTION:  
 *    *   Serial HEAT2D Example - C Version
 *     *   This example is based on a simplified 
 *      *   two-dimensional heat equation domain decomposition.  The initial 
 *       *   temperature is computed to be high in the middle of the domain and 
 *        *   zero at the boundaries.  The boundaries are held at zero throughout 
 *         *   the simulation.  During the time-stepping, an array containing two 
 *          *   domains is used; these domains alternate between old data and new data.
 *           * AUTHOR: D. Turner
 *            * Last Revised: 04/15/05 Blaise Barney
 *             ****************************************************************************/
#include <stdio.h>
#include <omp.h>
#include <stdlib.h>
#define NXPROB 1000
#define NYPROB 1000
void inidat(int, int, float *);
void prtdat(int, int, float *, char *);
void update(int,int,float *, float *);

struct Parms
{ 
  float cx;
  float cy;
  int nts;
} parms = {0.1, 0.1, 1000};

int main(int argc, char *argv[])
{
float u[2][NXPROB][NYPROB];
int ix, iy, iz, it;
double wtime = omp_get_wtime();
printf("Starting serial version of 2D heat example...\n");
printf("Using [%d][%d] grid.\n",NXPROB, NYPROB);

/* Initialize grid and create input file */
printf("Initializing grid and creating input file:");
inidat(NXPROB, NYPROB, &u[0][0][0]);
prtdat(NXPROB, NYPROB, &u[0][0][0], "initial.dat");
for (ix = 0; ix <= NXPROB-1; ix++) {
   u[1][ix][0] = u[0][ix][0];
   u[1][ix][NYPROB-1] = u[0][ix][NYPROB-1];
   }
for (iy = 0; iy <= NYPROB-1; iy++) {
   u[1][0][iy] = u[0][0][iy];
   u[1][NXPROB-1][iy] = u[0][NXPROB-1][iy];
   }

/* Iterate over all timesteps and create output file */
printf("Iterating over %d time steps...\n",parms.nts);
iz = 0;
for (it = 1; it <= parms.nts; it++) {
   update(NXPROB, NYPROB, &u[iz][0][0], &u[1-iz][0][0]);
   iz = 1 - iz;
   }
printf("Done. Created output file: ");
prtdat(NXPROB, NYPROB, &u[iz][0][0], "final.dat");
printf("Time taken by the serial program: %lf\n",omp_get_wtime()-wtime);
}


/****************************************************************************
 *  *  subroutine update
 *   ****************************************************************************/
void
update(int nx, int ny, float *u1, float *u2)
{
   int ix, iy;

   for (ix = 1; ix <= nx-2; ix++) {
      for (iy = 1; iy <= ny-2; iy++) {
         *(u2+ix*ny+iy) = *(u1+ix*ny+iy)  + 
         parms.cx * (*(u1+(ix+1)*ny+iy) + *(u1+(ix-1)*ny+iy) - 
         2.0 * *(u1+ix*ny+iy)) +
         parms.cy * (*(u1+ix*ny+iy+1) + *(u1+ix*ny+iy-1) - 
         2.0 * *(u1+ix*ny+iy));
         }
      }
}

/*****************************************************************************
 *  *  subroutine inidat
 *   *****************************************************************************/
void
inidat(int nx, int ny, float *u1)
{
   int ix, iy;

   for (ix = 0; ix <= nx-1; ix++) 
      for (iy = 0; iy <= ny-1; iy++) 
         *(u1+ix*ny+iy) = (float)(ix * (nx - ix - 1) * iy * (ny - iy - 1));
}

/**************************************************************************
 *  * subroutine prtdat
 *   **************************************************************************/
void
prtdat(int nx, int ny, float *u1, char *fnam)
{
   int ix, iy;
   FILE *fp;

   fp = fopen(fnam, "w");
   for (iy = ny-1; iy >= 0; iy--) {
      for (ix = 0; ix <= nx-1; ix++) {
        fprintf(fp, "%8.3f", *(u1+ix*ny+iy));
        if (ix != nx-1) {
           fprintf(fp, " ");
           }
       else {
          fprintf(fp, "\n");
          }
       }
    }
   fclose(fp);
   printf(" %s\n",fnam);
}
