Actual source code: ex54.c


  2: static char help[] = "Creates a matrix from quadrilateral finite elements in 2D, Laplacian \n\
  3:   -ne <size>       : problem size in number of elements (eg, -ne 31 gives 32^2 grid)\n\
  4:   -alpha <v>      : scaling of material coefficient in embedded circle\n\n";

  6: #include <petscksp.h>

  8: int main(int argc, char **args)
  9: {
 10:   Mat         Amat, Pmat;
 11:   PetscInt    i, m, M, its, Istart, Iend, j, Ii, ix, ne = 4;
 12:   PetscReal   x, y, h;
 13:   Vec         xx, bb;
 14:   KSP         ksp;
 15:   PetscReal   soft_alpha = 1.e-3;
 16:   MPI_Comm    comm;
 17:   PetscMPIInt npe, mype;
 18:   PetscScalar DD[4][4], DD2[4][4];
 19: #if defined(PETSC_USE_LOG)
 20:   PetscLogStage stage;
 21: #endif
 22: #define DIAG_S 0.0
 23:   PetscScalar DD1[4][4] = {
 24:     {5.0 + DIAG_S, -2.0,         -1.0,         -2.0        },
 25:     {-2.0,         5.0 + DIAG_S, -2.0,         -1.0        },
 26:     {-1.0,         -2.0,         5.0 + DIAG_S, -2.0        },
 27:     {-2.0,         -1.0,         -2.0,         5.0 + DIAG_S}
 28:   };

 31:   PetscInitialize(&argc, &args, (char *)0, help);
 32:   comm = PETSC_COMM_WORLD;
 33:   MPI_Comm_rank(comm, &mype);
 34:   MPI_Comm_size(comm, &npe);
 35:   PetscOptionsGetInt(NULL, NULL, "-ne", &ne, NULL);
 36:   h = 1. / ne;
 37:   /* ne*ne; number of global elements */
 38:   PetscOptionsGetReal(NULL, NULL, "-alpha", &soft_alpha, NULL);
 39:   M = (ne + 1) * (ne + 1); /* global number of nodes */

 41:   /* create stiffness matrix (2) */
 42:   MatCreate(comm, &Amat);
 43:   MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, M, M);
 44:   MatSetType(Amat, MATAIJ);
 45:   MatSetOption(Amat, MAT_SPD, PETSC_TRUE);
 46:   MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE);
 47:   MatSetFromOptions(Amat);
 48:   MatSeqAIJSetPreallocation(Amat, 81, NULL);
 49:   MatMPIAIJSetPreallocation(Amat, 81, NULL, 57, NULL);

 51:   MatCreate(comm, &Pmat);
 52:   MatSetSizes(Pmat, PETSC_DECIDE, PETSC_DECIDE, M, M);
 53:   MatSetType(Pmat, MATMPIAIJ);
 54:   MatSetFromOptions(Pmat);
 55:   MatSeqAIJSetPreallocation(Pmat, 81, NULL);
 56:   MatMPIAIJSetPreallocation(Pmat, 81, NULL, 57, NULL);

 58:   /* vectors */
 59:   MatCreateVecs(Amat, &bb, &xx);
 60:   VecSet(bb, .0);
 61:   /* generate element matrices -- see ex56.c on how to use different data set */
 62:   {
 63:     DD1[0][0] = 0.66666666666666663;
 64:     DD1[0][1] = -0.16666666666666669;
 65:     DD1[0][2] = -0.33333333333333343;
 66:     DD1[0][3] = -0.16666666666666666;
 67:     DD1[1][0] = -0.16666666666666669;
 68:     DD1[1][1] = 0.66666666666666663;
 69:     DD1[1][2] = -0.16666666666666666;
 70:     DD1[1][3] = -0.33333333333333343;
 71:     DD1[2][0] = -0.33333333333333343;
 72:     DD1[2][1] = -0.16666666666666666;
 73:     DD1[2][2] = 0.66666666666666663;
 74:     DD1[2][3] = -0.16666666666666663;
 75:     DD1[3][0] = -0.16666666666666666;
 76:     DD1[3][1] = -0.33333333333333343;
 77:     DD1[3][2] = -0.16666666666666663;
 78:     DD1[3][3] = 0.66666666666666663;

 80:     /* BC version of element */
 81:     for (i = 0; i < 4; i++) {
 82:       for (j = 0; j < 4; j++) {
 83:         if (i < 2 || j < 2) {
 84:           if (i == j) DD2[i][j] = .1 * DD1[i][j];
 85:           else DD2[i][j] = 0.0;
 86:         } else DD2[i][j] = DD1[i][j];
 87:       }
 88:     }
 89:   }
 90:   {
 91:     PetscReal *coords;
 92:     PC         pc;
 93:     /* forms the element stiffness for the Laplacian and coordinates */
 94:     MatGetOwnershipRange(Amat, &Istart, &Iend);
 95:     m = Iend - Istart;
 96:     PetscMalloc1(2 * m, &coords);
 97:     for (Ii = Istart, ix = 0; Ii < Iend; Ii++, ix++) {
 98:       j = Ii / (ne + 1);
 99:       i = Ii % (ne + 1);
100:       /* coords */
101:       x                  = h * (Ii % (ne + 1));
102:       y                  = h * (Ii / (ne + 1));
103:       coords[2 * ix]     = x;
104:       coords[2 * ix + 1] = y;
105:       if (i < ne && j < ne) {
106:         PetscInt jj, ii, idx[4];
107:         /* radius */
108:         PetscReal radius = PetscSqrtReal((x - .5 + h / 2) * (x - .5 + h / 2) + (y - .5 + h / 2) * (y - .5 + h / 2));
109:         PetscReal alpha  = 1.0;
110:         idx[0]           = Ii;
111:         idx[1]           = Ii + 1;
112:         idx[2]           = Ii + (ne + 1) + 1;
113:         idx[3]           = Ii + (ne + 1);
114:         if (radius < 0.25) alpha = soft_alpha;
115:         for (ii = 0; ii < 4; ii++) {
116:           for (jj = 0; jj < 4; jj++) DD[ii][jj] = alpha * DD1[ii][jj];
117:         }
118:         MatSetValues(Pmat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES);
119:         if (j > 0) {
120:           MatSetValues(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES);
121:         } else {
122:           /* a BC */
123:           for (ii = 0; ii < 4; ii++) {
124:             for (jj = 0; jj < 4; jj++) DD[ii][jj] = alpha * DD2[ii][jj];
125:           }
126:           MatSetValues(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES);
127:         }
128:       }
129:       if (j > 0) {
130:         PetscScalar v  = h * h;
131:         PetscInt    jj = Ii;
132:         VecSetValues(bb, 1, &jj, &v, INSERT_VALUES);
133:       }
134:     }
135:     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
136:     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
137:     MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY);
138:     MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY);
139:     VecAssemblyBegin(bb);
140:     VecAssemblyEnd(bb);

