为何设置PETSc向量值后DMDAVecGetArray无法读取修改值？

问题 #

以下程序试图设置PETSc向量的值，期望调用DMDAVecGetArray能读取到设置的内容，但在调用DMDAVecRestoreArray后，设置的值并未被保留：

#include <petscdmda.h>
#include <iostream>

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, (char*)0, NULL);

    DM da;
    Vec vec;
    PetscScalar *array1, *array2;

    DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
                 2, 2, 2, // grid dimensions
                 1, 1, PETSC_DECIDE, // number of dof, stencil width, number of processors in each dimension
                 9, 0, // stencil type, boundary type
                 NULL, NULL, NULL, // number of nodes in each dimension on each processor
                 &da);

    DMSetFromOptions(da);
    DMSetUp(da);
    DMCreateGlobalVector(da, &vec);

    DMDAVecGetArray(da, vec, &array1);
    array1[0] = 12345;
    // returns 12345:
    std::cout << "First value after setting: " << array1[0] << std::endl;
    DMDAVecRestoreArray(da, vec, &array1);

    DMDAVecGetArray(da, vec, &array2);
    // should also return 12345, but returns 4.63557e-310:
    std::cout << "First value after restoring and getting again: " << array2[0] << std::endl;
    DMDAVecRestoreArray(da, vec, &array2);

    VecDestroy(&vec);
    DMDestroy(&da);
    PetscFinalize();

    return 0;
}

编译命令：

g++ -o test test.cpp -I$PETSC_DIR/include -I$PETSC_DIR/$PETSC_ARCH/include -L$PETSC_DIR/$PETSC_ARCH/lib -lpetsc

原因分析 #

核心问题是**DMDACreate3d的参数传递顺序完全错误**，导致DMDA对象初始化异常，进而使得后续的向量内存管理逻辑失效。

DMDACreate3d的正确参数顺序为：

PetscErrorCode DMDACreate3d(
    MPI_Comm comm,
    DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz,
    DMDAStencilType stencil_type,
    PetscInt M, PetscInt N, PetscInt P, // 全局网格的x/y/z维度
    PetscInt m, PetscInt n, PetscInt p, // 每个维度的处理器数量
    PetscInt dof,                       // 每个网格点的自由度
    PetscInt s,                         // 模板宽度
    PetscInt *lx, PetscInt *ly, PetscInt *lz, // 本地网格维度（可选）
    DM *da
);

你的代码中错误地将处理器数量和自由度/模板宽度的位置颠倒，还传入了无效的9作为自由度参数，导致DMDA创建的向量存储结构异常，修改的数据无法被正确持久化。

修复方案 #

修正DMDACreate3d的参数顺序，根据需求传入正确的参数。例如，若需要每个网格点1个自由度、模板宽度为1，且让PETSc自动分配处理器数量，修改后的调用如下：

DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
             2, 2, 2, // 全局网格x/y/z维度
             PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, // 自动分配各维度处理器数
             1, // 每个网格点的自由度
             1, // 模板宽度
             NULL, NULL, NULL, // 自动分配本地网格维度
             &da);

多进程场景下，修改向量数据后需调用VecAssemblyBegin和VecAssemblyEnd完成跨进程数据组装，单进程场景下可省略这两步。

修正后的完整代码：

#include <petscdmda.h>
#include <iostream>

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, (char*)0, NULL);

    DM da;
    Vec vec;
    PetscScalar *array1, *array2;

    // 修正参数顺序的DMDACreate3d调用
    DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
                 2, 2, 2, // 全局网格x/y/z维度
                 PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, // 自动分配各维度处理器数
                 1, // 每个网格点的自由度
                 1, // 模板宽度
                 NULL, NULL, NULL, // 自动分配本地网格维度
                 &da);

    DMSetFromOptions(da);
    DMSetUp(da);
    DMCreateGlobalVector(da, &vec);

    DMDAVecGetArray(da, vec, &array1);
    array1[0] = 12345;
    std::cout << "First value after setting: " << array1[0] << std::endl;
    DMDAVecRestoreArray(da, vec, &array1);

    // 多进程场景下需添加以下两行，单进程可省略
    // VecAssemblyBegin(vec);
    // VecAssemblyEnd(vec);

    DMDAVecGetArray(da, vec, &array2);
    std::cout << "First value after restoring and getting again: " << array2[0] << std::endl;
    DMDAVecRestoreArray(da, vec, &array2);

    VecDestroy(&vec);
    DMDestroy(&da);
    PetscFinalize();

    return 0;
}

编译运行后，两次输出的数值都会是12345，符合预期。

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