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: }

 54: int main(int argc,char **argv)
 55: {
 57:   KSP            ksp;
 58:   DM             da,shell;
 59:   PetscInt       levels;

 61:   PetscInitialize(&argc,&argv,(char*)0,help);

 63:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 64:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-129,1,1,0,&da);
 65:   MyDMShellCreate(PETSC_COMM_WORLD,da,&shell);
 66:   /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */
 67:   DMGetRefineLevel(da,&levels);
 68:   DMSetRefineLevel(shell,levels);

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

 76:   KSPDestroy(&ksp);
 77:   DMDestroy(&shell);
 78:   DMDestroy(&da);
 79:   PetscFinalize();

 81:   return 0;
 82: }

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

 91:   DMShellGetContext(shell,(void**)&da);
 92:   DMCreateMatrix(da,A);
 93:   return 0;
 94: }

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

103:   DMShellGetContext(dm1,(void**)&da1);
104:   DMShellGetContext(dm2,(void**)&da2);
105:   DMCreateInterpolation(da1,da2,mat,vec);
106:   return 0;
107: }

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

116:   DMShellGetContext(dm1,(void**)&da1);
117:   DMShellGetContext(dm2,(void**)&da2);
118:   DMCreateInterpolation(da1,da2,mat,NULL);
119:   return 0;
120: }

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

129:   DMShellGetContext(shell,(void**)&da);
130:   DMCreateGlobalVector(da,x);
131:   VecSetDM(*x,shell);
132:   return 0;
133: }

137: static PetscErrorCode CreateLocalVector(DM shell,Vec *x)
138: {
140:   DM             da;

142:   DMShellGetContext(shell,(void**)&da);
143:   DMCreateLocalVector(da,x);
144:   VecSetDM(*x,shell);
145:   return 0;
146: }

150: static PetscErrorCode Refine(DM shell,MPI_Comm comm,DM *dmnew)
151: {
153:   DM             da,dafine;

155:   DMShellGetContext(shell,(void**)&da);
156:   DMRefine(da,comm,&dafine);
157:   MyDMShellCreate(PetscObjectComm((PetscObject)shell),dafine,dmnew);
158:   return 0;
159: }

163: static PetscErrorCode Coarsen(DM shell,MPI_Comm comm,DM *dmnew)
164: {
166:   DM             da,dacoarse;

168:   DMShellGetContext(shell,(void**)&da);
169:   DMCoarsen(da,comm,&dacoarse);
170:   MyDMShellCreate(PetscObjectComm((PetscObject)shell),dacoarse,dmnew);
171:   /* discard an "extra" reference count to dacoarse */
172:   DMDestroy(&dacoarse);
173:   return 0;
174: }

178: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
179: {
181:   PetscInt       mx,idx[2];
182:   PetscScalar    h,v[2];
183:   DM             da,shell;

186:   KSPGetDM(ksp,&shell);
187:   DMShellGetContext(shell,(void**)&da);
188:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
189:   h      = 1.0/((mx-1));
190:   VecSet(b,h);
191:   idx[0] = 0; idx[1] = mx -1;
192:   v[0]   = v[1] = 0.0;
193:   VecSetValues(b,2,idx,v,INSERT_VALUES);
194:   VecAssemblyBegin(b);
195:   VecAssemblyEnd(b);
196:   return(0);
197: }

201: static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
202: {
204:   PetscInt       i,mx,xm,xs;
205:   PetscScalar    v[3],h;
206:   MatStencil     row,col[3];
207:   DM             da,shell;

210:   KSPGetDM(ksp,&shell);
211:   DMShellGetContext(shell,(void**)&da);
212:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
213:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
214:   h    = 1.0/(mx-1);

216:   for (i=xs; i<xs+xm; i++) {
217:     row.i = i;
218:     if (i==0 || i==mx-1) {
219:       v[0] = 2.0/h;
220:       MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
221:     } else {
222:       v[0]  = (-1.0)/h;col[0].i = i-1;
223:       v[1]  = (2.0)/h;col[1].i = row.i;
224:       v[2]  = (-1.0)/h;col[2].i = i+1;
225:       MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
226:     }
227:   }
228:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
229:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
230:   return(0);
231: }