1 /*
2 Partial differential equation
3
4 d d u = 1, 0 < x < 1,
5 -- --
6 dx dx
7 with boundary conditions
8
9 u = 0 for x = 0, x = 1
10
11 This uses multigrid to solve the linear system
12
13 Demonstrates how to build a DMSHELL for managing multigrid. The DMSHELL simply creates a
14 DMDA1d to construct all the needed PETSc objects.
15
16 */
17
18 static char help[] = "Solves 1D constant coefficient Laplacian using DMSHELL and multigrid.\n\n";
19
20 #include <petscdm.h>
21 #include <petscdmda.h>
22 #include <petscdmshell.h>
23 #include <petscksp.h>
24
25 static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, PetscCtx);
26 static PetscErrorCode ComputeRHS(KSP, Vec, PetscCtx);
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 PetscCtxDestroyFn DestroyCtx;
35
MyDMShellCreate(MPI_Comm comm,DM da,DM * shell)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, DestroyCtx));
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
main(int argc,char ** argv)52 int main(int argc, char **argv)
53 {
54 KSP ksp;
55 DM da, shell;
56 PetscInt levels;
57
58 PetscFunctionBeginUser;
59 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
68
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));
74
75 PetscCall(KSPDestroy(&ksp));
76 PetscCall(DMDestroy(&shell));
77 PetscCall(PetscFinalize());
78 return 0;
79 }
80
DestroyCtx(PetscCtxRt ctx)81 static PetscErrorCode DestroyCtx(PetscCtxRt ctx)
82 {
83 PetscFunctionBeginUser;
84 PetscCall(DMDestroy((DM *)ctx));
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
CreateMatrix(DM shell,Mat * A)88 static PetscErrorCode CreateMatrix(DM shell, Mat *A)
89 {
90 DM da;
91
92 PetscFunctionBeginUser;
93 PetscCall(DMShellGetContext(shell, &da));
94 PetscCall(DMCreateMatrix(da, A));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
CreateInterpolation(DM dm1,DM dm2,Mat * mat,Vec * vec)98 static PetscErrorCode CreateInterpolation(DM dm1, DM dm2, Mat *mat, Vec *vec)
99 {
100 DM da1, da2;
101
102 PetscFunctionBeginUser;
103 PetscCall(DMShellGetContext(dm1, &da1));
104 PetscCall(DMShellGetContext(dm2, &da2));
105 PetscCall(DMCreateInterpolation(da1, da2, mat, vec));
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
CreateRestriction(DM dm1,DM dm2,Mat * mat)109 static PetscErrorCode CreateRestriction(DM dm1, DM dm2, Mat *mat)
110 {
111 DM da1, da2;
112 Mat tmat;
113
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 }
122
CreateGlobalVector(DM shell,Vec * x)123 static PetscErrorCode CreateGlobalVector(DM shell, Vec *x)
124 {
125 DM da;
126
127 PetscFunctionBeginUser;
128 PetscCall(DMShellGetContext(shell, &da));
129 PetscCall(DMCreateGlobalVector(da, x));
130 PetscCall(VecSetDM(*x, shell));
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
CreateLocalVector(DM shell,Vec * x)134 static PetscErrorCode CreateLocalVector(DM shell, Vec *x)
135 {
136 DM da;
137
138 PetscFunctionBeginUser;
139 PetscCall(DMShellGetContext(shell, &da));
140 PetscCall(DMCreateLocalVector(da, x));
141 PetscCall(VecSetDM(*x, shell));
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
Refine(DM shell,MPI_Comm comm,DM * dmnew)145 static PetscErrorCode Refine(DM shell, MPI_Comm comm, DM *dmnew)
146 {
147 DM da, dafine;
148
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 }
155
Coarsen(DM shell,MPI_Comm comm,DM * dmnew)156 static PetscErrorCode Coarsen(DM shell, MPI_Comm comm, DM *dmnew)
157 {
158 DM da, dacoarse;
159
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 }
166
ComputeRHS(KSP ksp,Vec b,PetscCtx ctx)167 static PetscErrorCode ComputeRHS(KSP ksp, Vec b, PetscCtx ctx)
168 {
169 PetscInt mx, idx[2];
170 PetscScalar h, v[2];
171 DM da, shell;
172
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 }
187
ComputeMatrix(KSP ksp,Mat J,Mat jac,PetscCtx ctx)188 static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, PetscCtx ctx)
189 {
190 PetscInt i, mx, xm, xs;
191 PetscScalar v[3], h;
192 MatStencil row, col[3];
193 DM da, shell;
194
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);
201
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 }
221
222 /*TEST
223
224 test:
225 nsize: 4
226 args: -ksp_monitor -pc_type mg -da_refine 3
227
228 TEST*/
229