Actual source code: ex65.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  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,(void**)&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,(void**)&da1);
 98:   DMShellGetContext(dm2,(void**)&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,(void**)&da1);
110:   DMShellGetContext(dm2,(void**)&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,(void**)&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,(void**)&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,(void**)&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,(void**)&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,(void**)&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,(void**)&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: }


217: /*TEST

219:    test:
220:       nsize: 4
221:       args: -ksp_monitor -pc_type mg -da_refine 3

223: TEST*/