This is incorrect and will not work.
for (j = 0; j < daInfo.my; ++j) {
dirichlet_rows[n_dirichlet_rows++] = j * daInfo.mx
<https://urldefense.us/v3/__http://dainfo.mx/__;!!G_uCfscf7eWS!faAU5SYXVdbeL4y_c26LP7pbZnC1p3U_6-f4qNMpQMw2dXi4plc1QMmEh9_6S1xCYurJg4RIuXOrnKz85IkV6m3YcLROpRRf$>
* ncomp;
dirichlet_rows[n_dirichlet_rows++] = (j+1) * daInfo.mx
<https://urldefense.us/v3/__http://dainfo.mx/__;!!G_uCfscf7eWS!faAU5SYXVdbeL4y_c26LP7pbZnC1p3U_6-f4qNMpQMw2dXi4plc1QMmEh9_6S1xCYurJg4RIuXOrnKz85IkV6m3YcLROpRRf$>
* ncomp - ncomp;
}
You are assuming that the PETSc global numbering of the matrix rows/columns is
the same as the natural ordering (on a 2d mesh) across the entire mesh. It is
not, rather all nodes are numbered on the first MPI process, followed by all
the on second etc. Thus mapping between the PETSc ordering and the natural
ordering is cumbersome.
But, no worries, MatZeroRowsColumnsStencil() allows you to indicate the
rows/columns to zero using the same stencil information you use to fill the
matrix, so you never need to worry about the mapping between the PETSc ordering
and the global natural ordering.
Barry
> On Aug 9, 2024, at 4:02 PM, Junming Duan <junming.duan.m...@gmail.com> wrote:
>
> Dear all,
>
> I tried to use MatZeroRowsColumns to eliminate Dirichlet boundary nodes.
> However, it cannot eliminate correctly in parallel.
> Please see the attached code which uses DMDA to create the matrix.
> When I used one process, it works as expected.
> For two processes, the domain is split in the x direction. But the 10th row,
> 20th column is not eliminated as observed when using one process. The results
> for two processes are also attached.
> I have input the same rows to be eliminated for both processes.
> Thank you for any help.
>
> #include <stdio.h>
> #include <petscdmda.h>
>
> int main(int argc, char **argv)
> {
> PetscInt M = 5, N = 3, m = PETSC_DECIDE, n = PETSC_DECIDE, ncomp = 2;
> PetscInt i, j;
> DMDALocalInfo daInfo;
> DM da;
> Mat A;
> Vec x, b;
> MatStencil row, col[5];
> PetscScalar v[5];
> PetscInt n_dirichlet_rows = 0, dirichlet_rows[2*(M+N)];
>
> PetscFunctionBeginUser;
> PetscCall(PetscInitialize(&argc, &argv, (char *)0, NULL));
> PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED,
> DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, M, N, m, n, ncomp, 1, NULL, NULL,
> &da));
> PetscCall(DMSetFromOptions(da));
> PetscCall(DMSetUp(da));
>
> PetscCall(DMView(da, PETSC_VIEWER_STDOUT_WORLD));
>
> PetscCall(DMDAGetLocalInfo(da, &daInfo));
> PetscCall(DMSetMatrixPreallocateOnly(da, PETSC_TRUE));
>
> PetscCall(DMCreateMatrix(da, &A));
> PetscCall(MatCreateVecs(A, &x, &b));
> PetscCall(MatZeroEntries(A));
> PetscCall(VecZeroEntries(x));
> PetscCall(VecZeroEntries(b));
>
> for (j = daInfo.ys; j < daInfo.ys + daInfo.ym; ++j) {
> for (i = daInfo.xs; i < daInfo.xs + daInfo.xm; ++i) {
> row.j = j;
> row.i = i;
> row.c = 0;
>
> col[0].j = j;
> col[0].i = i;
> col[0].c = 0;
> v[0] = row.i + col[0].i + 2;
> col[1].j = j;
> col[1].i = i - 1;
> col[1].c = 0;
> v[1] = row.i + col[1].i + 2;
> col[2].j = j;
> col[2].i = i + 1;
> col[2].c = 0;
> v[2] = row.i + col[2].i + 2;
> col[3].j = j - 1;
> col[3].i = i;
> col[3].c = 0;
> v[3] = row.i + col[2].i + 2;
> col[4].j = j + 1;
> col[4].i = i;
> col[4].c = 0;
> v[4] = row.i + col[2].i + 2;
>
> PetscCall(MatSetValuesStencil(A, 1, &row, 5, col, v, ADD_VALUES));
>
> row.j = j;
> row.i = i;
> row.c = 1;
> col[0].j = j;
> col[0].i = i;
> col[0].c = 1;
> v[1] = row.j + col[1].j + 2;
> col[1].j = j - 1;
> col[1].i = i;
> col[1].c = 1;
> v[1] = row.j + col[1].j + 2;
> col[2].j = j + 1;
> col[2].i = i;
> col[2].c = 1;
> v[2] = row.j + col[2].j + 2;
> col[3].j = j;
> col[3].i = i - 1;
> col[3].c = 1;
> v[3] = row.j + col[1].j + 2;
> col[4].j = j;
> col[4].i = i + 1;
> col[4].c = 1;
> v[4] = row.j + col[2].j + 2;
>
> PetscCall(MatSetValuesStencil(A, 1, &row, 5, col, v, ADD_VALUES));
>
> }
> }
> PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
> PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
> MatView(A, 0);
>
> for (j = 0; j < daInfo.my; ++j) {
> dirichlet_rows[n_dirichlet_rows++] = j * daInfo.mx
> <https://urldefense.us/v3/__http://dainfo.mx/__;!!G_uCfscf7eWS!faAU5SYXVdbeL4y_c26LP7pbZnC1p3U_6-f4qNMpQMw2dXi4plc1QMmEh9_6S1xCYurJg4RIuXOrnKz85IkV6m3YcLROpRRf$>
> * ncomp;
> dirichlet_rows[n_dirichlet_rows++] = (j+1) * daInfo.mx
> <https://urldefense.us/v3/__http://dainfo.mx/__;!!G_uCfscf7eWS!faAU5SYXVdbeL4y_c26LP7pbZnC1p3U_6-f4qNMpQMw2dXi4plc1QMmEh9_6S1xCYurJg4RIuXOrnKz85IkV6m3YcLROpRRf$>
> * ncomp - ncomp;
> }
> PetscCall(PetscPrintf(PETSC_COMM_SELF, "n_dirichlet_rows: %d\n",
> n_dirichlet_rows));
> for (j = 0; j < n_dirichlet_rows; ++j) {
> PetscCall(PetscPrintf(PETSC_COMM_SELF, "%d, ", dirichlet_rows[j]));
> }
> PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
> PetscCall(MatZeroRowsColumns(A, n_dirichlet_rows, dirichlet_rows, 1, NULL,
> NULL));
> MatView(A, 0);
>
> PetscCall(VecDestroy(&x));
> PetscCall(VecDestroy(&b));
> PetscCall(DMDestroy(&da));
> PetscCall(PetscFinalize());
> return 0;
> }
>
> ———————————————————
>
> n_dirichlet_rows: 6
> 0, 8, 10, 18, 20, 28,
> Mat Object: 2 MPI processes
> type: mpiaij
> row 0: (0, 1.) (2, 0.) (10, 0.)
> row 1: (1, 2.) (3, 3.) (11, 3.)
> row 2: (0, 0.) (2, 4.) (4, 5.) (12, 0.)
> row 3: (1, 1.) (3, 4.) (5, 3.) (13, 3.)
> row 4: (2, 5.) (4, 6.) (6, 0.) (14, 0.)
> row 5: (3, 1.) (5, 6.) (7, 3.) (15, 3.)
> row 6: (4, 0.) (6, 1.) (8, 0.) (16, 0.)
> row 7: (5, 1.) (7, 8.) (9, 3.) (17, 3.)
> row 8: (6, 0.) (8, 1.) (18, 0.)
> row 9: (7, 1.) (9, 10.) (19, 3.)
> row 10: (0, 0.) (10, 2.) (12, 0.) (20, 3.)
> row 11: (1, 3.) (11, 2.) (13, 5.) (21, 5.)
> row 12: (2, 0.) (10, 0.) (12, 1.) (14, 0.) (22, 0.)
> row 13: (3, 3.) (11, 3.) (13, 4.) (15, 5.) (23, 5.)
> row 14: (4, 0.) (12, 0.) (14, 1.) (16, 0.) (24, 0.)
> row 15: (5, 3.) (13, 3.) (15, 6.) (17, 5.) (25, 5.)
> row 16: (6, 0.) (14, 0.) (16, 8.) (18, 9.) (26, 9.)
> row 17: (7, 3.) (15, 3.) (17, 8.) (19, 5.) (27, 5.)
> row 18: (8, 0.) (16, 9.) (18, 10.) (28, 0.)
> row 19: (9, 3.) (17, 3.) (19, 10.) (29, 5.)
> row 20: (10, 3.) (20, 2.) (22, 3.)
> row 21: (11, 5.) (21, 2.) (23, 7.)
> row 22: (12, 0.) (20, 3.) (22, 4.) (24, 5.)
> row 23: (13, 5.) (21, 5.) (23, 4.) (25, 7.)
> row 24: (14, 0.) (22, 5.) (24, 6.) (26, 7.)
> row 25: (15, 5.) (23, 5.) (25, 6.) (27, 7.)
> row 26: (16, 9.) (24, 7.) (26, 8.) (28, 0.)
> row 27: (17, 5.) (25, 5.) (27, 8.) (29, 7.)
> row 28: (18, 0.) (26, 0.) (28, 1.)
> row 29: (19, 5.) (27, 5.) (29, 10.)
>
> Junming
