# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <string.h>
# include <time.h>

int main ( );
double he_double_product_integral ( int i, int j );
double he_triple_product_integral ( int i, int j, int k );
double r8_factorial ( int n );
void timestamp ( );

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for PCE_BURGERS.

  Discussion:

    The time-dependent viscous Burgers equation to be solved is:

      du/dt = - d ( u*(1/2-u)) /dx + nu d2u/dx2

    with boundary conditions

      u(-3.0) = 0.0, u(+3.0) = 1.0.

    The viscosity nu is assumed to be an uncertain quantity with
    normal distribution of known mean and variance.

    A polynomial chaos expansion is to be used, with Hermite polynomial
    basis functions h(i,x), 0 <= i <= n.

    Because the first two Hermite polynomials are simply 1 and x, 
    we have that 

      nu = nu_mean * h(0,x) + nu_variance * h(1,x).

    We replace the time derivative by an explicit Euler approximation,
    so that the equation now describes the value of U(x,t+dt) in terms
    of known data at time t.

    Now assume that the solution U(x,t) can be approximated
    by the truncated expansion:

      U(x,t) = sum ( 0 <= i <= n ) c(i,t) * h(i,x)

    In the equation, we replace U by its expansion, and then multiply
    successively by each of the basis functions h(*,x) to get a set of
    n+1 equations that can be used to determine the values of c(i,t+dt).

    This process is repeated until the desired final time is reached.

    At any time, the coefficients c(0,t) contain information definining
    the expected value of u(x,t) at that time, while the higher order coefficients
    can be used to deterimine higher moments.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    15 September 2012

  Author:

    Original FORTRAN90 version by Gianluca Iaccarino.
    C version by John Burkardt.

  Local parameters:

    Local, double DT, the timestep.

    Local, double DX, the spacing between grid points.

    Local, int N, the number of intervals in the spatial domain.

    Local, double NUMEAN, the mean of viscosity.

    Local, double NUVARIANCE, the variance of viscosity.

    Local, int P, the order of the PC expansion.

    Local, double T, the current time.

    Local, double TF, the final integration time.

    Local, double U1[(N+1)*(P+1)], the PCE representation at the current time.

    Local, double U2[(N+1)*(P+1)], the PCE representation for the next time.

    Local, double X[N+1], the grid points.
