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*/