/* pde2_4th.c   solution of d^2u(x,y)/dx^2 + d^2u(x,y)/dy^2 = f(x,y)
 *              on  0 < x < 1, 0 < y < 1.0, f(x,y) = -2cos(x)sin(y)
 *              boundaries and solution u(x,y) = cos(x)sin(y)
 *
 * find solution for h=hx=0.1  h=hx=hy
 * then h=0.05 1/20, 0.025 1/40, 0.0125 1/80
 * order(4,4) compared to order (2,2) 
 * u(x,y)=1/20 (u(x+h,y+h)+u(x+h,y-h)+u(x-h,y+h)+u(x-h,y-h)) +
 *        1/5  (u(x+h,y)  +u(x-h,y)  +u(x,y+h)  +u(x,y-h)) -
 *        h^2/12 (f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)+8*f(x,y))
 *
 * compared to order (2,2)
 * u(x,y) = 1/4 (u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)-h^2 f(x,y)) - 
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#define max(a,b) ((a)>(b)?(a):(b))
#define pi 3.141592653589793

static double u(double xx, double yy);
static double f(double xx, double yy);
static double ff(double xx, double yy, double hh);
static double us[100][100]; /* 4th order solution */
static double ut[100][100]; /* temporary */
static double ws[100][100]; /* second order solution */
static double wt[100][100]; /* temporary */

