/* abc_eq.c  discretize, build linear equations, solve
 * solve  a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
 * d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y) = c(x,y)
 *
 * where, for this test case
 *  a1(x,y) = exp(x/2)*exp(y)/2
 *  b1(x,y) = 0.7/(x*x*y*y + 0.5)
 *  c1(x,y) = (4 - exp(x) - exp(y/2))*2
 *  d1(x,y) = x*x+y
 *  e1(x,y) = x*y*y
 *  f1(x,y) = 3.0*x + 2.0*y
 *
 *  c(x,y) = 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)            +
 *     0.7*(6.0*x + 8.0*y + 5.0)/(x*x*y*y+0.5)              +
 *     (8.0 - 2.0*exp(x) - 2.0*exp(y/2.0))*(12.0*y + 8.0*x) +
 *     (x*x+y)*(3.0*x*x + 6.0*x*y + 4.0*y*y + 5.0*y + 6.0)  +
 *     x*y*y*(6.0*y*y + 3.0*x*x + 8.0*x*y + 5.0*x +7.0)     +
 *     3.0*x + 2.0*y;
 *
 * The solution (normally only the boundary values are given)
 * u(x,y) = x^3 + 2y^3 + 3x^2 y + 4xy^2 + 5xy + 6x + 7y + 8
 *
 * The solver uses  #include "abc.h" to gain access to the users problem
 * Compile with abc.c to be able to call functions defined above
 *
 * replace continuous derivatives with discrete derivatives
 * solve for U(i,j) numerical approximation U(x_i,y_j)
 * A * U = F, solve simultaneous equations for U
 *
 * subscripts in code use x[i], y[ii]
 * derivatives are computed with b,hx,a included as in
 * cx[j] for first derivative of x at point j
 * cxx[j] for second derivative of x at point j
 * cx[j]*cy[jj] for d^2u/dy dx at j, jj
 * 
 * The boundaries xmin..xmax, ymin..ymax  are known
 * The number of grid points in each dimension is known nx, ny
 * hx=(xmax-xmin)/(nx-1), hy=(ymax-ymin)/(ny-1),
 * thus xg(i)=xmin+i*hx, yg(j)=ymin+j*hy for this test
 * non uniform grids could have non uniformly spaced values in xg, yg
 * then the linear system of equations is built,
 * nxy=(nx-2)*(ny-2) equations because the boundary value
 *                   is known on each end.
 * The matrix A is nxy rows by nxy columns (nxy+1) columns including F
 * The forcing function vector F is nxy elements adjunct to A
 * The solution vector U is nxy elements.
 * The building of the A matix requires 4 nested loops for rows
 * i  in 1..nx-1, j  in 1..ny-1,   for rows of A
 * then each term in the differential equation fills in that row
 * and cs = (nx-2)*(ny-2) the RHS, fills the last column of A
 * subscripting functions are used to compute linear subscripts of A, U
 * row=(i-1)*(ny-2)+(ii-1)
 * row*(nxy+1)+col for A
 *
 * care must be taken to move the negative of boundary elements to
 * the right hand size of the equation to be solved. These are subtracted
 * from the computed value of the forcing function, cs, for that row.
 * tests are j==0 || j==nx-1    jj==0 || jj==ny-1 for boundaries
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "deriv.h"
#include "simeq.h"
#include "abc.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

/* nx, ny, xmin, xmax, ymin, ymax defined in abc.h */
#define nxy (nx-2)*(ny-2)
#define nn max(nx,ny)

int s(int i, int ii) /* for x,y */
{
  return (i-1)*(ny-2) + (ii-1);
} /* end s */

int sk(int k, int i, int ii) /* for x,y */
{
  return k*(nxy+1) + (i-1)*(ny-2) + (ii-1);

} /* end sk */

int ss(int i, int ii, int j, int jj)
{
  return s(i,ii)*(nxy+1) + s(j,jj);
} /* end ss */

static double A[nxy*(nxy+1)]; /* last column is RHS */
static double U[nxy];
static double xg[nx];
static double yg[ny];
static double Ua[nxy];
static double hx, hy;

void printA()
{
  int i, ii, j, jj, cs, k;
  cs = (nx-2)*(ny-2); /* constant RHS column */
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      k = s(i,ii); /* row index */
      for(j=1; j<nx-1; j++)
      {
        for(jj=1; jj<ny-1; jj++)
        {
          printf("A[row(%d,%d),col(%d,%d)] =%f\n",
                 i,ii,j,jj, A[ss(i,ii,j,jj)]);
	} /* jj */
      } /* j */
      printf("RHS A[%d,%d]=%f\n\n", k, cs, A[k*(nxy+1)+cs]);
             /* equivalent A[sk(k,i,ii+1)] */
    } /* ii */
  } /* i */
  printf("\n");
} /* end printA */

