当前位置：首页 > news >正文

使用Hypre BoomerAMG求解大规模泊松方程示例

news 来源：原创 2025/4/22 3:44:05

使用Hypre BoomerAMG求解大规模泊松方程示例

下面是一个使用Hypre库中的BoomerAMG预处理器求解大规模泊松方程的C++示例程序。这个示例假设你已经在系统上安装了Hypre库。

程序说明

使用有限体积法离散二维泊松方程
网格大小为7000x7000左右（约5000万自由度）
使用BoomerAMG作为预条件子的共轭梯度法求解

#include <iostream>
#include <cmath>
#include "HYPRE_struct_ls.h"
#include "HYPRE_krylov.h"// 定义网格大小
const int NX = 7000;  // x方向网格数
const int NY = 7000;  // y方向网格数
const int N = NX * NY; // 总自由度int main(int argc, char *argv[]) {// 初始化HypreHYPRE_Init();// 创建矩阵、向量和线性系统HYPRE_StructMatrix A;HYPRE_StructVector b, x;HYPRE_StructSolver solver, precond;// 1. 设置网格和离散化// 这里我们使用二维结构化网格int ndim = 2;int ilower[2] = {0, 0};int iupper[2] = {NX-1, NY-1};// 创建网格HYPRE_StructGrid grid;HYPRE_StructGridCreate(MPI_COMM_WORLD, ndim, &grid);HYPRE_StructGridSetExtents(grid, ilower, iupper);HYPRE_StructGridAssemble(grid);// 2. 创建矩阵结构// 定义5点模板（中心、上、下、左、右）int stencil_size = 5;int stencil_indices[5] = {0, 1, 2, 3, 4}; // 中心, 上, 下, 左, 右int offsets[5][2] = {{0,0}, {0,1}, {0,-1}, {-1,0}, {1,0}};HYPRE_StructStencil stencil;HYPRE_StructStencilCreate(ndim, stencil_size, &stencil);for (int i = 0; i < stencil_size; i++) {HYPRE_StructStencilSetElement(stencil, i, offsets[i]);}// 创建矩阵对象HYPRE_StructMatrixCreate(MPI_COMM_WORLD, grid, stencil, &A);HYPRE_StructMatrixInitialize(A);// 3. 组装矩阵（有限体积法离散泊松方程）// 假设均匀网格，Δx = Δy = hdouble h = 1.0 / (NX - 1);double h2 = h * h;double values[stencil_size];// 内部点使用标准5点模板values[0] = 4.0 / h2;  // 中心values[1] = -1.0 / h2;  // 上values[2] = -1.0 / h2;  // 下values[3] = -1.0 / h2;  // 左values[4] = -1.0 / h2;  // 右// 设置矩阵值（这里简化处理，实际应用中需要考虑边界条件）int nentries = stencil_size;int nvalues = N * nentries;double *all_values = new double[nvalues];int *rows = new int[N];for (int i = 0; i < N; i++) {rows[i] = i;for (int j = 0; j < nentries; j++) {all_values[i * nentries + j] = values[j];}}HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries, stencil_indices, all_values);HYPRE_StructMatrixAssemble(A);// 4. 创建和组装右端向量bHYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &b);HYPRE_StructVectorInitialize(b);double *b_values = new double[N];for (int i = 0; i < N; i++) {b_values[i] = 1.0; // 假设单位源项}HYPRE_StructVectorSetBoxValues(b, ilower, iupper, b_values);HYPRE_StructVectorAssemble(b);// 5. 创建初始猜测解向量xHYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &x);HYPRE_StructVectorInitialize(x);double *x_values = new double[N];for (int i = 0; i < N; i++) {x_values[i] = 0.0; // 初始猜测为0}HYPRE_StructVectorSetBoxValues(x, ilower, iupper, x_values);HYPRE_StructVectorAssemble(x);// 6. 设置求解器（PCG + BoomerAMG）HYPRE_StructPCGCreate(MPI_COMM_WORLD, &solver);// 设置BoomerAMG作为预条件子HYPRE_StructBoomerAMGCreate(&precond);// 设置BoomerAMG参数（这些参数可以根据问题调整）HYPRE_StructBoomerAMGSetMaxIter(precond, 1); // 作为预条件子，1次迭代足够HYPRE_StructBoomerAMGSetTol(precond, 0.0);HYPRE_StructBoomerAMGSetPrintLevel(precond, 1); // 打印基本信息HYPRE_StructBoomerAMGSetCoarsenType(precond, 6); // Falgout coarseningHYPRE_StructBoomerAMGSetRelaxType(precond, 3); // GS/Jacobi混合松弛HYPRE_StructBoomerAMGSetNumSweeps(precond, 1);HYPRE_StructBoomerAMGSetMaxLevels(precond, 20);// 设置PCG参数HYPRE_StructPCGSetMaxIter(solver, 200);HYPRE_StructPCGSetTol(solver, 1.0e-6);HYPRE_StructPCGSetTwoNorm(solver, 1);HYPRE_StructPCGSetPrintLevel(solver, 2); // 打印迭代信息HYPRE_StructPCGSetLogging(solver, 1);// 设置预条件子HYPRE_StructPCGSetPrecond(solver, (HYPRE_StructSolverPrecond) HYPRE_StructBoomerAMGSolve,(HYPRE_StructSolverPrecond) HYPRE_StructBoomerAMGSetup,precond);// 7. 求解线性系统HYPRE_StructPCGSetup(solver, A, b, x);HYPRE_StructPCGSolve(solver, A, b, x);// 8. 获取求解信息int num_iterations;double final_res_norm;HYPRE_StructPCGGetNumIterations(solver, &num_iterations);HYPRE_StructPCGGetFinalRelativeResidualNorm(solver, &final_res_norm);std::cout << "Number of iterations = " << num_iterations << std::endl;std::cout << "Final relative residual norm = " << final_res_norm << std::endl;// 9. 清理内存delete[] all_values;delete[] rows;delete[] b_values;delete[] x_values;HYPRE_StructGridDestroy(grid);HYPRE_StructStencilDestroy(stencil);HYPRE_StructMatrixDestroy(A);HYPRE_StructVectorDestroy(b);HYPRE_StructVectorDestroy(x);HYPRE_StructPCGDestroy(solver);HYPRE_StructBoomerAMGDestroy(precond);// 结束HypreHYPRE_Finalize();return 0;
}

编译和运行

编译命令（假设Hypre已正确安装）：

mpicxx -o hypre_poisson hypre_poisson.cpp -lHYPRE

运行命令（使用MPI并行）：

mpirun -np 4 ./hypre_poisson

关键点说明

网格离散：使用有限体积法离散二维泊松方程，生成5点模板的稀疏矩阵。
BoomerAMG设置：
- 使用Falgout coarsening (类型6)
- 混合GS/Jacobi松弛 (类型3)
- 设置适当的层数和松弛次数
求解器配置：
- 使用PCG（预条件共轭梯度法）作为外层迭代
- BoomerAMG作为预条件子
性能考虑：
- 对于5000万自由度的问题，需要足够的MPI进程
- 可能需要调整BoomerAMG参数以获得最佳性能
- 实际应用中应考虑更复杂边界条件的处理