5.6 Решение задачи Дирихле для уравнения Пуассона

Как и в главе 3.4 приведем пример параллельного решения [8], [4] уравнения Пуассона методом конечных разностей на стандартном семиточечном шаблоне:

Δu = - f.

Пусть f = xyz(1 - x)(1 - y)(1 - z), граничные условия нулевые, размерность сетки n3.

В начале включаем заголовочные файлы необходимых пакетов из Trilinos.

#include "AztecOO.h" 
#include "Epetra_Vector.h" 
#include "Epetra_CrsMatrix.h" 
#include "Epetra_LinearProblem.h" 
#include "ml_MultiLevelPreconditioner.h"

Объявляем переменные описывающие вычислительную сетку и задаем размерность.

int main(int argc, char *argv[]) { 
    double h = 0.001; 
    int nx = 100; 
    int ny = 100; 
    int nz = 100; 
    int n3 = nx * ny * nz; 
    double Lx = nx * h; 
    double Ly = ny * h; 
    double Lz = nz * h;

Инициализируем комуникаторы MPI и Epetra.

    MPI_Init(&argc, &argv); 
    Epetra_MpiComm Comm(MPI_COMM_WORLD);

Инициализируем схемы распрeделения данных.

    Epetra_Map *Map = new Epetra_Map(n3, 0, Comm); 
    unsigned int NumMyElements = Map->NumMyElements(); 
    int *MyGlobalElements = new int[NumMyElements]; 
    Map->MyGlobalElements(MyGlobalElements);

Создаем матрицу описывающую линейный оператор решаемой системы уравнений.

    int *NumEntriesPerRow = new int[n3]; 
    for (unsigned int i = 0; i < NumMyElements; i++) { 
        unsigned int v = MyGlobalElements[i]; 
        unsigned int iz = v / nx / ny; 
        unsigned int iy = (v - nx * ny * iz) / nx; 
        unsigned int ix = v - nx * ny * iz - nx * iy; 
        unsigned int count = 7; 
        if ((ix == 0) || (ix == nx - 1)) count--; 
        if ((iy == 0) || (iy == ny - 1)) count--; 
        if ((iz == 0) || (iz == nz - 1)) count--; 
        NumEntriesPerRow[i] = count; 
    } 
    Epetra_CrsMatrix *A = new Epetra_CrsMatrix(Copy, *Map, *NumEntriesPerRow); 
    int indices[7]; 
    double values[7]; 
    values[0] = -6.0; 
    values[1] =  1.0; 
    values[2] =  1.0; 
    values[3] =  1.0; 
    values[4] =  1.0; 
    values[5] =  1.0; 
    values[6] =  1.0; 
    for (unsigned int i = 0; i < NumMyElements; i++) 
    { 
        unsigned int v = MyGlobalElements[i]; 
        unsigned int iz = v / nx / ny; 
        unsigned int iy = (v - nx * ny * iz) / nx; 
        unsigned int ix = v - nx * ny * iz - nx * iy; 
        unsigned int count = 0; 
        indices[count++] = v; 
        if (ix !=      0) indices[count++] = v - 1; 
        if (ix != nx - 1) indices[count++] = v + 1; 
        if (iy !=      0) indices[count++] = v - nx; 
        if (iy != ny - 1) indices[count++] = v + nx; 
        if (iz !=      0) indices[count++] = v - nx * ny; 
        if (iz != nz - 1) indices[count++] = v + nx * ny; 
        A->InsertGlobalValues(v, count, values, indices); 
    } 
    A->FillComplete();

Создаем вектор, в котором будет записываться решение системы.

    Epetra_Vector x(*Map, false); 
    x.PutScalar(0.0);

Создаем вектор для правой части и инициализируем его.

    Epetra_Vector b(*Map, false); 
    for (unsigned int v = 0; v < NumMyElements; v++) 
    { 
        unsigned int iz = v / nx / ny; 
        unsigned int iy = (v - nx * ny * iz) / nx; 
        unsigned int ix = v - nx * ny * iz - nx * iy; 
        int index = MyGlobalElements[v]; 
        double z = iz * h; 
        double y = iy * h; 
        double x = ix * h; 
        double value = h * h * x * y * z * (Lx - x) * (Ly - y) * (Lz - z); 
        b.ReplaceGlobalValues(1, &value, &index); 
    }

Настраиваем и задаем предобуславливатель.

    Teuchos::ParameterList MLList; 
    ML_Epetra::SetDefaults("SA", MLList); 
    MLList.set("ML output", 10); 
    MLList.set("max levels", 5); 
    MLList.set("increasing or decreasing", "increasing"); 
    MLList.set("aggregation: type", "Uncoupled"); 
    MLList.set("smoother: type", "Chebyshev"); 
    MLList.set("smoother: sweeps", 3); 
    MLList.set("smoother: pre or post", "both"); 
    MLList.set("coarse: type", "Jacobi"); 
    ML_Epetra::MultiLevelPreconditioner* MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList, true); 
    MLPrec->PrintUnused(0);

Инициализируем линейный решатель и находим решение для составленной системы.

    Epetra_LinearProblem Problem(A, &x, &b); 
    AztecOO solver(Problem); 
    solver.SetPrecOperator(MLPrec); 
    solver.SetAztecOption(AZ_solver, AZ_cg); 
    //solver.SetAztecOption(AZ_output, 32); 
    solver.Iterate(500, 1e-7);

Удаляем объекты и завершаем работу программы.

    delete Map; 
    delete MyGlobalElements; 
    delete NumEntriesPerRow; 
    delete A; 
    delete MLPrec; 
    MPI_Finalize(); 
    return(EXIT_SUCCESS); 
}