/* pde_bihar44tl_eq.c  discretization solution biharmonic 4th order 4D 
 *                     using pthreads
 * solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + utttt(x,y,z,t) + 
 *        2*uxxyy(x,y,z,t) + 2*uxxzz(x,y,z,t)  + 2*uxxtt(x,y,z,t)  + 
 *        2*uyyzz(x,y,z,t) + 2*uyytt(x,y,z,t)  + 2*uzztt(x,y,z,t)  - 
 *        16*u(x,y,z,t) = f(x,y,z,t)
 *
 *  f(x,y,z,t) := 0.0
 *
 *  in the hyper cube  xmin< x <xmax and  ymin< y <ymax  ,
 *                     zmin< z <zmax and  tmin< t <tmax 
 *
 *     subject to boundary conditions of the first
 *     kind, Dirichlet
 *
 *  Dirichlet boundary conditions are computed using:
 *  ub(x,y,z,t) := sin(x+y+z+t)
 *
 *  uses   tsimeq.c (threads)  nderiv.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include "tsimeq.h"
#include "lderiv.h"

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


#define nx  12
#define ny  12
#define nz  12
#define nt  12
/* DOF not including boundary */
#define nxyzt (nx-2)*(ny-2)*(nz-2)*(nt-2)

static int debug = 0;
static double xmin = 0.0; /* problem parameters */
static double xmax = 1.57; 
static double ymin = 0.0;
static double ymax = 1.57;
static double zmin = 0.0;
static double zmax = 1.57;
static double tmin = 0.0;
static double tmax = 1.57;

/* global nx, hx, xg, ny, hy, yg, nz, hz, zg, nt, ht, tg  needed by functions*/
static double xg[nx]; /* X grid */
static double yg[ny]; /* Y grid */
static double zg[ny]; /* z grid */
static double tg[ny]; /* t grid */
static double hx, hy, hz, ht;
static double ug[nxyzt];          /* solution, does not include boundary */
static double kg[nxyzt][nxyzt+1]; /* simultaneous equations, incl RHS */ 
static double maxerr;

static double cxxxx[nx][nx];
static double cyyyy[ny][ny];
static double czzzz[nz][nz];
static double ctttt[nt][nt];
static double cxx[nx][nx];
static double cyy[ny][ny];
static double czz[nz][nz];
static double ctt[nt][nt];
/* cxxyy etc. built dynamically */

/* function prototypes */
static double F(double x, double y, double z, double t); /* RHS */
static double ub(double x, double y, double z, double t);
static void discretize(int i, int ii, int iii, int iiii);
static void dsetup();
static void write_soln();
static void check_soln();
static int s(int i, int ii, int iii, int iiii) /* index 1..nx-1  etc */
{
  return (i-1) + (ii-1)*(nx-2) + (iii-1)*(nx-2)*(ny-2) +
         (iiii-1)*(nx-2)*(ny-2)*(nz-2);
} /* end s */

