# include "poisson.h" //****************************************************************************80 int main(int argc, char* argv[]) //****************************************************************************80 // // Purpose: // // MAIN is the main program for POISSON_SERIAL. // // Discussion: // // POISSON_SERIAL is a program for solving the Poisson problem. // // This program runs serially. Its output is used as a benchmark for // comparison with similar programs run in a parallel environment. // // The Poisson equation // // - DEL^2 U(x,y) = F(x,y) // // is solved on the unit square [0,1] x [0,1] using a grid of NX by // NX evenly spaced points. The first and last points in each direction // are boundary points. // // The boundary conditions and F are set so that the exact solution is // // U(x,y) = sin ( pi * x * y ) // // so that // // - DEL^2 U(x,y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y ) // // The Jacobi iteration is repeatedly applied until convergence is detected. // // For convenience in writing the discretized equations, we assume that NX = NY. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 26 October 2011 // // Author: // // John Burkardt // { bool converged; double diff; double dx; double dy; double error; double** f; int i; int it; int it_max = 100000; int j; int nx = 256; int ny = 128; double tolerance = 0.000001; double** u; double u_norm; double** udiff; double** uexact; double** unew; double unew_norm; double x; double y; //Check args if (argc < 2) { printf("USAGE:\tpoisson 0 (tiny)\n\tpoisson 1 (small)\n\tpoisson 2 (large)\n"); return 0; } // Read Input int type = atoi(argv[1]); //tiny if (type == 0) { nx = 16; ny = 16; it_max = 1000; } //small if (type == 1) { nx = 64; ny = 64; it_max = 100000; } //large if (type == 2) { nx = 256; ny = 256; it_max = 100000; } //Alloc memory f = new double* [nx]; u = new double* [nx]; udiff = new double* [nx]; uexact = new double* [nx]; unew = new double* [nx]; for (i = 0; i < nx; i++) { f[i] = new double[ny]; u[i] = new double[ny]; udiff[i] = new double[ny]; uexact[i] = new double[ny]; unew[i] = new double[ny]; } dx = 1.0 / (double)(nx - 1); dy = 1.0 / (double)(ny - 1); // // Print a message. // cout << "\n"; cout << "POISSON_SERIAL:\n"; cout << " C++ version\n"; cout << " A program for solving the Poisson equation.\n"; cout << "\n"; cout << " -DEL^2 U = F(X,Y)\n"; cout << "\n"; cout << " on the rectangle 0 <= X <= 1, 0 <= Y <= 1.\n"; cout << "\n"; cout << " F(X,Y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )\n"; cout << "\n"; cout << " The number of interior X grid points is " << nx << "\n"; cout << " The number of interior Y grid points is " << ny << "\n"; cout << " The X grid spacing is " << dx << "\n"; cout << " The Y grid spacing is " << dy << "\n"; // // Initialize the data. // rhs(nx, ny, f); // // Set the initial solution estimate. // We are "allowed" to pick up the boundary conditions exactly. // for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { if (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) { unew[i][j] = f[i][j]; } else { unew[i][j] = 0.0; } } } unew_norm = r8mat_rms(nx, ny, unew); // // Set up the exact solution. // for (j = 0; j < ny; j++) { y = (double)(j) / (double)(ny - 1); for (i = 0; i < nx; i++) { x = (double)(i) / (double)(nx - 1); uexact[i][j] = u_exact(x, y); } } u_norm = r8mat_rms(nx, ny, uexact); cout << " RMS of exact solution = " << u_norm << "\n"; // // Do the iteration. // converged = false; cout << "\n"; cout << " Step ||Unew|| ||Unew-U|| ||Unew-Exact||\n"; cout << "\n"; for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { udiff[i][j] = unew[i][j] - uexact[i][j]; } } error = r8mat_rms(nx, ny, udiff); cout << " " << setw(4) << 0 << " " << setw(14) << unew_norm << " " << " " << " " << setw(14) << error << "\n"; for (it = 1; it <= it_max; it++) { for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { u[i][j] = unew[i][j]; } } // // UNEW is derived from U by one Jacobi step. // sweep(nx, ny, dx, dy, f, u, unew); // // Check for convergence. // u_norm = unew_norm; unew_norm = r8mat_rms(nx, ny, unew); for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { udiff[i][j] = unew[i][j] - u[i][j]; } } diff = r8mat_rms(nx, ny, udiff); for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { udiff[i][j] = unew[i][j] - uexact[i][j]; } } error = r8mat_rms(nx, ny, udiff); if (diff <= tolerance) { cout << " " << setw(4) << it << " " << setw(14) << unew_norm << " " << setw(14) << diff << " " << setw(14) << error << "\n"; converged = true; break; } } if (converged) { cout << " The iteration has converged.\n"; } else { cout << " The iteration has NOT converged.\n"; } // // Terminate. // cout << "\n"; cout << "POISSON_SERIAL:\n"; cout << " Normal end of execution.\n"; cout << "\n"; //Free memory for (i = 0; i < nx; i++) { delete[] f[i]; delete[] u[i]; delete[] udiff[i]; delete[] uexact[i]; delete[] unew[i]; } delete[] f; delete[] u; delete[] udiff; delete[] uexact; delete[] unew; return 0; } //****************************************************************************80