int main(int argc, char *argv[])
{
  double xmin, xmax, ymin, ymax;
  double h;
  double sum4, sum4prev, sum2, sum2prev;
  double maxdif, maxerr, dif4, dif4prev, dif2, dif2prev;
  int i, j, hval, n, m;

  printf("pde2_4th.c running \n");
  printf("sum2,4 are L2 norm of errors \n");
  printf("dif2,4 are average difference from analytic solution \n");
  xmin = 0.0;      /* values from problem statement */
  xmax = 1.0;
  ymin = 0.0;
  ymax = 1.0;
  sum4 = 0.0;
  sum2 = 0.0;
  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  h = 0.1;         /*  space x, h = step size */
  for(hval=1; hval<=4; hval++) /* halve h each pass */
  {
    n = (int)((xmax-xmin)/h)+1;  /* intervals, subscripts 0 thru n */
                               /* [0][j] [n][j] [i][0] [i][n] are boundaries */
    printf("hval=%d, h=%g, n=%d \n", hval, h, n);

    /* set boundary conditions */
    for(i=0; i<=n; i++)
    {
      us[i][0] = u(xmin+(double)(i)*h, ymin);
      ut[i][0] = u(xmin+(double)(i)*h, ymin);
      us[i][n] = u(xmin+(double)(i)*h, ymax);
      ut[i][n] = u(xmin+(double)(i)*h, ymax);
      ws[i][0] = u(xmin+(double)(i)*h, ymin);
      wt[i][0] = u(xmin+(double)(i)*h, ymin);
      ws[i][n] = u(xmin+(double)(i)*h, ymax);
      wt[i][n] = u(xmin+(double)(i)*h, ymax);
    }
    for(j=0; j<=n; j++)
    {
      us[0][j] = u(xmin, ymin+(double)(j)*h);
      ut[0][j] = u(xmin, ymin+(double)(j)*h);
      us[n][j] = u(xmax, ymin+(double)(j)*h);
      ut[n][j] = u(xmax, ymin+(double)(j)*h);
      ws[0][j] = u(xmin, ymin+(double)(j)*h);
      wt[0][j] = u(xmin, ymin+(double)(j)*h);
      ws[n][j] = u(xmax, ymin+(double)(j)*h);
      wt[n][j] = u(xmax, ymin+(double)(j)*h);
    }
    for(i=1; i<n; i++)
    {
      for(j=1; j<n; j++)
      {
        us[i][j] = 0.0;
        ws[i][j] = 0.0;
      }
    }

    /* iterate solution until change small enough */
    maxdif = 999.9;
    m = 1;
    while(maxdif>0.0001)
    {
      for(i=1; i<n; i++)
      {
        for(j=1; j<n; j++)
	{
          ut[i][j] = (us[i+1][j+1]+us[i-1][j+1]+us[i+1][j-1]+us[i-1][j-1])/20.0 +
	             (us[i+1][j]  +us[i-1][j]  +us[i][j+1]  +us[i][j-1])/5.0 -
	             ff(xmin+(double)(i)*h, ymin+(double)(j)*h, h);
          wt[i][j] = (ws[i-1][j]+ws[i+1][j]+ws[i][j+1]+ws[i][j-1] -
		      h*h*f(xmin+(double)(i)*h, ymin+(double)(j)*h))/4.0;
	}
      }
      if(m<=n*n) /* "precondition" */
      {
        for(i=1; i<n; i++)
        {
          for(j=1; j<n; j++)
	  {
            us[i][j] = ut[i][j];
            ws[i][j] = wt[i][j];
	  }
        }
      }
      else
      {
        dif4prev = dif4;
        dif2prev = dif2;
        dif4 = 0.0;
        dif2 = 0.0;
        maxdif = 0.0;
        maxerr = 0.0;
        for(i=1; i<n; i++)
        {
          for(j=1; j<n; j++)
	  {
            maxdif = max(max(abs(us[i][j]-ut[i][j]),
                             abs(ws[i][j]-wt[i][j])),maxdif);
            maxerr = max(max(abs(us[i][j]-
                          u(xmin+(double)(i)*h, ymin+(double)(j)*h)),
                          abs(ws[i][j]-
                          u(xmin+(double)(i)*h, ymin+(double)(j)*h))),maxerr);
            dif4 = dif4 + abs(us[i][j]-ut[i][j]);
            us[i][j] = ut[i][j];
            dif2 = dif2 + abs(ws[i][j]-wt[i][j]);
            ws[i][j] = wt[i][j];
	  }
        }
        dif4 = dif4/(double)((n-1)*(n-1));
        dif2 = dif2/(double)((n-1)*(n-1));
        printf("m=%d, dif4=%g, dif2=%g, maxdif=%g, maxerr=%g \n",
             m, dif4, dif2, maxdif, maxerr);
      }
      m = m+1;
    }

    /* compute L2 norm of each error and ratio */
    sum4 = 0.0;
    sum2 = 0.0;
    for(i=0; i<=n; i++) /* checking boundary for printout */
    {
      for(j=0; j<=n; j++)
      {
        dif4 = us[i][j]-u(xmin+(double)(i)*h, ymin+(double)(j)*h);
        dif2 = ws[i][j]-u(xmin+(double)(i)*h, ymin+(double)(j)*h);
        sum4 = sum4 + dif4*dif4;
        sum2 = sum2 + dif2*dif2;
        if(hval==1) printf("us[%d][%d]=%g, u=%g, dif4=%g \n",
                           i, j, us[i][j],
                  u(xmin+(double)(i)*h, ymin+(double)(j)*h), dif4);
      }
    }
    sum4 = sqrt(sum4/(double)((n-1)*(n-1)));
    sum2 = sqrt(sum2/(double)((n-1)*(n-1)));
    printf("sum4=%g, sum2=%g \n",sum4, sum2);
    if(hval>1)
    {
      printf("ratio4=%g, ratio2=%g \n \n",
	     sum4/sum4prev, sum2/sum2prev);
    }
    sum4prev = sum4;
    sum2prev = sum2; 
    h=h/2.0;
  } /* loop on size of h */
  return 0;
}

static double u(double xx, double yy) /* boundary condition and solution */
{
  double val;
 
  val = cos(xx)*sin(yy);
  return val;
} /* end  u */

static double f(double xx, double yy) /* forcing function */
{
  double val;
 
  val = -2.0*cos(xx)*sin(yy);
  return val;
} /* end  f */

static double ff(double xx, double yy, double hh) /* initial condition u(x,0) */
{
  double val;
  /* 3/10 * 1/12 = 1/40 */
  val = (hh*hh/40.0)*(f(xx+hh,yy)+f(xx-hh,yy)+f(xx,yy+hh)+f(xx,yy-hh)+8.0*f(xx,yy));
  return val;
} /* end  ff */
/* pde2_4th.c */