*/
{
  double conv;
  double dp;
  double dt;
  double dx;
  int i;
  int it;
  int ix;
  int j;
  int k;
  int n;
  int nt;
  double numean;
  double nuvariance;
  FILE *output;
  char output_filename[81];
  int p;
  double t1;
  double t2;
  double term1;
  double term2;
  double tf;
  double ti;
  double tp;
  double *u1;
  double *u2;
  double *umean;
  double *uvariance;
  double visc[2];
  double *x;

  p = 5 ;
  n = 32;
  nt = 2000;
  ti = 0.0;
  tf = 2.0;
  dt = ( tf - ti ) / ( double ) ( nt );
  numean = 0.25;
  nuvariance = 0.08;

  timestamp ( );
  printf ( "\n" );
  printf ( "PCE_BURGERS:\n" );
  printf ( "  C version\n" );
  printf ( "\n" );
  printf ( "  Polynomial Chaos Solution\n" );
  printf ( "  1D Burgers equation\n" );
  printf ( "  Original version by Gianluca Iaccarino\n" );
  printf ( "\n" );
  printf ( "  PCE order       = %d\n", p );
  printf ( "  Number of cells = %d\n", n );
  printf ( "  Time step       = %g\n", dt );
  printf ( "  Initial time    = %g\n", ti );
  printf ( "  Final time      = %g\n", tf );
  printf ( "  Viscosity Mean  = %g\n", numean );
  printf ( "  Viscosity Var   = %g\n", nuvariance );
  printf ( "\n" );

  u1 = ( double * ) malloc ( (n+1)*(p+1) * sizeof ( double ) );
  u2 = ( double * ) malloc ( (n+1)*(p+1) * sizeof ( double ) );
  x = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
/*
  Define some numerical parameters
*/
  dx = 6.0 / ( double ) ( n );
  conv = dt / ( 2.0 * dx );
/*
  The expansion for viscosity stops at the linear term.
*/
  visc[0] = numean * dt / ( dx * dx );
  visc[1] = nuvariance * dt / ( dx * dx );
/*
  Define a uniform grid
*/
  for ( i = 0; i <= n; i++ )
  {
    x[i] = ( ( double ) ( n - i ) * ( -3.0 )   
           + ( double ) (     i ) * ( +3.0 ) ) 
           / ( double ) ( n     );
  }
/*
  Set the initial conditions
*/
  for ( j = 0; j <= p; j++ )
  {
    for ( i = 0; i <= n; i++ )
    {
      u1[i+j*(n+1)] = 0.0;
    }
  }

  for ( i = 0; i <= n; i++ )
  {
    u1[i+0*(n+1)] = 0.5 + x[i] / 6.0;
  }
/*
  Write the current solution.
*/
  strcpy ( output_filename, "burgers.history.txt" );

  output = fopen ( output_filename, "wt" );

  fprintf ( output, "----------\n" );
  fprintf ( output, "T = %g\n", ti );
  for ( i = 0; i <= n; i++ )
  {
    fprintf ( output, "  %14g", x[i] );
    for ( j = 0; j <= p; j++ )
    {
      fprintf ( output, "  %14g", u1[i+j*(n+1)] );
    }
    fprintf ( output, "\n" );
  }
/*
  Time integration
*/
  t1 = ti;

  for ( it = 1; it <= nt; it++ )
  {
    t2 = ( ( double ) ( nt - it ) * ti   
         + ( double ) (      it ) * tf ) 
         / ( double ) ( nt      );
/*
  Boundary conditions.
*/
    for ( j = 0; j <= p; j++ )
    {
      u2[0+j*(n+1)] = 0.0;
    }
    u2[n+0*(n+1)] = 1.0;
    for ( j = 1; j <= p; j++ )
    {
      u2[n+j*(n+1)] = 0.0;
    }

    for ( k = 0; k <= p; k++ )
    {
      dp = he_double_product_integral ( k, k );

      for ( ix = 1; ix < n; ix++ )
      {
/*
  Viscous term.
*/
        term1 = visc[0] * ( u1[ix+1+k*(n+1)] - 2.0 * u1[ix+k*(n+1)] + u1[ix-1+k*(n+1)] );
        i = 1;
        for ( j = 0; j <= p; j++ )
        {
          tp = he_triple_product_integral ( i, j, k );
          term1 = term1 + visc[i] 
            * ( u1[ix+1+j*(n+1)] - 2.0 * u1[ix+j*(n+1)] + u1[ix-1+j*(n+1)] ) * tp / dp;
        }
/*
  Convective term.
*/
        term2 = - conv * 0.5 * ( u1[ix+1+k*(n+1)] - u1[ix-1+k*(n+1)] );
        for ( j = 0; j <= p; j++ )
        {
          for ( i = 0; i <= p; i++ )
          {
            tp = he_triple_product_integral ( i, j, k );
            term2 = term2 + ( conv 
              * u1[ix+i*(n+1)] * ( u1[ix+1+j*(n+1)] - u1[ix-1+j*(n+1)] ) * tp ) / dp;
          }
        }
        u2[ix+k*(n+1)] = u1[ix+k*(n+1)] + term1 + term2;
      }
    }

    t1 = t2;
    for ( j = 0; j <= p; j++ )
    {
      for ( i = 0; i <= n; i++ )
      {
        u1[i+j*(n+1)] = u2[i+j*(n+1)];
      }
    }
/*
  Print solution every 100 time steps.
*/
    if ( ( it % 100 ) == 0 )
    {
      fprintf ( output, "----------\n" );
      fprintf ( output, "T = %g\n", t1 );
      for ( i = 0; i <= n; i++ )
      {
        fprintf ( output, "  %14g", x[i] );
        for ( j = 0; j <= p; j++ )
        {
          fprintf ( output, "  %14g", u1[i+j*(n+1)] );
        }
        fprintf ( output, "\n" );
      }
    }
  }
  fclose ( output );
  printf ( "  Time history in \"%s\".\n", output_filename );
/*
  Compute the mean and variance.
*/
  umean = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
  uvariance = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );

  for ( i = 0; i <= n; i++ )
  {
    umean[i] = u1[i+0*(n+1)];
  }

  for ( i = 0; i <= n; i++ )
  {
    uvariance[i] = 0.0;
    for ( j = 1; j <= p; j++ )
    {
      dp = he_double_product_integral ( j, j );
      uvariance[i] = uvariance[i] + pow ( u1[i+j*(n+1)], 2 ) * dp;
    }
  }
