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

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

Δu = - f.
Для решения системы линейных алгебраических уравнений применим метод PCG с предобуславливателем AMG и используем IJ.

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

В начале программы подключаются заголовочные файлы библиотеки Hypre, чтобы использовать необходимые решатели и интерфейс представления данных.

#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include "_hypre_utilities.h" 
#include "HYPRE_krylov.h" 
#include "HYPRE.h" 
#include "HYPRE_parcsr_ls.h"

Зададим точность приближенного решения

 #define tol 1e-7

Определим необходимые переменные сеточных данных; переменные для разделения на процессы; матрицу линейной системы, A, вектор правой части, b, вектор приближенного решения, x, линейный решатель, solver, предобуславливатель, precond и инициализируем MPI

 
int main (int argc, char *argv[]){ 
   int i; 
   int N, n; 
   double h, h2; 
   int myid, num_procs; 
   int ilower, iupper; 
   int local_size; 
 
   HYPRE_IJMatrix A; 
   HYPRE_ParCSRMatrix parcsr_A; 
   HYPRE_IJVector b; 
   HYPRE_ParVector par_b; 
   HYPRE_IJVector x; 
   HYPRE_ParVector par_x; 
   HYPRE_Solver solver, precond; 
 
   MPI_Init(&argc, &argv); 
   MPI_Comm_rank(MPI_COMM_WORLD, &myid); 
   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

Зададим параметры расчетной сетки

   n = 10; 
   N = n*n*n; 
   h = 1.0+(n+1); 
   h2 = h*h;

Каждый процесс решает только свой блок неизвестных, который задается параметрами ilower и iuppper. Количество неизвестных для текущего процесса задано в переменной local_size.

   local_size = N+num_procs; 
 
   ilower = local_size*myid; 
   iupper = local_size*(myid+1) - 1;

Создадим матрицу c использованием интерфейса IJ, где в качестве формата хранения выбран parallel csr format storage.

   HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A); 
   HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR); 
   HYPRE_IJMatrixInitialize(A);

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

      int nnz; 
      double values[7]; 
      int cols[7]; 
 
      for (i = ilower; i <= iupper; i++) 
      { 
          nnz = 0; 
          if ( (i-n*n)>=0 ){ 
              cols[nnz] = i-n*n; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          if ( (i-n)>=0 ) { 
              cols[nnz] = i-n; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          if ( (i-1)>=0 ) { 
              cols[nnz] = i-1; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          cols[nnz] = i; 
          values[nnz] = 6.0; 
          nnz++; 
          if ( (i+1)<n ) { 
              cols[nnz] = i+1; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          if ( (i+n)<n ) { 
              cols[nnz] = i+n; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          if ( (i+n*n)<n ) { 
              cols[nnz] = i+n*n; 
              values[nnz] = -1.0; 
              nnz++; 
          } 
          HYPRE_IJMatrixSetValues(A, 1, &nnz, &i, cols, values); 
          HYPRE_IJMatrixAssemble(A); 
          HYPRE_IJMatrixGetObject(A, (void**) &parcsr_A); 
      }

Зададим правую часть, b и начальное значение решения, x

      HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&b); 
      HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR); 
      HYPRE_IJVectorInitialize(b); 
 
      HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&x); 
      HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR); 
      HYPRE_IJVectorInitialize(x); 
 
      double *rhs_values, *x_values; 
      int    *rows; 
 
      rhs_values = calloc(local_size, sizeof(double)); 
      x_values = calloc(local_size, sizeof(double)); 
      rows = calloc(local_size, sizeof(int)); 
 
      int ii = 0; 
      for (i = ilower; i <= iupper; i++){ 
         int iz = i + n + n; 
         int iy = (i - n * n * iz) + n; 
         int ix = i - n * n * iz - n * iy; 
         double vz = iz * h; 
         double vy = iy * h; 
         double vx = ix * h; 
         double vals = - h2 * vx * vy * vz * (1.0 - vx) * (1.0 - vy) * (1.0 - vz); 
         rhs_values[ii] = vals; 
 
         x_values[ii] = 0.0; 
         rows[ii] = i; 
         ii++; 
      } 
 
      HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values); 
      HYPRE_IJVectorSetValues(x, local_size, rows, x_values); 
 
      free(x_values); 
      free(rhs_values); 
      free(rows); 
 
      HYPRE_IJVectorAssemble(b); 
      HYPRE_IJVectorGetObject(b, (void **) &par_b); 
 
      HYPRE_IJVectorAssemble(x); 
      HYPRE_IJVectorGetObject(x, (void **) &par_x);

Определим решатель PCG с предобуславливателем AMG:

      int num_iterations; 
      double final_res_norm; 
 
      HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver); 
 
      HYPRE_PCGSetMaxIter(solver, 20); /* max iterations */ 
      HYPRE_PCGSetTol(solver, tol); /* conv. tolerance */ 
      HYPRE_PCGSetTwoNorm(solver, 1); /* use the two norm as the stopping criteria */ 
      HYPRE_PCGSetPrintLevel(solver, 2); /* print solve info */ 
      HYPRE_PCGSetLogging(solver, 1); /* needed to get run info later */ 
 
      HYPRE_BoomerAMGCreate(&precond); 
      HYPRE_BoomerAMGSetPrintLevel(precond, 1); /* print amg solution info */ 
      HYPRE_BoomerAMGSetStrongThreshold(precond,0.5); 
      HYPRE_BoomerAMGSetCoarsenType(precond, 6); 
      HYPRE_BoomerAMGSetRelaxType(precond, 0); /* Sym G.S.+Jacobi hybrid */ 
      HYPRE_BoomerAMGSetNumSweeps(precond, 1); 
      HYPRE_BoomerAMGSetTol(precond, 0.0); /* conv. tolerance zero */ 
      HYPRE_BoomerAMGSetMaxIter(precond, 1); /* do only one iteration! */ 
      HYPRE_BoomerAMGSetCycleType(precond,1); 
 
      HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve, 
                          (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);

Решаем уравнение

      HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x); 
      HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);

Выводим количество итераций и норму невязки

      HYPRE_PCGGetNumIterations(solver, &num_iterations); 
      HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm); 
      if (myid == 0){ 
         printf("\n"); 
         printf("Iterations = %d\n", num_iterations); 
         printf("Final Relative Residual Norm = %e\n", final_res_norm); 
         printf("\n"); 
      }

Сохраним решение в файле

      HYPRE_IJVectorPrint(x, "ij.out.x");

Очистим память

      HYPRE_ParCSRPCGDestroy(solver); 
      HYPRE_BoomerAMGDestroy(precond); 
 
      HYPRE_IJMatrixDestroy(A); 
      HYPRE_IJVectorDestroy(b); 
      HYPRE_IJVectorDestroy(x); 
 
      MPI_Finalize(); 
 
      return(0); 
}