142:     /* Setup solver */
143:     KSPCreate(PETSC_COMM_WORLD, &ksp);
144:     KSPSetFromOptions(ksp);

146:     /* finish KSP/PC setup */
147:     KSPSetOperators(ksp, Amat, Amat);

149:     KSPGetPC(ksp, &pc);
150:     PCSetCoordinates(pc, 2, m, coords);
151:     PetscFree(coords);
152:   }

154:   if (!PETSC_TRUE) {
155:     PetscViewer viewer;
156:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
157:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
158:     MatView(Amat, viewer);
159:     PetscViewerPopFormat(viewer);
160:     PetscViewerDestroy(&viewer);
161:   }

163:   /* solve */
164: #if defined(PETSC_USE_LOG)
165:   PetscLogStageRegister("Solve", &stage);
166:   PetscLogStagePush(stage);
167: #endif
168:   VecSet(xx, .0);

170:   KSPSetUp(ksp);

172:   KSPSolve(ksp, bb, xx);

174: #if defined(PETSC_USE_LOG)
175:   PetscLogStagePop();
176: #endif

178:   KSPGetIterationNumber(ksp, &its);

180:   if (!PETSC_TRUE) {
181:     PetscReal   norm, norm2;
182:     PetscViewer viewer;
183:     Vec         res;
184:     PetscViewerASCIIOpen(comm, "rhs.m", &viewer);
185:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
186:     VecView(bb, viewer);
187:     PetscViewerPopFormat(viewer);
188:     PetscViewerDestroy(&viewer);
189:     VecNorm(bb, NORM_2, &norm2);

191:     PetscViewerASCIIOpen(comm, "solution.m", &viewer);
192:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
193:     VecView(xx, viewer);
194:     PetscViewerPopFormat(viewer);
195:     PetscViewerDestroy(&viewer);

197:     VecDuplicate(xx, &res);
198:     MatMult(Amat, xx, res);
199:     VecAXPY(bb, -1.0, res);
200:     VecDestroy(&res);
201:     VecNorm(bb, NORM_2, &norm);
202:     PetscPrintf(PETSC_COMM_WORLD, "[%d]%s |b-Ax|/|b|=%e, |b|=%e\n", 0, PETSC_FUNCTION_NAME, (double)(norm / norm2), (double)norm2);

204:     PetscViewerASCIIOpen(comm, "residual.m", &viewer);
205:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
206:     VecView(bb, viewer);
207:     PetscViewerPopFormat(viewer);
208:     PetscViewerDestroy(&viewer);
209:   }

211:   /* Free work space */
212:   KSPDestroy(&ksp);
213:   VecDestroy(&xx);
214:   VecDestroy(&bb);
215:   MatDestroy(&Amat);
216:   MatDestroy(&Pmat);

218:   PetscFinalize();
219:   return 0;
220: }

222: /*TEST

224:    build:
225:       requires: !complex

227:    test:
228:       nsize: 4
229:       args: -ne 19 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -ksp_monitor -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 5 -pc_gamg_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -pc_gamg_aggressive_coarsening 0

231:    test:
232:       suffix: seqaijmkl
233:       nsize: 4
234:       requires: mkl_sparse
235:       args: -ne 19 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -ksp_monitor -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 5 -pc_gamg_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -mat_seqaij_type seqaijmkl -pc_gamg_aggressive_coarsening 0

237:    test:
238:       suffix: Classical
239:       args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -pc_gamg_type classical -ksp_monitor -ksp_converged_reason -mg_levels_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -mat_coarsen_type mis
240:       output_file: output/ex54_classical.out

242:    test:
243:       suffix: geo
244:       nsize: 4
245:       args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 4 -pc_gamg_type geo -pc_gamg_coarse_eq_limit 200 -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.05 -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1e-3 -ksp_norm_type unpreconditioned
246:       requires: triangle
247:       output_file: output/ex54_0.out

249: TEST*/