Actual source code: noflux_check.c

  1: static char help[] = "Check to see of DM_BOUNDARY_MIRROR works in 3D for DMDA with star stencil\n";

  3: #include "petscdmda.h"

  5: /* Contributed by Gourav Kumbhojkar */

  7: PetscErrorCode globalKMat_3d(Mat K, DMDALocalInfo info)
  8: {
  9:   MatStencil  row, col[7];
 10:   PetscScalar vals[7];
 11:   PetscInt    ncols;

 13:   PetscFunctionBeginUser;
 14:   for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
 15:     for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
 16:       for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
 17:         ncols = 0;
 18:         row.i = i;
 19:         row.j = j;
 20:         row.k = k;

 22:         col[0].i      = i;
 23:         col[0].j      = j;
 24:         col[0].k      = k;
 25:         vals[ncols++] = -6.; //ncols=1

 27:         col[ncols].i  = i - 1;
 28:         col[ncols].j  = j;
 29:         col[ncols].k  = k;
 30:         vals[ncols++] = 1.; //ncols=2

 32:         col[ncols].i  = i + 1;
 33:         col[ncols].j  = j;
 34:         col[ncols].k  = k;
 35:         vals[ncols++] = 1.; //ncols=3

 37:         col[ncols].i  = i;
 38:         col[ncols].j  = j - 1;
 39:         col[ncols].k  = k;
 40:         vals[ncols++] = 1.; //ncols=4

 42:         col[ncols].i  = i;
 43:         col[ncols].j  = j + 1;
 44:         col[ncols].k  = k;
 45:         vals[ncols++] = 1.; //ncols=5

 47:         col[ncols].i  = i;
 48:         col[ncols].j  = j;
 49:         col[ncols].k  = k + 1;
 50:         vals[ncols++] = 1.; //ncols=6

 52:         col[ncols].i  = i;
 53:         col[ncols].j  = j;
 54:         col[ncols].k  = k - 1;
 55:         vals[ncols++] = 1.; //ncols=7

 57:         PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
 58:       }
 59:     }
 60:   }
 61:   PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
 62:   PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: PetscErrorCode globalKMat_2d(Mat K, DMDALocalInfo info)
 67: {
 68:   MatStencil  row, col[5];
 69:   PetscScalar vals[5];
 70:   PetscInt    ncols;

 72:   PetscFunctionBeginUser;
 73:   for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
 74:     for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
 75:       ncols = 0;
 76:       row.i = i;
 77:       row.j = j;

 79:       col[0].i      = i;
 80:       col[0].j      = j;
 81:       vals[ncols++] = -4.; //ncols=1

 83:       col[ncols].i  = i - 1;
 84:       col[ncols].j  = j;
 85:       vals[ncols++] = 1.; //ncols=2

 87:       col[ncols].i  = i;
 88:       col[ncols].j  = j - 1;
 89:       vals[ncols++] = 1.; //ncols=3

 91:       col[ncols].i  = i + 1;
 92:       col[ncols].j  = j;
 93:       vals[ncols++] = 1.; //ncols=4

 95:       col[ncols].i  = i;
 96:       col[ncols].j  = j + 1;
 97:       vals[ncols++] = 1.; //ncols=5

 99:       PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
100:     }
101:   }
102:   PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
103:   PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: int main(int argc, char **argv)
108: {
109:   DM                     da3d, da2d;
110:   DMDALocalInfo          info3d, info2d;
111:   Mat                    K3d, K2d;
112:   PetscInt               ne, num_pts;
113:   ISLocalToGlobalMapping ltgm3d, ltgm2d;
114:   Vec                    row2d, row3d;
115:   PetscReal              norm2d, norm3d;

117:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
118:   ne      = 8;
119:   num_pts = ne + 1;

121:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2d));
122:   PetscCall(DMSetUp(da2d));
123:   PetscCall(DMDAGetLocalInfo(da2d, &info2d));
124:   PetscCall(DMCreateMatrix(da2d, &K2d));
125:   PetscCall(DMGetLocalToGlobalMapping(da2d, &ltgm2d));
126:   PetscCall(ISLocalToGlobalMappingView(ltgm2d, PETSC_VIEWER_STDOUT_WORLD));
127:   //PetscFinalize();
128:   PetscCall(globalKMat_2d(K2d, info2d));
129:   PetscCall(MatView(K2d, PETSC_VIEWER_STDOUT_WORLD));
130:   PetscCall(MatCreateVecs(K2d, &row2d, NULL));

132:   PetscCall(MatGetRowSum(K2d, row2d));
133:   PetscCall(VecNorm(row2d, NORM_2, &norm2d));

135:   PetscCheck(norm2d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "2D atrix row sum should be zero");
136:   PetscCall(VecDestroy(&row2d));
137:   PetscCall(MatDestroy(&K2d));
138:   PetscCall(DMDestroy(&da2d));

140:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, num_pts + 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &da3d));
141:   PetscCall(DMSetUp(da3d));
142:   PetscCall(DMCreateMatrix(da3d, &K3d));
143:   PetscCall(DMDAGetLocalInfo(da3d, &info3d));
144:   PetscCall(DMGetLocalToGlobalMapping(da3d, &ltgm3d));
145:   PetscCall(ISLocalToGlobalMappingView(ltgm3d, PETSC_VIEWER_STDOUT_WORLD));
146:   PetscCall(globalKMat_3d(K3d, info3d));
147:   PetscCall(MatView(K3d, PETSC_VIEWER_STDOUT_WORLD));
148:   PetscCall(MatCreateVecs(K3d, &row3d, NULL));
149:   PetscCall(MatGetRowSum(K3d, row3d));
150:   PetscCall(VecNorm(row3d, NORM_2, &norm3d));
151:   PetscCheck(norm3d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "3D atrix row sum should be zero");
152:   PetscCall(VecDestroy(&row3d));

154:   PetscCall(DMDestroy(&da3d));
155:   PetscCall(MatDestroy(&K3d));
156:   return PetscFinalize();
157: }

159: /*TEST

161:    test:

163:    test:
164:       suffix: 2
165:       nsize: 2

167:    test:
168:       suffix: 4
169:       nsize: 4

171:    test:
172:       suffix: 8
173:       nsize: 8

175: TEST*/