fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <string.h>
  5.  
  6. #define NX 100                  // Number of spatial points
  7. #define NT_TOTAL 500            // Total time steps for T = 0.5
  8. #define X_MAX 1.0               // Maximum x value
  9. #define PI 3.141592653589793
  10.  
  11. // Function prototypes
  12. double erf_sol(double x, double t, double Re);
  13. void explicit_scheme(double *u, double Re, int Nt, double dt, double dx);
  14. void crank_nicolson_scheme(double *u, double Re, int Nt, double dt, double dx);
  15.  
  16. // Main function
  17. int main() {
  18.     double x[NX], u_initial[NX];
  19.     double T_values[] = {0.36, 0.47};      // Times to evaluate the solution
  20.     double Re_values[] = {10, 50};          // Reynolds numbers
  21.     double u_exact[NX], u_explicit[NX], u_cn[NX];
  22.     double dt = 0.5 / NT_TOTAL;             // Time step size
  23.     double dx = 2 * X_MAX / (NX - 1);      // Spatial step size
  24.  
  25.     // Generate grid and initialize conditions
  26.     for (int i = 0; i < NX; ++i) {
  27.         x[i] = -X_MAX + i * dx;
  28.         u_initial[i] = (x[i] < 0) ? 1.0 : 0.0; // Initial condition
  29.     }
  30.  
  31.     // Open gnuplot
  32.     FILE *gp = popen("gnuplot -persistent", "w");
  33.     fprintf(gp, "set multiplot layout 2,2\n");
  34.  
  35.     for (int r = 0; r < 2; ++r) {
  36.         double Re = Re_values[r];
  37.  
  38.         for (int t_idx = 0; t_idx < 2; ++t_idx) {
  39.             double T = T_values[t_idx];
  40.             int Nt = (int)(T / dt);  // Number of time steps
  41.  
  42.             // Exact solution
  43.             for (int i = 0; i < NX; ++i) {
  44.                 u_exact[i] = erf_sol(x[i], T, Re);
  45.             }
  46.  
  47.             // Explicit solution
  48.             memcpy(u_explicit, u_initial, NX * sizeof(double));
  49.             explicit_scheme(u_explicit, Re, Nt, dt, dx);
  50.  
  51.             // Crank-Nicolson solution
  52.             memcpy(u_cn, u_initial, NX * sizeof(double));
  53.             crank_nicolson_scheme(u_cn, Re, Nt, dt, dx);
  54.  
  55.             // Plot results for explicit scheme vs exact solution
  56.             fprintf(gp, "plot '-' title 'Explicit Re=%.1f, T=%.2f' with lines, '-' title 'Exact' with lines\n", Re, T);
  57.             for (int i = 0; i < NX; ++i) {
  58.                 fprintf(gp, "%f %f\n", x[i], u_explicit[i]);
  59.             }
  60.             fprintf(gp, "e\n");
  61.             for (int i = 0; i < NX; ++i) {
  62.                 fprintf(gp, "%f %f\n", x[i], u_exact[i]);
  63.             }
  64.             fprintf(gp, "e\n");
  65.  
  66.             // Plot error for explicit scheme
  67.             fprintf(gp, "plot '-' title 'Error Explicit' with lines\n");
  68.             for (int i = 0; i < NX; ++i) {
  69.                 fprintf(gp, "%f %f\n", x[i], fabs(u_explicit[i] - u_exact[i]));
  70.             }
  71.             fprintf(gp, "e\n");
  72.  
  73.             // Plot results for Crank-Nicolson scheme vs exact solution
  74.             fprintf(gp, "plot '-' title 'Crank-Nicolson Re=%.1f, T=%.2f' with lines, '-' title 'Exact' with lines\n", Re, T);
  75.             for (int i = 0; i < NX; ++i) {
  76.                 fprintf(gp, "%f %f\n", x[i], u_cn[i]);
  77.             }
  78.             fprintf(gp, "e\n");
  79.             for (int i = 0; i < NX; ++i) {
  80.                 fprintf(gp, "%f %f\n", x[i], u_exact[i]);
  81.             }
  82.             fprintf(gp, "e\n");
  83.  
  84.             // Plot error for Crank-Nicolson scheme
  85.             fprintf(gp, "plot '-' title 'Error Crank-Nicolson' with lines\n");
  86.             for (int i = 0; i < NX; ++i) {
  87.                 fprintf(gp, "%f %f\n", x[i], fabs(u_cn[i] - u_exact[i]));
  88.             }
  89.             fprintf(gp, "e\n");
  90.         }
  91.     }
  92.  
  93.     fprintf(gp, "unset multiplot\n");
  94.     pclose(gp);
  95.  
  96.     return 0;
  97. }
  98.  
  99. // Error function approximation
  100. double erf_sol(double x, double t, double Re) {
  101.     return 0.5 * (1 - erf(x / (2 * sqrt((1 / Re) * t))));
  102. }
  103.  
  104. // Explicit scheme function
  105. void explicit_scheme(double *u, double Re, int Nt, double dt, double dx) {
  106.     double u_new[NX];
  107.  
  108.     for (int n = 0; n < Nt; ++n) {
  109.         memcpy(u_new, u, NX * sizeof(double));
  110.         for (int i = 1; i < NX - 1; ++i) {
  111.             u_new[i] = u[i] - dt * u[i] * (u[i + 1] - u[i - 1]) / (2 * dx)
  112.                         + (dt / Re) * (u[i + 1] - 2 * u[i] + u[i - 1]) / (dx * dx);
  113.         }
  114.         memcpy(u, u_new, NX * sizeof(double));
  115.     }
  116. }
  117.  
  118. // Crank-Nicolson scheme function
  119. void crank_nicolson_scheme(double *u, double Re, int Nt, double dt, double dx) {
  120.     double alpha = dt / (2 * dx * dx * Re);
  121.     double B[NX];
  122.     double main_diag[NX], upper_diag[NX - 1], lower_diag[NX - 1];
  123.  
  124.     // Initialize the diagonals
  125.     for (int i = 0; i < NX; ++i) {
  126.         main_diag[i] = 1 + 2 * alpha;
  127.         if (i < NX - 1) {
  128.             upper_diag[i] = -alpha;
  129.             lower_diag[i] = -alpha;
  130.         }
  131.     }
  132.  
  133.     for (int n = 0; n < Nt; ++n) {
  134.         for (int i = 1; i < NX - 1; ++i) {
  135.             B[i] = u[i] - dt * u[i] * (u[i + 1] - u[i - 1]) / (4 * dx)
  136.                    + alpha * (u[i + 1] - 2 * u[i] + u[i - 1]);
  137.         }
  138.         B[0] = 1.0;  // Left boundary condition
  139.         B[NX - 1] = 0.0;  // Right boundary condition
  140.  
  141.         // Thomas Algorithm (tridiagonal matrix solver)
  142.         // Forward elimination
  143.         for (int i = 1; i < NX; ++i) {
  144.             double m = upper_diag[i - 1] / main_diag[i - 1];
  145.             main_diag[i] -= m * lower_diag[i - 1];
  146.             B[i] -= m * B[i - 1];
  147.         }
  148.  
  149.         // Back substitution
  150.         u[NX - 1] = B[NX - 1] / main_diag[NX - 1];
  151.         for (int i = NX - 2; i >= 0; --i) {
  152.             u[i] = (B[i] - lower_diag[i] * u[i + 1]) / main_diag[i];
  153.         }
  154.     }
  155. }
  156.  
Success #stdin #stdout #stderr 0.01s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
sh: 1: gnuplot: not found
https://ideone.com/06lsuj
language:
C (gcc 8.3)
created:
3 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!