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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpf_nuderivd.h"
#define  digits 100
#include "simeq.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

#define nx 27  /* including boundary at 0 and nx-1 */
static double xmin = 0.0;
static double pi = 3.1415926535897932384626433832795028841971;
static double xmax; /*  = pi * 9.0/2.0; */
static double hx;

#define neqn (nx-2) /* DOF number of equations */
#define cs neqn            /* RHS */

static double xg[nx];
/* nx-2   internal, non boundary cells */
static double cxx[nx];
static double cx[nx];
static double us[neqn]; /* solution vector */
static double ut[neqn*(neqn+1)]; /* includes RHS */

static double U(double x)
{               /* solution for(boundary and optional check */
  return sin(x)/4.0 + cos(x)/4.0;
} /* end U */

static double c(double x)
{               /* RHS of differential equation to solve */
  return -sin(x)/2.0 + cos(x);
} /*  end c */

static int S(int i) /* zero based index includes boundary */
{                           /* this is row of ut[] */
  return (i-1); /* index into us and ut, RHS, no boundary */
} /* end S */

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

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

static void init_matrix()
{
  int i, ii, k;

  /* cs is RHS, constant column subscript */
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      ut[S(i)*(neqn+1)+S(ii)] = 0.0;
    }
    ut[S(i)*(neqn+1)+cs] = 0.0;
  }

  printf("internal cells zeroed \n");
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    ii = S(i); /* equation index */
    /* for each term */

    /* d^2 u(x)/dx^2 */
    mpf_nuderivd(digits, 2,nx,i,xg,cxx);
    for(k=0; k<nx; k++)
    {
      if(k==0 || k==nx-1)
      {
        ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] - cxx[k]*U(xg[k]);
      }	 
      else
      {
        ut[ii*(neqn+1)+S(k)] = ut[ii*(neqn+1)+S(k)] + cxx[k];
      }
    } /* k */

    /* 3 d u(x)/dx */
    mpf_nuderivd(digits, 1,nx,i,xg,cx);
    for(k=0; k<nx; k++)
    {
      if(k==0 || k==nx-1)
      {
        ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] - 3.0*cx[k]*U(xg[k]);
      }
      else
      {
        ut[ii*(neqn+1)+S(k)] = ut[ii*(neqn+1)+S(k)] + 3.0*cx[k];
      }
    } /* k */
    /* 2 u(x) */
    ut[ii*(neqn+1)+ii] = ut[ii*(neqn+1)+ii] + 2.0; /* ii=S(i) */

    /* RHS constant */
    ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] + c(xg[i]);
    /* end terms */
  } /* i */
} /* end init_matrix */

static void print_B()
{
  int i, ii;

  printf("matrix B includes RHS\n");
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      printf("i=%2d, ii=%2d, ut=%g\n",
	     i, ii ,ut[S(i)*(neqn+1)+S(ii)]);
    } /* ii */
    printf("RHS=          %g\n", ut[S(i)*(neqn+1)+cs]);
  } /* i */
  printf("\n");
  printf("\n");
} /*  end print_B */

static void check_rows() /* only for(checking */
{	/* row is in ut solution coordinates */
  int i, ii, row, col;
  double sum;

  for(i=1; i<nx-1; i++)
  {
    row = S(i);
    sum = 0.0;
    for(ii=1; ii<nx-1; ii++)
    {
      col = S(ii);
      sum = sum + ut[row*(neqn+1)+col]*U(xg[ii]); /* u is check solution */
    } /* ii */
    sum = sum - ut[row*(neqn+1)+cs];
    if(abs(sum)>0.001)
    {
      printf("row i=%2d, is bad err=%g\n", i, sum);
    }
  } /*i */
} /* end check_rows */

static void check_soln()  /* check when solution unknown */
{
  double smaxerr, err, x, uxx, ux;
  int i, j, ii, jj, k;

  /*  uxx(x) + 3 ux(x) + 2 u(x) = c(x) */
  printf("check_soln against PDE\n");
  smaxerr = 0.0;

  for(i=1; i<nx-1; i++)  /* through dof equations */
  {
    mpf_nuderivd(digits, 2,nx,i,xg,cxx);
    x = xg[i];
    mpf_nuderivd(digits, 1,nx,i,xg,cx);
    err = -c(x); /* add other side of equation terms */
    uxx = 0.0;
    for(k=0; k<nx; k++) /* compute numerical derivative at (i) */
    {
      if(k==0 || k==nx-1)
      {
        uxx = uxx + cxx[k]*U(xg[k]); /* boundary */
      }
      else
      {
        uxx = uxx + cxx[k]*us[S(k)];
      }
    }
    ux = 0.0;
    for(k=0; k<nx; k++)
    {
      if(k==0 || k==nx-1)
      {
        ux = ux + cx[k]*U(xg[i]); /* boundary */
      }
      else
      {
        ux = ux + cx[k]*us[S(i)];
      }
    }
    err = err + uxx + 3.0*ux + 2.0*U(x);
    if(abs(err) > smaxerr) smaxerr = abs(err);
  } /* i x */

  printf("check soln against PDE max error=%g\n", smaxerr);
} /* end Check_Soln */

/* write_soln2  for gnuplot */
static void  write_soln2(double Ux[], char file_name[])
{
  FILE * Fout;
  int i;

  Fout = fopen(file_name, "w");
  if(Fout==NULL)
  {
    printf("ERROR unable to open %s for writing \n", file_name);
    exit(1);
  }
  printf("writing %s\n", file_name);
  for(i=0; i<nx; i++)
  {
    if(i==0 || i==nx-1)
    {
      fprintf(Fout,"%f %f \n",
              xg[i], U(xg[i])); /* boundary */
    }
    else
    { 
      fprintf(Fout,"%f %f \n",
              xg[i], Ux[S(i)]); /* computed solution */
    }
  } /* i */
  fprintf(Fout,"\n");
  fclose(Fout);
  printf("finished writing %s\n", file_name);
}  /* end write_soln2 */

int main(int argc, char * argv[])
{
  int i;

  printf("pde12_eq.c running\n");
  printf("differential equation to solve\n");
  printf("d^2u/dx^2 + 3 du/dx + 2 u(x) = c(x,y) = -sin(x)/2 + cos(x)\n");
  printf("uniform grid on  0 to 9/2 Pi\n");
  printf("known Solution, for testing method\n");
  printf("u(x) = sin(x)/4 + cos(x)/4 \n");
  printf("\n");
  xmax = pi * 9.0/2.0;
  hx = (xmax-xmin)/(double)(nx-1);

  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    printf("xg(%2d)=%f\n", i, xg[i]);
  }
  printf("xmin=%f, xmax=%f, hx=%f, nx=%2d\n", xmin, xmax, hx, nx);
  printf("u(0.5)=%f\n", U(0.5));
  printf("c(0.5)=%f\n", c(0.5));
  printf("\n");

  init_matrix();
  printf("matrix initialized\n");
  /* print_B(); */
  check_rows();
  simeqb(neqn, ut, us);
  check_soln();
  printf("\n");
  print();
  write_soln2(us, "pde12_eq_c.dat");
  printf("finished pde12_eq.c\n");
  return 0;
} /* end pde12_eq.c */