Fortran可分配数组与PETSc向量/矩阵高效转换技术问询

核心思路是将PETSc的核心对象（向量、矩阵、求解器）封装在独立模块中，仅初始化一次，后续时间步仅更新数据而非重建对象，同时通过内存直接访问或批量赋值的方式实现与Fortran可分配数组的高效交互，完全隔离现有代码与PETSc的细节。

1. 模块封装与资源复用 #

创建一个Fortran模块，将PETSc的Vec、Mat、KSP（求解器）等对象声明为模块全局变量（或封装在派生类型中），在模拟初始化阶段完成这些对象的创建与配置，时间步循环中仅更新数据，模拟结束时统一销毁。这彻底避免了每步创建/销毁大型对象的开销。

关键要点： #

• 若你的CFD问题网格固定，稀疏矩阵的非零元结构（行/列索引）通常不会变化，仅矩阵元素值随时间更新——这种场景下，PETSc矩阵只需初始化一次结构，后续仅更新数值即可，性能最优。
• 即使网格有动态调整，也可通过判断结构是否变化，选择性重建矩阵，而非每步重建。

2. 稀疏矩阵高效转换 #

假设你的现有代码用**COO格式（三个可分配数组：row_idx(:)、col_idx(:)、val(:)）**存储稀疏矩阵（多数CFD稀疏矩阵实现会采用这种格式），可按以下步骤处理：

初始化阶段（仅执行一次）： #

1. 调用MatCreate()创建PETSc矩阵对象，指定矩阵大小（全局行数/列数）。
2. 调用MatSetType()设置矩阵类型（串行用MATSEQAIJ，并行用MATMPIAIJ）。
3. 调用MatSetPreallocation()预分配非零元空间（可传入每行预估非零元数，或直接用COO数组的长度）。
4. 通过MatSetValues()将COO格式的行、列、值数组批量导入PETSc矩阵，再调用MatAssemblyBegin()和MatAssemblyEnd()完成矩阵组装。

时间步更新阶段（每步执行）： #

若矩阵结构未变，直接调用MatSetValues()传入更新后的val(:)数组（行、列索引复用初始化时的数组），再执行MatAssemblyBegin/End(MAT_FLUSH_ASSEMBLY)完成数值更新——此操作仅覆盖现有值，无需重新分配内存。

3. 向量数据高效传递 #

对于右端项B和解向量X，避免频繁调用VecSetValues()，而是通过内存直接映射的方式实现数据交互：

• 调用VecGetArrayF90()获取PETSc向量对应的Fortran数组指针，直接将现有可分配数组的数据拷贝到该指针指向的内存（或反向拷贝解向量数据）。
• 操作完成后调用VecRestoreArrayF90()释放指针。

这种方式跳过了PETSc的高层接口开销，性能接近原生数组操作。

4. 求解器调用封装 #

在初始化阶段完成KSP求解器的配置：

1. 调用KSPCreate()创建KSP对象。
2. 调用KSPSetOperators()关联PETSc矩阵。
3. 设置求解器类型为LSQR：KSPSetType(ksp, KSPLSQR)。
4. 配置收敛准则、最大迭代数、精度等参数（如KSPSetTolerances()）。

时间步循环中，仅需更新矩阵和右端项向量，然后调用KSPSolve(ksp, b_vec, x_vec)完成求解，最后将解向量的数据拷贝回你的Fortran可分配数组。

示例代码片段 #

module petsc_wrapper
  use petscvec
  use petscmat
  use petscsnes
  use petscksp
  implicit none

  ! PETSc核心对象（模块全局，仅初始化一次）
  type(Vec) :: b_petsc, x_petsc
  type(Mat) :: a_petsc
  type(KSP) :: ksp_solver

  ! 原代码的稀疏矩阵COO格式数组（假设已在其他模块定义，或作为参数传入）
  integer, allocatable :: row_idx(:), col_idx(:)
  real(8), allocatable :: mat_val(:), rhs(:), sol(:)
  integer :: n_global, nnz ! 全局矩阵大小、非零元数量