int main(int argc, char * argv[])
{
  double x, y, z, t, u;
  int i, ii, iii, iiii, j, jj, jjj, jjjj;
  int i4; /* s(i,ii,iii,iiii) */
  int j4; /* s(j,jj,jjj,jjjj) */
  double val, err, avgerr, tstart, tnext, tnow;
  int stat;

  tstart = (double)time(NULL);
  printf("pde_bihar44tl_eq.c running.\n");
  printf("solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) +\n");
  printf("       uzzzz(x,y,z,t) + utttt(x,y,z,t) +\n"); 
  printf("     2*uxxyy(x,y,z,t) + 2*uxxzz(x,y,z,t)  + 2*uxxtt(x,y,z,t)  +\n"); 
  printf("     2*uyyzz(x,y,z,t) + 2*uyytt(x,y,z,t)  + 2*uzztt(x,y,z,t)  -\n"); 
  printf("     16*u(x,y,z,t) = f(x,y,z,t) \n");
  printf(" \n");
  printf("f(x,y,z,t) := 0.0\n");
  printf(" \n");
  printf("in the hyper cube   xmin< x <xmax  and  ymin< y <ymax,\n");
  printf("                    zmin< z <zmax  and  tmin< t <tmax,\n");
  printf(" \n");
  printf("   subject to boundary conditions of the first kind, Dirichlet\n");
  printf("   ub(x,y,z,t) := sin(x+y+z+t) \n");
  printf("xmin=%7.3f, xmax=%7.3f, ymin=%7.3f, ymax=%7.3f\n",
         xmin, xmax, ymin, ymax);
  printf("zmin=%7.3f, zmax=%7.3f, tmin=%7.3f, tmax=%7.3f\n",
         zmin, zmax, tmin, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n",nx, ny, nz, nt);
  hx = (xmax-xmin)/(double)(nx-1);
  printf("x grid  hx=%f \n", hx);
  for(i=0; i<nx; i++) /* does not need to be uniform */
  {
    xg[i] = xmin+i*hx;
    printf("xg[%d]=%f\n", i, xg[i]);
  }  
  printf(" \n");
  hy = (ymax-ymin)/(double)(ny-1);
  printf("y grid hy=%f \n", hy);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
    printf("yg[%d]=%f\n", ii, yg[ii]);
  }  
  printf(" \n");
  hz = (zmax-zmin)/(double)(nz-1);
  printf("z grid hz=%f \n", hz);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+iii*hz;
    printf("zg[%d]=%f\n", iii, zg[iii]);
  }  
  printf(" \n");
  ht = (tmax-tmin)/(double)(nt-1);
  printf("t grid  ht=%f \n", ht);
  for(iiii=0; iiii<nt; iiii++)
  {
    tg[iiii] = tmin+iiii*ht;
    printf("tg[%d]=%f\n", iiii, tg[iiii]);
  }  
  printf(" \n");

  /* initialize kg  */ 
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<nx-1; j++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(jj=1; jj<ny-1; jj++)
        {
          for(iii=1; iii<nz-1; iii++)
          {
            for(jjj=1; jjj<nz-1; jjj++)
            {
              for(iiii=1; iiii<nt-1; iiii++)
              {
                for(jjjj=1; jjjj<nt-1; jjjj++)
                {
                  kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjj)] = 0.0;
		}
	      }
	    }
	  }
        }
      }
    }
  }
  
  dsetup(); /* setup discrete coefficients */
 
  /* compute entries in stiffness, equation matrix */ 
  printf("compute stiffness, equation matrix \n");
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          /* compute discrete derivatives for s(i,ii,iii,iiii) equation */  
          discretize(i, ii, iii, iiii);
                        /* sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] */
        } /* end iiii loop */
      } /* end iii loop */
    } /* end ii loop */
  } /* end i loop */

  printf("k computed discrete matrix \n");
  printf("f computed forcing function, RHS \n");
  tnext = (double)time(NULL);
  printf("time for building equations = %f seconds \n", tnext-tstart);

  /* threaded parallel simeq, nx-2 threads */
  tsimeqb(nx-2, nxyzt, (double*)kg, ug); /* B is kg fg */
  tnow = (double)time(NULL);
  printf("time for solve system of equations= %f seconds \n", tnow-tnext);
  tnext = tnow;

  write_soln(); /* to file x y z t u=computed solution */

  check_soln(); /* against PDE */
  /* independent of knowing analytic solution */

  avgerr = 0.0;
  maxerr = 0.0;
  printf("ug computed by discretization, u analytic, error \n");
  for(i=1; i<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      for(iii=1; iii<nz-1; iii++)
      {
        z = zg[iii];
        for(iiii=1; iiii<nt-1; iiii++)
        {
          t = tg[iiii];
          u = ub(x,y,z,t);
          i4 = s(i,ii,iii,iiii);
          err = abs(ug[i4]-u);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          printf("ug[%d,%d,%d,%d]=%f, u=%f, err=%g \n",
                 i, ii, iii, iiii, ug[i4], u, ug[i4]-u);
        }
      }
    }
  }
  printf("xmax=%f, ymax=%f, xmax=%f, tmax=%f \n", xmax, ymax, zmax, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  printf("                  maxerr=%g, avgerr=%g \n",
         maxerr,avgerr/(double)(nxyzt));
  printf("\n");
  fflush(stdout);
  printf("use plot4d_gl to see results \n");

  tnow = (double)time(NULL);
  printf("time for checking and writing = %f seconds \n", tnow-tnext);
  printf("total time = %f seconds \n", tnow-tstart);
  return 0;
} /* end pde_bihar44tl_eq main */

