/* pde2sin_sparse.c    2d uniform grid boundary value PDE */
/*                 known solution for testing method  */
/*                 u(x,y) = sin(x-y)                  */
/*                                                    */
/*                 differential equation to solve     */
/*                 du/dx + du/dy=0                    */
/* uniform grid on rectangle 0 to 2Pi, 0 to 2Pi       */
/* set up and solve system of linear equations        */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "deriv.h"
#include "sparse.h"
#undef abs
#define abs(a) ((a)<0?(-(a)):(a))

/* do NOT make nx=ny, this creates a singular matrix  */
/* including boundary at 0 and nx-1, 0 and ny-1    */
#define nx 13 
#define ny 12
#define neqn ((nx-2)*(ny-2))
#define cs neqn
static double pi = 3.1415926535897932384626433832795028841971;
static double xmin, xmax, ymin, ymax, hx, hy;

static double xg[nx];
static double yg[ny];
/* nx-2 by ny-2  internal, boundary cells */
static double cx[nx];
static double cy[ny];
static double coefdx, coefdy;
static double us[neqn];  /* solution vector */
/* static double ut[neqn][neqn+1];  includes RHS, now B,eqn,col */
static struct sparse B;

static double u(double x, double y)
{
  /* solution for boundary and optionalcheck */
    return sin(x-y);
} /* end u */

static double c(double x, double y)
{
  /* RHS of differential equation to solve */
  return 0.0;
} /* end c */

static int S(int i, int j)
{
  if(i<1 || i>nx-1 || j<1 || j>ny-1) printf("BAD S i=%2d, j=%2d\n", i, j);
  return (i-1)*(ny-2)+(j-1); /* equation index into us and B */
} /* end S */

static void print()
{
  int i, j;
  double error = 0.0;
  double max_error = 0.0;
  double avg_error = 0.0;
  double val;

  printf("\n");
  printf("exact solution u, computed us, error\n");
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      printf("xg[%2d]=%7.3f, ", i, xg[i]);
      printf("yg[%2d]=%7.3f, ", j, yg[j]);
      printf("u=%7.3f, ", u(xg[i],yg[j]));
      if(i==0 || i==nx-1 || j==0 || j==ny-1) /* boundary */
      {
        val = u(xg[i],yg[j]);
        error = 0.0;
      }
      else
      {
        val = us[S(i,j)];
        error = abs(val-u(xg[i],yg[j]));
        avg_error = avg_error + error;
      }
      if(error > max_error) max_error = error;
      printf("us=%g, ", val);
      printf("err=%g\n", error);
    }
  }
  printf("avg_error=%g, max_error=%g\n", avg_error/(double)(neqn),
         max_error);
} /* end print */

static void check_row(int row) /* only for checking */
{
  /* row is in B solution coordinates */
  int i, j, ii;
  double sum;

  sum = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      ii = S(i,j);
      sum = sum + sparseGet(B,row,ii)*u(xg[i],yg[j]); /* u is check solution */
           /* was ut[row][ii] */
    }
  }
  sum = sum - sparseGet(B,row,cs); /* was ut[row][cs]; */
  if(abs(sum)>0.001)
    printf("BAD row=%4d, err=%g\n", row, sum);
} /* end check_row */

static void init_matrix()
{
  int i, j, ii, jj;
  int k, kj, ik;
  /* cs is RHS, constant column subscript */
  /* zero internal cells, not with sparse */

  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = S(i,j); /* equation */
      /* for each term */

      /* coefdx * d u(x,y)/dx */
      coefdx = 1.0;
      rderiv(1,nx,i,hx,cx);
      for(k=0; k<nx; k++)
      {
        cx[k] = cx[k] * coefdx;
        if(k==0 || k==nx-1)
	{
          sparseAdd(B, ii, cs, - cx[k]*u(xg[k],yg[j]));
	}
        else
	{
          kj = S(k,j);
          sparseAdd(B, ii, kj, cx[k]);
        }
      } /* k */
      /* coefdy * d u(x,y)/dy */
      coefdy = 1.0;
      rderiv(1,ny,j,hy,cy);
      for(k=0; k<ny; k++)
      {
        cy[k] = cy[k] * coefdy;
        if(k==0 || k==ny-1)
	{
          sparseAdd(B, ii, cs, - cy[k]*u(xg[i],yg[k]));
	}
        else
	{
          ik = S(i,k);
          sparseAdd(B, ii, ik, cy[k]);
        }
      } /* k */
      /* RHS constant */
      /* ut[ii][cs] = ut[ii][cs] + c(xg[i],yg[j]); */
      /* end terms */
      check_row(ii);
    } /* j */
  } /* i */
} /* end init_matrix */

