/* pde3.c  just checking second order approximation         */
/*      u(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  */
/*       du/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  */
/*   d^2u/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  */
/*       du/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  */
/*   d^2u/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  */
/*                                                          */
/* solve d^2u/dx^2 + d^2u/dy^2 = 12x + 10y = c(x,y)         */
/* by second order method on square -1,-1 to 1,1            */

/* second order numerical difference                        */
/* d^2u/dx^2 + d^2u/dy^2 =                                  */
/* 1/h^2(u(x-h,y)-2u(x,y)+u(x+h,y))+                        */
/* 1/h^2(u(x,y-h)-2u(x,y)+u(x,y+h)) + O(h^2)                */
/*  = c(x,y)                                                */
/* solve for new u(x,y) =                                   */
/* (u(x-h,y)+u(x+h,y)+u(x,y-h)+u(x,y+h)-c(x,y)*h^2)/4       */

#include <stdio.h>
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))

static int nx = 11;
static int ny = 11;
static double us[100][100]; /* solution being iterated        */
static double ut[100][100]; /* temporary for partial results  */
                 /* nx-2 by ny-2 internal, non boundary cells */
static const double xmin = -1.0;
static const double ymin = -1.0;
static const double xmax =  1.0;
static const double ymax =  1.0;
static double h;   /* x  and  y step = (xmax-xmin)/(n-1) = (ymax-ymin) */
static double diff=1.0e100;  /* sum abs difference form last iteration */
static double error;         /* sum of absolute exact-computed         */
static double avgerror;      /* average error                          */
static double maxerror;      /* maximum cell error                     */

static double u(double x, double y) /* for boundary and check */
{
  return x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0;
}

static double c(double x, double y) /* from solve */
{
  return 12.0*x + 10.0*y;
}

static double u00(int i, int j) /* new iterated value */
{
  return ( us[i-1][j] + us[i+1][j] + us[i][j-1] + us[i][j+1] -
           c(xmin+i*h,ymin+j*h)*h*h )/4.0;
}

static void init_boundary()
{
  int i, j;
  
  /* set boundary on four sides */
  for(i=0; i<nx; i++)
  {
    us[i][0]    = u(xmin+i*h, ymin);
    us[i][ny-1] = u(xmin+i*h, ymax);
  }
  for(j=0; j<ny; j++)
  {
    us[0][j]    = u(xmin, ymin+j*h);
    us[nx-1][j] = u(xmax, ymin+j*h);
  }
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      us[i][j] = 0.0;
} /* end init_boundary */

static void iterate()
{
  int i, j;
  
  diff = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      ut[i][j] = u00(i,j); 
      diff = diff + abs(ut[i][j]-us[i][j]);
    }
  }
  /* copy temp to solution */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      us[i][j] = ut[i][j];
} /* end iterate */

static void check() /* typically not known, instrumentation */
{
  int i, j;
  double uexact, err;
  
  error = 0.0;
  maxerror = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      uexact = u(xmin+i*h,ymin+j*h);
      err = abs(us[i][j]-uexact);
      error = error + err;
      if(err>maxerror) maxerror=err;
    }
  }
} /* end check */

static void printus()
{
  int i, j;
  
  printf("\n computed solution\n");
  for(i=0; i<min(nx,11); i++)
  {
    for(j=0; j<min(ny,11); j++)
    {
      printf("%6.3f ", us[i][j]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printus */

static void printexact()
{
  int i, j;
  
  printf("\n exact solution\n");
  for(i=0; i<min(nx,11); i++)
  {
    for(j=0; j<min(ny,11); j++)
    {
      printf("%6.3f ", u(xmin+i*h,ymin+j*h));
    }
    printf("\n");
  }
  printf("\n");
} /* end printexact */

static void printdiff()
{
  int i, j;
  
  printf("\n difference exact-computed solution\n");
  printf("\n"); 
  for(i=1; i<min(nx-1,10); i++)
  {
    printf("       ");
    for(j=1; j<min(ny-1,10); j++)
    {
      printf("%7.4f", u(xmin+i*h,ymin+j*h)-us[i][j]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printdiff */

int main(int argc, char * argv[])
{
  int iter=0;
  
  printf("pde3.c running\n");
  h = (xmax-xmin)/(nx-1); /* step size in x */
  /*= (ymax-ymin)/(ny-1);    step size in y */
  printf("xmin=%g, xmax=%g, nx=%d, h=%g \n", xmin, xmax, nx, h);
  printf("ymin=%g, ymax=%g, ny=%d, h=%g \n", ymin, ymax, ny, h);
  
  printexact();
  init_boundary();
  printf(" initial condition\n");
  printus();
  
  while(diff>0.0001 && iter<200)
  {
    iterate();
    check();
    avgerror = error/((nx-2)*(ny-2));
    printf("iter=%d, diff=%g, error=%g, avg=%g, max=%g\n",
           iter, diff, error, avgerror, maxerror);
    iter++;
    if(iter%50==0) { printus(); printdiff(); }
  }
  return 0;
} /* end main of pde3.c */