/* PDE problem definition functions  f and ub */

static double F(double x, double y, double z, double t)
{
  return 0.0;

}

static double ub(double x, double y, double z, double t)
{
  return sin(x+y+z+t); /* boundary, and used here for checking */
}

static void discretize(int i, int ii, int iii, int iiii)
                   /* sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)...] */
                   /* this code unique to this specific PDE */
{
  int cs = nxyzt;        /* RHS index */
  int i4;
  int j, jj, jjj, jjjj;
  double val;

  i4 = s(i,ii,iii,iiii); /* equation index */
  /* for each term in each equation */

  /* for uxxxx  d^4u/dx^4 */
  for(j=0; j<nx; j++)
  {
    val = cxxxx[i][j];
    if(j==0 || j==nx-1)
      kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[iii],tg[iiii]);
    else kg[i4][s(j,ii,iii,iiii)] += val;
  }

  /* for uyyyy  d^4u/dy^4 */
  for(jj=0; jj<ny; jj++)
  {
    val = cyyyy[ii][jj];
    if(jj==0 || jj==ny-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[jj],zg[iii],tg[iiii]);
    else kg[i4][s(i,jj,iii,iiii)] += val;
  }

  /* for uzzzz  d^4u/dz^4 */
  for(jjj=0; jjj<nz; jjj++)
  {
    val = cyyyy[iii][jjj];
    if(jjj==0 || jjj==nz-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[ii],zg[jjj],tg[iiii]);
    else kg[i4][s(i,ii,jjj,iiii)] += val;
  }

  /* for utttt  d^4u/dt^4 */
  for(jjjj=0; jjjj<nt; jjjj++)
  {
    val = cyyyy[iiii][jjjj];
    if(jjjj==0 || jjjj==nt-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[ii],zg[iii],tg[jjjj]);
    else kg[i4][s(i,ii,iii,jjjj)] += val;
  }

  /* for uxxyy  d^4u/dx^2dy^2 */
  for(j=0; j<nx; j++)
  {
    for(jj=0; jj<ny; jj++)
    {
      val = 2.0*cxx[i][j]*cyy[ii][jj];
      if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
	kg[i4][cs] -= val*ub(xg[j],yg[jj],zg[iii],tg[iiii]);
      else kg[i4][s(j,jj,iii,iiii)] += val;
    } /* end jj */
  } /* end j */

  /* for uxxzz  d^4u/dx^2dz^2 */
  for(j=0; j<nx; j++)
  {
    for(jjj=0; jjj<nz; jjj++)
    {
      val = 2.0*cxx[i][j]*czz[iii][jjj];
      if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[jjj],tg[iiii]);
      else kg[i4][s(j,ii,jjj,iiii)] += val;
    } /* end jjj */
  } /* end j */

  /* for uxxtt  d^4u/dx^2dt^2 */
  for(j=0; j<nx; j++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*cxx[i][j]*ctt[iiii][jjjj];
      if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[iii],tg[jjjj]);
      else kg[i4][s(j,ii,iii,jjjj)] += val;
    } /* end jjjj */
  } /* end j */

  /* for uyyzz  d^4u/dy^2dz^2 */
  for(jj=0; jj<ny; jj++)
  {
    for(jjj=0; jjj<nz; jjj++)
    {
      val = 2.0*cyy[ii][jj]*czz[iii][jjj];
      if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
	kg[i4][cs] -= val*ub(xg[i],yg[jj],zg[jjj],tg[iiii]);
      else kg[i4][s(i,jj,jjj,iiii)] += val;
    } /* end jjj */
  } /* end jj */

  /* for uyytt  d^4u/dy^2dt^2 */
  for(jj=0; jj<ny; jj++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*cyy[ii][jj]*ctt[iiii][jjjj];
      if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[i],yg[jj],zg[iii],tg[jjjj]);
      else kg[i4][s(i,jj,iii,jjjj)] += val;
    } /* end jjjj */
  } /* end jj */

  /* for uzztt  d^4u/dz^2dt^2 */
  for(jjj=0; jjj<nz; jjj++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*czz[iii][jjj]*ctt[iiii][jjjj];
      if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[i],yg[ii],zg[jjj],tg[jjjj]);
      else kg[i4][s(i,ii,jjj,jjjj)] += val;
    } /* end jjjj */
  } /* end jjj */

  kg[i4][i4] -= 16.0; /* U */ 
} /* end discretize */