static void check_soln()  /* check when solution unknown */
{
  int i, ii, j, jj;
  double smaxerr, err, x, y, ux, uy;
  /*  ux(x,y) + uy(x,y) = c(x,y) */
  printf("check_soln against PDE\n");
  smaxerr = 0.0;

  for(i=1; i<nx-1; i++) /* through dof equations */
  {
    rderiv(1,nx,i,hx,cx);
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      rderiv(1,ny,ii,hy,cy);
      y = yg[ii];
      err = -c(x,y); /* add other side of equation terms */
      ux = 0.0;
      for(j=0; j<nx; j++) /* compute numerical derivative at (i,ii) */
      {
        if(j==0 || j==nx-1)
	  ux = ux + cx[j]*u(xg[j],yg[ii]); /* boundary */
        else
          ux = ux + cx[j]*us[S(j,ii)];
      }
      uy = 0.0;
      for(jj=0; jj<ny; jj++)
      {
        if(jj==0 || jj==ny-1)
	  uy = uy + cy[jj]*u(xg[i],yg[jj]); /* boundary */
        else
          uy = uy + cy[jj]*us[S(i,jj)];
      }
      err = err + ux + uy;
      if(abs(err) > smaxerr) smaxerr = abs(err);
    } /* ii Y */
  } /* i X */
  printf("check_soln max error=%g\n", smaxerr);
} /* end Check_Soln */

int main(int argc, char *argv[])
{
  int i, j;
  /* pde2sin_sparse  */
  printf("pde2sin_sparse.c   running\n");
  printf("differential equation to solve\n");
  printf("du/dx + du/dy = c(x,y)\n");
  printf("c(x,y)=0\n");
  printf("uniform grid on rectangle 0,2Pi to 0,2Pi\n");
  printf("known Solution, for testing method\n");
  printf("u(x,y) = sin(x-y)\n");
  printf("do not make nx=ny, this creates a singular matrix\n");
  xmin = 0.0;
  xmax = 2.0*pi;
  ymin = 0.0;
  ymax = 2.0*pi;
  hx = (xmax-xmin)/(double)(nx-1);
  hy = (ymax-ymin)/(double)(ny-1);

  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)(i)*hx;
    printf("xg[%2d]=%f\n", i, xg[i]);
  }
  for(j=0; j<ny; j++)
  {
    yg[j] = ymin+(double)(j)*hy;
    printf("yg[%2d]=%f\n", j, yg[j]);
  }
  printf("xmin=%7.3f, ", xmin);
  printf("xmax=%7.3f, ", xmax);
  printf("hx=%7.5f, ", hx);
  printf("nx=%3d\n", nx);
  printf("ymin=%7.3f, ", ymin);
  printf("ymax=%7.3f, ", ymax);
  printf("hy=%7.5f, ",hy);
  printf("ny=%3d\n", ny);
  printf("u(xg[0],yg[0])=%f\n", u(xg[0],yg[0]));
  printf("c(xg[0],yg[0])=%f\n", c(xg[0],yg[0]));

  /* debug stuff */
  printf("debug print\n");
  rderiv(1,nx,2,hx,cx);
  for(i=0; i<nx; i++) printf("rderiv 1,2,%2d cx[%2d]=%g\n", nx, i, cx[i]);
  rderiv(1,ny,2,hy,cy);
  for(i=0; i<ny; i++) printf("rderiv 1,2,%2d cy[%2d]=%g\n", ny, i, cy[i]);

  sparseInit(&B, neqn);

  init_matrix();
  printf("matrix initialized\n");

  /* debug stuff */
  printf("matrix B includes RHS\n");
  sparseWrite_All(B);

  sparseSimeq(B, us); /* returned us */
  check_soln();
  print();
  printf("pde2sin_sparse.c finished\n");
  return 0;
} /* end pde2sin_sparse.c */