Actual source code: ex32.c

  1: /*T
  2:    Concepts: KSP^solving a system of linear equations
  3:    Concepts: KSP^Laplacian, 2d
  4:    Processors: n
  5: T*/

  7: /*
  8: Laplacian in 2D. Modeled by the partial differential equation

 10:    div  grad u = f,  0 < x,y < 1,

 12: with forcing function

 14:    f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}

 16: with pure Neumann boundary conditions

 18: The functions are cell-centered

 20: This uses multigrid to solve the linear system

 22:        Contributed by Andrei Draganescu <aidraga@sandia.gov>

 24: Note the nice multigrid convergence despite the fact it is only using
 25: peicewise constant interpolation/restriction. This is because cell-centered multigrid
 26: does not need the same rule: 

 28:     polynomial degree(interpolation) + polynomial degree(restriction) + 2 > degree of PDE

 30: that vertex based multigrid needs.
 31: */

 33: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

 35:  #include petscda.h
 36:  #include petscksp.h
 37:  #include petscmg.h
 38:  #include petscdmmg.h


 43: typedef enum {DIRICHLET, NEUMANN} BCType;

 45: typedef struct {
 46:   PetscScalar   nu;
 47:   BCType        bcType;
 48: } UserContext;

 52: int main(int argc,char **argv)
 53: {
 54:   DMMG           *dmmg;
 55:   DA             da;
 56:   UserContext    user;
 57:   PetscReal      norm;
 58:   const char     *bcTypes[2] = {"dirichlet","neumann"};
 60:   PetscInt       l,bc;

 62:   PetscInitialize(&argc,&argv,(char *)0,help);
 63: 
 64:   DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
 65:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 66:   DASetInterpolationType(da, DA_Q0);
 67: 
 68:   DMMGSetDM(dmmg,(DM)da);
 69:   DADestroy(da);
 70:   for (l = 0; l < DMMGGetLevels(dmmg); l++) {
 71:     DMMGSetUser(dmmg,l,&user);
 72:   }
 73: 
 74:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
 75:   user.nu     = 0.1;
 76:   PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", 0.1, &user.nu, PETSC_NULL);
 77:   bc          = (PetscInt)NEUMANN;
 78:   PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
 79:   user.bcType = (BCType)bc;
 80:   PetscOptionsEnd();
 81: 
 82:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
 83:   if (user.bcType == NEUMANN) {
 84:     DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
 85:   }
 86: 
 87:   DMMGSolve(dmmg);
 88: 
 89:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
 90:   VecAXPY(DMMGGetRHS(dmmg),-1.0,DMMGGetr(dmmg));
 91:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
 92:   /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",norm); */

 94:   DMMGDestroy(dmmg);
 95:   PetscFinalize();
 96:   return 0;
 97: }

101: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
102: {
103:   DA             da = (DA)dmmg->dm;
104:   UserContext    *user = (UserContext *) dmmg->user;
106:   PetscInt       i,j,mx,my,xm,ym,xs,ys;
107:   PetscScalar    Hx,Hy;
108:   PetscScalar    **array;

111:   DAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0);
112:   Hx   = 1.0 / (PetscReal)(mx);
113:   Hy   = 1.0 / (PetscReal)(my);
114:   DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
115:   DAVecGetArray(da, b, &array);
116:   for (j=ys; j<ys+ym; j++){
117:     for(i=xs; i<xs+xm; i++){
118:       array[j][i] = PetscExpScalar(-(((PetscReal)i+0.5)*Hx)*(((PetscReal)i+0.5)*Hx)/user->nu)*PetscExpScalar(-(((PetscReal)j+0.5)*Hy)*(((PetscReal)j+0.5)*Hy)/user->nu)*Hx*Hy;
119:     }
120:   }
121:   DAVecRestoreArray(da, b, &array);
122:   VecAssemblyBegin(b);
123:   VecAssemblyEnd(b);

125:   /* force right hand side to be consistent for singular matrix */
126:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
127:   if (user->bcType == NEUMANN)
128:     {
129:       MatNullSpace nullspace;
130: 
131:       KSPGetNullSpace(dmmg->ksp,&nullspace);
132:       MatNullSpaceRemove(nullspace,b,PETSC_NULL);
133:     }
134:   return(0);
135: }

137: 
140: PetscErrorCode ComputeJacobian(DMMG dmmg, Mat J,Mat jac)
141: {
142:   DA             da = (DA) dmmg->dm;
143:   UserContext    *user = (UserContext *) dmmg->user;
145:   PetscInt       i,j,mx,my,xm,ym,xs,ys,num, numi, numj;
146:   PetscScalar    v[5],Hx,Hy,HydHx,HxdHy;
147:   MatStencil     row, col[5];

150:   DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
151:   Hx    = 1.0 / (PetscReal)(mx);
152:   Hy    = 1.0 / (PetscReal)(my);
153:   HxdHy = Hx/Hy;
154:   HydHx = Hy/Hx;
155:   DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
156:   for (j=ys; j<ys+ym; j++)
157:     {
158:       for(i=xs; i<xs+xm; i++)
159:         {
160:           row.i = i; row.j = j;
161:           if (i==0 || j==0 || i==mx-1 || j==my-1)
162:             {
163:               if (user->bcType == DIRICHLET)
164:                 {
165:                   SETERRQ(PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
166:                   v[0] = 2.0*(HxdHy + HydHx);
167:                   MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
168: 
169:                 }
170:               else if (user->bcType == NEUMANN)
171:                 {
172:                   num = 0; numi=0; numj=0;
173:                   if (j!=0)
174:                     {
175:                       v[num] = -HxdHy;
176:                       col[num].i = i;
177:                       col[num].j = j-1;
178:                       num++; numj++;
179:                     }
180:                   if (i!=0)
181:                     {
182:                       v[num] = -HydHx;
183:                       col[num].i = i-1;
184:                       col[num].j = j;
185:                       num++; numi++;
186:                     }
187:                   if (i!=mx-1)
188:                     {
189:                       v[num] = -HydHx;
190:                       col[num].i = i+1;
191:                       col[num].j = j;
192:                       num++; numi++;
193:                     }
194:                   if (j!=my-1)
195:                     {
196:                       v[num] = -HxdHy;
197:                       col[num].i = i;
198:                       col[num].j = j+1;
199:                       num++; numj++;
200:                     }
201:                   v[num]   = (PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx; col[num].i = i;   col[num].j = j;
202:                   num++;
203:                   MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
204:                 }
205:             }
206:           else
207:             {
208:               v[0] = -HxdHy;              col[0].i = i;   col[0].j = j-1;
209:               v[1] = -HydHx;              col[1].i = i-1; col[1].j = j;
210:               v[2] = 2.0*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
211:               v[3] = -HydHx;              col[3].i = i+1; col[3].j = j;
212:               v[4] = -HxdHy;              col[4].i = i;   col[4].j = j+1;
213:               MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
214:             }
215:         }
216:     }
217:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
218:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
219:   return(0);
220: }