// pde_sparse14.java  high order discretization,
// solve a4(x)*uxxxx(x) + a0(x)*u(x) = F(x)
//     F(x) = a4(x)*96  a0(x)*(4 x^4 + 2 x  + 1)
//  a4(x) = x+1      a0(x) = x^2 -1
// boundary conditions computed using u(x)
// with only two boundary conditions, see pde_14a_eq.java for 4
// analytic solution is u(x) = 4 x^4 + 2 x  + 1
//
// solve for Ui numerical approximation U(x_i)
// matrix solve simultaneous equations for U
// needs simeq and nderiv

class pde_sparse14
{
  final int nx = 10; // number of X coordinates
  // initialize k and f  including boundary 
  sparse kf = new sparse(nx);  
  double xg[] = new double[nx]; // uniform spacing not required, with nuderiv
  double cuxxxx[][] = new double[nx][nx]; // for each partial derivative in
  nderiv D = new nderiv(); //  has rderiv
  double  hx; // could use nuderiv
  double xmin = 0.2; // problem parameters 
  double xmax = 3.3;
  boolean debug = false;

  pde_sparse14()
  {
    double x;
    double val, err, avgerr, maxerr;

    double ug[] = new double[nx];
    double Ua[] = new double[nx];
    double cx[] = new double[nx];
    double t_start, t_init, t_kf, t_simeq, t_end;

    System.out.println("pde_sparse14.c running");
    System.out.println("Given a4(x)*uxxxx(x) + a0(x)*u(x) = F(x)");
    System.out.println("  F(x) = a4(x)*96  a0(x)*(4 x^4 + 2 x  + 1) ");
    System.out.println("  a4(x) = x+1      a0(x) = x^2 -1");
    System.out.println("with 2 boundary conditions, at least one Dirichlet");
    System.out.println("xmin+hx<x<xmax-hx  inside Boundaries");
    System.out.println("Analytic solution u(x) = 4 x^4 + 2 x + 1");
    System.out.println(" ");

    System.out.println("xmin="+xmin+", xmax="+xmax);
    System.out.println("nx="+nx);

    t_start = System.currentTimeMillis();
    System.out.println("x grid ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      System.out.println("xg["+i+"]="+xg[i]+", soln="+uana(xg[i]));
    } // i  
    System.out.println(" ");
    build_discrete(); // one time generation of cuxxxx

    for(int i=0; i<nx; i++) // for checking
    {
      Ua[i] = uana(xg[i]); // including boundary
    }

    // kf initially all zero (empty)

    kf.add(0,0,1.0);
    kf.add(0,nx,uana(xg[0])); // Dirichlet boundary
    kf.add(nx-1,nx-1,1.0);
    kf.add(nx-1,nx,uana(xg[nx-1])); // Dirichlet boundary

    System.out.println(" ");
    t_init = System.currentTimeMillis();

    System.out.println("discreteize stiffness matrix and RHS");
    for(int i=1; i<nx-1; i++) // not including boundary
    {
      // compute discrete term for kf[row=i]
      coef(i); // build one equation at a time  
    } // end i loop 


    System.out.println("k computed stiffness matrix");
    System.out.println("f computed forcing function");
    t_kf = System.currentTimeMillis();
    System.out.println("kf build took "+((t_kf-t_init)/1000.0)+" seconds");

    if(debug)
    {
      kf.write_all();
      System.out.println(" ");
    } 
    System.out.println("solve "+nx+" equations in "+nx+" unknowns");
    kf.simeq(ug);

    System.out.println("equations solved");
    t_simeq = System.currentTimeMillis();
    System.out.println("simeq took "+((t_simeq-t_kf)/1000.0)+" seconds");

    System.out.println(" ");
    System.out.println("check_soln on known solution");
    check_soln(Ua);
    System.out.println("check_soln based on numeric approximation of PDE");
    check_soln(ug);
    System.out.println(" ");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed,  Ua analytic,  error  ");
    for(int i=0; i<nx; i++)
    {
       err = Math.abs(ug[i]-Ua[i]);
       if(err>maxerr) maxerr = err;
       avgerr = avgerr + err;
       System.out.println("ug["+i+"]="+
                                ug[i]+
			 ",  Ua="+Ua[i]+
			 ",  err="+(ug[i]-Ua[i]));
    }
    System.out.println("                         maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx));
    t_end = System.currentTimeMillis();
    System.out.println("total took "+((t_end-t_start)/1000.0)+" seconds");

    System.out.println(" ");
  } // end pde_sparse14 constructor

  void  build_discrete()
  {
    // build cuxxxx[point][term]
    for(int i=0; i<nx; i++)
    {
      D.rderiv(4, nx, i, hx, cuxxxx[i]);
      if(debug)
      {
        for(int j=0; j<nx; j++)
	  System.out.println("cuxxxx["+i+"]["+j+"]= "+cuxxxx[i][j]);
        System.out.println(" ");
      }
    } 
    System.out.println(" ");
  } // end build_discrete

  void coef(int i) // row is equation, zero based, boundaries elsewhere
  {
    double val = 0.0;
    double valrhs = 0.0;
    int row = i;
    int rhs = nx;

    // RHS may need to be upadted  kf[row][rhs] = valrhs
    // evaluate function coefficients at xg[k]

    // basic RHS
    kf.add(row,rhs,F(xg[row])); // right hand side

    // for PDE with a0(x)*u(x) on diagonal, set a0(x)
    kf.add(row,row,a0(xg[row])) ; //  U(x) coefficient for this PDE

    for(int j=0; j<nx; j++)
    {
      if(j==0 || j==nx) // use boundary
      {
	kf.add(row,rhs,-a4(xg[row])*cuxxxx[row][j]*uana(xg[j]));
                                                   // known boundary
	                                           // or U[j] or Ua[j]
      }
      else
      {
	kf.add(row,j,a4(xg[row])*cuxxxx[row][j]);
      }
    } // end j
  } // end coef

  // PDE problem definition functions
  double F(double x)
                                           // forcing function, F(x)
  {
      return a4(x)*96.0 + a0(x)*(4.0*x*x*x*x + 2.0*x + 1.0);
  }

  double uana(double x)
                                      // analytic solution and boundaries
  {
    return 4.0*x*x*x*x + 2.0*x +1.0;
  }

    double u_prime(double x) // for Neumann boundaries
  {
    return 16.0*x*x*x + 2.0;
  }

  double a4(double x)
  {
      // return 1.0;
    return x+1.0;
  }

  double a0(double x)
  {
      // return 1.0;
    return x*x-1.0;
  }

  // check solution against PDE
  void check_soln(double U[])
  {
    // check a4(x)*uxxxx(x) + a0(x)*u(x) = F(x)
    double maxerr, err, sumerr, sumsqerr, rmserr, avgerr;
    double eval_deriv, x;
    int ii, cnt;

    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;

    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      // compute derivative coefficients in x direction

      eval_deriv = 0.0;

      for(int j=0; j<nx; j++)
      {
        eval_deriv += a4(x)*cuxxxx[i][j]*U[j];
      }
      eval_deriv += a0(x)*U[i];
      err = eval_deriv-F(x);  
      sumerr += Math.abs(err);
      sumsqerr += err*err;
      maxerr = Math.max(Math.abs(err),maxerr);
      cnt++;
    }
    // extra check on Neumann boundary ?

    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln

  public static void main (String[] args)
  {
    new pde_sparse14();
  }
} // end class pde_sparse14