Actual source code: ex54.c

petsc-master 2019-11-16
Report Typos and Errors

  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 coeficient in embedded circle\n\n";

  6:  #include <petscksp.h>

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

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

 39:   /* create stiffness matrix (2) */
 40:   MatCreate(comm,&Amat);
 41:   MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,M,M);
 42:   MatSetType(Amat,MATMPIAIJ);
 43:   MatSetFromOptions(Amat);
 44:   MatSeqAIJSetPreallocation(Amat,81,NULL);
 45:   MatMPIAIJSetPreallocation(Amat,81,NULL,57,NULL);

 47:   MatCreate(comm,&Pmat);
 48:   MatSetSizes(Pmat,PETSC_DECIDE,PETSC_DECIDE,M,M);
 49:   MatSetType(Pmat,MATMPIAIJ);
 50:   MatSetFromOptions(Pmat);
 51:   MatSeqAIJSetPreallocation(Pmat,81,NULL);
 52:   MatMPIAIJSetPreallocation(Pmat,81,NULL,57,NULL);

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

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

132:     /* Setup solver */
133:     KSPCreate(PETSC_COMM_WORLD,&ksp);
134:     KSPSetFromOptions(ksp);

136:     /* finish KSP/PC setup */
137:     KSPSetOperators(ksp, Amat, Amat);

139:     KSPGetPC(ksp,&pc);
140:     PCSetCoordinates(pc, 2, m, coords);
141:     PetscFree(coords);
142:   }

144:   if (!PETSC_TRUE) {
145:     PetscViewer viewer;
146:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
147:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
148:     MatView(Amat,viewer);
149:     PetscViewerPopFormat(viewer);
150:     PetscViewerDestroy(&viewer);
151:   }

153:   /* solve */
154: #if defined(PETSC_USE_LOG)
155:   PetscLogStageRegister("Solve", &stage);
156:   PetscLogStagePush(stage);
157: #endif
158:   VecSet(xx,.0);

160:   KSPSetUp(ksp);

162:   KSPSolve(ksp,bb,xx);

164: #if defined(PETSC_USE_LOG)
165:   PetscLogStagePop();
166: #endif

168:   KSPGetIterationNumber(ksp,&its);

170:   if (!PETSC_TRUE) {
171:     PetscReal   norm,norm2;
172:     PetscViewer viewer;
173:     Vec         res;
174:     PetscViewerASCIIOpen(comm, "rhs.m", &viewer);
175:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
176:     VecView(bb,viewer);
177:     PetscViewerPopFormat(viewer);
178:     PetscViewerDestroy(&viewer);
179:     VecNorm(bb, NORM_2, &norm2);

181:     PetscViewerASCIIOpen(comm, "solution.m", &viewer);
182:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
183:     VecView(xx,viewer);
184:     PetscViewerPopFormat(viewer);
185:     PetscViewerDestroy(&viewer);

187:     VecDuplicate(xx, &res);
188:     MatMult(Amat, xx, res);
189:     VecAXPY(bb, -1.0, res);
190:     VecDestroy(&res);
191:     VecNorm(bb,NORM_2,&norm);
192:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,PETSC_FUNCTION_NAME,norm/norm2,norm2);

194:     PetscViewerASCIIOpen(comm, "residual.m", &viewer);
195:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
196:     VecView(bb,viewer);
197:     PetscViewerPopFormat(viewer);
198:     PetscViewerDestroy(&viewer);
199:   }

201:   /* Free work space */
202:   KSPDestroy(&ksp);
203:   VecDestroy(&xx);
204:   VecDestroy(&bb);
205:   MatDestroy(&Amat);
206:   MatDestroy(&Pmat);

208:   PetscFinalize();
209:   return ierr;
210: }



214: /*TEST

216:    build:
217:       requires: !complex

219:    test:
220:       nsize: 4
221:       args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -ksp_converged_reason -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -mg_levels_pc_type sor -pc_gamg_use_sa_esteig

223:    test:
224:       suffix: seqaijmkl
225:       nsize: 4
226:       requires: mkl_sparse
227:       args: -ne 19 -alpha 1.e-3 -pc_type gamg -pc_gamg_agg_nsmooths 1  -mg_levels_ksp_max_it 3 -ksp_monitor -ksp_converged_reason -ksp_type cg -mat_seqaij_type seqaijmkl

229:    test:
230:       suffix: Classical
231:       args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type classical -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -ksp_converged_reason -mg_levels_esteig_ksp_type cg
232:       output_file: output/ex54_classical.out

234:    test:
235:       suffix: geo
236:       nsize: 4
237:       args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type geo -pc_gamg_coarse_eq_limit 200 -mg_levels_pc_type jacobi -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -ksp_monitor_short -mg_levels_esteig_ksp_type cg  -mg_levels_ksp_max_it 3
238:       requires: triangle
239:       output_file: output/ex54_0.out

241: TEST*/