contains

  ! 初始化PETSc对象（模拟开始时调用）
  subroutine petsc_init()
    integer :: ierr

    ! 初始化PETSc环境
    call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
    CHKERRQ(ierr)

    ! 创建向量
    call VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n_global, b_petsc, ierr)
    call VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n_global, x_petsc, ierr)
    CHKERRQ(ierr)

    ! 创建稀疏矩阵
    call MatCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n_global, n_global, &
                      nnz, PETSC_NULL_INTEGER, a_petsc, ierr)
    call MatSetType(a_petsc, MATMPIAIJ, ierr)
    CHKERRQ(ierr)

    ! 导入初始矩阵数据
    call MatSetValues(a_petsc, nnz, row_idx, nnz, col_idx, mat_val, INSERT_VALUES, ierr)
    call MatAssemblyBegin(a_petsc, MAT_FINAL_ASSEMBLY, ierr)
    call MatAssemblyEnd(a_petsc, MAT_FINAL_ASSEMBLY, ierr)
    CHKERRQ(ierr)

    ! 配置KSP求解器（LSQR）
    call KSPCreate(PETSC_COMM_WORLD, ksp_solver, ierr)
    call KSPSetOperators(ksp_solver, a_petsc, a_petsc, ierr)
    call KSPSetType(ksp_solver, KSPLSQR, ierr)
    call KSPSetTolerances(ksp_solver, 1d-8, PETSC_DEFAULT, PETSC_DEFAULT, 1000, ierr)
    CHKERRQ(ierr)
  end subroutine petsc_init

  ! 时间步更新与求解（每时间步调用）
  subroutine petsc_solve_step()
    integer :: ierr
    real(8), pointer :: b_ptr(:), x_ptr(:)

    ! 更新矩阵数值（结构不变时）
    call MatSetValues(a_petsc, nnz, row_idx, nnz, col_idx, mat_val, INSERT_VALUES, ierr)
    call MatAssemblyBegin(a_petsc, MAT_FLUSH_ASSEMBLY, ierr)
    call MatAssemblyEnd(a_petsc, MAT_FLUSH_ASSEMBLY, ierr)
    CHKERRQ(ierr)

    ! 将原代码的右端项rhs拷贝到PETSc向量
    call VecGetArrayF90(b_petsc, b_ptr, ierr)
    b_ptr(:) = rhs(:)
    call VecRestoreArrayF90(b_petsc, b_ptr, ierr)
    call VecAssemblyBegin(b_petsc, ierr)
    call VecAssemblyEnd(b_petsc, ierr)
    CHKERRQ(ierr)

    ! 求解
    call KSPSolve(ksp_solver, b_petsc, x_petsc, ierr)
    CHKERRQ(ierr)

    ! 将解拷贝回原代码的sol数组
    call VecGetArrayF90(x_petsc, x_ptr, ierr)
    sol(:) = x_ptr(:)
    call VecRestoreArrayF90(x_petsc, x_ptr, ierr)
    CHKERRQ(ierr)
  end subroutine petsc_solve_step

  ! 销毁PETSc对象（模拟结束时调用）
  subroutine petsc_finalize()
    integer :: ierr

    call VecDestroy(b_petsc, ierr)
    call VecDestroy(x_petsc, ierr)
    call MatDestroy(a_petsc, ierr)
    call KSPDestroy(ksp_solver, ierr)
    call PetscFinalize(ierr)
    CHKERRQ(ierr)
  end subroutine petsc_finalize

end module petsc_wrapper

注意事项 #

• 确保PETSc的Fortran接口正确链接（编译时需指定PETSc的mpif90或petscfc编译器）。
• 并行场景下，需注意全局索引与本地索引的对应关系——若你的原代码是串行的，可改用VecCreateSeq()和MatCreateSeqAIJ()。
• 若矩阵结构随时间变化（如动态网格），需在结构变化时调用MatDestroy()重建矩阵，而非仅更新数值。

内容的提问来源于stack exchange，提问作者cutus_low