# include # include # include # include # include "ppc_array.h" # include "ppc_porous_medium.h" # include "ppc_xmalloc.h" int main ( int argc, char **argv ); double barenblatt ( double x, double t, double m, double c, double delta ); void pme1_solve ( int m, int n, double T ); double pme1_bcL ( double t ); double pme1_bcR ( double t ); double pme1_exact ( double x, double t ); double pme1_ic ( double x ); void pme2_solve ( int m, int n, double T ); double pme2_bcL ( double t ); double pme2_bcR ( double t ); double pme2_ic ( double x ); void show_usage ( char *progname ); void timestamp ( ); /******************************************************************************/ int main ( int argc, char **argv ) /******************************************************************************/ /* Purpose: ppc_porous_medium_test() tests ppc_porous_medium(). Licensing: This code is distributed under the MIT license. Modified: 20 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 */ { char *endptr; int m; int n; double T; timestamp ( ); printf ( "\n" ); printf ( "ppc_porous_medium_test ( ):\n" ); printf ( " C version.\n" ); printf ( " Test ppc_mesh().\n" ); if ( argc != 4 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } T = strtod ( argv[1], &endptr ); if ( *endptr != '\0' || T <= 0.0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } n = strtol ( argv[2], &endptr, 10 ); if ( *endptr != '\0' || n < 0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } m = strtol ( argv[3], &endptr, 10 ); if ( *endptr != '\0' || m < 0 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } /* Do PME1 */ pme1_solve ( m, n, T ); /* Do PME2 */ pme2_solve ( m, n, T ); /* Terminate. */ printf ( "\n" ); printf ( "ppc_porous_medium_test ( ):\n" ); printf ( " Normal end of execution.\n" ); timestamp ( ); return EXIT_SUCCESS; } /******************************************************************************/ double barenblatt ( double x, double t, double m, double c, double delta ) /******************************************************************************/ /* Purpose: barenblatt() evaluates Barenblatt's porous medium exact solution. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: 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. double m, c, delta: parameters. Output: double value: the solution value at (x,t). */ { double alpha; double beta; double bot; double factor; double gamma; double value; 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 ) { value = 1.0 / pow ( t + delta, beta ) * pow ( factor, alpha ); } else { value = 0.0; } return value; } /******************************************************************************/ void pme1_solve ( int m, int n, double T ) /******************************************************************************/ /* Purpose: pme1_solve() solves porous medium problem #1. Licensing: This code is distributed under the MIT license. Modified: 21 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 */ { double a; double b; double diff_max; double dt; double dx; double r; double s; printf ( "\n" ); printf ( "pme1_solve():\n" ); printf ( " Solve porous medium problem #1:\n" ); /* Set the problem specification. */ struct pme_spec pme1 = { .a = -1.0, .b = +1.0, .T = T, .n = n, .m = m, .ic = pme1_ic, .bcL = pme1_bcL, .bcR = pme1_bcR, .exact_sol = pme1_exact, .maple_out = "pme1.mpl", .matlab_out = "pme1.m", .geomview_out = "pme1.gv", }; /* Report the parameters. */ a = pme1.a; b = pme1.b; dt = T / m; dx = ( b - a ) / n; r = dt / 2.0 / dx / dx; s = sqrt ( r ); printf ( "\n" ); printf ( " Integrate over time t for 0 <= t <= T = %g\n", T ); printf ( " Take m = %d equal steps in time\n", m ); printf ( " dt = %g\n", dt ); printf ( " a = %g <= x <= b = %g\n", a, b ); printf ( " Use n = %d subintervals in space \n", n ); printf ( " dx = %g\n", dx ); printf ( " r = dt / 2 / dx^2 = %g\n", r ); printf ( " s = sqrt(r) = %g\n", s ); /* Prepare storage for solution U. */ make_matrix ( pme1.u, m + 1, n + 1 ); /* Solve the system. */ pme_solve ( &pme1 ); /* Report maximum difference between computed and exact solutions. */ if ( pme1.exact_sol ) { diff_max = error_vs_exact ( &pme1 ); printf ( "\n" ); printf ( " Maximum difference between computed and exact: %g\n", diff_max ); } /* Free the memory. */ free_matrix ( pme1.u ); return; } /******************************************************************************/ double pme1_bcL ( double t ) /******************************************************************************/ /* Purpose: pme1_bcL(): left boundary condition for porous medium problem #1. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double t: the time. Output: double pme1_bcL: the value of solution #1 at (-1,t). */ { double c; double delta; double m; double xL; double value; xL = -1.0; m = 3.0; c = sqrt ( 3.0 ) / 15.0; delta = 1.0 / 75.0; value = barenblatt ( xL, t, m, c, delta ); return value; } /******************************************************************************/ double pme1_bcR ( double t ) /******************************************************************************/ /* Purpose: pme1_bcR(): right boundary condition for porous medium problem #1. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double t: the time. Output: double pme1_bcR: the value of solution #1 at (+1,t). */ { double c; double delta; double m; double xR; double value; xR = +1.0; m = 3.0; c = sqrt ( 3.0 ) / 15.0; delta = 1.0 / 75.0; value = barenblatt ( xR, t, m, c, delta ); return value; } /******************************************************************************/ double pme1_exact ( double x, double t ) /******************************************************************************/ /* Purpose: pme1_exact() returns exact solution #1 of porous medium problem #1. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: 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 pme1_exact: the value of the exact solution #1 at (x,t). */ { double c; double delta; double m; double value; m = 3.0; c = sqrt ( 3.0 ) / 15.0; delta = 1.0 / 75.0; value = barenblatt ( x, t, m, c, delta ); return value; } /******************************************************************************/ double pme1_ic ( double x ) /******************************************************************************/ /* Purpose: pme1_ic(): initial condition for porous medium problem #1. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double x: the position. Output: double pme1_ic: the value of solution #1 at (x,0). */ { double c; double delta; double m; double t; double value; t = 0.0; m = 3.0; c = sqrt ( 3.0 ) / 15.0; delta = 1.0 / 75.0; value = barenblatt ( x, t, m, c, delta ); return value; } /******************************************************************************/ void pme2_solve ( int m, int n, double T ) /******************************************************************************/ /* Purpose: pme2_solve() solves porous medium problem #2. Licensing: This code is distributed under the MIT license. Modified: 20 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 */ { double a; double b; double diff_max; double dt; double dx; double r; double s; printf ( "\n" ); printf ( "pme2_solve():\n" ); printf ( " Solve porous medium problem #2:\n" ); /* Set the problem specification. */ struct pme_spec pme2 = { .a = -1.0, .b = +1.0, .T = T, .n = n, .m = m, .ic = pme2_ic, .bcL = pme2_bcL, .bcR = pme2_bcR, .exact_sol = NULL, .maple_out = "pme2.mpl", .matlab_out = "pme2.m", .geomview_out = "pme2.gv", }; /* Report the parameters. */ a = pme2.a; b = pme2.b; dt = T / m; dx = ( b - a ) / n; r = dt / 2.0 / dx / dx; s = sqrt ( r ); printf ( "\n" ); printf ( " Integrate over time t for 0 <= t <= T = %g\n", T ); printf ( " Take m = %d equal steps in time\n", m ); printf ( " dt = %g\n", dt ); printf ( " a = %g <= x <= b = %g\n", a, b ); printf ( " Use n = %d subintervals in space \n", n ); printf ( " dx = %g\n", dx ); printf ( " r = dt / 2 / dx^2 = %g\n", r ); printf ( " s = sqrt(r) = %g\n", s ); /* Prepare storage for solution U. */ make_matrix ( pme2.u, m + 1, n + 1 ); /* Solve the system. */ pme_solve ( &pme2 ); /* Report maximum difference between computed and exact solutions. (Don't, because we don't have the exact solution...) */ if ( pme2.exact_sol ) { diff_max = error_vs_exact ( &pme2 ); printf ( "\n" ); printf ( " Maximum difference between computed and exact: %g\n", diff_max ); } /* Free the memory. */ free_matrix ( pme2.u ); return; } /******************************************************************************/ double pme2_bcL ( double t ) /******************************************************************************/ /* Purpose: pme2_bcL(): left boundary condition for porous medium problem #2. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double t: the time. Output: double pme1_bcL: the value of the solution at (-1,t). */ { double value; value = 0.0; return value; } /******************************************************************************/ double pme2_bcR ( double t ) /******************************************************************************/ /* Purpose: pme2_bcR(): right boundary condition for porous medium problem #2.. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double t: the time. Output: double pme2_bcR: the value of solution #2 at (+1,t). */ { double value; value = 0.0; return value; } /******************************************************************************/ double pme2_ic ( double x ) /******************************************************************************/ /* Purpose: pme2_ic(): initial condition for porous medium problem #2. Licensing: This code is distributed under the MIT license. Modified: 28 May 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double x: the position. Output: double pme2_ic: the value of solution #2 at (x,0). */ { double value; if ( -0.5 <= x && x < 0.0 ) { value = -0.5; } else if ( 0.0 <= x && x <= 0.5 ) { value = +0.5; } else { value = 0.0; } return value; } /******************************************************************************/ void show_usage ( char *progname ) /******************************************************************************/ /* Purpose: show_usage_and_exit() prints usage instructions. Licensing: This code is distributed under the MIT license. Modified: 18 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: char *progname: the name of the program being executed. */ { printf ( "\n" ); printf ( "show_usage():\n" ); printf ( " %s T n m\n", progname ); printf ( " T = solve equation over 0 <= t <= T.\n" ); printf ( " n = number of x subintervals in [a,b].\n" ); printf ( " m = number of t subintervals in [0,T].\n" ); return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 Author: John Burkardt */ { # 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 ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }