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: {
40: DMShellCreate(comm,shell);
41: DMShellSetContext(*shell,da);
42: DMShellSetCreateMatrix(*shell,CreateMatrix);
43: DMShellSetCreateGlobalVector(*shell,CreateGlobalVector);
44: DMShellSetCreateLocalVector(*shell,CreateLocalVector);
45: DMShellSetRefine(*shell,Refine);
46: DMShellSetCoarsen(*shell,Coarsen);
47: DMShellSetCreateInterpolation(*shell,CreateInterpolation);
48: DMShellSetCreateRestriction(*shell,CreateRestriction);
49: return 0;
50: }
52: int main(int argc,char **argv)
53: {
55: KSP ksp;
56: DM da,shell;
57: PetscInt levels;
59: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
60: KSPCreate(PETSC_COMM_WORLD,&ksp);
61: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,129,1,1,0,&da);
62: DMSetFromOptions(da);
63: DMSetUp(da);
64: 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: DMGetRefineLevel(da,&levels);
67: DMSetRefineLevel(shell,levels);
69: KSPSetDM(ksp,shell);
70: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
71: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
72: KSPSetFromOptions(ksp);
73: KSPSolve(ksp,NULL,NULL);
75: KSPDestroy(&ksp);
76: DMDestroy(&shell);
77: DMDestroy(&da);
78: PetscFinalize();
79: return ierr;
80: }
82: static PetscErrorCode CreateMatrix(DM shell,Mat *A)
83: {
85: DM da;
87: DMShellGetContext(shell,&da);
88: DMCreateMatrix(da,A);
89: return 0;
90: }
92: static PetscErrorCode CreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
93: {
94: DM da1,da2;
97: DMShellGetContext(dm1,&da1);
98: DMShellGetContext(dm2,&da2);
99: DMCreateInterpolation(da1,da2,mat,vec);
100: return 0;
101: }
103: static PetscErrorCode CreateRestriction(DM dm1,DM dm2,Mat *mat)
104: {
105: DM da1,da2;
107: Mat tmat;
109: DMShellGetContext(dm1,&da1);
110: DMShellGetContext(dm2,&da2);
111: DMCreateInterpolation(da1,da2,&tmat,NULL);
112: MatTranspose(tmat,MAT_INITIAL_MATRIX,mat);
113: MatDestroy(&tmat);
114: return 0;
115: }
117: static PetscErrorCode CreateGlobalVector(DM shell,Vec *x)
118: {
120: DM da;
122: DMShellGetContext(shell,&da);
123: DMCreateGlobalVector(da,x);
124: VecSetDM(*x,shell);
125: return 0;
126: }
128: static PetscErrorCode CreateLocalVector(DM shell,Vec *x)
129: {
131: DM da;
133: DMShellGetContext(shell,&da);
134: DMCreateLocalVector(da,x);
135: VecSetDM(*x,shell);
136: return 0;
137: }
139: static PetscErrorCode Refine(DM shell,MPI_Comm comm,DM *dmnew)
140: {
142: DM da,dafine;
144: DMShellGetContext(shell,&da);
145: DMRefine(da,comm,&dafine);
146: MyDMShellCreate(PetscObjectComm((PetscObject)shell),dafine,dmnew);
147: return 0;
148: }
150: static PetscErrorCode Coarsen(DM shell,MPI_Comm comm,DM *dmnew)
151: {
153: DM da,dacoarse;
155: DMShellGetContext(shell,&da);
156: DMCoarsen(da,comm,&dacoarse);
157: MyDMShellCreate(PetscObjectComm((PetscObject)shell),dacoarse,dmnew);
158: /* discard an "extra" reference count to dacoarse */
159: DMDestroy(&dacoarse);
160: return 0;
161: }
163: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
164: {
166: PetscInt mx,idx[2];
167: PetscScalar h,v[2];
168: DM da,shell;
171: KSPGetDM(ksp,&shell);
172: DMShellGetContext(shell,&da);
173: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
174: h = 1.0/((mx-1));
175: VecSet(b,h);
176: idx[0] = 0; idx[1] = mx -1;
177: v[0] = v[1] = 0.0;
178: VecSetValues(b,2,idx,v,INSERT_VALUES);
179: VecAssemblyBegin(b);
180: VecAssemblyEnd(b);
181: return(0);
182: }
184: static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
185: {
187: PetscInt i,mx,xm,xs;
188: PetscScalar v[3],h;
189: MatStencil row,col[3];
190: DM da,shell;
193: KSPGetDM(ksp,&shell);
194: DMShellGetContext(shell,&da);
195: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
196: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
197: h = 1.0/(mx-1);
199: for (i=xs; i<xs+xm; i++) {
200: row.i = i;
201: if (i==0 || i==mx-1) {
202: v[0] = 2.0/h;
203: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
204: } else {
205: v[0] = (-1.0)/h;col[0].i = i-1;
206: v[1] = (2.0)/h;col[1].i = row.i;
207: v[2] = (-1.0)/h;col[2].i = i+1;
208: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
209: }
210: }
211: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
212: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
213: return(0);
214: }
216: /*TEST
218: test:
219: nsize: 4
220: args: -ksp_monitor -pc_type mg -da_refine 3
222: TEST*/