static void dsetup()
{
  int i, ii, iii, iiii, j;

  /* each such as cxxxx[i][j] is j_th coefficient at DOF i */
  for(i=0; i<nx; i++) rderiv(4, nx, i, hx, cxxxx[i]);
  if(debug)
  {
    printf("debug print of cxxxx[i][j] \n");
    for(i=0; i<nx; i++)
    {
      printf("i=%d ",i);
      for(j=0; j<nx; j++) printf("%f ", cxxxx[i][j]);
      printf(" \n");
    }
  }
  for(i=0; i<nx; i++) rderiv(2, nx, i, hx, cxx[i]);
  if(debug)
  {
    printf("debug print of cxx[i][j] \n");
    for(i=0; i<nx; i++)
    {
      printf("i=%d ",i);
      for(j=0; j<nx; j++) printf("%f ", cxx[i][j]);
      printf(" \n");
    }
  }
  for(ii=0; ii<ny; ii++) rderiv(4, ny, ii, hy, cyyyy[ii]);
  for(ii=0; ii<ny; ii++) rderiv(2, ny, ii, hy, cyy[ii]);
  for(iii=0; iii<nz; iii++) rderiv(4, nz, iii, hz, czzzz[iii]);
  for(iii=0; iii<nz; iii++) rderiv(2, nz, iii, hz, czz[iii]);
  for(iiii=0; iiii<nt; iiii++) rderiv(4, nt, iiii, ht, ctttt[iiii]);
  for(iiii=0; iiii<nt; iiii++) rderiv(2, nt, iiii, ht, ctt[iiii]);
} /* end dsetup */

static void write_soln()
{
  FILE * fout;
  double x, y, z, t, u;
  int i4;
  int i, ii, iii, iiii;

  fout = fopen("pde_bihar44tl_eq_c.dat", "w");
  if(fout==NULL)
  {
    printf("file pde_bihar44tl_eq_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing pde_bihar44tl_eq_c.dat \n"); /* including boundary */
  fprintf(fout,"4 %d %d %d %d \n", nx, ny, nz, nt);
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          if(i==0 || i==nx-1 || ii==0 || ii==ny-1 ||
             iii==0 || iii==nz-1 || iiii==0 || iiii==nt-1)
	  {
            u = ub(x,y,z,t);
	  }
          else
	  {
            i4 = s(i,ii,iii,iiii);
            u = ug[i4];
	  }
          fprintf(fout," %f %f %f %f %f \n",x,y,z,t,u);
	} /* end t */
      } /* end z */
    } /* end y */
  } /* end x */
  fprintf(fout," \n");
  fflush(fout);
  fclose(fout);
  printf("finished writing pde_bihar44tl_eq_c.dat \n");
} /* end write_soln */