double eval_derivative(int xord, int yord,
                       int i, int ii, double UU[])
{
  int j, jj;
  double cx[nx];
  double cy[ny];
  double p1[nn];
  double discrete;
  double x, y;

  rderiv(xord, nx, i, hx, cx);
  x = xg[i];
  rderiv(yord, ny, ii, hy, cy);
  y = yg[ii];
  discrete = 0.0;

  /* cases of one partial x, y  to any order */
  if(xord!=0 && yord==0)
    for(j=0; j<nx; j++)
      if(j==0 || j==nx-1) discrete += cx[j]*u(xg[j],y);
      else discrete += cx[j]*UU[s(j,ii)];
  else if(xord==0 && yord!=0)
    for(jj=0; jj<ny; jj++)
      if(jj==0 || jj==ny-1) discrete += cy[jj]*u(x,yg[jj]);
      else discrete += cy[jj]*UU[s(i,jj)];

  /*  cases of two partials xy  to any order */
  else if(xord!=0 && yord!=0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jj=0; jj<ny; jj++)
        if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
             p1[j] += cy[jj]*u(xg[j],yg[jj]);
        else p1[j] += cy[jj]*UU[s(j,jj)];
      discrete += cx[j]*p1[j];
    }
  }
  return discrete;
} /* end eval_deriv */

void check_soln(double UU[])
{
  /* check  a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
   * d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y) = c(x,y)
   */
  double val, err, smaxerr;
  int i, ii;
  double x, y;

  smaxerr = 0.0;
  for(i=1; i<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      val = 0.0;
      /* a1(x,y)*uxx(x,y) */
      val += a1(x,y)*eval_derivative(2, 0, i, ii, UU);
      /* + b1(x,y)*uxy(x,y) */
      val += b1(x,y)*eval_derivative(1, 1, i, ii, UU);
      /* + c1(x,y)*uyy(x,y) */
      val += c1(x,y)*eval_derivative(0, 2, i, ii, UU);
      /* + d1(x,y)*ux(x,y) */
      val += d1(x,y)*eval_derivative(1, 0, i, ii, UU);
      /* + e1(x,y)*uy(x,y) */
      val += e1(x,y)*eval_derivative(0, 1, i, ii, UU);
      /* + f1(x,y)*u(x,y) */
      val += f1(x,y)*u(xg[i],yg[ii]);
      /* -  c(x,y) should be zero */
      err = abs(val - c(xg[i],yg[ii]));
      if(err>smaxerr) smaxerr = err;
    } /* end ii */
  } /* end i */
  printf("check_soln maxerr=%g \n",smaxerr);
} /* end check_soln */
 

