1 static char help[] = "3D, tri-quadratic hexahedra (Q1), displacement finite element formulation\n\ 2 of linear elasticity. E=1.0, nu=1/3.\n\ 3 Unit cube domain with Dirichlet boundary\n\n"; 4 5 #include <petscdmplex.h> 6 #include <petscsnes.h> 7 #include <petscds.h> 8 #include <petscdmforest.h> 9 10 static PetscReal s_soft_alpha=1.e-3; 11 static PetscReal s_mu=0.4; 12 static PetscReal s_lambda=0.4; 13 14 static void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 15 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 16 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 17 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 18 { 19 f0[0] = 1; /* x direction pull */ 20 f0[1] = -x[2]; /* add a twist around x-axis */ 21 f0[2] = x[1]; 22 } 23 24 static void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 25 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 26 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 27 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 28 { 29 const PetscInt Ncomp = dim; 30 PetscInt comp, d; 31 for (comp = 0; comp < Ncomp; ++comp) { 32 for (d = 0; d < dim; ++d) { 33 f1[comp*dim+d] = 0.0; 34 } 35 } 36 } 37 38 /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */ 39 static void f1_u_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux, 40 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 41 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 42 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 43 { 44 PetscReal trace,mu=s_mu,lambda=s_lambda,rad; 45 PetscInt i,j; 46 for (i=0,rad=0.;i<dim;i++) { 47 PetscReal t=x[i]; 48 rad += t*t; 49 } 50 rad = PetscSqrtReal(rad); 51 if (rad>0.25) { 52 mu *= s_soft_alpha; 53 lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */ 54 } 55 for (i=0,trace=0; i < dim; ++i) { 56 trace += PetscRealPart(u_x[i*dim+i]); 57 } 58 for (i=0; i < dim; ++i) { 59 for (j=0; j < dim; ++j) { 60 f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]); 61 } 62 f1[i*dim+i] += lambda*trace; 63 } 64 } 65 66 /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */ 67 static void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 68 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 69 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 70 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 71 { 72 PetscReal trace,mu=s_mu,lambda=s_lambda; 73 PetscInt i,j; 74 for (i=0,trace=0; i < dim; ++i) { 75 trace += PetscRealPart(u_x[i*dim+i]); 76 } 77 for (i=0; i < dim; ++i) { 78 for (j=0; j < dim; ++j) { 79 f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]); 80 } 81 f1[i*dim+i] += lambda*trace; 82 } 83 } 84 85 /* 3D elasticity */ 86 #define IDX(ii,jj,kk,ll) (27*ii+9*jj+3*kk+ll) 87 88 void g3_uu_3d_private( PetscScalar g3[], const PetscReal mu, const PetscReal lambda) 89 { 90 if (1) { 91 g3[0] += lambda; 92 g3[0] += mu; 93 g3[0] += mu; 94 g3[4] += lambda; 95 g3[8] += lambda; 96 g3[10] += mu; 97 g3[12] += mu; 98 g3[20] += mu; 99 g3[24] += mu; 100 g3[28] += mu; 101 g3[30] += mu; 102 g3[36] += lambda; 103 g3[40] += lambda; 104 g3[40] += mu; 105 g3[40] += mu; 106 g3[44] += lambda; 107 g3[50] += mu; 108 g3[52] += mu; 109 g3[56] += mu; 110 g3[60] += mu; 111 g3[68] += mu; 112 g3[70] += mu; 113 g3[72] += lambda; 114 g3[76] += lambda; 115 g3[80] += lambda; 116 g3[80] += mu; 117 g3[80] += mu; 118 } else { 119 int i,j,k,l; 120 static int cc=-1; 121 cc++; 122 for (i = 0; i < 3; ++i) { 123 for (j = 0; j < 3; ++j) { 124 for (k = 0; k < 3; ++k) { 125 for (l = 0; l < 3; ++l) { 126 if (k==l && i==j) g3[IDX(i,j,k,l)] += lambda; 127 if (i==k && j==l) g3[IDX(i,j,k,l)] += mu; 128 if (i==l && j==k) g3[IDX(i,j,k,l)] += mu; 129 if (k==l && i==j && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += lambda;\n",IDX(i,j,k,l)); 130 if (i==k && j==l && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l)); 131 if (i==l && j==k && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l)); 132 } 133 } 134 } 135 } 136 } 137 } 138 139 static void g3_uu_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux, 140 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 141 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 142 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 143 { 144 PetscReal mu=s_mu, lambda=s_lambda,rad; 145 PetscInt i; 146 for (i=0,rad=0.;i<dim;i++) { 147 PetscReal t=x[i]; 148 rad += t*t; 149 } 150 rad = PetscSqrtReal(rad); 151 if (rad>0.25) { 152 mu *= s_soft_alpha; 153 lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */ 154 } 155 g3_uu_3d_private(g3,mu,lambda); 156 } 157 158 static void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 159 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 160 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 161 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 162 { 163 g3_uu_3d_private(g3,s_mu,s_lambda); 164 } 165 166 static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 167 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 168 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 169 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 170 { 171 const PetscInt Ncomp = dim; 172 PetscInt comp; 173 174 for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0; 175 } 176 177 /* PI_i (x_i^4 - x_i^2) */ 178 static void f0_u_x4(PetscInt dim, PetscInt Nf, PetscInt NfAux, 179 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 180 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 181 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 182 { 183 const PetscInt Ncomp = dim; 184 PetscInt comp,i; 185 186 for (comp = 0; comp < Ncomp; ++comp) { 187 f0[comp] = 1e5; 188 for (i = 0; i < Ncomp; ++i) { 189 f0[comp] *= /* (comp+1)* */(x[i]*x[i]*x[i]*x[i] - x[i]*x[i]); /* assumes (0,1]^D domain */ 190 } 191 } 192 } 193 194 PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 195 { 196 const PetscInt Ncomp = dim; 197 PetscInt comp; 198 199 for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0; 200 return 0; 201 } 202 203 int main(int argc,char **args) 204 { 205 Mat Amat; 206 PetscErrorCode ierr; 207 SNES snes; 208 KSP ksp; 209 MPI_Comm comm; 210 PetscMPIInt rank; 211 #if defined(PETSC_USE_LOG) 212 PetscLogStage stage[17]; 213 #endif 214 PetscBool test_nonzero_cols = PETSC_FALSE,use_nearnullspace = PETSC_TRUE,attach_nearnullspace = PETSC_FALSE; 215 Vec xx,bb; 216 PetscInt iter,i,N,dim = 3,cells[3] = {1,1,1},max_conv_its,local_sizes[7],run_type = 1; 217 DM dm,distdm,basedm; 218 PetscBool flg; 219 char convType[256]; 220 PetscReal Lx,mdisp[10],err[10]; 221 const char * const options[10] = {"-ex56_dm_refine 0", 222 "-ex56_dm_refine 1", 223 "-ex56_dm_refine 2", 224 "-ex56_dm_refine 3", 225 "-ex56_dm_refine 4", 226 "-ex56_dm_refine 5", 227 "-ex56_dm_refine 6", 228 "-ex56_dm_refine 7", 229 "-ex56_dm_refine 8", 230 "-ex56_dm_refine 9"}; 231 PetscFunctionBeginUser; 232 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 233 comm = PETSC_COMM_WORLD; 234 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 235 /* options */ 236 ierr = PetscOptionsBegin(comm,NULL,"3D bilinear Q1 elasticity options","");CHKERRQ(ierr); 237 { 238 i = 3; 239 ierr = PetscOptionsIntArray("-cells", "Number of (flux tube) processor in each dimension", "ex56.c", cells, &i, NULL);CHKERRQ(ierr); 240 241 Lx = 1.; /* or ne for rod */ 242 max_conv_its = 3; 243 ierr = PetscOptionsInt("-max_conv_its","Number of iterations in convergence study","",max_conv_its,&max_conv_its,NULL);CHKERRQ(ierr); 244 if (max_conv_its<=0 || max_conv_its>7) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_USER, "Bad number of iterations for convergence test (%D)",max_conv_its); 245 ierr = PetscOptionsReal("-lx","Length of domain","",Lx,&Lx,NULL);CHKERRQ(ierr); 246 ierr = PetscOptionsReal("-alpha","material coefficient inside circle","",s_soft_alpha,&s_soft_alpha,NULL);CHKERRQ(ierr); 247 ierr = PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);CHKERRQ(ierr); 248 ierr = PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);CHKERRQ(ierr); 249 ierr = PetscOptionsBool("-attach_mat_nearnullspace","MatNearNullSpace API test (via MatSetNearNullSpace)","",attach_nearnullspace,&attach_nearnullspace,NULL);CHKERRQ(ierr); 250 ierr = PetscOptionsInt("-run_type","0: twisting load on cantalever, 1: 3rd order accurate convergence test","",run_type,&run_type,NULL);CHKERRQ(ierr); 251 } 252 ierr = PetscOptionsEnd();CHKERRQ(ierr); 253 ierr = PetscLogStageRegister("Mesh Setup", &stage[16]);CHKERRQ(ierr); 254 for (iter=0 ; iter<max_conv_its ; iter++) { 255 char str[] = "Solve 0"; 256 str[6] += iter; 257 ierr = PetscLogStageRegister(str, &stage[iter]);CHKERRQ(ierr); 258 } 259 /* create DM, Plex calls DMSetup */ 260 ierr = PetscLogStagePush(stage[16]);CHKERRQ(ierr); 261 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr); 262 { 263 DMLabel label; 264 IS is; 265 ierr = DMCreateLabel(dm, "boundary");CHKERRQ(ierr); 266 ierr = DMGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 267 ierr = DMPlexMarkBoundaryFaces(dm, 1, label);CHKERRQ(ierr); 268 if (!run_type) { 269 ierr = DMGetStratumIS(dm, "boundary", 1, &is);CHKERRQ(ierr); 270 ierr = DMCreateLabel(dm,"Faces");CHKERRQ(ierr); 271 if (is) { 272 PetscInt d, f, Nf; 273 const PetscInt *faces; 274 PetscInt csize; 275 PetscSection cs; 276 Vec coordinates ; 277 DM cdm; 278 ierr = ISGetLocalSize(is, &Nf);CHKERRQ(ierr); 279 ierr = ISGetIndices(is, &faces);CHKERRQ(ierr); 280 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 281 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 282 ierr = DMGetLocalSection(cdm, &cs);CHKERRQ(ierr); 283 /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */ 284 for (f = 0; f < Nf; ++f) { 285 PetscReal faceCoord; 286 PetscInt b,v; 287 PetscScalar *coords = NULL; 288 PetscInt Nv; 289 ierr = DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);CHKERRQ(ierr); 290 Nv = csize/dim; /* Calculate mean coordinate vector */ 291 for (d = 0; d < dim; ++d) { 292 faceCoord = 0.0; 293 for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]); 294 faceCoord /= Nv; 295 for (b = 0; b < 2; ++b) { 296 if (PetscAbs(faceCoord - b) < PETSC_SMALL) { /* domain have not been set yet, still [0,1]^3 */ 297 ierr = DMSetLabelValue(dm, "Faces", faces[f], d*2+b+1);CHKERRQ(ierr); 298 } 299 } 300 } 301 ierr = DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);CHKERRQ(ierr); 302 } 303 ierr = ISRestoreIndices(is, &faces);CHKERRQ(ierr); 304 } 305 ierr = ISDestroy(&is);CHKERRQ(ierr); 306 ierr = DMGetLabel(dm, "Faces", &label);CHKERRQ(ierr); 307 ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr); 308 } 309 } 310 { 311 PetscInt dimEmbed, i; 312 PetscInt nCoords; 313 PetscScalar *coords,bounds[] = {0,1,-.5,.5,-.5,.5,}; /* x_min,x_max,y_min,y_max */ 314 Vec coordinates; 315 bounds[1] = Lx; 316 if (run_type==1) { 317 for (i = 0; i < 2*dim; i++) bounds[i] = (i%2) ? 1 : 0; 318 } 319 ierr = DMGetCoordinatesLocal(dm,&coordinates);CHKERRQ(ierr); 320 ierr = DMGetCoordinateDim(dm,&dimEmbed);CHKERRQ(ierr); 321 if (dimEmbed != dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"dimEmbed != dim %D",dimEmbed); 322 ierr = VecGetLocalSize(coordinates,&nCoords);CHKERRQ(ierr); 323 if (nCoords % dimEmbed) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Coordinate vector the wrong size"); 324 ierr = VecGetArray(coordinates,&coords);CHKERRQ(ierr); 325 for (i = 0; i < nCoords; i += dimEmbed) { 326 PetscInt j; 327 PetscScalar *coord = &coords[i]; 328 for (j = 0; j < dimEmbed; j++) { 329 coord[j] = bounds[2 * j] + coord[j] * (bounds[2 * j + 1] - bounds[2 * j]); 330 } 331 } 332 ierr = VecRestoreArray(coordinates,&coords);CHKERRQ(ierr); 333 ierr = DMSetCoordinatesLocal(dm,coordinates);CHKERRQ(ierr); 334 } 335 336 /* convert to p4est, and distribute */ 337 338 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr); 339 ierr = PetscOptionsFList("-dm_type","Convert DMPlex to another format (should not be Plex!)","ex56.c",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr); 340 ierr = PetscOptionsEnd(); 341 if (flg) { 342 DM newdm; 343 ierr = DMConvert(dm,convType,&newdm);CHKERRQ(ierr); 344 if (newdm) { 345 const char *prefix; 346 PetscBool isForest; 347 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);CHKERRQ(ierr); 348 ierr = PetscObjectSetOptionsPrefix((PetscObject)newdm,prefix);CHKERRQ(ierr); 349 ierr = DMIsForest(newdm,&isForest);CHKERRQ(ierr); 350 if (!isForest) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Converted to non Forest?"); 351 ierr = DMDestroy(&dm);CHKERRQ(ierr); 352 dm = newdm; 353 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Convert failed?"); 354 } else { 355 PetscPartitioner part; 356 /* Plex Distribute mesh over processes */ 357 ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr); 358 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 359 ierr = DMPlexDistribute(dm, 0, NULL, &distdm);CHKERRQ(ierr); 360 if (distdm) { 361 const char *prefix; 362 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);CHKERRQ(ierr); 363 ierr = PetscObjectSetOptionsPrefix((PetscObject)distdm,prefix);CHKERRQ(ierr); 364 ierr = DMDestroy(&dm);CHKERRQ(ierr); 365 dm = distdm; 366 } 367 } 368 ierr = PetscLogStagePop();CHKERRQ(ierr); 369 basedm = dm; dm = NULL; 370 371 for (iter=0 ; iter<max_conv_its ; iter++) { 372 ierr = PetscLogStagePush(stage[16]);CHKERRQ(ierr); 373 /* make new DM */ 374 ierr = DMClone(basedm, &dm);CHKERRQ(ierr); 375 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, "ex56_");CHKERRQ(ierr); 376 ierr = PetscObjectSetName( (PetscObject)dm,"Mesh");CHKERRQ(ierr); 377 if (max_conv_its > 1) { 378 /* If max_conv_its == 1, then we are not doing a convergence study. */ 379 ierr = PetscOptionsInsertString(NULL,options[iter]);CHKERRQ(ierr); 380 } 381 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); /* refinement done here in Plex, p4est */ 382 /* snes */ 383 ierr = SNESCreate(comm, &snes);CHKERRQ(ierr); 384 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 385 /* fem */ 386 { 387 const PetscInt Ncomp = dim; 388 const PetscInt components[] = {0,1,2}; 389 const PetscInt Nfid = 1, Npid = 1; 390 const PetscInt fid[] = {1}; /* The fixed faces (x=0) */ 391 const PetscInt pid[] = {2}; /* The faces with loading (x=L_x) */ 392 PetscFE fe; 393 PetscDS prob; 394 DMLabel label; 395 DM cdm = dm; 396 397 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, PETSC_FALSE, NULL, PETSC_DECIDE, &fe);CHKERRQ(ierr); /* elasticity */ 398 ierr = PetscObjectSetName((PetscObject) fe, "deformation");CHKERRQ(ierr); 399 /* FEM prob */ 400 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 401 ierr = DMCreateDS(dm);CHKERRQ(ierr); 402 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 403 /* setup problem */ 404 if (run_type==1) { 405 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d);CHKERRQ(ierr); 406 ierr = PetscDSSetResidual(prob, 0, f0_u_x4, f1_u_3d);CHKERRQ(ierr); 407 } else { 408 PetscWeakForm wf; 409 PetscInt bd, i; 410 411 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d_alpha);CHKERRQ(ierr); 412 ierr = PetscDSSetResidual(prob, 0, f0_u, f1_u_3d_alpha);CHKERRQ(ierr); 413 414 ierr = DMGetLabel(dm, "Faces", &label);CHKERRQ(ierr); 415 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "traction", label, Npid, pid, 0, Ncomp, components, NULL, NULL, NULL, &bd);CHKERRQ(ierr); 416 ierr = PetscDSGetBoundary(prob, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 417 for (i = 0; i < Npid; ++i) {ierr = PetscWeakFormSetIndexBdResidual(wf, label, pid[i], 0, 0, 0, f0_bd_u_3d, 0, f1_bd_u);CHKERRQ(ierr);} 418 } 419 /* bcs */ 420 if (run_type==1) { 421 PetscInt id = 1; 422 ierr = DMGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 423 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) zero, NULL, NULL, NULL);CHKERRQ(ierr); 424 } else { 425 ierr = DMGetLabel(dm, "Faces", &label);CHKERRQ(ierr); 426 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, Nfid, fid, 0, Ncomp, components, (void (*)(void)) zero, NULL, NULL, NULL);CHKERRQ(ierr); 427 } 428 while (cdm) { 429 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 430 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 431 } 432 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 433 } 434 /* vecs & mat */ 435 ierr = DMCreateGlobalVector(dm,&xx);CHKERRQ(ierr); 436 ierr = VecDuplicate(xx, &bb);CHKERRQ(ierr); 437 ierr = PetscObjectSetName((PetscObject) bb, "b");CHKERRQ(ierr); 438 ierr = PetscObjectSetName((PetscObject) xx, "u");CHKERRQ(ierr); 439 ierr = DMCreateMatrix(dm, &Amat);CHKERRQ(ierr); 440 ierr = MatSetOption(Amat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); /* Some matrix kernels can take advantage of symmetry if we set this. */ 441 ierr = MatSetOption(Amat,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); /* Inform PETSc that Amat is always symmetric, so info set above isn't lost. */ 442 ierr = MatSetBlockSize(Amat,3);CHKERRQ(ierr); 443 ierr = MatSetOption(Amat,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); 444 ierr = VecGetSize(bb,&N);CHKERRQ(ierr); 445 local_sizes[iter] = N; 446 ierr = PetscInfo2(snes,"%D global equations, %D vertices\n",N,N/dim);CHKERRQ(ierr); 447 if ((use_nearnullspace || attach_nearnullspace) && N/dim > 1) { 448 /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */ 449 DM subdm; 450 MatNullSpace nearNullSpace; 451 PetscInt fields = 0; 452 PetscObject deformation; 453 ierr = DMCreateSubDM(dm, 1, &fields, NULL, &subdm);CHKERRQ(ierr); 454 ierr = DMPlexCreateRigidBody(subdm, 0, &nearNullSpace);CHKERRQ(ierr); 455 ierr = DMGetField(dm, 0, NULL, &deformation);CHKERRQ(ierr); 456 ierr = PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);CHKERRQ(ierr); 457 ierr = DMDestroy(&subdm);CHKERRQ(ierr); 458 if (attach_nearnullspace) { 459 ierr = MatSetNearNullSpace(Amat,nearNullSpace);CHKERRQ(ierr); 460 } 461 ierr = MatNullSpaceDestroy(&nearNullSpace);CHKERRQ(ierr); /* created by DM and destroyed by Mat */ 462 } 463 ierr = DMPlexSetSNESLocalFEM(dm,NULL,NULL,NULL);CHKERRQ(ierr); 464 ierr = SNESSetJacobian(snes, Amat, Amat, NULL, NULL);CHKERRQ(ierr); 465 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 466 ierr = DMSetUp(dm);CHKERRQ(ierr); 467 ierr = PetscLogStagePop();CHKERRQ(ierr); 468 ierr = PetscLogStagePush(stage[16]);CHKERRQ(ierr); 469 /* ksp */ 470 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 471 ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE);CHKERRQ(ierr); 472 /* test BCs */ 473 ierr = VecZeroEntries(xx);CHKERRQ(ierr); 474 if (test_nonzero_cols) { 475 if (!rank) { 476 ierr = VecSetValue(xx,0,1.0,INSERT_VALUES);CHKERRQ(ierr); 477 } 478 ierr = VecAssemblyBegin(xx);CHKERRQ(ierr); 479 ierr = VecAssemblyEnd(xx);CHKERRQ(ierr); 480 } 481 ierr = VecZeroEntries(bb);CHKERRQ(ierr); 482 ierr = VecGetSize(bb,&i);CHKERRQ(ierr); 483 local_sizes[iter] = i; 484 ierr = PetscInfo2(snes,"%D equations in vector, %D vertices\n",i,i/dim);CHKERRQ(ierr); 485 ierr = PetscLogStagePop();CHKERRQ(ierr); 486 /* solve */ 487 ierr = PetscLogStagePush(stage[iter]);CHKERRQ(ierr); 488 ierr = SNESSolve(snes, bb, xx);CHKERRQ(ierr); 489 ierr = PetscLogStagePop();CHKERRQ(ierr); 490 ierr = VecNorm(xx,NORM_INFINITY,&mdisp[iter]);CHKERRQ(ierr); 491 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 492 { 493 PetscViewer viewer = NULL; 494 PetscViewerFormat fmt; 495 ierr = PetscOptionsGetViewer(comm,NULL,"ex56_","-vec_view",&viewer,&fmt,&flg);CHKERRQ(ierr); 496 if (flg) { 497 ierr = PetscViewerPushFormat(viewer,fmt);CHKERRQ(ierr); 498 ierr = VecView(xx,viewer);CHKERRQ(ierr); 499 ierr = VecView(bb,viewer);CHKERRQ(ierr); 500 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 501 } 502 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 503 } 504 /* Free work space */ 505 ierr = DMDestroy(&dm);CHKERRQ(ierr); 506 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 507 ierr = VecDestroy(&xx);CHKERRQ(ierr); 508 ierr = VecDestroy(&bb);CHKERRQ(ierr); 509 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 510 } 511 ierr = DMDestroy(&basedm);CHKERRQ(ierr); 512 if (run_type==1) err[0] = 59.975208 - mdisp[0]; /* error with what I think is the exact solution */ 513 else err[0] = 171.038 - mdisp[0]; 514 for (iter=1 ; iter<max_conv_its ; iter++) { 515 if (run_type==1) err[iter] = 59.975208 - mdisp[iter]; 516 else err[iter] = 171.038 - mdisp[iter]; 517 ierr = PetscPrintf(PETSC_COMM_WORLD,"[%d] %D) N=%12D, max displ=%9.7e, disp diff=%9.2e, error=%4.3e, rate=%3.2g\n",rank,iter,local_sizes[iter],(double)mdisp[iter], 518 (double)(mdisp[iter]-mdisp[iter-1]),(double)err[iter],(double)(PetscLogReal(err[iter-1]/err[iter])/PetscLogReal(2.)));CHKERRQ(ierr); 519 } 520 521 ierr = PetscFinalize(); 522 return ierr; 523 } 524 525 /*TEST 526 527 test: 528 suffix: 0 529 nsize: 4 530 requires: !single 531 args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-10 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -matptap_via scalable -ex56_dm_view 532 timeoutfactor: 2 533 534 # HYPRE PtAP broken with complex numbers 535 test: 536 suffix: hypre 537 requires: hypre !single !complex 538 nsize: 4 539 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -ksp_converged_reason -use_mat_nearnullspace true -petscpartitioner_type simple 540 541 test: 542 suffix: ml 543 requires: ml !single 544 nsize: 4 545 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1.e-8 -pc_type ml -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 3 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type sor -petscpartitioner_type simple -use_mat_nearnullspace 546 547 test: 548 suffix: hpddm 549 requires: hpddm slepc !single define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES) 550 nsize: 4 551 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fgmres -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1.e-8 -pc_type hpddm -petscpartitioner_type simple -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 6 -pc_hpddm_coarse_p 1 -pc_hpddm_coarse_pc_type svd 552 553 test: 554 nsize: 4 555 requires: parmetis !single 556 args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-10 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -pc_gamg_mat_partitioning_type parmetis -pc_gamg_repartition true -snes_converged_reason 557 558 test: 559 suffix: bddc 560 nsize: 4 561 requires: !single 562 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}} -pc_type bddc 563 564 testset: 565 nsize: 4 566 requires: !single 567 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-10 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type aij -pc_type bddc -attach_mat_nearnullspace {{0 1}separate output} 568 test: 569 suffix: bddc_approx_gamg 570 args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop 571 # HYPRE PtAP broken with complex numbers 572 test: 573 requires: hypre !complex 574 suffix: bddc_approx_hypre 575 args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop -prefix_push pc_bddc_neumann_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop 576 test: 577 requires: ml 578 suffix: bddc_approx_ml 579 args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop 580 581 test: 582 suffix: fetidp 583 nsize: 4 584 requires: !single 585 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}} 586 587 test: 588 suffix: bddc_elast 589 nsize: 4 590 requires: !single 591 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -pc_type bddc -pc_bddc_monolithic -attach_mat_nearnullspace 592 593 test: 594 suffix: fetidp_elast 595 nsize: 4 596 requires: !single 597 args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -fetidp_bddc_pc_bddc_monolithic -attach_mat_nearnullspace 598 599 600 testset: 601 nsize: 4 602 requires: !single 603 args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-10 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_converged_reason -snes_monitor_short -ex56_dm_view -petscpartitioner_type simple -pc_gamg_process_eq_limit 20 604 output_file: output/ex56_cuda.out 605 606 test: 607 suffix: cuda 608 requires: cuda 609 args: -ex56_dm_mat_type aijcusparse -ex56_dm_vec_type cuda 610 611 test: 612 suffix: viennacl 613 requires: viennacl 614 args: -ex56_dm_mat_type aijviennacl -ex56_dm_vec_type viennacl 615 616 test: 617 suffix: kokkos 618 requires: kokkos_kernels 619 args: -ex56_dm_mat_type aijkokkos -ex56_dm_vec_type kokkos 620 # Don't run AIJMKL caes with complex scalars because of convergence issues. 621 # Note that we need to test both single and multiple MPI rank cases, because these use different sparse MKL routines to implement the PtAP operation. 622 test: 623 suffix: seqaijmkl 624 nsize: 1 625 requires: define(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) !single !complex 626 args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -mat_block_size 3 -ex56_dm_view -run_type 1 -mat_seqaij_type seqaijmkl 627 timeoutfactor: 2 628 629 test: 630 suffix: mpiaijmkl 631 nsize: 2 632 requires: define(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) !single !complex 633 args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -mat_block_size 3 -ex56_dm_view -run_type 1 -mat_seqaij_type seqaijmkl 634 timeoutfactor: 2 635 636 TEST*/ 637