# include <cmath>
# include <cstdlib>

using namespace std;

# include "porous_medium_exact.hpp"

//****************************************************************************80

void porous_medium_exact ( double x, double t, double &u, double &ut, 
  double &ux, double &uxx )

//****************************************************************************80
//
//  Purpose:
//
//    porous_medium_exact() evaluates an exact solution of the porous medium equation.
//
//  Discussion:
//
//    du/dt = Del^2 ( u^m )
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    20 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Grigory Barenblatt,
//    On some unsteady fluid and gas motions in a porous medium,
//    Prikladnaya Matematika i Mekhanika (Applied Mathematics and Mechanics, 
//    Volume 16, Number 1, pages 67-78, 1952.
//
//    Rouben Rostamian,
//    Programming Projects in C 
//    for Students of Engineering, Science, and Mathematics,
//    SIAM, 2014,
//    ISBN: 978-1-611973-49-5
//
//  Input:
//
//    double x, t: the position and time.
//
//  Output:
//
//    double &u, &ut, &ux, &uxx: the values of the exact solution, its time
//    derivative, and its first and second spatial derivatives at (x,t).
//
{
  double alpha;
  double beta;
  double bot;
  double c;
  double delta;
  double factor;
  double gamma;
  double m;

  porous_medium_parameters ( 
    NULL, NULL,   NULL, NULL, NULL,
    &c,   &delta, &m,   NULL, NULL );

  alpha =         1.0 / ( m - 1.0 );
  beta =          1.0 / ( m + 1.0 );
  gamma = ( m - 1.0 ) / ( 2.0 * m * ( m + 1.0 ) );

  bot = pow ( t + delta, beta );
  factor = c - gamma * pow ( x / bot, 2 );

  if ( 0.0 < factor )
  {
    u = 1.0 / pow ( t + delta, beta ) * pow ( factor, alpha );

    ut = 2.0 * alpha * beta * gamma * pow ( t + delta, -1.0-3.0*beta ) 
      * x * x * pow ( factor, alpha - 1.0 ) 
      - beta * pow ( t + delta, -1.0-beta) * pow ( factor, alpha );

    ux = - 2.0 * alpha * gamma 
      * pow ( t + delta, -3.0 * beta ) * x * pow ( factor, alpha-1.0 );

    uxx = 4.0 * ( alpha - 1.0 ) * alpha * gamma * gamma
      * pow ( t + delta, -5.0*beta ) * x * x * pow ( factor, alpha-2.0 ) 
      - 2.0 * alpha * gamma * pow ( t + delta, -3.0 * beta ) 
      * pow ( factor, alpha-1.0 );
  }
  else
  {
    u = 0.0;
    ut = 0.0;
    ux = 0.0;
    uxx = 0.0;
  }
  return;
}
//****************************************************************************80

void porous_medium_parameters ( double *c_in, double *delta_in, double *m_in, 
  double *t0_in, double *tstop_in, double *c_out, double *delta_out, 
  double *m_out, double *t0_out, double *tstop_out )

//****************************************************************************80
//
//  Purpose:
//
//    porous_medium_parameters() returns parameters for the porous medium equation.
//
//  Discussion:
//
//    du/dt = Del^2 ( u^m )
//
//    If input values are specified, this resets the default parameters.
//    Otherwise, the output will be the current defaults.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double *c_in: a linear factor for the volume of the solution at any time. 
//
//    double *delta_in: an offset to the time t.
//
//    double *m_in: the power of u in the equation.  m must be greater than 1.
//
//    double *t0_in: the initial time.
//
//    double *tstop_in: the final time.
//
//  Output:
//
//    double *c_out: a linear factor for the volume of the solution at any time. 
//
//    double *delta_out: an offset to the time t.
//
//    double *m_out: the power of u in the equation.  m must be greater than 1.
//
//    double *t0_out: the initial time.
//
//    double *tstop_out: the final time.
//
{
  static double c_default = sqrt ( 3.0 ) / 15.0;
  static double delta_default = 1.0 / 75.0;
  static double m_default = 3.0;
  static double t0_default = 0.0;
  static double tstop_default = 4.0;
//
//  New values, if supplied on input, overwrite the current values.
//
  if ( c_in )
  {
    c_default = *c_in;
  }

  if ( delta_in )
  {
    delta_default = *delta_in;
  }

  if ( m_in )
  {
    m_default = *m_in;
  }

  if ( t0_in )
  {
    t0_default = *t0_in;
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
//
//  The current values are copied to the output.
//
  if ( c_out )
  {
    *c_out = c_default;
  }

  if ( delta_out )
  {
    *delta_out = delta_default;
  }

  if ( m_out )
  {
    *m_out = m_default;
  }

  if ( t0_out )
  {
    *t0_out = t0_default;
  }

  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }

  return;
}
//****************************************************************************80

double porous_medium_residual ( double t, double x )

//****************************************************************************80
//
//  Purpose:
//
//    porous_medium_residual() computes the residual of the porous medium equation.
//
//  Discussion:
//
//    ut = Del^2 ( u^m )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double t, x: the time and position where the solution is evaluated.
//
//  Output:
//
//    double r: the residual at that time and position.
//
{
  double m;
  double r;
  double u;
  double ut;
  double ux;
  double uxx;

  porous_medium_parameters ( 
    NULL, NULL, NULL, NULL, NULL, 
    NULL, NULL, &m,   NULL, NULL );

  porous_medium_exact ( t, x, u, ut, ux, uxx );

  r = ut - m * ( m - 1 ) * pow ( u, m-2 ) * ux * ux - m * pow ( u, m-1 ) * uxx;

  return r;
}