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

Чтобы понять как работает PETSc [5], [12], рассмотрим простой не параллельный пример, который решает одномерную задачу Лапласа. Этот пример иллюстрирует решение системы линейных уравнений, полученной после стандартной конечно-разностной аппроксимации, с помощью методов подпространств Крылова.

В начале мы должны подключить заголовочный файл petscksp.h, чтобы использовать решатель KSP.

static char help[] = "Решает  трехдиагональную    систему    линейных    уравнений    с    использованием    KSP -   решателя.\n\n"; 
#include "petscksp.h"

Отметим, что этот файл автоматически подключает следующие заголовочные файлы:

Объявляем вектора приближенного решения x, правой части b, и точного решения u; матрицу линейной системы A; линейный решатель ksp и предобуславливатель pc; норму ошибки norm.

int main(int argc, char **args) 
{ 
    Vec            x, b, u; 
    Mat            A; 
    KSP            ksp; 
    PC             pc; 
    PetscReal      norm; 
    PetscErrorCode ierr; 
    PetscInt       i, n = 10, col[3], its; 
    PetscMPIInt    size; 
    PetscScalar    neg_one = -1.0, one = 1.0, value[3]; 
    PetscInitialize(&argc, &args, (char *)0, help); 
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr); 
    if (size != 1) 
        SETERRQ(1, "This is a uniprocessor example only!"); 
    ierr = PetscOptionsGetInt(PETSC_NULL, "-n", &n, PETSC_NULL);CHKERRQ(ierr);

Создаем векторы. Первый вектор создается с нуля, а остальные создаются дублированием.

    ierr = VecCreate(PETSC_COMM_WORLD, &x);CHKERRQ(ierr); 
    ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr); 
    ierr = VecSetSizes(x, PETSC_DECIDE, n);CHKERRQ(ierr); 
    ierr = VecSetFromOptions(x);CHKERRQ(ierr); 
 
    ierr = VecDuplicate(x, &b);CHKERRQ(ierr); 
    ierr = VecDuplicate(x, &u);CHKERRQ(ierr);

Создаем матрицу.

    ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr); 
    ierr = MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);CHKERRQ(ierr); 
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);

Заполняем и собираем матрицу

    value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 
    for (i=1; i<n-1; i++) { 
        col[0] = i-1; col[1] = i; col[2] = i+1; 
        ierr = MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);CHKERRQ(ierr); 
    } 
    i = n - 1; col[0] = n - 2; col[1] = n - 1; 
    ierr = MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);CHKERRQ(ierr); 
    i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0; 
    ierr = MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);CHKERRQ(ierr); 
 
    ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 
    ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

Задаем точное решение и затем вычисляем вектор правой части.

    ierr = VecSet(u, one);CHKERRQ(ierr); 
    ierr = MatMult(A, u, b);CHKERRQ(ierr);

Создаем линейный решатель

    ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr);

Задаем оператор, который определяет систему линейных уравнений.

    ierr = KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

Задаем предобуславливатель типа Якоби и точность

    ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 
    ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr); 
    ierr = KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);CHKERRQ(ierr); 
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

Решаем систему линейных уравнений

    ierr = KSPSolve(ksp, b, x);CHKERRQ(ierr);

Выводим информацию о решателе на экран

    ierr = KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

Находим норму ошибки вычисленного значения с точным

    ierr = VecAXPY(x, neg_one, u);CHKERRQ(ierr); 
    ierr = VecNorm(x, NORM_2, &norm);CHKERRQ(ierr);

Выводим на экран количество проделанных итераций

    ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %A, Iterations %D\n", 
                       norm, its);CHKERRQ(ierr);

Освобождаем память

    ierr = VecDestroy(x);CHKERRQ(ierr); 
    ierr = VecDestroy(u);CHKERRQ(ierr); 
    ierr = VecDestroy(b);CHKERRQ(ierr); 
    ierr = MatDestroy(A);CHKERRQ(ierr); 
    ierr = KSPDestroy(ksp);CHKERRQ(ierr); 
    ierr = PetscFinalize();CHKERRQ(ierr); 
    return 0; 
}