Actual source code: ex65.c
2: /*
3: Partial differential equation
5: d d u = 1, 0 < x < 1,
6: -- --
7: dx dx
8: with boundary conditions
10: u = 0 for x = 0, x = 1
12: This uses multigrid to solve the linear system
14: Demonstrates how to build a DMSHELL for managing multigrid. The DMSHELL simply creates a
15: DMDA1d to construct all the needed PETSc objects.
17: */
19: static char help[] = "Solves 1D constant coefficient Laplacian using DMSHELL and multigrid.\n\n";
21: #include <petscdm.h>
22: #include <petscdmda.h>
23: #include <petscdmshell.h>
24: #include <petscksp.h>
26: static PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
27: static PetscErrorCode ComputeRHS(KSP,Vec,void*);
28: static PetscErrorCode CreateMatrix(DM,Mat*);
29: static PetscErrorCode CreateGlobalVector(DM,Vec*);
30: static PetscErrorCode CreateLocalVector(DM,Vec*);
31: static PetscErrorCode Refine(DM,MPI_Comm,DM*);
32: static PetscErrorCode Coarsen(DM,MPI_Comm,DM*);
33: static PetscErrorCode CreateInterpolation(DM,DM,Mat*,Vec*);
34: static PetscErrorCode CreateRestriction(DM,DM,Mat*);
36: static PetscErrorCode MyDMShellCreate(MPI_Comm comm,DM da,DM *shell)
37: {
39: DMShellCreate(comm,shell);
40: DMShellSetContext(*shell,da);
41: DMShellSetCreateMatrix(*shell,CreateMatrix);
42: DMShellSetCreateGlobalVector(*shell,CreateGlobalVector);
43: DMShellSetCreateLocalVector(*shell,CreateLocalVector);
44: DMShellSetRefine(*shell,Refine);
45: DMShellSetCoarsen(*shell,Coarsen);
46: DMShellSetCreateInterpolation(*shell,CreateInterpolation);
47: DMShellSetCreateRestriction(*shell,CreateRestriction);
48: return 0;
49: }
51: int main(int argc,char **argv)
52: {
53: KSP ksp;
54: DM da,shell;
55: PetscInt levels;
57: PetscInitialize(&argc,&argv,(char*)0,help);
58: KSPCreate(PETSC_COMM_WORLD,&ksp);
59: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,129,1,1,0,&da);
60: DMSetFromOptions(da);
61: DMSetUp(da);
62: MyDMShellCreate(PETSC_COMM_WORLD,da,&shell);
63: /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */
64: DMGetRefineLevel(da,&levels);
65: DMSetRefineLevel(shell,levels);
67: KSPSetDM(ksp,shell);
68: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
69: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
70: KSPSetFromOptions(ksp);
71: KSPSolve(ksp,NULL,NULL);
73: KSPDestroy(&ksp);
74: DMDestroy(&shell);
75: DMDestroy(&da);
76: PetscFinalize();
77: return 0;
78: }
80: static PetscErrorCode CreateMatrix(DM shell,Mat *A)
81: {
82: DM da;
84: DMShellGetContext(shell,&da);
85: DMCreateMatrix(da,A);
86: return 0;
87: }
89: static PetscErrorCode CreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
90: {
91: DM da1,da2;
93: DMShellGetContext(dm1,&da1);
94: DMShellGetContext(dm2,&da2);
95: DMCreateInterpolation(da1,da2,mat,vec);
96: return 0;
97: }
99: static PetscErrorCode CreateRestriction(DM dm1,DM dm2,Mat *mat)
100: {
101: DM da1,da2;
102: Mat tmat;
104: DMShellGetContext(dm1,&da1);
105: DMShellGetContext(dm2,&da2);
106: DMCreateInterpolation(da1,da2,&tmat,NULL);
107: MatTranspose(tmat,MAT_INITIAL_MATRIX,mat);
108: MatDestroy(&tmat);
109: return 0;
110: }
112: static PetscErrorCode CreateGlobalVector(DM shell,Vec *x)
113: {
114: DM da;
116: DMShellGetContext(shell,&da);
117: DMCreateGlobalVector(da,x);
118: VecSetDM(*x,shell);
119: return 0;
120: }
122: static PetscErrorCode CreateLocalVector(DM shell,Vec *x)
123: {
124: DM da;
126: DMShellGetContext(shell,&da);
127: DMCreateLocalVector(da,x);
128: VecSetDM(*x,shell);
129: return 0;
130: }
132: static PetscErrorCode Refine(DM shell,MPI_Comm comm,DM *dmnew)
133: {
134: DM da,dafine;
136: DMShellGetContext(shell,&da);
137: DMRefine(da,comm,&dafine);
138: MyDMShellCreate(PetscObjectComm((PetscObject)shell),dafine,dmnew);
139: return 0;
140: }
142: static PetscErrorCode Coarsen(DM shell,MPI_Comm comm,DM *dmnew)
143: {
144: DM da,dacoarse;
146: DMShellGetContext(shell,&da);
147: DMCoarsen(da,comm,&dacoarse);
148: MyDMShellCreate(PetscObjectComm((PetscObject)shell),dacoarse,dmnew);
149: /* discard an "extra" reference count to dacoarse */
150: DMDestroy(&dacoarse);
151: return 0;
152: }
154: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
155: {
156: PetscInt mx,idx[2];
157: PetscScalar h,v[2];
158: DM da,shell;
161: KSPGetDM(ksp,&shell);
162: DMShellGetContext(shell,&da);
163: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
164: h = 1.0/((mx-1));
165: VecSet(b,h);
166: idx[0] = 0; idx[1] = mx -1;
167: v[0] = v[1] = 0.0;
168: VecSetValues(b,2,idx,v,INSERT_VALUES);
169: VecAssemblyBegin(b);
170: VecAssemblyEnd(b);
171: return 0;
172: }
174: static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
175: {
176: PetscInt i,mx,xm,xs;
177: PetscScalar v[3],h;
178: MatStencil row,col[3];
179: DM da,shell;
182: KSPGetDM(ksp,&shell);
183: DMShellGetContext(shell,&da);
184: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
185: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
186: h = 1.0/(mx-1);
188: for (i=xs; i<xs+xm; i++) {
189: row.i = i;
190: if (i==0 || i==mx-1) {
191: v[0] = 2.0/h;
192: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
193: } else {
194: v[0] = (-1.0)/h;col[0].i = i-1;
195: v[1] = (2.0)/h;col[1].i = row.i;
196: v[2] = (-1.0)/h;col[2].i = i+1;
197: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
198: }
199: }
200: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
201: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
202: return 0;
203: }
205: /*TEST
207: test:
208: nsize: 4
209: args: -ksp_monitor -pc_type mg -da_refine 3
211: TEST*/