static void check_soln() /* against PDE */
{
  /* for all interior points */
  double x, y, z, t, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr, eval_deriv;
  int cnt = 0;
  int i, ii, iii, iiii, j, jj, jjj, jjjj;

  printf("check_soln values against PDE \n");
  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* check DOF */
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      for(iii=1; iii<nz-1; iii++)
      {
        z = zg[iii];
        for(iiii=1; iiii<nt-1; iiii++)
        {
          t = tg[iiii];
          eval_deriv = 0.0;

          /* PDE uxxxx term */          
          for(j=0; j<nx; j++)
          {
            if(j==0 || j==nx-1)
	      eval_deriv += cxxxx[i][j]*ub(xg[j],y,z,t);
            else
	      eval_deriv += cxxxx[i][j]*ug[s(j,ii,iii,iiii)];
          } 
          /* PDE uyyyy term */
          for(jj=0; jj<ny; jj++)
          {
            if(jj==0 || jj==ny-1)
	      eval_deriv += cyyyy[ii][jj]*ub(x,yg[jj],z,t);
	    else
	      eval_deriv += cyyyy[ii][jj]*ug[s(i,jj,iii,iiii)];
          }
          /* PDE uzzzz term */
          for(jjj=0; jjj<nz; jjj++)
          {
            if(jjj==0 || jjj==nz-1)
	      eval_deriv += czzzz[iii][jjj]*ub(x,y,zg[jjj],t);
	    else
	      eval_deriv += czzzz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
          } 
          /* PDE utttt term */
          for(jjjj=0; jjjj<nt; jjjj++)
          {
            if(jjjj==0 || jjjj==nt-1)
	      eval_deriv += ctttt[iiii][jjjj]*ub(x,y,z,tg[jjjj]);
	    else
	      eval_deriv += ctttt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
          }

          /* PDE  2*uxxyy term */
          for(j=0; j<nx; j++)
          {
            for(jj=0; jj<ny; jj++)
            {
              if(j==0 || j==nx-1 || jj==0 || jj==ny-1) /* boundary, not DOF */
                eval_deriv += 2.0*cxx[i][j]*cyy[ii][jj]*ub(xg[j],yg[jj],z,t);
	      else
		eval_deriv += 2.0*cxx[i][j]*cyy[ii][jj]*ug[s(j,jj,iii,iiii)];
            } /* end jj */
          } /* end j */

          /* PDE  2*uxxzz term */
          for(j=0; j<nx; j++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
                eval_deriv += 2.0*cxx[i][j]*czz[iii][jjj]*ub(xg[j],y,zg[jjj],t);
	      else
                eval_deriv += 2.0*cxx[i][j]*czz[iii][jjj]*ug[s(j,ii,jjj,iiii)];
            } /* end jjj */
          } /* end j */

          /* PDE  2*uxxtt term */
          for(j=0; j<nx; j++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
                eval_deriv += 2.0*cxx[i][j]*ctt[iiii][jjjj]*ub(xg[j],y,z,tg[jjjj]);
	      else
                eval_deriv += 2.0*cxx[i][j]*ctt[iiii][jjjj]*ug[s(j,ii,iii,jjjj)];
            } /* end jjjj */
          } /* end j */

          /* PDE  2*uyyzz term */
          for(jj=0; jj<ny; jj++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
                eval_deriv += 2.0*cyy[ii][jj]*czz[iii][jjj]*ub(x,yg[jj],zg[jjj],t);
	      else
                eval_deriv += 2.0*cyy[ii][jj]*czz[iii][jjj]*ug[s(i,jj,jjj,iiii)];
            } /* end jjj */
          } /* end jj */

          /* PDE  2*uyytt term */
          for(jj=0; jj<ny; jj++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
                eval_deriv += 2.0*cyy[ii][jj]*ctt[iiii][jjjj]*ub(x,yg[jj],z,tg[jjjj]);
	      else
                eval_deriv += 2.0*cyy[ii][jj]*ctt[iiii][jjjj]*ug[s(i,jj,iii,jjjj)];
            } /* end jjjj */
          } /* end jj */

          /* PDE  2*uzztt term */
          for(jjj=0; jjj<nz; jjj++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
                eval_deriv += 2.0*czz[iii][jjj]*ctt[iiii][jjjj]*ub(x,y,zg[jjj],tg[jjjj]);
	      else
                eval_deriv += 2.0*czz[iii][jjj]*ctt[iiii][jjjj]*ug[s(i,ii,jjj,jjjj)];
            } /* end jjjj */
          } /* end jjj */

          /* PDE -16 * u(x,y,z,t) term */
          eval_deriv -= 16.0*ug[s(i,ii,iii,iiii)];

          err = eval_deriv-F(x,y,z,t); /* subtract RHS */
          err = abs(err);
          sumerr += err;
          sumsqerr += err*err;
          maxerr = max(err,maxerr);
          cnt++;
	} /* end t */
      } /* end z */ 
    } /* end y */
  } /* end x */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("check solution against PDE \n");
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 
/* end pde_bihar44tl_eq.c */