Actual source code: ex32.c
petsc-3.12.0 2019-09-29
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*/