Actual source code: ex32.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: /*
  2:   Laplacian in 3D. Use for testing BAIJ matrix.
  3:   Modeled by the partial differential equation

  5:    - Laplacian u = 1,0 < x,y,z < 1,

  7:    with boundary conditions
  8:    u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
  9: */

 11: static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";

 13:  #include <petscdm.h>
 14:  #include <petscdmda.h>
 15:  #include <petscksp.h>

 17: extern PetscErrorCode ComputeMatrix(DM,Mat);
 18: extern PetscErrorCode ComputeRHS(DM,Vec);

 20: int main(int argc,char **argv)
 21: {
 23:   KSP            ksp;
 24:   PC             pc;
 25:   Vec            x,b;
 26:   DM             da;
 27:   Mat            A,Atrans;
 28:   PetscInt       dof=1,M=8;
 29:   PetscBool      flg,trans=PETSC_FALSE;

 31:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 32:   PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
 33:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 34:   PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);

 36:   DMDACreate(PETSC_COMM_WORLD,&da);
 37:   DMSetDimension(da,3);
 38:   DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
 39:   DMDASetStencilType(da,DMDA_STENCIL_STAR);
 40:   DMDASetSizes(da,M,M,M);
 41:   DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 42:   DMDASetDof(da,dof);
 43:   DMDASetStencilWidth(da,1);
 44:   DMDASetOwnershipRanges(da,NULL,NULL,NULL);
 45:   DMSetFromOptions(da);
 46:   DMSetUp(da);

 48:   DMCreateGlobalVector(da,&x);
 49:   DMCreateGlobalVector(da,&b);
 50:   ComputeRHS(da,b);
 51:   DMSetMatType(da,MATBAIJ);
 52:   DMSetFromOptions(da);
 53:   DMCreateMatrix(da,&A);
 54:   ComputeMatrix(da,A);

 56:   /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
 57:   MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);
 58:   MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);
 59:   MatScale(A,0.5);
 60:   MatDestroy(&Atrans);

 62:   /* Test sbaij matrix */
 63:   flg  = PETSC_FALSE;
 64:   PetscOptionsGetBool(NULL,NULL, "-test_sbaij1", &flg,NULL);
 65:   if (flg) {
 66:     Mat       sA;
 67:     PetscBool issymm;
 68:     MatIsTranspose(A,A,0.0,&issymm);
 69:     if (issymm) {
 70:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 71:     } else {PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric\n");}
 72:     MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
 73:     MatDestroy(&A);
 74:     A    = sA;
 75:   }

 77:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 78:   KSPSetFromOptions(ksp);
 79:   KSPSetOperators(ksp,A,A);
 80:   KSPGetPC(ksp,&pc);
 81:   PCSetDM(pc,(DM)da);

 83:   if (trans) {
 84:     KSPSolveTranspose(ksp,b,x);
 85:   } else {
 86:     KSPSolve(ksp,b,x);
 87:   }

 89:   /* check final residual */
 90:   flg  = PETSC_FALSE;
 91:   PetscOptionsGetBool(NULL,NULL, "-check_final_residual", &flg,NULL);
 92:   if (flg) {
 93:     Vec       b1;
 94:     PetscReal norm;
 95:     KSPGetSolution(ksp,&x);
 96:     VecDuplicate(b,&b1);
 97:     MatMult(A,x,b1);
 98:     VecAXPY(b1,-1.0,b);
 99:     VecNorm(b1,NORM_2,&norm);
100:     PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
101:     VecDestroy(&b1);
102:   }

104:   KSPDestroy(&ksp);
105:   VecDestroy(&x);
106:   VecDestroy(&b);
107:   MatDestroy(&A);
108:   DMDestroy(&da);
109:   PetscFinalize();
110:   return ierr;
111: }

113: PetscErrorCode ComputeRHS(DM da,Vec b)
114: {
116:   PetscInt       mx,my,mz;
117:   PetscScalar    h;

120:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
121:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
122:   VecSet(b,h);
123:   return(0);
124: }

126: PetscErrorCode ComputeMatrix(DM da,Mat B)
127: {
129:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
130:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
131:   MatStencil     row,col;

134:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
135:   /* For simplicity, this example only works on mx=my=mz */
136:   if (mx != my || mx != mz) SETERRQ3(PETSC_COMM_SELF,1,"This example only works with mx %D = my %D = mz %D\n",mx,my,mz);

138:   Hx      = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
139:   HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;

141:   PetscMalloc1(2*dof*dof+1,&v);
142:   v_neighbor = v + dof*dof;
143:   PetscArrayzero(v,2*dof*dof+1);
144:   k3         = 0;
145:   for (k1=0; k1<dof; k1++) {
146:     for (k2=0; k2<dof; k2++) {
147:       if (k1 == k2) {
148:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
149:         v_neighbor[k3] = -HxHydHz;
150:       } else {
151:         v[k3]          = k1/(dof*dof);
152:         v_neighbor[k3] = k2/(dof*dof);
153:       }
154:       k3++;
155:     }
156:   }
157:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

159:   for (k=zs; k<zs+zm; k++) {
160:     for (j=ys; j<ys+ym; j++) {
161:       for (i=xs; i<xs+xm; i++) {
162:         row.i = i; row.j = j; row.k = k;
163:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boudary points */
164:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
165:         } else { /* interior points */
166:           /* center */
167:           col.i = i; col.j = j; col.k = k;
168:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);

170:           /* x neighbors */
171:           col.i = i-1; col.j = j; col.k = k;
172:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
173:           col.i = i+1; col.j = j; col.k = k;
174:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

176:           /* y neighbors */
177:           col.i = i; col.j = j-1; col.k = k;
178:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
179:           col.i = i; col.j = j+1; col.k = k;
180:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

182:           /* z neighbors */
183:           col.i = i; col.j = j; col.k = k-1;
184:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
185:           col.i = i; col.j = j; col.k = k+1;
186:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
187:         }
188:       }
189:     }
190:   }
191:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
192:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
193:   PetscFree(v);
194:   return(0);
195: }



199: /*TEST

201:    test:
202:       args: -ksp_monitor_short -dm_mat_type sbaij -ksp_monitor_short -pc_type cholesky -ksp_view

204: TEST*/