main.cpp

# 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