int main(int argc, char *argv[])
{
  double x, y;
  int i, ii, j, jj;
  int k, cs;
  /* coeficients of terms in first equation */
  double coefdxx; /* computed below if not constant */
  double coefdxy;
  double coefdyy;
  double coefdx;
  double coefdy;
  double coefu;
  /* derivative point arrays, dynamic in this version */
  double cx[nx], cy[ny]; 
  double cxx[nx], cyy[ny]; 
  double val, err, avgerr, maxerr;
  double time_start;
  double time_now;
  double time_last;

  printf("abc_eq.c running \n");
  printf("Given  a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +\n");
  printf("       d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y) = c(x,y) \n");
  printf("where, for this test case \n");
  printf("  a1(x,y) = exp(x/2)*exp(y)/2 \n");
  printf("  b1(x,y) = 0.7/(x*x*y*y + 0.5) \n");
  printf("  c1(x,y) = (4 - exp(x) - exp(y/2))*2 \n");
  printf("  d1(x,y) = x*x+y \n");
  printf("  e1(x,y) = x*y*y \n");
  printf("  f1(x,y) = 3.0*x + 2.0*y \n");
  printf(" \n");
  printf("  c(x,y) = 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)            + \n");
  printf("     0.7*(6.0*x + 8.0*y + 5.0)/(x*x*y*y+0.5)              + \n");
  printf("     (8.0 - 2.0*exp(x) - 2.0*exp(y/2.0))*(12.0*y + 8.0*x) + \n");
  printf("     (x*x+y)*(3.0*x*x + 6.0*x*y + 4.0*y*y + 5.0*y + 6.0)  + \n");
  printf("     x*y*y*(6.0*y*y + 3.0*x*x + 8.0*x*y + 5.0*x +7.0)     + \n");
  printf("     3.0*x + 2.0*y; \n");
  printf(" \n");
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("Analytic solution \n");
  printf(" u(x,y) = x^3 + 2y^3 + 3x^2 y + 4xy^2 + 5xy + 6x + 7y + 8 \n");
  printf(" \n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  printf("nx=%d, ny=%d \n", nx, ny);
  time_start = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time start = %f seconds \n", time_start);
  time_last = time_start;

  printf("x grid and analytic solution at ymin \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    Ua[i] = u(xg[i],ymin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f)=%7.4f \n",
            i, xg[i], ymin, Ua[i]);
  } /* i */  
  printf("\n");
  printf("y grid and analytic solution at xmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+(double)ii*hy;
    Ua[ii] = u(xmin, yg[ii]); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f)=%7.4f \n",
            ii, xmin, yg[ii], Ua[ii]);
  } /* ii */
  printf("\n");

  printf("check rderiv(2,nx,0,hx,cxx \n");
  rderiv(2,nx,0,hx,cxx);
  for(i=0; i<nx; i++) printf("cxx[%d]=%f \n", i, cxx[i]);
  printf("check rderiv(2,nx,nx-1,hx,cxx \n");
  rderiv(2,nx,nx-1,hx,cxx);
  for(i=0; i<nx; i++) printf("cxx[%d]=%f \n", i, cxx[i]);
  printf("\n");

  printf("put solution in Ua vector \n");
  /* put solution in Ua vector */
  for(i=1; i<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      Ua[s(i,ii)] = u(x,y);
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("initialize A \n");
  /* initialize A */
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(j=1; j<nx-1; j++)
      {
        for(jj=1; jj<ny-1; jj++)
        {
          A[ss(i,ii,j,jj)] = 0.0;
	} /* jj */
      } /* j */
    } /* ii */
  } /* i */

  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  /* compute entries in A matrix, inside boundary */
  printf("compute A matrix \n");
  cs = (nx-2)*(ny-2); /* constant RHS column */
  val = 0.0;
  for(i=1; i<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      k = s(i,ii); /* row index */

      /* for each term in first equation */

      /* for a1(x,y)*d^2u/dx^2 */
      coefdxx = a1(x,y);
      rderiv(2,nx,i,hx,cxx);
      for(j=0; j<nx; j++)
      {
        val = coefdxx*cxx[j];
        if(j==0 || j==nx-1) A[k*(nxy+1)+cs] -= val*u(xg[j],yg[ii]);
        else A[sk(k,j,ii)] += val;
      }

      /* for b1(x,y)*d^2u/dx dy */
      coefdxy = b1(x,y);
      rderiv(1,nx,i,hx,cx);
      rderiv(1,ny,ii,hy,cy);
      for(j=0; j<nx; j++)
      {
        for(jj=0; jj<ny; jj++)
        {
          val = coefdxy*cx[j]*cy[jj];
          if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
	       A[k*(nxy+1)+cs] -= val*u(xg[j],yg[jj]);
          else A[sk(k,j,jj)] += val;
	} /* end jj */
      } /* end j */

      /* for c1(x,y)*d^2u/dy^2 */
      coefdyy = c1(x,y);
      rderiv(2,ny,ii,hy,cyy);
      for(jj=0; jj<ny; jj++)
      {
        val = coefdyy*cyy[jj];
        if(jj==0 || jj==ny-1) A[k*(nxy+1)+cs] -= val*u(xg[i],yg[jj]);
        else A[sk(k,i,jj)] += val;
      }

      /* for d1(x,y)*du/dx */
      coefdx = d1(x,y);
      rderiv(1,nx,i,hx,cx);
      for(j=0; j<nx; j++)
      {
        val = coefdx*cx[j];
        if(j==0 || j==nx-1) A[k*(nxy+1)+cs] -= val*u(xg[j],yg[ii]);
        else A[sk(k,j,ii)] += val;
      }

      /* for e1(x,y)*du/dy */
      coefdy = e1(x,y);
      rderiv(1,ny,ii,hy,cy);
      for(jj=0; jj<ny; jj++)
      {
        val = coefdy*cy[jj];
        if(jj==0 || jj==ny-1) A[k*(nxy+1)+cs] -= val*u(xg[i],yg[jj]);
        else A[sk(k,i,jj)] += val;
      }

      /* for f1(x,y)*u(x,y) */
      coefu = f1(x,y);
      A[sk(k,i,ii)] += coefu;

      /* for f(x,y) RHS constant */
      A[k*(nxy+1)+cs] += c(xg[i],yg[ii]);

      /* end terms */

    } /* end ii loop */
  } /* end i loop */

  /* printA(); lots of output */

  printf("A computed stiffness matrix \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("solve %d equations in %d unknowns \n", nxy, nxy);

  simeqb(nxy, A, U);

  printf("equations solved \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");
  
  printf("check PDE against known solution\n");
  check_soln(Ua);
  printf("check PDE against computer solution\n");
  check_soln(U);
  printf("\n");

  avgerr = 0.0;
  maxerr = 0.0;
  printf("          U computed, Ua analytic, error \n");
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      err = abs(U[s(i,ii)]-Ua[s(i,ii)]);
      if(err>maxerr) maxerr = err;
      avgerr = avgerr + err;
      printf("ug[%d,%d]=%9.4f, Ua=%9.4f, err=%g \n",
             i, ii, U[s(i,ii)], Ua[s(i,ii)], err);
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  printf("\n");


  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nx*ny));
  printf("\n");
  printf("end abc_eq.c \n");
  return 0;
} /* end main */

/* end abc_eq.c */