# include # include # include # include "ppc_array.h" # include "ppc_plot3d.h" # include "ppc_porous_medium.h" # include "ppc_xmalloc.h" double croot ( double k ); double error_versus_exact ( struct pme_spec *prob ); double phi ( double x ); void plot_with_geomview ( struct plot3d *data ); void plot_with_maple ( struct plot3d *data ); void plot_with_matlab ( struct plot3d *data ); void plot3d ( struct plot3d *data ); void pme_solve ( struct pme_spec *prob ); void write_plotting_script ( struct pme_spec *prob ); /******************************************************************************/ double croot ( double k ) /******************************************************************************/ /* Purpose: croot() solves the cubic equation eta + eta^3 = k 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 k: the right hand side of the equation. Output: double value: the solution to the cubic equation. */ { double s = sqrt ( 12.0 + 81.0 * k * k ); double gamma = pow ( 108.0 * k + 12.0 * s, 1.0 / 3.0 ); double value = gamma / 6.0 - 2.0 / gamma; return value; } /******************************************************************************/ double error_vs_exact ( struct pme_spec *prob ) /******************************************************************************/ /* Purpose: error_vs_exact() compares computed solution to formula for exact solution. Licensing: This code is distributed under the MIT license. Modified: 19 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: struct pme_spec *prob: data defining the problem. Output: double diff_max: the maximum value of the absolute value of the difference between the computed and exact solutions at grid points. */ { double diff_max; int i; int j; double t; double u_comp; double u_exact; double x; diff_max = 0.0; /* For each time value t: */ for ( i = 0; i <= prob->m; i++ ) { t = i * prob->T / prob->m; /* For each space value x: */ for ( j = 0; j <= prob->n; j++ ) { x = ( ( prob->n - j ) * prob->a + j * prob->b ) / prob->n; u_exact = prob->exact_sol ( x, t ); u_comp = prob->u[i][j]; diff_max = fmax ( diff_max, fabs ( u_exact - u_comp ) ); } } return diff_max; } /******************************************************************************/ double phi ( double u ) /******************************************************************************/ /* Purpose: phi() evaluates phi(u) in the porous medium equation du/dt = d2 phi(u)/dx^2. Licensing: This code is distributed under the MIT license. Modified: 29 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 u: the argument. Output: double value: the function value. */ { double value; value = pow ( u, 3 ); return value; } /******************************************************************************/ void plot_with_geomview ( struct plot3d *data ) /******************************************************************************/ /* Purpose: plot_with_geomview() plots porous medium solution data using geomview(). Licensing: This code is distributed under the MIT license. Modified: 19 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 Input: struct plot3d *data: the plot data. */ { double dt; double dx; FILE *fp; int m; int n; m = data->m; n = data->n; dx = ( data->b - data->a ) / n; dt = data->T / m; fp = fopen ( data->geomview_out, "w" ); if ( fp == NULL ) { fprintf ( stderr, "*** Cannot open script file %s for writing\n", data->geomview_out ); return; } fprintf ( fp, "# Geomview script for porous medium solver\n" ); fprintf ( fp, "{ appearance { +edge }\n" ); fprintf ( fp, "MESH %d %d\n", n+1, m+1 ); for ( int i = 0; i <= m; i++ ) { double t = i * dt; for ( int j = 0; j <= n; j++ ) { double x = data->a + j * dx; fprintf ( fp, "%g %g %g\n", x, t, data->u[i][j] ); } } fprintf ( fp, "}\n" ); fclose ( fp ); printf ( " Geomview script written to %s\n", data->geomview_out ); return; } /******************************************************************************/ void plot_with_maple ( struct plot3d *data ) /******************************************************************************/ /* Purpose: plot_with_maple() plots porous medium solution data using maple(). Licensing: This code is distributed under the MIT license. Modified: 19 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 Input: struct plot3d *data: the plot data. */ { double dt; double dx; FILE *fp; int m; int n; m = data->m; n = data->n; dx = ( data->b - data->a ) / n; dt = data->T / m; fp = fopen ( data->maple_out, "w" ); if ( fp == NULL ) { fprintf ( stderr, "*** Cannot open file %s for writing\n", data->maple_out ); return; } fprintf ( fp, "# Maple script for porous medium solver\n"); fprintf ( fp, "plots:-display(MESH([\n" ); for ( int i = 0; i <= m; i++ ) { double t = i * dt; fprintf ( fp, "\t[" ); for ( int j = 0; j <= n; j++ ) { double x = data->a + j * dx; fprintf ( fp, "[%g,%g,%g],", x, t, data->u[i][j] ); } fprintf ( fp, "NULL],\n" ); } fprintf ( fp, "NULL]), labels=[x,t,u]);" ); fclose ( fp ); printf ( " Maple script written to %s\n", data->maple_out ); return; } /******************************************************************************/ void plot_with_matlab ( struct plot3d *data ) /******************************************************************************/ /* Purpose: plot_with_matlab() plots porous medium solution data using matlab(). Licensing: This code is distributed under the MIT license. Modified: 19 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 Input: struct plot3d *data: the plot data. */ { // names of matlab's plotting variables # define XVEC "PLOT_WITH_MATLAB_FD_XVEC" # define TVEC "PLOT_WITH_MATLAB_FD_TVEC" # define ZMAT "PLOT_WITH_MATLAB_FD_ZMAT" double dt; double dx; FILE *fp; int m; int n; m = data->m; n = data->n; dx = ( data->b - data->a ) / n; dt = data->T / m; fp = fopen ( data->matlab_out, "w" ); if ( fp == NULL ) { fprintf ( stderr, "*** Cannot open file %s for writing\n", data->matlab_out ); return; } fprintf ( fp, "%% Matlab script for porous medium solver\n" ); fprintf ( fp, "%s = [\n", XVEC ); for ( int j = 0; j <= n; j++ ) { double x = data->a + j * dx; fprintf ( fp, "%g ", x ); } fprintf ( fp, "\n];\n" ); fprintf ( fp, "%s = [\n", TVEC ); for ( int i = 0; i <= m; i++ ) { double t = i * dt; fprintf ( fp, "%g ", t ); } fprintf ( fp, "\n];\n" ); fprintf ( fp, "%s = [\n", ZMAT ); for ( int i = 0; i <= m; i++ ) { for ( int j = 0; j <= n; j++ ) { fprintf ( fp, "%g ", data->u[i][j] ); } fprintf ( fp, "\n" ); } fprintf ( fp, "];\n" ); fprintf ( fp, "surf(%s, %s, %s)\n", XVEC, TVEC, ZMAT ); fprintf ( fp, "xlabel('x')\n" ); fprintf ( fp, "ylabel('t')\n" ); fprintf ( fp, "zlabel('u')\n" ); fprintf ( fp, "colormap jet\n" ); fclose ( fp ); printf ( " Matlab script written to %s\n", data->matlab_out ); return; # undef XVEC # undef TVEC # undef ZMAT } /******************************************************************************/ void plot3d ( struct plot3d *data ) /******************************************************************************/ /* Purpose: plot3d() plots ppc_fd1() data using geomview(), maple(), or matlab(). Licensing: This code is distributed under the MIT license. Modified: 18 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 Input: struct plot3d *data: the plot data. */ { if ( data->geomview_out != NULL ) { plot_with_geomview ( data ); } if ( data->maple_out != NULL ) { plot_with_maple ( data ); } if ( data->matlab_out != NULL ) { plot_with_matlab ( data ); } return; } /******************************************************************************/ void pme_solve ( struct pme_spec *prob ) /******************************************************************************/ /* Purpose: pme_solve() solves the porous medium equation. Licensing: This code is distributed under the MIT license. Modified: 29 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: struct pme_spec *prob: a problem specification. Output: struct pme_spec *prob: the problem specification has been updated to include the solution. */ { double a; double b; double c; double dt; double dx; int i; int j; int m; int n; double r; double s; double t; double T; double **u; double *v; double x; m = prob->m; n = prob->n; a = prob->a; b = prob->b; dx = ( b - a ) / n; T = prob->T; dt = T / m; u = prob->u; r = dt / ( 2.0 * dx * dx ); s = sqrt ( r ); make_vector ( v, n + 1 ); for ( j = 0; j <= n; j++ ) { x = a + j * dx; u[0][j] = prob->ic ( x ); } for ( i = 1; i <= m; i++ ) { t = i * dt; /* The left half timestep. */ v[0] = prob->bcL ( t - dt / 2 ); for ( j = 1; j <= n - 1; j++ ) { c = r * phi ( v[j-1] ) + u[i-1][j] - r * phi ( u[i-1][j] ) + r * phi ( u[i-1][j+1] ); v[j] = croot ( s * c ) / s; } v[n] = prob->bcR ( t - dt / 2 ); /* The right half timestep. */ u[i][n] = prob->bcR ( t ); for ( j = n - 1; 1 <= j; j-- ) { c = r * phi ( v[j-1] ) + v[j] - r * phi ( v[j] ) + r * phi ( u[i][j+1] ); u[i][j] = croot ( s * c ) / s; } u[i][0] = prob->bcL ( t ); } free_vector ( v ); write_plotting_script ( prob ); if ( prob->exact_sol != NULL ) { prob->error = error_vs_exact ( prob ); } return; } /******************************************************************************/ void write_plotting_script ( struct pme_spec *prob ) /******************************************************************************/ /* Purpose: write_plotting_script() prepares solution data to be plotted. 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: struct pme_spec *prob: the porous medium problem structure. */ { struct plot3d data; data.a = prob->a; data.b = prob->b; data.T = prob->T; data.n = prob->n; data.m = prob->m; data.u = prob->u; data.maple_out = prob->maple_out; data.matlab_out = prob->matlab_out; data.geomview_out = prob->geomview_out; plot3d ( &data ); }