Actual source code: ex65.c

```  1: /*
2:  Partial differential equation

4:    d   d u = 1, 0 < x < 1,
5:    --   --
6:    dx   dx
7: with boundary conditions

9:    u = 0 for x = 0, x = 1

11:    This uses multigrid to solve the linear system

13:    Demonstrates how to build a DMSHELL for managing multigrid. The DMSHELL simply creates a
14:    DMDA1d to construct all the needed PETSc objects.

16: */

18: static char help[] = "Solves 1D constant coefficient Laplacian using DMSHELL and multigrid.\n\n";

20: #include <petscdm.h>
21: #include <petscdmda.h>
22: #include <petscdmshell.h>
23: #include <petscksp.h>

25: static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
26: static PetscErrorCode ComputeRHS(KSP, Vec, void *);
27: static PetscErrorCode CreateMatrix(DM, Mat *);
28: static PetscErrorCode CreateGlobalVector(DM, Vec *);
29: static PetscErrorCode CreateLocalVector(DM, Vec *);
30: static PetscErrorCode Refine(DM, MPI_Comm, DM *);
31: static PetscErrorCode Coarsen(DM, MPI_Comm, DM *);
32: static PetscErrorCode CreateInterpolation(DM, DM, Mat *, Vec *);
33: static PetscErrorCode CreateRestriction(DM, DM, Mat *);
34: static PetscErrorCode Destroy(void *);

36: static PetscErrorCode MyDMShellCreate(MPI_Comm comm, DM da, DM *shell)
37: {
38:   PetscFunctionBeginUser;
39:   PetscCall(DMShellCreate(comm, shell));
40:   PetscCall(DMShellSetContext(*shell, da));
41:   PetscCall(DMShellSetCreateMatrix(*shell, CreateMatrix));
42:   PetscCall(DMShellSetCreateGlobalVector(*shell, CreateGlobalVector));
43:   PetscCall(DMShellSetCreateLocalVector(*shell, CreateLocalVector));
44:   PetscCall(DMShellSetRefine(*shell, Refine));
45:   PetscCall(DMShellSetCoarsen(*shell, Coarsen));
46:   PetscCall(DMShellSetCreateInterpolation(*shell, CreateInterpolation));
47:   PetscCall(DMShellSetCreateRestriction(*shell, CreateRestriction));
48:   PetscCall(DMShellSetDestroyContext(*shell, Destroy));
49:   PetscFunctionReturn(PETSC_SUCCESS);
50: }

52: int main(int argc, char **argv)
53: {
54:   KSP      ksp;
55:   DM       da, shell;
56:   PetscInt levels;

58:   PetscFunctionBeginUser;
59:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
60:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
61:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 129, 1, 1, 0, &da));
62:   PetscCall(DMSetFromOptions(da));
63:   PetscCall(DMSetUp(da));
64:   PetscCall(MyDMShellCreate(PETSC_COMM_WORLD, da, &shell));
65:   /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */
66:   PetscCall(DMGetRefineLevel(da, &levels));
67:   PetscCall(DMSetRefineLevel(shell, levels));

69:   PetscCall(KSPSetDM(ksp, shell));
70:   PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, NULL));
71:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, NULL));
72:   PetscCall(KSPSetFromOptions(ksp));
73:   PetscCall(KSPSolve(ksp, NULL, NULL));

75:   PetscCall(KSPDestroy(&ksp));
76:   PetscCall(DMDestroy(&shell));
77:   PetscCall(PetscFinalize());
78:   return 0;
79: }

81: static PetscErrorCode Destroy(void *ctx)
82: {
83:   PetscFunctionBeginUser;
84:   PetscCall(DMDestroy((DM *)&ctx));
85:   PetscFunctionReturn(PETSC_SUCCESS);
86: }

88: static PetscErrorCode CreateMatrix(DM shell, Mat *A)
89: {
90:   DM da;

92:   PetscFunctionBeginUser;
93:   PetscCall(DMShellGetContext(shell, &da));
94:   PetscCall(DMCreateMatrix(da, A));
95:   PetscFunctionReturn(PETSC_SUCCESS);
96: }

98: static PetscErrorCode CreateInterpolation(DM dm1, DM dm2, Mat *mat, Vec *vec)
99: {
100:   DM da1, da2;

102:   PetscFunctionBeginUser;
103:   PetscCall(DMShellGetContext(dm1, &da1));
104:   PetscCall(DMShellGetContext(dm2, &da2));
105:   PetscCall(DMCreateInterpolation(da1, da2, mat, vec));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: static PetscErrorCode CreateRestriction(DM dm1, DM dm2, Mat *mat)
110: {
111:   DM  da1, da2;
112:   Mat tmat;

114:   PetscFunctionBeginUser;
115:   PetscCall(DMShellGetContext(dm1, &da1));
116:   PetscCall(DMShellGetContext(dm2, &da2));
117:   PetscCall(DMCreateInterpolation(da1, da2, &tmat, NULL));
118:   PetscCall(MatTranspose(tmat, MAT_INITIAL_MATRIX, mat));
119:   PetscCall(MatDestroy(&tmat));
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: static PetscErrorCode CreateGlobalVector(DM shell, Vec *x)
124: {
125:   DM da;

127:   PetscFunctionBeginUser;
128:   PetscCall(DMShellGetContext(shell, &da));
129:   PetscCall(DMCreateGlobalVector(da, x));
130:   PetscCall(VecSetDM(*x, shell));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: static PetscErrorCode CreateLocalVector(DM shell, Vec *x)
135: {
136:   DM da;

138:   PetscFunctionBeginUser;
139:   PetscCall(DMShellGetContext(shell, &da));
140:   PetscCall(DMCreateLocalVector(da, x));
141:   PetscCall(VecSetDM(*x, shell));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: static PetscErrorCode Refine(DM shell, MPI_Comm comm, DM *dmnew)
146: {
147:   DM da, dafine;

149:   PetscFunctionBeginUser;
150:   PetscCall(DMShellGetContext(shell, &da));
151:   PetscCall(DMRefine(da, comm, &dafine));
152:   PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dafine, dmnew));
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

156: static PetscErrorCode Coarsen(DM shell, MPI_Comm comm, DM *dmnew)
157: {
158:   DM da, dacoarse;

160:   PetscFunctionBeginUser;
161:   PetscCall(DMShellGetContext(shell, &da));
162:   PetscCall(DMCoarsen(da, comm, &dacoarse));
163:   PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dacoarse, dmnew));
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

167: static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
168: {
169:   PetscInt    mx, idx[2];
170:   PetscScalar h, v[2];
171:   DM          da, shell;

173:   PetscFunctionBeginUser;
174:   PetscCall(KSPGetDM(ksp, &shell));
175:   PetscCall(DMShellGetContext(shell, &da));
176:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
177:   h = 1.0 / (mx - 1);
178:   PetscCall(VecSet(b, h));
179:   idx[0] = 0;
180:   idx[1] = mx - 1;
181:   v[0] = v[1] = 0.0;
182:   PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES));
183:   PetscCall(VecAssemblyBegin(b));
184:   PetscCall(VecAssemblyEnd(b));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
189: {
190:   PetscInt    i, mx, xm, xs;
191:   PetscScalar v[3], h;
192:   MatStencil  row, col[3];
193:   DM          da, shell;

195:   PetscFunctionBeginUser;
196:   PetscCall(KSPGetDM(ksp, &shell));
197:   PetscCall(DMShellGetContext(shell, &da));
198:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
199:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
200:   h = 1.0 / (mx - 1);

202:   for (i = xs; i < xs + xm; i++) {
203:     row.i = i;
204:     if (i == 0 || i == mx - 1) {
205:       v[0] = 2.0 / h;
206:       PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
207:     } else {
208:       v[0]     = (-1.0) / h;
209:       col[0].i = i - 1;
210:       v[1]     = (2.0) / h;
211:       col[1].i = row.i;
212:       v[2]     = (-1.0) / h;
213:       col[2].i = i + 1;
214:       PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES));
215:     }
216:   }
217:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
218:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*TEST

224:    test:
225:       nsize: 4
226:       args: -ksp_monitor -pc_type mg -da_refine 3

228: TEST*/
```