dp2d_mpi.h

Flens库-一个在C++的矩阵运算库

#include <flens/flens.h>
#include <fftw3.h>
#include <poisson_flens/poisson_flens.h>


namespace flens {

//== problem set ===============================================================

MpiCart mpiCart;
int N;
DistributedGridVector2D u, f, sol, r, error;

void
setProblemSize(int n)
{
    N = n;
    u = f = sol = r = DistributedGridVector2D(mpiCart, N+1);
}

// problem 1: -u_xx -u_yy = 5*pi^2*sin(pi*x)*pi^2*sin(2*pi*y), BC: u = 0
// solution:   u(x) = sin(pi*x) * sin(2*pi*y)
void
problem1()
{
    double h = 1./f.rh;
    DistributedGridVector2D::LocalGrid F = f.localGrid();
    DistributedGridVector2D::LocalGrid SOL = sol.localGrid();

    for (int i=F.firstRow(); i<=F.lastRow(); ++i) {
        for (int j=F.firstCol(); j<=F.lastCol(); ++j) {
            double x = (i+f.i0)*h;
            double y = (j+f.j0)*h;
            F(i,j) = 5*M_PI*M_PI * sin(M_PI*x) * sin(2*M_PI*y);
            SOL(i,j) = sin(M_PI*x) * sin(2*M_PI*y);
        }
    }
}

// problem 2: -u''(x) = 0, BC: u = 1
// solution:   u(x) = 1
void
problem2()
{
    assert(0);
    // needs to be modified ... but time is short!
    /*
    f.grid = 0;
    u.grid(0,_) = 1;
    u.grid(N+1,_) = 1;
    u.grid(_,0) = 1;
    u.grid(_,N+1) = 1;
    sol.grid = 1;
    */
}

// problem 3: -u_xx -u_yy = 2*pi^2*sin(pi*x) * pi^2*sin(pi*y), BC: u = 0
// solution:   u(x) = sin(pi*x) * sin(pi*y)
void
problem3()
{
    double h = 1./f.rh;
    DistributedGridVector2D::LocalGrid F = f.localGrid();
    DistributedGridVector2D::LocalGrid SOL = sol.localGrid();

    for (int i=F.firstRow(); i<=F.lastRow(); ++i) {
        for (int j=F.firstCol(); j<=F.lastCol(); ++j) {
            double x = (i+f.i0)*h;
            double y = (j+f.j0)*h;
            F(i,j) = 2*M_PI*M_PI * sin(M_PI*x) * sin(M_PI*y);
            SOL(i,j) = sin(M_PI*x) * sin(M_PI*y);
        }
    }
}

//== error statistic ===========================================================

void
errorStat(int it=-1)
{
    DirichletPoisson2D A(N+1);
    r = f - A*u;
    error = sol - u;

    double rNormInf = normInf(r);
    double rNormL2 = normL2(r);
    double errorNormInf = normInf(error);
    double errorNormL2 = normL2(error);

    if ((mpiCart.row==0) && (mpiCart.col==0)) {
        if (it>=0) {
            std::cout.width(3);
            std::cout << it << ") | ";
        }

        std::cout.precision(12);
        std::cout.setf(std::ios::fixed);
        std::cout.width(17);
        std::cout << rNormL2 << " | ";

        std::cout.precision(12);
        std::cout.setf(std::ios::fixed);
        std::cout.width(17);
        std::cout << rNormInf << " | ";

        std::cout.precision(12);
        std::cout.setf(std::ios::fixed);
        std::cout.width(16);
        std::cout << errorNormL2 << " | ";

        std::cout.precision(12);
        std::cout.setf(std::ios::fixed);
        std::cout.width(16);
        std::cout << errorNormInf << " | " << std::endl;
    }
}

} // namespace flens