# include <math.h>
# include <stdlib.h>

# include "porous_medium_exact.h"

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

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

/******************************************************************************/
/*
  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:

    19 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:

    real x, t: the position and time.

  Output:

    real 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;
}
/******************************************************************************/

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 )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double porous_medium_residual ( double t, double x )

/******************************************************************************/
/*
  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;
}