/*
  Write data about the solution at the final time.
*/
  strcpy ( output_filename, "burgers.moments.txt" );
  output = fopen ( output_filename, "wt" );
  fprintf ( output, "X E[U] Var[U]\n" );
  for ( i = 0; i <= n; i++ )
  {
    fprintf ( output, "  %18g  %18g  %18g\n", x[i], umean[i], uvariance[i] );
  }
  fclose ( output );
  printf ( "  Moments written to \"%s\".\n", output_filename );

  strcpy ( output_filename, "burgers.modes.txt" );
  output = fopen ( output_filename, "wt" );
  fprintf ( output, "X U_0 ... U_P\n" );
  for ( i = 0; i <= n; i++ )
  {
    fprintf ( output, "  %20g", x[i] );
    for ( j = 0; j <= p; j++ )
    {
      fprintf ( output, "  %20g", u1[i+j*(n+1)] );
    }
    fprintf ( output, "\n" );
  }
  fclose ( output );
  printf ( "  Final modes written to \"%s\".\n", output_filename );
/*
  Free memory.
*/
  free ( u1 );
  free ( u2 );
  free ( umean );
  free ( uvariance );
  free ( x );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "PCE_BURGERS:\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

double he_double_product_integral ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    HE_DOUBLE_PRODUCT_INTEGRAL: integral of He(i,x)*He(j,x)*e^(-x^2/2).

  Discussion:

    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x) exp(-x^2/2) dx

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    16 March 2012

  Author:

    John Burkardt

  Reference:

    Dongbin Xiu,
    Numerical Methods for Stochastic Computations: A Spectral Method Approach,
    Princeton, 2010,
    ISBN13: 978-0-691-14212-8,
    LC: QA274.23.X58.

  Parameters:

    Input, int I, J, the polynomial indices.

    Output, double HE_DOUBLE_PRODUCT_INTEGRAL, the value of the integral.
*/
{
  double value;

  if ( i == j )
  {
    value = r8_factorial ( i );
  }
  else
  {
    value = 0.0;
  }
  return value;
}
/******************************************************************************/

double he_triple_product_integral ( int i, int j, int k )

/******************************************************************************/
/*
  Purpose:

    HE_TRIPLE_PRODUCT_INTEGRAL: integral of He(i,x)*He(j,x)*He(k,x)*e^(-x^2/2).

  Discussion:

    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x)*He(k,x) exp(-x^2/2) dx

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 March 2012

  Author:

    John Burkardt

  Reference:

    Dongbin Xiu,
    Numerical Methods for Stochastic Computations: A Spectral Method Approach,
    Princeton, 2010,
    ISBN13: 978-0-691-14212-8,
    LC: QA274.23.X58.

  Parameters:

    Input, int I, J, K, the polynomial indices.

    Output, double HE_TRIPLE_PRODUCT_INTEGRAL, the value of the integral.
*/
{
  int s;
  double value;

  s = ( i + j + k ) / 2;

  if ( s < i || s < j || s < k )
  {
    value = 0.0;
  }
  else if ( ( ( i + j + k ) % 2 ) != 0 )
  {
    value = 0.0;
  }
  else
  {
    value = r8_factorial ( i ) / r8_factorial ( s - i ) 
          * r8_factorial ( j ) / r8_factorial ( s - j ) 
          * r8_factorial ( k ) / r8_factorial ( s - k );
  }

  return value;
}
/******************************************************************************/

double r8_factorial ( int n )

/******************************************************************************/
/*
  Purpose:

    R8_FACTORIAL computes the factorial of N.

  Discussion:

    factorial ( N ) = product ( 1 <= I <= N ) I

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    16 January 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the argument of the factorial function.
    If N is less than 1, the function value is returned as 1.

    Output, double R8_FACTORIAL, the factorial of N.
*/
{
  int i;
  double value;

  value = 1.0;

  for ( i = 1; i <= n; i++ )
  {
    value = value * ( double ) ( i );
  }

  return value;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}