1 static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements.\n\ 2 We solve the Poisson problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports discretized auxiliary fields (conductivity) as well as\n\ 5 multilevel nonlinear solvers.\n\n\n"; 6 7 /* 8 A visualization of the adaptation can be accomplished using: 9 10 -dm_adapt_view hdf5:$PWD/adapt.h5 -sol_adapt_view hdf5:$PWD/adapt.h5::append -dm_adapt_pre_view hdf5:$PWD/orig.h5 -sol_adapt_pre_view hdf5:$PWD/orig.h5::append 11 12 Information on refinement: 13 14 -info :~sys,vec,is,mat,ksp,snes,ts 15 */ 16 17 #include <petscdmplex.h> 18 #include <petscdmadaptor.h> 19 #include <petscsnes.h> 20 #include <petscds.h> 21 #include <petscviewerhdf5.h> 22 23 typedef enum {NEUMANN, DIRICHLET, NONE} BCType; 24 typedef enum {RUN_FULL, RUN_EXACT, RUN_TEST, RUN_PERF} RunType; 25 typedef enum {COEFF_NONE, COEFF_ANALYTIC, COEFF_FIELD, COEFF_NONLINEAR, COEFF_CIRCLE, COEFF_CROSS, COEFF_CHECKERBOARD_0, COEFF_CHECKERBOARD_1} CoeffType; 26 27 typedef struct { 28 PetscInt debug; /* The debugging level */ 29 RunType runType; /* Whether to run tests, or solve the full problem */ 30 PetscBool jacobianMF; /* Whether to calculate the Jacobian action on the fly */ 31 PetscLogEvent createMeshEvent; 32 PetscBool showInitial, showSolution, restart, quiet, nonzInit; 33 /* Domain and mesh definition */ 34 PetscInt dim; /* The topological mesh dimension */ 35 DMBoundaryType periodicity[3]; /* The domain periodicity */ 36 PetscInt cells[3]; /* The initial domain division */ 37 char filename[2048]; /* The optional mesh file */ 38 PetscBool interpolate; /* Generate intermediate mesh elements */ 39 PetscReal refinementLimit; /* The largest allowable cell volume */ 40 PetscBool viewHierarchy; /* Whether to view the hierarchy */ 41 PetscBool simplex; /* Simplicial mesh */ 42 /* Problem definition */ 43 BCType bcType; 44 CoeffType variableCoefficient; 45 PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx); 46 PetscBool fieldBC; 47 void (**exactFields)(PetscInt, PetscInt, PetscInt, 48 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 49 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 50 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]); 51 PetscBool bdIntegral; /* Compute the integral of the solution on the boundary */ 52 /* Reproducing tests from SISC 40(3), pp. A1473-A1493, 2018 */ 53 PetscInt div; /* Number of divisions */ 54 PetscInt k; /* Parameter for checkerboard coefficient */ 55 PetscInt *kgrid; /* Random parameter grid */ 56 /* Solver */ 57 PC pcmg; /* This is needed for error monitoring */ 58 PetscBool checkksp; /* Whether to check the KSPSolve for runType == RUN_TEST */ 59 } AppCtx; 60 61 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 62 { 63 u[0] = 0.0; 64 return 0; 65 } 66 67 static PetscErrorCode ecks(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 68 { 69 u[0] = x[0]; 70 return 0; 71 } 72 73 /* 74 In 2D for Dirichlet conditions, we use exact solution: 75 76 u = x^2 + y^2 77 f = 4 78 79 so that 80 81 -\Delta u + f = -4 + 4 = 0 82 83 For Neumann conditions, we have 84 85 -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (bottom) 86 -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top) 87 -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left) 88 -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right) 89 90 Which we can express as 91 92 \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y) 93 94 The boundary integral of this solution is (assuming we are not orienting the edges) 95 96 \int^1_0 x^2 dx + \int^1_0 (1 + y^2) dy + \int^1_0 (x^2 + 1) dx + \int^1_0 y^2 dy = 1/3 + 4/3 + 4/3 + 1/3 = 3 1/3 97 */ 98 static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 99 { 100 *u = x[0]*x[0] + x[1]*x[1]; 101 return 0; 102 } 103 104 static void quadratic_u_field_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 105 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 106 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 107 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uexact[]) 108 { 109 uexact[0] = a[0]; 110 } 111 112 static PetscErrorCode circle_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 113 { 114 const PetscReal alpha = 500.; 115 const PetscReal radius2 = PetscSqr(0.15); 116 const PetscReal r2 = PetscSqr(x[0] - 0.5) + PetscSqr(x[1] - 0.5); 117 const PetscReal xi = alpha*(radius2 - r2); 118 119 *u = PetscTanhScalar(xi) + 1.0; 120 return 0; 121 } 122 123 static PetscErrorCode cross_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 124 { 125 const PetscReal alpha = 50*4; 126 const PetscReal xy = (x[0]-0.5)*(x[1]-0.5); 127 128 *u = PetscSinReal(alpha*xy) * (alpha*PetscAbsReal(xy) < 2*PETSC_PI ? 1 : 0.01); 129 return 0; 130 } 131 132 static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 133 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 134 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 135 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 136 { 137 f0[0] = 4.0; 138 } 139 140 static void f0_circle_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 141 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 142 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 143 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 144 { 145 const PetscReal alpha = 500.; 146 const PetscReal radius2 = PetscSqr(0.15); 147 const PetscReal r2 = PetscSqr(x[0] - 0.5) + PetscSqr(x[1] - 0.5); 148 const PetscReal xi = alpha*(radius2 - r2); 149 150 f0[0] = (-4.0*alpha - 8.0*PetscSqr(alpha)*r2*PetscTanhReal(xi)) * PetscSqr(1.0/PetscCoshReal(xi)); 151 } 152 153 static void f0_cross_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 154 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 155 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 156 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 157 { 158 const PetscReal alpha = 50*4; 159 const PetscReal xy = (x[0]-0.5)*(x[1]-0.5); 160 161 f0[0] = PetscSinReal(alpha*xy) * (alpha*PetscAbsReal(xy) < 2*PETSC_PI ? 1 : 0.01); 162 } 163 164 static void f0_checkerboard_0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 165 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 166 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 167 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 168 { 169 f0[0] = -20.0*PetscExpReal(-(PetscSqr(x[0] - 0.5) + PetscSqr(x[1] - 0.5))); 170 } 171 172 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 173 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 174 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 175 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 176 { 177 PetscInt d; 178 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d]; 179 } 180 181 static void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux, 182 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 183 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 184 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 185 { 186 PetscInt comp; 187 for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0; 188 } 189 190 /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */ 191 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 192 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 193 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 194 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 195 { 196 PetscInt d; 197 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 198 } 199 200 /* < \nabla v, \nabla u + {\nabla u}^T > 201 This just gives \nabla u, give the perdiagonal for the transpose */ 202 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 203 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 204 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 205 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 206 { 207 PetscInt d; 208 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 209 } 210 211 /* 212 In 2D for x periodicity and y Dirichlet conditions, we use exact solution: 213 214 u = sin(2 pi x) 215 f = -4 pi^2 sin(2 pi x) 216 217 so that 218 219 -\Delta u + f = 4 pi^2 sin(2 pi x) - 4 pi^2 sin(2 pi x) = 0 220 */ 221 static PetscErrorCode xtrig_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 222 { 223 *u = PetscSinReal(2.0*PETSC_PI*x[0]); 224 return 0; 225 } 226 227 static void f0_xtrig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 228 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 229 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 230 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 231 { 232 f0[0] = -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[0]); 233 } 234 235 /* 236 In 2D for x-y periodicity, we use exact solution: 237 238 u = sin(2 pi x) sin(2 pi y) 239 f = -8 pi^2 sin(2 pi x) 240 241 so that 242 243 -\Delta u + f = 4 pi^2 sin(2 pi x) sin(2 pi y) + 4 pi^2 sin(2 pi x) sin(2 pi y) - 8 pi^2 sin(2 pi x) = 0 244 */ 245 static PetscErrorCode xytrig_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 246 { 247 *u = PetscSinReal(2.0*PETSC_PI*x[0])*PetscSinReal(2.0*PETSC_PI*x[1]); 248 return 0; 249 } 250 251 static void f0_xytrig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 252 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 253 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 254 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 255 { 256 f0[0] = -8.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[0]); 257 } 258 259 /* 260 In 2D for Dirichlet conditions with a variable coefficient, we use exact solution: 261 262 u = x^2 + y^2 263 f = 6 (x + y) 264 nu = (x + y) 265 266 so that 267 268 -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0 269 */ 270 static PetscErrorCode nu_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 271 { 272 *u = x[0] + x[1]; 273 return 0; 274 } 275 276 static PetscErrorCode checkerboardCoeff(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 277 { 278 AppCtx *user = (AppCtx *) ctx; 279 PetscInt div = user->div; 280 PetscInt k = user->k; 281 PetscInt mask = 0, ind = 0, d; 282 283 PetscFunctionBeginUser; 284 for (d = 0; d < dim; ++d) mask = (mask + (PetscInt) (x[d]*div)) % 2; 285 if (user->kgrid) { 286 for (d = 0; d < dim; ++d) { 287 if (d > 0) ind *= dim; 288 ind += (PetscInt) (x[d]*div); 289 } 290 k = user->kgrid[ind]; 291 } 292 u[0] = mask ? 1.0 : PetscPowRealInt(10.0, -k); 293 PetscFunctionReturn(0); 294 } 295 296 void f0_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 297 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 298 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 299 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 300 { 301 f0[0] = 6.0*(x[0] + x[1]); 302 } 303 304 /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */ 305 void f1_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 306 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 307 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 308 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 309 { 310 PetscInt d; 311 for (d = 0; d < dim; ++d) f1[d] = (x[0] + x[1])*u_x[d]; 312 } 313 314 void f1_field_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 315 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 316 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 317 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 318 { 319 PetscInt d; 320 for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d]; 321 } 322 323 /* < \nabla v, \nabla u + {\nabla u}^T > 324 This just gives \nabla u, give the perdiagonal for the transpose */ 325 void g3_analytic_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 326 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 327 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 328 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 329 { 330 PetscInt d; 331 for (d = 0; d < dim; ++d) g3[d*dim+d] = x[0] + x[1]; 332 } 333 334 void g3_field_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 335 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 336 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 337 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 338 { 339 PetscInt d; 340 for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0]; 341 } 342 343 /* 344 In 2D for Dirichlet conditions with a nonlinear coefficient (p-Laplacian with p = 4), we use exact solution: 345 346 u = x^2 + y^2 347 f = 16 (x^2 + y^2) 348 nu = 1/2 |grad u|^2 349 350 so that 351 352 -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0 353 */ 354 void f0_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 355 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 356 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 357 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 358 { 359 f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]); 360 } 361 362 /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */ 363 void f1_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 364 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 365 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 366 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 367 { 368 PetscScalar nu = 0.0; 369 PetscInt d; 370 for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d]; 371 for (d = 0; d < dim; ++d) f1[d] = 0.5*nu*u_x[d]; 372 } 373 374 /* 375 grad (u + eps w) - grad u = eps grad w 376 377 1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u 378 = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u 379 = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u) 380 = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>) 381 */ 382 void g3_analytic_nonlinear_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 383 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 384 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 385 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 386 { 387 PetscScalar nu = 0.0; 388 PetscInt d, e; 389 for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d]; 390 for (d = 0; d < dim; ++d) { 391 g3[d*dim+d] = 0.5*nu; 392 for (e = 0; e < dim; ++e) { 393 g3[d*dim+e] += u_x[d]*u_x[e]; 394 } 395 } 396 } 397 398 /* 399 In 3D for Dirichlet conditions we use exact solution: 400 401 u = 2/3 (x^2 + y^2 + z^2) 402 f = 4 403 404 so that 405 406 -\Delta u + f = -2/3 * 6 + 4 = 0 407 408 For Neumann conditions, we have 409 410 -\nabla u \cdot -\hat z |_{z=0} = (2z)|_{z=0} = 0 (bottom) 411 -\nabla u \cdot \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top) 412 -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (front) 413 -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back) 414 -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left) 415 -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right) 416 417 Which we can express as 418 419 \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z) 420 */ 421 static PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 422 { 423 *u = 2.0*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/3.0; 424 return 0; 425 } 426 427 static void quadratic_u_field_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 428 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 429 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 430 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uexact[]) 431 { 432 uexact[0] = a[0]; 433 } 434 435 static void bd_integral_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux, 436 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 437 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 438 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *uint) 439 { 440 uint[0] = u[0]; 441 } 442 443 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 444 { 445 const char *bcTypes[3] = {"neumann", "dirichlet", "none"}; 446 const char *runTypes[4] = {"full", "exact", "test", "perf"}; 447 const char *coeffTypes[8] = {"none", "analytic", "field", "nonlinear", "circle", "cross", "checkerboard_0", "checkerboard_1"}; 448 PetscInt bd, bc, run, coeff, n; 449 PetscBool rand = PETSC_FALSE, flg; 450 PetscErrorCode ierr; 451 452 PetscFunctionBeginUser; 453 options->debug = 0; 454 options->runType = RUN_FULL; 455 options->dim = 2; 456 options->periodicity[0] = DM_BOUNDARY_NONE; 457 options->periodicity[1] = DM_BOUNDARY_NONE; 458 options->periodicity[2] = DM_BOUNDARY_NONE; 459 options->cells[0] = 2; 460 options->cells[1] = 2; 461 options->cells[2] = 2; 462 options->filename[0] = '\0'; 463 options->interpolate = PETSC_TRUE; 464 options->refinementLimit = 0.0; 465 options->bcType = DIRICHLET; 466 options->variableCoefficient = COEFF_NONE; 467 options->fieldBC = PETSC_FALSE; 468 options->jacobianMF = PETSC_FALSE; 469 options->showInitial = PETSC_FALSE; 470 options->showSolution = PETSC_FALSE; 471 options->restart = PETSC_FALSE; 472 options->viewHierarchy = PETSC_FALSE; 473 options->simplex = PETSC_TRUE; 474 options->quiet = PETSC_FALSE; 475 options->nonzInit = PETSC_FALSE; 476 options->bdIntegral = PETSC_FALSE; 477 options->checkksp = PETSC_FALSE; 478 options->div = 4; 479 options->k = 1; 480 options->kgrid = NULL; 481 482 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr); 483 ierr = PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);CHKERRQ(ierr); 484 run = options->runType; 485 ierr = PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 4, runTypes[options->runType], &run, NULL);CHKERRQ(ierr); 486 487 options->runType = (RunType) run; 488 489 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 490 bd = options->periodicity[0]; 491 ierr = PetscOptionsEList("-x_periodicity", "The x-boundary periodicity", "ex12.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);CHKERRQ(ierr); 492 options->periodicity[0] = (DMBoundaryType) bd; 493 bd = options->periodicity[1]; 494 ierr = PetscOptionsEList("-y_periodicity", "The y-boundary periodicity", "ex12.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);CHKERRQ(ierr); 495 options->periodicity[1] = (DMBoundaryType) bd; 496 bd = options->periodicity[2]; 497 ierr = PetscOptionsEList("-z_periodicity", "The z-boundary periodicity", "ex12.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);CHKERRQ(ierr); 498 options->periodicity[2] = (DMBoundaryType) bd; 499 n = 3; 500 ierr = PetscOptionsIntArray("-cells", "The initial mesh division", "ex12.c", options->cells, &n, NULL);CHKERRQ(ierr); 501 ierr = PetscOptionsString("-f", "Mesh filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);CHKERRQ(ierr); 502 ierr = PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr); 503 ierr = PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);CHKERRQ(ierr); 504 bc = options->bcType; 505 ierr = PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);CHKERRQ(ierr); 506 options->bcType = (BCType) bc; 507 coeff = options->variableCoefficient; 508 ierr = PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,8,coeffTypes[options->variableCoefficient],&coeff,NULL);CHKERRQ(ierr); 509 options->variableCoefficient = (CoeffType) coeff; 510 511 ierr = PetscOptionsBool("-field_bc", "Use a field representation for the BC", "ex12.c", options->fieldBC, &options->fieldBC, NULL);CHKERRQ(ierr); 512 ierr = PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);CHKERRQ(ierr); 513 ierr = PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);CHKERRQ(ierr); 514 ierr = PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);CHKERRQ(ierr); 515 ierr = PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);CHKERRQ(ierr); 516 ierr = PetscOptionsBool("-dm_view_hierarchy", "View the coarsened hierarchy", "ex12.c", options->viewHierarchy, &options->viewHierarchy, NULL);CHKERRQ(ierr); 517 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex12.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 518 ierr = PetscOptionsBool("-quiet", "Don't print any vecs", "ex12.c", options->quiet, &options->quiet, NULL);CHKERRQ(ierr); 519 ierr = PetscOptionsBool("-nonzero_initial_guess", "nonzero initial guess", "ex12.c", options->nonzInit, &options->nonzInit, NULL);CHKERRQ(ierr); 520 ierr = PetscOptionsBool("-bd_integral", "Compute the integral of the solution on the boundary", "ex12.c", options->bdIntegral, &options->bdIntegral, NULL);CHKERRQ(ierr); 521 if (options->runType == RUN_TEST) { 522 ierr = PetscOptionsBool("-run_test_check_ksp", "Check solution of KSP", "ex12.c", options->checkksp, &options->checkksp, NULL);CHKERRQ(ierr); 523 } 524 ierr = PetscOptionsInt("-div", "The number of division for the checkerboard coefficient", "ex12.c", options->div, &options->div, NULL);CHKERRQ(ierr); 525 ierr = PetscOptionsInt("-k", "The exponent for the checkerboard coefficient", "ex12.c", options->k, &options->k, NULL);CHKERRQ(ierr); 526 ierr = PetscOptionsBool("-k_random", "Assign random k values to checkerboard", "ex12.c", rand, &rand, NULL);CHKERRQ(ierr); 527 ierr = PetscOptionsEnd(); 528 ierr = PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);CHKERRQ(ierr); 529 530 if (rand) { 531 PetscRandom r; 532 PetscReal val; 533 PetscInt N = PetscPowInt(options->div, options->dim), i; 534 535 ierr = PetscMalloc1(N, &options->kgrid);CHKERRQ(ierr); 536 ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr); 537 ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr); 538 ierr = PetscRandomSetInterval(r, 0.0, options->k);CHKERRQ(ierr); 539 ierr = PetscRandomSetSeed(r, 1973);CHKERRQ(ierr); 540 ierr = PetscRandomSeed(r);CHKERRQ(ierr); 541 for (i = 0; i < N; ++i) { 542 ierr = PetscRandomGetValueReal(r, &val);CHKERRQ(ierr); 543 options->kgrid[i] = 1 + (PetscInt) val; 544 } 545 ierr = PetscRandomDestroy(&r);CHKERRQ(ierr); 546 } 547 PetscFunctionReturn(0); 548 } 549 550 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 551 { 552 DM plex; 553 DMLabel label; 554 PetscErrorCode ierr; 555 556 PetscFunctionBeginUser; 557 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 558 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 559 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 560 ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr); 561 ierr = DMDestroy(&plex);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 566 { 567 PetscInt dim = user->dim; 568 const char *filename = user->filename; 569 PetscBool interpolate = user->interpolate; 570 PetscReal refinementLimit = user->refinementLimit; 571 size_t len; 572 PetscErrorCode ierr; 573 574 PetscFunctionBeginUser; 575 ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr); 576 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 577 if (!len) { 578 PetscInt d; 579 580 if (user->periodicity[0] || user->periodicity[1] || user->periodicity[2]) for (d = 0; d < dim; ++d) user->cells[d] = PetscMax(user->cells[d], 3); 581 ierr = DMPlexCreateBoxMesh(comm, dim, user->simplex, user->cells, NULL, NULL, user->periodicity, interpolate, dm);CHKERRQ(ierr); 582 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 583 } else { 584 ierr = DMPlexCreateFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 585 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 586 } 587 { 588 PetscPartitioner part; 589 DM refinedMesh = NULL; 590 DM distributedMesh = NULL; 591 592 /* Refine mesh using a volume constraint */ 593 if (refinementLimit > 0.0) { 594 ierr = DMPlexSetRefinementLimit(*dm, refinementLimit);CHKERRQ(ierr); 595 ierr = DMRefine(*dm, comm, &refinedMesh);CHKERRQ(ierr); 596 if (refinedMesh) { 597 const char *name; 598 599 ierr = PetscObjectGetName((PetscObject) *dm, &name);CHKERRQ(ierr); 600 ierr = PetscObjectSetName((PetscObject) refinedMesh, name);CHKERRQ(ierr); 601 ierr = DMDestroy(dm);CHKERRQ(ierr); 602 *dm = refinedMesh; 603 } 604 } 605 /* Distribute mesh over processes */ 606 ierr = DMPlexGetPartitioner(*dm,&part);CHKERRQ(ierr); 607 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 608 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 609 if (distributedMesh) { 610 ierr = DMDestroy(dm);CHKERRQ(ierr); 611 *dm = distributedMesh; 612 } 613 } 614 if (interpolate) { 615 if (user->bcType == NEUMANN) { 616 DMLabel label; 617 618 ierr = DMCreateLabel(*dm, "boundary");CHKERRQ(ierr); 619 ierr = DMGetLabel(*dm, "boundary", &label);CHKERRQ(ierr); 620 ierr = DMPlexMarkBoundaryFaces(*dm, 1, label);CHKERRQ(ierr); 621 } else if (user->bcType == DIRICHLET) { 622 PetscBool hasLabel; 623 624 ierr = DMHasLabel(*dm,"marker",&hasLabel);CHKERRQ(ierr); 625 if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);} 626 } 627 } 628 { 629 char convType[256]; 630 PetscBool flg; 631 632 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr); 633 ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr); 634 ierr = PetscOptionsEnd(); 635 if (flg) { 636 DM dmConv; 637 638 ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr); 639 if (dmConv) { 640 ierr = DMDestroy(dm);CHKERRQ(ierr); 641 *dm = dmConv; 642 } 643 } 644 } 645 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */ 646 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 647 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 648 if (user->viewHierarchy) { 649 DM cdm = *dm; 650 PetscInt i = 0; 651 char buf[256]; 652 653 while (cdm) { 654 ierr = DMSetUp(cdm);CHKERRQ(ierr); 655 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 656 ++i; 657 } 658 cdm = *dm; 659 while (cdm) { 660 PetscViewer viewer; 661 PetscBool isHDF5, isVTK; 662 663 --i; 664 ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr); 665 ierr = PetscViewerSetType(viewer,PETSCVIEWERHDF5);CHKERRQ(ierr); 666 ierr = PetscViewerSetOptionsPrefix(viewer,"hierarchy_");CHKERRQ(ierr); 667 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 668 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);CHKERRQ(ierr); 669 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);CHKERRQ(ierr); 670 if (isHDF5) { 671 ierr = PetscSNPrintf(buf, 256, "ex12-%d.h5", i);CHKERRQ(ierr); 672 } else if (isVTK) { 673 ierr = PetscSNPrintf(buf, 256, "ex12-%d.vtu", i);CHKERRQ(ierr); 674 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr); 675 } else { 676 ierr = PetscSNPrintf(buf, 256, "ex12-%d", i);CHKERRQ(ierr); 677 } 678 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 679 ierr = PetscViewerFileSetName(viewer,buf);CHKERRQ(ierr); 680 ierr = DMView(cdm, viewer);CHKERRQ(ierr); 681 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 682 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 683 } 684 } 685 ierr = PetscLogEventEnd(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr); 686 PetscFunctionReturn(0); 687 } 688 689 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) 690 { 691 PetscDS prob; 692 const PetscInt id = 1; 693 PetscErrorCode ierr; 694 695 PetscFunctionBeginUser; 696 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 697 switch (user->variableCoefficient) { 698 case COEFF_NONE: 699 if (user->periodicity[0]) { 700 if (user->periodicity[1]) { 701 ierr = PetscDSSetResidual(prob, 0, f0_xytrig_u, f1_u);CHKERRQ(ierr); 702 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 703 } else { 704 ierr = PetscDSSetResidual(prob, 0, f0_xtrig_u, f1_u);CHKERRQ(ierr); 705 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 706 } 707 } else { 708 ierr = PetscDSSetResidual(prob, 0, f0_u, f1_u);CHKERRQ(ierr); 709 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 710 } 711 break; 712 case COEFF_ANALYTIC: 713 ierr = PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);CHKERRQ(ierr); 714 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);CHKERRQ(ierr); 715 break; 716 case COEFF_FIELD: 717 ierr = PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);CHKERRQ(ierr); 718 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);CHKERRQ(ierr); 719 break; 720 case COEFF_NONLINEAR: 721 ierr = PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);CHKERRQ(ierr); 722 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);CHKERRQ(ierr); 723 break; 724 case COEFF_CIRCLE: 725 ierr = PetscDSSetResidual(prob, 0, f0_circle_u, f1_u);CHKERRQ(ierr); 726 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 727 break; 728 case COEFF_CROSS: 729 ierr = PetscDSSetResidual(prob, 0, f0_cross_u, f1_u);CHKERRQ(ierr); 730 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 731 break; 732 case COEFF_CHECKERBOARD_0: 733 ierr = PetscDSSetResidual(prob, 0, f0_checkerboard_0_u, f1_field_u);CHKERRQ(ierr); 734 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);CHKERRQ(ierr); 735 break; 736 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient); 737 } 738 switch (user->dim) { 739 case 2: 740 switch (user->variableCoefficient) { 741 case COEFF_CIRCLE: 742 user->exactFuncs[0] = circle_u_2d;break; 743 case COEFF_CROSS: 744 user->exactFuncs[0] = cross_u_2d;break; 745 case COEFF_CHECKERBOARD_0: 746 user->exactFuncs[0] = zero;break; 747 default: 748 if (user->periodicity[0]) { 749 if (user->periodicity[1]) { 750 user->exactFuncs[0] = xytrig_u_2d; 751 } else { 752 user->exactFuncs[0] = xtrig_u_2d; 753 } 754 } else { 755 user->exactFuncs[0] = quadratic_u_2d; 756 user->exactFields[0] = quadratic_u_field_2d; 757 } 758 } 759 if (user->bcType == NEUMANN) {ierr = PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);CHKERRQ(ierr);} 760 break; 761 case 3: 762 user->exactFuncs[0] = quadratic_u_3d; 763 user->exactFields[0] = quadratic_u_field_3d; 764 if (user->bcType == NEUMANN) {ierr = PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);CHKERRQ(ierr);} 765 break; 766 default: 767 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim); 768 } 769 /* Setup constants */ 770 switch (user->variableCoefficient) { 771 case COEFF_CHECKERBOARD_0: 772 { 773 PetscScalar constants[2]; 774 775 constants[0] = user->div; 776 constants[1] = user->k; 777 ierr = PetscDSSetConstants(prob, 2, constants);CHKERRQ(ierr); 778 } 779 break; 780 default: break; 781 } 782 ierr = PetscDSSetExactSolution(prob, 0, user->exactFuncs[0], user);CHKERRQ(ierr); 783 /* Setup Boundary Conditions */ 784 if (user->bcType != NONE) { 785 ierr = DMAddBoundary(dm, user->bcType == DIRICHLET ? (user->fieldBC ? DM_BC_ESSENTIAL_FIELD : DM_BC_ESSENTIAL) : DM_BC_NATURAL, 786 "wall", user->bcType == DIRICHLET ? "marker" : "boundary", 0, 0, NULL, 787 user->fieldBC ? (void (*)(void)) user->exactFields[0] : (void (*)(void)) user->exactFuncs[0], NULL, 1, &id, user);CHKERRQ(ierr); 788 } 789 PetscFunctionReturn(0); 790 } 791 792 static PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user) 793 { 794 PetscErrorCode (*matFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx) = {nu_2d}; 795 void *ctx[1]; 796 Vec nu; 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 ctx[0] = user; 801 if (user->variableCoefficient == COEFF_CHECKERBOARD_0) {matFuncs[0] = checkerboardCoeff;} 802 ierr = DMCreateLocalVector(dmAux, &nu);CHKERRQ(ierr); 803 ierr = PetscObjectSetName((PetscObject) nu, "Coefficient");CHKERRQ(ierr); 804 ierr = DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctx, INSERT_ALL_VALUES, nu);CHKERRQ(ierr); 805 ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);CHKERRQ(ierr); 806 ierr = VecDestroy(&nu);CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810 static PetscErrorCode SetupBC(DM dm, DM dmAux, AppCtx *user) 811 { 812 PetscErrorCode (*bcFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 813 Vec uexact; 814 PetscInt dim; 815 PetscErrorCode ierr; 816 817 PetscFunctionBegin; 818 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 819 if (dim == 2) bcFuncs[0] = quadratic_u_2d; 820 else bcFuncs[0] = quadratic_u_3d; 821 ierr = DMCreateLocalVector(dmAux, &uexact);CHKERRQ(ierr); 822 ierr = DMProjectFunctionLocal(dmAux, 0.0, bcFuncs, NULL, INSERT_ALL_VALUES, uexact);CHKERRQ(ierr); 823 ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) uexact);CHKERRQ(ierr); 824 ierr = VecDestroy(&uexact);CHKERRQ(ierr); 825 PetscFunctionReturn(0); 826 } 827 828 static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user) 829 { 830 DM dmAux, coordDM; 831 PetscErrorCode ierr; 832 833 PetscFunctionBegin; 834 /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */ 835 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 836 if (!feAux) PetscFunctionReturn(0); 837 ierr = DMClone(dm, &dmAux);CHKERRQ(ierr); 838 ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr); 839 ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr); 840 ierr = DMSetField(dmAux, 0, NULL, (PetscObject) feAux);CHKERRQ(ierr); 841 ierr = DMCreateDS(dmAux);CHKERRQ(ierr); 842 if (user->fieldBC) {ierr = SetupBC(dm, dmAux, user);CHKERRQ(ierr);} 843 else {ierr = SetupMaterial(dm, dmAux, user);CHKERRQ(ierr);} 844 ierr = DMDestroy(&dmAux);CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 849 { 850 DM cdm = dm; 851 const PetscInt dim = user->dim; 852 PetscFE fe, feAux = NULL; 853 PetscBool simplex = user->simplex; 854 MPI_Comm comm; 855 PetscErrorCode ierr; 856 857 PetscFunctionBeginUser; 858 /* Create finite element for each field and auxiliary field */ 859 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 860 ierr = PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe);CHKERRQ(ierr); 861 ierr = PetscObjectSetName((PetscObject) fe, "potential");CHKERRQ(ierr); 862 if (user->variableCoefficient == COEFF_FIELD || user->variableCoefficient == COEFF_CHECKERBOARD_0) { 863 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "mat_", -1, &feAux);CHKERRQ(ierr); 864 ierr = PetscObjectSetName((PetscObject) feAux, "coefficient");CHKERRQ(ierr); 865 ierr = PetscFECopyQuadrature(fe, feAux);CHKERRQ(ierr); 866 } else if (user->fieldBC) { 867 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "bc_", -1, &feAux);CHKERRQ(ierr); 868 ierr = PetscFECopyQuadrature(fe, feAux);CHKERRQ(ierr); 869 } 870 /* Set discretization and boundary conditions for each mesh */ 871 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 872 ierr = DMCreateDS(dm);CHKERRQ(ierr); 873 ierr = SetupProblem(dm, user);CHKERRQ(ierr); 874 while (cdm) { 875 ierr = SetupAuxDM(cdm, feAux, user);CHKERRQ(ierr); 876 if (user->bcType == DIRICHLET && user->interpolate) { 877 PetscBool hasLabel; 878 879 ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr); 880 if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);} 881 } 882 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 883 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 884 } 885 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 886 ierr = PetscFEDestroy(&feAux);CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 #include "petsc/private/petscimpl.h" 891 892 /*@C 893 KSPMonitorError - Outputs the error at each iteration of an iterative solver. 894 895 Collective on KSP 896 897 Input Parameters: 898 + ksp - the KSP 899 . its - iteration number 900 . rnorm - 2-norm, preconditioned residual value (may be estimated). 901 - ctx - monitor context 902 903 Level: intermediate 904 905 .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault() 906 @*/ 907 static PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx) 908 { 909 AppCtx *user = (AppCtx *) ctx; 910 DM dm; 911 Vec du = NULL, r; 912 PetscInt level = 0; 913 PetscBool hasLevel; 914 #if defined(PETSC_HAVE_HDF5) 915 PetscViewer viewer; 916 char buf[256]; 917 #endif 918 PetscErrorCode ierr; 919 920 PetscFunctionBegin; 921 ierr = KSPGetDM(ksp, &dm);CHKERRQ(ierr); 922 /* Calculate solution */ 923 { 924 PC pc = user->pcmg; /* The MG PC */ 925 DM fdm = NULL, cdm = NULL; 926 KSP fksp, cksp; 927 Vec fu, cu = NULL; 928 PetscInt levels, l; 929 930 ierr = KSPBuildSolution(ksp, NULL, &du);CHKERRQ(ierr); 931 ierr = PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);CHKERRQ(ierr); 932 ierr = PCMGGetLevels(pc, &levels);CHKERRQ(ierr); 933 ierr = PCMGGetSmoother(pc, levels-1, &fksp);CHKERRQ(ierr); 934 ierr = KSPBuildSolution(fksp, NULL, &fu);CHKERRQ(ierr); 935 for (l = levels-1; l > level; --l) { 936 Mat R; 937 Vec s; 938 939 ierr = PCMGGetSmoother(pc, l-1, &cksp);CHKERRQ(ierr); 940 ierr = KSPGetDM(cksp, &cdm);CHKERRQ(ierr); 941 ierr = DMGetGlobalVector(cdm, &cu);CHKERRQ(ierr); 942 ierr = PCMGGetRestriction(pc, l, &R);CHKERRQ(ierr); 943 ierr = PCMGGetRScale(pc, l, &s);CHKERRQ(ierr); 944 ierr = MatRestrict(R, fu, cu);CHKERRQ(ierr); 945 ierr = VecPointwiseMult(cu, cu, s);CHKERRQ(ierr); 946 if (l < levels-1) {ierr = DMRestoreGlobalVector(fdm, &fu);CHKERRQ(ierr);} 947 fdm = cdm; 948 fu = cu; 949 } 950 if (levels-1 > level) { 951 ierr = VecAXPY(du, 1.0, cu);CHKERRQ(ierr); 952 ierr = DMRestoreGlobalVector(cdm, &cu);CHKERRQ(ierr); 953 } 954 } 955 /* Calculate error */ 956 ierr = DMGetGlobalVector(dm, &r);CHKERRQ(ierr); 957 ierr = DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);CHKERRQ(ierr); 958 ierr = VecAXPY(r,-1.0,du);CHKERRQ(ierr); 959 ierr = PetscObjectSetName((PetscObject) r, "solution error");CHKERRQ(ierr); 960 /* View error */ 961 #if defined(PETSC_HAVE_HDF5) 962 ierr = PetscSNPrintf(buf, 256, "ex12-%D.h5", level);CHKERRQ(ierr); 963 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);CHKERRQ(ierr); 964 ierr = VecView(r, viewer);CHKERRQ(ierr); 965 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 966 #endif 967 ierr = DMRestoreGlobalVector(dm, &r);CHKERRQ(ierr); 968 PetscFunctionReturn(0); 969 } 970 971 /*@C 972 SNESMonitorError - Outputs the error at each iteration of an iterative solver. 973 974 Collective on SNES 975 976 Input Parameters: 977 + snes - the SNES 978 . its - iteration number 979 . rnorm - 2-norm of residual 980 - ctx - user context 981 982 Level: intermediate 983 984 .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution() 985 @*/ 986 static PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx) 987 { 988 AppCtx *user = (AppCtx *) ctx; 989 DM dm; 990 Vec u, r; 991 PetscInt level = -1; 992 PetscBool hasLevel; 993 #if defined(PETSC_HAVE_HDF5) 994 PetscViewer viewer; 995 #endif 996 char buf[256]; 997 PetscErrorCode ierr; 998 999 PetscFunctionBegin; 1000 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1001 /* Calculate error */ 1002 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 1003 ierr = DMGetGlobalVector(dm, &r);CHKERRQ(ierr); 1004 ierr = PetscObjectSetName((PetscObject) r, "solution error");CHKERRQ(ierr); 1005 ierr = DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);CHKERRQ(ierr); 1006 ierr = VecAXPY(r, -1.0, u);CHKERRQ(ierr); 1007 /* View error */ 1008 ierr = PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);CHKERRQ(ierr); 1009 ierr = PetscSNPrintf(buf, 256, "ex12-%D.h5", level);CHKERRQ(ierr); 1010 #if defined(PETSC_HAVE_HDF5) 1011 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);CHKERRQ(ierr); 1012 ierr = VecView(r, viewer);CHKERRQ(ierr); 1013 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1014 /* Cleanup */ 1015 ierr = DMRestoreGlobalVector(dm, &r);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 #else 1018 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"You need to configure with --download-hdf5"); 1019 #endif 1020 } 1021 1022 int main(int argc, char **argv) 1023 { 1024 DM dm; /* Problem specification */ 1025 SNES snes; /* nonlinear solver */ 1026 Vec u; /* solution vector */ 1027 Mat A,J; /* Jacobian matrix */ 1028 MatNullSpace nullSpace; /* May be necessary for Neumann conditions */ 1029 AppCtx user; /* user-defined work context */ 1030 JacActionCtx userJ; /* context for Jacobian MF action */ 1031 PetscReal error = 0.0; /* L_2 error in the solution */ 1032 PetscBool isFAS; 1033 PetscErrorCode ierr; 1034 1035 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1036 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1037 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 1038 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1039 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 1040 ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr); 1041 1042 ierr = PetscMalloc2(1, &user.exactFuncs, 1, &user.exactFields);CHKERRQ(ierr); 1043 ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr); 1044 1045 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 1046 ierr = PetscObjectSetName((PetscObject) u, "potential");CHKERRQ(ierr); 1047 1048 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 1049 if (user.jacobianMF) { 1050 PetscInt M, m, N, n; 1051 1052 ierr = MatGetSize(J, &M, &N);CHKERRQ(ierr); 1053 ierr = MatGetLocalSize(J, &m, &n);CHKERRQ(ierr); 1054 ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr); 1055 ierr = MatSetSizes(A, m, n, M, N);CHKERRQ(ierr); 1056 ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr); 1057 ierr = MatSetUp(A);CHKERRQ(ierr); 1058 #if 0 1059 ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);CHKERRQ(ierr); 1060 #endif 1061 1062 userJ.dm = dm; 1063 userJ.J = J; 1064 userJ.user = &user; 1065 1066 ierr = DMCreateLocalVector(dm, &userJ.u);CHKERRQ(ierr); 1067 if (user.fieldBC) {ierr = DMProjectFieldLocal(dm, 0.0, userJ.u, user.exactFields, INSERT_BC_VALUES, userJ.u);CHKERRQ(ierr);} 1068 else {ierr = DMProjectFunctionLocal(dm, 0.0, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);CHKERRQ(ierr);} 1069 ierr = MatShellSetContext(A, &userJ);CHKERRQ(ierr); 1070 } else { 1071 A = J; 1072 } 1073 1074 nullSpace = NULL; 1075 if (user.bcType != DIRICHLET) { 1076 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);CHKERRQ(ierr); 1077 ierr = MatSetNullSpace(A, nullSpace);CHKERRQ(ierr); 1078 } 1079 1080 ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr); 1081 ierr = SNESSetJacobian(snes, A, J, NULL, NULL);CHKERRQ(ierr); 1082 1083 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1084 1085 if (user.fieldBC) {ierr = DMProjectField(dm, 0.0, u, user.exactFields, INSERT_ALL_VALUES, u);CHKERRQ(ierr);} 1086 else {ierr = DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);CHKERRQ(ierr);} 1087 if (user.restart) { 1088 #if defined(PETSC_HAVE_HDF5) 1089 PetscViewer viewer; 1090 1091 ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr); 1092 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 1093 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 1094 ierr = PetscViewerFileSetName(viewer, user.filename);CHKERRQ(ierr); 1095 ierr = PetscViewerHDF5PushGroup(viewer, "/fields");CHKERRQ(ierr); 1096 ierr = VecLoad(u, viewer);CHKERRQ(ierr); 1097 ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr); 1098 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1099 #endif 1100 } 1101 if (user.showInitial) { 1102 Vec lv; 1103 ierr = DMGetLocalVector(dm, &lv);CHKERRQ(ierr); 1104 ierr = DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr); 1105 ierr = DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr); 1106 ierr = DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);CHKERRQ(ierr); 1107 ierr = DMRestoreLocalVector(dm, &lv);CHKERRQ(ierr); 1108 } 1109 if (user.viewHierarchy) { 1110 SNES lsnes; 1111 KSP ksp; 1112 PC pc; 1113 PetscInt numLevels, l; 1114 PetscBool isMG; 1115 1116 ierr = PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);CHKERRQ(ierr); 1117 if (isFAS) { 1118 ierr = SNESFASGetLevels(snes, &numLevels);CHKERRQ(ierr); 1119 for (l = 0; l < numLevels; ++l) { 1120 ierr = SNESFASGetCycleSNES(snes, l, &lsnes);CHKERRQ(ierr); 1121 ierr = SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);CHKERRQ(ierr); 1122 } 1123 } else { 1124 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 1125 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 1126 ierr = PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);CHKERRQ(ierr); 1127 if (isMG) { 1128 user.pcmg = pc; 1129 ierr = PCMGGetLevels(pc, &numLevels);CHKERRQ(ierr); 1130 for (l = 0; l < numLevels; ++l) { 1131 ierr = PCMGGetSmootherDown(pc, l, &ksp);CHKERRQ(ierr); 1132 ierr = KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);CHKERRQ(ierr); 1133 } 1134 } 1135 } 1136 } 1137 if (user.runType == RUN_FULL || user.runType == RUN_EXACT) { 1138 PetscErrorCode (*initialGuess[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx) = {zero}; 1139 1140 if (user.nonzInit) initialGuess[0] = ecks; 1141 if (user.runType == RUN_FULL) { 1142 ierr = DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 1143 } 1144 if (user.debug) { 1145 ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");CHKERRQ(ierr); 1146 ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1147 } 1148 ierr = VecViewFromOptions(u, NULL, "-guess_vec_view");CHKERRQ(ierr); 1149 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 1150 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 1151 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1152 1153 if (user.showSolution) { 1154 ierr = PetscPrintf(PETSC_COMM_WORLD, "Solution\n");CHKERRQ(ierr); 1155 ierr = VecChop(u, 3.0e-9);CHKERRQ(ierr); 1156 ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1157 } 1158 } else if (user.runType == RUN_PERF) { 1159 Vec r; 1160 PetscReal res = 0.0; 1161 1162 ierr = SNESGetFunction(snes, &r, NULL, NULL);CHKERRQ(ierr); 1163 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr); 1164 ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");CHKERRQ(ierr); 1165 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 1166 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 1167 ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);CHKERRQ(ierr); 1168 } else { 1169 Vec r; 1170 PetscReal res = 0.0, tol = 1.0e-11; 1171 1172 /* Check discretization error */ 1173 ierr = SNESGetFunction(snes, &r, NULL, NULL);CHKERRQ(ierr); 1174 ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");CHKERRQ(ierr); 1175 if (!user.quiet) {ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 1176 ierr = DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);CHKERRQ(ierr); 1177 if (error < tol) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < %2.1e\n", (double)tol);CHKERRQ(ierr);} 1178 else {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)error);CHKERRQ(ierr);} 1179 /* Check residual */ 1180 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr); 1181 ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");CHKERRQ(ierr); 1182 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 1183 if (!user.quiet) {ierr = VecView(r, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 1184 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 1185 ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);CHKERRQ(ierr); 1186 /* Check Jacobian */ 1187 { 1188 Vec b; 1189 1190 ierr = SNESComputeJacobian(snes, u, A, A);CHKERRQ(ierr); 1191 ierr = VecDuplicate(u, &b);CHKERRQ(ierr); 1192 ierr = VecSet(r, 0.0);CHKERRQ(ierr); 1193 ierr = SNESComputeFunction(snes, r, b);CHKERRQ(ierr); 1194 ierr = MatMult(A, u, r);CHKERRQ(ierr); 1195 ierr = VecAXPY(r, 1.0, b);CHKERRQ(ierr); 1196 ierr = PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");CHKERRQ(ierr); 1197 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 1198 if (!user.quiet) {ierr = VecView(r, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 1199 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 1200 ierr = PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);CHKERRQ(ierr); 1201 /* check solver */ 1202 if (user.checkksp) { 1203 KSP ksp; 1204 1205 if (nullSpace) { 1206 ierr = MatNullSpaceRemove(nullSpace, u);CHKERRQ(ierr); 1207 } 1208 ierr = SNESComputeJacobian(snes, u, A, J);CHKERRQ(ierr); 1209 ierr = MatMult(A, u, b);CHKERRQ(ierr); 1210 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 1211 ierr = KSPSetOperators(ksp, A, J);CHKERRQ(ierr); 1212 ierr = KSPSolve(ksp, b, r);CHKERRQ(ierr); 1213 ierr = VecAXPY(r, -1.0, u);CHKERRQ(ierr); 1214 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 1215 ierr = PetscPrintf(PETSC_COMM_WORLD, "KSP Error: %g\n", (double)res);CHKERRQ(ierr); 1216 } 1217 ierr = VecDestroy(&b);CHKERRQ(ierr); 1218 } 1219 } 1220 ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr); 1221 { 1222 Vec nu; 1223 1224 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &nu);CHKERRQ(ierr); 1225 if (nu) {ierr = VecViewFromOptions(nu, NULL, "-coeff_view");CHKERRQ(ierr);} 1226 } 1227 1228 if (user.bdIntegral) { 1229 DMLabel label; 1230 PetscInt id = 1; 1231 PetscScalar bdInt = 0.0; 1232 PetscReal exact = 3.3333333333; 1233 1234 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 1235 ierr = DMPlexComputeBdIntegral(dm, u, label, 1, &id, bd_integral_2d, &bdInt, NULL);CHKERRQ(ierr); 1236 ierr = PetscPrintf(PETSC_COMM_WORLD, "Solution boundary integral: %.4g\n", (double) PetscAbsScalar(bdInt));CHKERRQ(ierr); 1237 if (PetscAbsReal(PetscAbsScalar(bdInt) - exact) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Invalid boundary integral %g != %g", (double) PetscAbsScalar(bdInt), (double)exact); 1238 } 1239 1240 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 1241 if (user.jacobianMF) {ierr = VecDestroy(&userJ.u);CHKERRQ(ierr);} 1242 if (A != J) {ierr = MatDestroy(&A);CHKERRQ(ierr);} 1243 ierr = MatDestroy(&J);CHKERRQ(ierr); 1244 ierr = VecDestroy(&u);CHKERRQ(ierr); 1245 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 1246 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1247 ierr = PetscFree2(user.exactFuncs, user.exactFields);CHKERRQ(ierr); 1248 ierr = PetscFree(user.kgrid);CHKERRQ(ierr); 1249 ierr = PetscFinalize(); 1250 return ierr; 1251 } 1252 1253 /*TEST 1254 # 2D serial P1 test 0-4 1255 test: 1256 suffix: 2d_p1_0 1257 requires: triangle 1258 args: -run_type test -refinement_limit 0.0 -bc_type dirichlet -interpolate 0 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1259 1260 test: 1261 suffix: 2d_p1_1 1262 requires: triangle 1263 args: -run_type test -refinement_limit 0.0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1264 1265 test: 1266 suffix: 2d_p1_2 1267 requires: triangle 1268 args: -run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1269 1270 test: 1271 suffix: 2d_p1_neumann_0 1272 requires: triangle 1273 args: -run_type test -refinement_limit 0.0 -bc_type neumann -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -dm_view ascii::ascii_info_detail 1274 1275 test: 1276 suffix: 2d_p1_neumann_1 1277 requires: triangle 1278 args: -run_type test -refinement_limit 0.0625 -bc_type neumann -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1279 1280 # 2D serial P2 test 5-8 1281 test: 1282 suffix: 2d_p2_0 1283 requires: triangle 1284 args: -run_type test -refinement_limit 0.0 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1285 1286 test: 1287 suffix: 2d_p2_1 1288 requires: triangle 1289 args: -run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1290 1291 test: 1292 suffix: 2d_p2_neumann_0 1293 requires: triangle 1294 args: -run_type test -refinement_limit 0.0 -bc_type neumann -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -dm_view ascii::ascii_info_detail 1295 1296 test: 1297 suffix: 2d_p2_neumann_1 1298 requires: triangle 1299 args: -run_type test -refinement_limit 0.0625 -bc_type neumann -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -dm_view ascii::ascii_info_detail 1300 1301 test: 1302 suffix: bd_int_0 1303 requires: triangle 1304 args: -run_type test -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -bd_integral -dm_view -quiet 1305 1306 test: 1307 suffix: bd_int_1 1308 requires: triangle 1309 args: -run_type test -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -bd_integral -dm_view -quiet 1310 1311 # 3D serial P1 test 9-12 1312 test: 1313 suffix: 3d_p1_0 1314 requires: ctetgen 1315 args: -run_type test -dim 3 -refinement_limit 0.0 -bc_type dirichlet -interpolate 0 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -dm_view -cells 1,1,1 1316 1317 test: 1318 suffix: 3d_p1_1 1319 requires: ctetgen 1320 args: -run_type test -dim 3 -refinement_limit 0.0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -dm_view -cells 1,1,1 1321 1322 test: 1323 suffix: 3d_p1_2 1324 requires: ctetgen 1325 args: -run_type test -dim 3 -refinement_limit 0.0125 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -dm_view -cells 1,1,1 1326 1327 test: 1328 suffix: 3d_p1_neumann_0 1329 requires: ctetgen 1330 args: -run_type test -dim 3 -bc_type neumann -interpolate 1 -petscspace_degree 1 -snes_fd -show_initial -dm_plex_print_fem 1 -dm_view -cells 1,1,1 1331 1332 # Analytic variable coefficient 13-20 1333 test: 1334 suffix: 13 1335 requires: triangle 1336 args: -run_type test -refinement_limit 0.0 -variable_coefficient analytic -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1337 test: 1338 suffix: 14 1339 requires: triangle 1340 args: -run_type test -refinement_limit 0.0625 -variable_coefficient analytic -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1341 test: 1342 suffix: 15 1343 requires: triangle 1344 args: -run_type test -refinement_limit 0.0 -variable_coefficient analytic -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1345 test: 1346 suffix: 16 1347 requires: triangle 1348 args: -run_type test -refinement_limit 0.0625 -variable_coefficient analytic -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1349 test: 1350 suffix: 17 1351 requires: ctetgen 1352 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient analytic -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1353 1354 test: 1355 suffix: 18 1356 requires: ctetgen 1357 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient analytic -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1358 1359 test: 1360 suffix: 19 1361 requires: ctetgen 1362 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient analytic -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1363 1364 test: 1365 suffix: 20 1366 requires: ctetgen 1367 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient analytic -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1368 1369 # P1 variable coefficient 21-28 1370 test: 1371 suffix: 21 1372 requires: triangle 1373 args: -run_type test -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1374 1375 test: 1376 suffix: 22 1377 requires: triangle 1378 args: -run_type test -refinement_limit 0.0625 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1379 1380 test: 1381 suffix: 23 1382 requires: triangle 1383 args: -run_type test -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1384 1385 test: 1386 suffix: 24 1387 requires: triangle 1388 args: -run_type test -refinement_limit 0.0625 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1389 1390 test: 1391 suffix: 25 1392 requires: ctetgen 1393 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1394 1395 test: 1396 suffix: 26 1397 requires: ctetgen 1398 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1399 1400 test: 1401 suffix: 27 1402 requires: ctetgen 1403 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1404 1405 test: 1406 suffix: 28 1407 requires: ctetgen 1408 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -mat_petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1409 1410 # P0 variable coefficient 29-36 1411 test: 1412 suffix: 29 1413 requires: triangle 1414 args: -run_type test -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1415 1416 test: 1417 suffix: 30 1418 requires: triangle 1419 args: -run_type test -refinement_limit 0.0625 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 1420 1421 test: 1422 suffix: 31 1423 requires: triangle 1424 args: -run_type test -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1425 1426 test: 1427 requires: triangle 1428 suffix: 32 1429 args: -run_type test -refinement_limit 0.0625 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1430 1431 test: 1432 requires: ctetgen 1433 suffix: 33 1434 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1435 1436 test: 1437 suffix: 34 1438 requires: ctetgen 1439 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient field -interpolate 1 -petscspace_degree 1 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1440 1441 test: 1442 suffix: 35 1443 requires: ctetgen 1444 args: -run_type test -dim 3 -refinement_limit 0.0 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1445 1446 test: 1447 suffix: 36 1448 requires: ctetgen 1449 args: -run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient field -interpolate 1 -petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1450 1451 # Full solve 39-44 1452 test: 1453 suffix: 39 1454 requires: triangle !single 1455 args: -run_type full -refinement_limit 0.015625 -interpolate 1 -petscspace_degree 2 -pc_type gamg -ksp_rtol 1.0e-10 -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason ::ascii_info_detail 1456 test: 1457 suffix: 40 1458 requires: triangle !single 1459 args: -run_type full -refinement_limit 0.015625 -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 2 -pc_type svd -ksp_rtol 1.0e-10 -snes_monitor_short -snes_converged_reason ::ascii_info_detail 1460 test: 1461 suffix: 41 1462 requires: triangle !single 1463 args: -run_type full -refinement_limit 0.03125 -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 1 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short 1464 test: 1465 suffix: 42 1466 requires: triangle !single 1467 args: -run_type full -refinement_limit 0.0625 -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 2 -dm_plex_print_fem 0 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short -fas_levels_2_snes_type newtonls -fas_levels_2_pc_type svd -fas_levels_2_ksp_rtol 1.0e-10 -fas_levels_2_snes_atol 1.0e-11 -fas_levels_2_snes_monitor_short 1468 test: 1469 suffix: 43 1470 requires: triangle !single 1471 nsize: 2 1472 args: -run_type full -refinement_limit 0.03125 -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 1 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short 1473 1474 test: 1475 suffix: 44 1476 requires: triangle !single 1477 nsize: 2 1478 args: -run_type full -refinement_limit 0.0625 -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 2 -dm_plex_print_fem 0 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short -fas_levels_2_snes_type newtonls -fas_levels_2_pc_type svd -fas_levels_2_ksp_rtol 1.0e-10 -fas_levels_2_snes_atol 1.0e-11 -fas_levels_2_snes_monitor_short 1479 1480 # These tests use a loose tolerance just to exercise the PtAP operations for MATIS and multiple PCBDDC setup calls inside PCMG 1481 testset: 1482 requires: triangle !single 1483 nsize: 3 1484 args: -interpolate -run_type full -petscspace_degree 1 -dm_mat_type is -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type bddc -pc_mg_galerkin pmat -ksp_rtol 1.0e-2 -snes_converged_reason -dm_refine_hierarchy 2 -snes_max_it 4 1485 test: 1486 suffix: gmg_bddc 1487 filter: sed -e "s/CONVERGED_FNORM_RELATIVE iterations 3/CONVERGED_FNORM_RELATIVE iterations 4/g" 1488 args: -mg_levels_pc_type jacobi 1489 test: 1490 filter: sed -e "s/iterations [0-4]/iterations 4/g" 1491 suffix: gmg_bddc_lev 1492 args: -mg_levels_pc_type bddc 1493 1494 # Restarting 1495 testset: 1496 suffix: restart 1497 requires: hdf5 triangle !complex 1498 args: -run_type test -refinement_limit 0.0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 1499 test: 1500 args: -dm_view hdf5:sol.h5 -vec_view hdf5:sol.h5::append 1501 test: 1502 args: -f sol.h5 -restart 1503 1504 # Periodicity 1505 test: 1506 suffix: periodic_0 1507 requires: triangle 1508 args: -run_type full -refinement_limit 0.0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -snes_converged_reason ::ascii_info_detail 1509 1510 test: 1511 requires: !complex 1512 suffix: periodic_1 1513 args: -quiet -run_type test -simplex 0 -x_periodicity periodic -y_periodicity periodic -vec_view vtk:test.vtu:vtk_vtu -interpolate 1 -petscspace_degree 1 -dm_refine 1 1514 1515 # 2D serial P1 test with field bc 1516 test: 1517 suffix: field_bc_2d_p1_0 1518 requires: triangle 1519 args: -run_type test -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1520 1521 test: 1522 suffix: field_bc_2d_p1_1 1523 requires: triangle 1524 args: -run_type test -dm_refine 1 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1525 1526 test: 1527 suffix: field_bc_2d_p1_neumann_0 1528 requires: triangle 1529 args: -run_type test -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1530 1531 test: 1532 suffix: field_bc_2d_p1_neumann_1 1533 requires: triangle 1534 args: -run_type test -dm_refine 1 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1535 1536 # 3D serial P1 test with field bc 1537 test: 1538 suffix: field_bc_3d_p1_0 1539 requires: ctetgen 1540 args: -run_type test -dim 3 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1541 1542 test: 1543 suffix: field_bc_3d_p1_1 1544 requires: ctetgen 1545 args: -run_type test -dim 3 -dm_refine 1 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1546 1547 test: 1548 suffix: field_bc_3d_p1_neumann_0 1549 requires: ctetgen 1550 args: -run_type test -dim 3 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1551 1552 test: 1553 suffix: field_bc_3d_p1_neumann_1 1554 requires: ctetgen 1555 args: -run_type test -dim 3 -dm_refine 1 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 1 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1556 1557 # 2D serial P2 test with field bc 1558 test: 1559 suffix: field_bc_2d_p2_0 1560 requires: triangle 1561 args: -run_type test -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1562 1563 test: 1564 suffix: field_bc_2d_p2_1 1565 requires: triangle 1566 args: -run_type test -dm_refine 1 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1567 1568 test: 1569 suffix: field_bc_2d_p2_neumann_0 1570 requires: triangle 1571 args: -run_type test -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1572 1573 test: 1574 suffix: field_bc_2d_p2_neumann_1 1575 requires: triangle 1576 args: -run_type test -dm_refine 1 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 1577 1578 # 3D serial P2 test with field bc 1579 test: 1580 suffix: field_bc_3d_p2_0 1581 requires: ctetgen 1582 args: -run_type test -dim 3 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1583 1584 test: 1585 suffix: field_bc_3d_p2_1 1586 requires: ctetgen 1587 args: -run_type test -dim 3 -dm_refine 1 -interpolate 1 -bc_type dirichlet -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1588 1589 test: 1590 suffix: field_bc_3d_p2_neumann_0 1591 requires: ctetgen 1592 args: -run_type test -dim 3 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1593 1594 test: 1595 suffix: field_bc_3d_p2_neumann_1 1596 requires: ctetgen 1597 args: -run_type test -dim 3 -dm_refine 1 -interpolate 1 -bc_type neumann -field_bc -petscspace_degree 2 -bc_petscspace_degree 2 -show_initial -dm_plex_print_fem 1 -cells 1,1,1 1598 1599 # Full solve simplex: Convergence 1600 test: 1601 suffix: tet_conv_p1_r0 1602 requires: ctetgen 1603 args: -run_type full -dim 3 -dm_refine 0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -dm_view -snes_converged_reason ::ascii_info_detail -pc_type lu -cells 1,1,1 1604 test: 1605 suffix: tet_conv_p1_r2 1606 requires: ctetgen 1607 args: -run_type full -dim 3 -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -dm_view -snes_converged_reason ::ascii_info_detail -pc_type lu -cells 1,1,1 1608 test: 1609 suffix: tet_conv_p1_r3 1610 requires: ctetgen 1611 args: -run_type full -dim 3 -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -dm_view -snes_converged_reason ::ascii_info_detail -pc_type lu -cells 1,1,1 1612 test: 1613 suffix: tet_conv_p2_r0 1614 requires: ctetgen 1615 args: -run_type full -dim 3 -dm_refine 0 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -dm_view -snes_converged_reason ::ascii_info_detail -pc_type lu -cells 1,1,1 1616 test: 1617 suffix: tet_conv_p2_r2 1618 requires: ctetgen 1619 args: -run_type full -dim 3 -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -dm_view -snes_converged_reason ::ascii_info_detail -pc_type lu -cells 1,1,1 1620 1621 # Full solve simplex: PCBDDC 1622 test: 1623 suffix: tri_bddc 1624 requires: triangle !single 1625 nsize: 5 1626 args: -run_type full -petscpartitioner_type simple -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -dm_mat_type is -pc_type bddc -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1627 1628 # Full solve simplex: PCBDDC 1629 test: 1630 suffix: tri_parmetis_bddc 1631 requires: triangle !single parmetis 1632 nsize: 4 1633 args: -run_type full -petscpartitioner_type parmetis -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -dm_mat_type is -pc_type bddc -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1634 1635 testset: 1636 args: -run_type full -petscpartitioner_type simple -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -dm_mat_type is -pc_type bddc -ksp_type gmres -snes_monitor_short -ksp_monitor_short -snes_view -simplex 0 -petscspace_poly_tensor -pc_bddc_corner_selection -cells 3,3 -ksp_rtol 1.e-9 -pc_bddc_use_edges 0 1637 nsize: 5 1638 output_file: output/ex12_quad_bddc.out 1639 filter: sed -e "s/aijcusparse/aij/g" -e "s/aijviennacl/aij/g" -e "s/factorization: cusparse/factorization: petsc/g" 1640 test: 1641 requires: !single 1642 suffix: quad_bddc 1643 test: 1644 requires: !single cuda 1645 suffix: quad_bddc_cuda 1646 args: -matis_localmat_type aijcusparse -pc_bddc_dirichlet_pc_factor_mat_solver_type cusparse -pc_bddc_neumann_pc_factor_mat_solver_type cusparse 1647 test: 1648 requires: !single viennacl 1649 suffix: quad_bddc_viennacl 1650 args: -matis_localmat_type aijviennacl 1651 1652 # Full solve simplex: ASM 1653 test: 1654 suffix: tri_q2q1_asm_lu 1655 requires: triangle !single 1656 args: -run_type full -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type asm -pc_asm_type restrict -pc_asm_blocks 4 -sub_pc_type lu -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1657 1658 test: 1659 suffix: tri_q2q1_msm_lu 1660 requires: triangle !single 1661 args: -run_type full -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type asm -pc_asm_type restrict -pc_asm_local_type multiplicative -pc_asm_blocks 4 -sub_pc_type lu -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1662 1663 test: 1664 suffix: tri_q2q1_asm_sor 1665 requires: triangle !single 1666 args: -run_type full -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type asm -pc_asm_type restrict -pc_asm_blocks 4 -sub_pc_type sor -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1667 1668 test: 1669 suffix: tri_q2q1_msm_sor 1670 requires: triangle !single 1671 args: -run_type full -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type asm -pc_asm_type restrict -pc_asm_local_type multiplicative -pc_asm_blocks 4 -sub_pc_type sor -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 1672 1673 # Full solve simplex: FAS 1674 test: 1675 suffix: fas_newton_0 1676 requires: triangle !single 1677 args: -run_type full -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 1 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short 1678 1679 test: 1680 suffix: fas_newton_1 1681 requires: triangle !single 1682 args: -run_type full -dm_refine_hierarchy 3 -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type lu -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_snes_linesearch_type basic -fas_levels_ksp_rtol 1.0e-10 -fas_levels_snes_monitor_short 1683 filter: sed -e "s/total number of linear solver iterations=14/total number of linear solver iterations=15/g" 1684 1685 test: 1686 suffix: fas_ngs_0 1687 requires: triangle !single 1688 args: -run_type full -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_refine_hierarchy 1 -snes_view -fas_levels_1_snes_type ngs -fas_levels_1_snes_monitor_short 1689 1690 test: 1691 suffix: fas_newton_coarse_0 1692 requires: pragmatic triangle 1693 TODO: broken 1694 args: -run_type full -dm_refine 2 -dm_plex_hash_location -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -pc_type svd -ksp_rtol 1.0e-10 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -dm_coarsen_hierarchy 1 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short 1695 1696 test: 1697 suffix: mg_newton_coarse_0 1698 requires: triangle pragmatic 1699 TODO: broken 1700 args: -run_type full -dm_refine 3 -interpolate 1 -petscspace_degree 1 -snes_monitor_short -ksp_monitor_true_residual -snes_converged_reason ::ascii_info_detail -dm_coarsen_hierarchy 3 -dm_plex_hash_location -snes_view -dm_view -ksp_type richardson -pc_type mg -pc_mg_levels 4 -snes_atol 1.0e-8 -ksp_atol 1.0e-8 -snes_rtol 0.0 -ksp_rtol 0.0 -ksp_norm_type unpreconditioned -mg_levels_ksp_type gmres -mg_levels_pc_type ilu -mg_levels_ksp_max_it 10 1701 1702 test: 1703 suffix: mg_newton_coarse_1 1704 requires: triangle pragmatic 1705 TODO: broken 1706 args: -run_type full -dm_refine 5 -interpolate 1 -petscspace_degree 1 -dm_coarsen_hierarchy 5 -dm_plex_hash_location -dm_plex_separate_marker -dm_plex_coarsen_bd_label marker -dm_plex_remesh_bd -ksp_type richardson -ksp_rtol 1.0e-12 -pc_type mg -pc_mg_levels 3 -mg_levels_ksp_max_it 2 -snes_converged_reason ::ascii_info_detail -snes_monitor -ksp_monitor_true_residual -mg_levels_ksp_monitor_true_residual -dm_view -ksp_view 1707 1708 test: 1709 suffix: mg_newton_coarse_2 1710 requires: triangle pragmatic 1711 TODO: broken 1712 args: -run_type full -dm_refine 5 -interpolate 1 -petscspace_degree 1 -dm_coarsen_hierarchy 5 -dm_plex_hash_location -dm_plex_separate_marker -dm_plex_remesh_bd -ksp_type richardson -ksp_rtol 1.0e-12 -pc_type mg -pc_mg_levels 3 -mg_levels_ksp_max_it 2 -snes_converged_reason ::ascii_info_detail -snes_monitor -ksp_monitor_true_residual -mg_levels_ksp_monitor_true_residual -dm_view -ksp_view 1713 1714 # Full solve tensor 1715 test: 1716 suffix: tensor_plex_2d 1717 args: -run_type test -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dm_refine_hierarchy 2 -cells 2,2 1718 1719 test: 1720 suffix: tensor_p4est_2d 1721 requires: p4est 1722 args: -run_type test -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dm_forest_initial_refinement 2 -dm_forest_minimum_refinement 0 -dm_plex_convert_type p4est -cells 2,2 1723 1724 test: 1725 suffix: tensor_plex_3d 1726 args: -run_type test -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dim 3 -dm_refine_hierarchy 1 -cells 2,2,2 1727 1728 test: 1729 suffix: tensor_p4est_3d 1730 requires: p4est 1731 args: -run_type test -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dm_forest_initial_refinement 1 -dm_forest_minimum_refinement 0 -dim 3 -dm_plex_convert_type p8est -cells 2,2,2 1732 1733 test: 1734 suffix: p4est_test_q2_conformal_serial 1735 requires: p4est 1736 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -cells 2,2 1737 1738 test: 1739 suffix: p4est_test_q2_conformal_parallel 1740 requires: p4est 1741 nsize: 7 1742 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -petscpartitioner_type simple -cells 2,2 1743 1744 test: 1745 suffix: p4est_test_q2_conformal_parallel_parmetis 1746 requires: parmetis p4est 1747 nsize: 4 1748 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -petscpartitioner_type parmetis -cells 2,2 1749 1750 test: 1751 suffix: p4est_test_q2_nonconformal_serial 1752 requires: p4est 1753 filter: grep -v "CG or CGNE: variant" 1754 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 1755 1756 test: 1757 suffix: p4est_test_q2_nonconformal_parallel 1758 requires: p4est 1759 filter: grep -v "CG or CGNE: variant" 1760 nsize: 7 1761 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -cells 2,2 1762 1763 test: 1764 suffix: p4est_test_q2_nonconformal_parallel_parmetis 1765 requires: parmetis p4est 1766 nsize: 4 1767 args: -run_type test -interpolate 1 -petscspace_degree 2 -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type parmetis -cells 2,2 1768 1769 test: 1770 suffix: p4est_exact_q2_conformal_serial 1771 requires: p4est !single !complex !__float128 1772 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -cells 2,2 1773 1774 test: 1775 suffix: p4est_exact_q2_conformal_parallel 1776 requires: p4est !single !complex !__float128 1777 nsize: 4 1778 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -cells 2,2 1779 1780 test: 1781 suffix: p4est_exact_q2_conformal_parallel_parmetis 1782 requires: parmetis p4est !single 1783 nsize: 4 1784 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -petscpartitioner_type parmetis -cells 2,2 1785 1786 test: 1787 suffix: p4est_exact_q2_nonconformal_serial 1788 requires: p4est 1789 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 1790 1791 test: 1792 suffix: p4est_exact_q2_nonconformal_parallel 1793 requires: p4est 1794 nsize: 7 1795 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -cells 2,2 1796 1797 test: 1798 suffix: p4est_exact_q2_nonconformal_parallel_parmetis 1799 requires: parmetis p4est 1800 nsize: 4 1801 args: -run_type exact -interpolate 1 -petscspace_degree 2 -fas_levels_snes_atol 1.e-10 -snes_max_it 1 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type none -fas_coarse_ksp_type preonly -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type none -fas_levels_ksp_type preonly -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type parmetis -cells 2,2 1802 1803 test: 1804 suffix: p4est_full_q2_nonconformal_serial 1805 requires: p4est !single 1806 filter: grep -v "variant HERMITIAN" 1807 args: -run_type full -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type jacobi -fas_coarse_ksp_type cg -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type jacobi -fas_levels_ksp_type cg -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 1808 1809 test: 1810 suffix: p4est_full_q2_nonconformal_parallel 1811 requires: p4est !single 1812 filter: grep -v "variant HERMITIAN" 1813 nsize: 7 1814 args: -run_type full -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type fas -snes_fas_levels 3 -fas_coarse_pc_type jacobi -fas_coarse_ksp_type cg -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type jacobi -fas_levels_ksp_type cg -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -cells 2,2 1815 1816 test: 1817 suffix: p4est_full_q2_nonconformal_parallel_bddcfas 1818 requires: p4est !single 1819 filter: grep -v "variant HERMITIAN" 1820 nsize: 7 1821 args: -run_type full -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type fas -snes_fas_levels 3 -dm_mat_type is -fas_coarse_pc_type bddc -fas_coarse_ksp_type cg -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type bddc -fas_levels_ksp_type cg -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -cells 2,2 1822 1823 test: 1824 suffix: p4est_full_q2_nonconformal_parallel_bddc 1825 requires: p4est !single 1826 filter: grep -v "variant HERMITIAN" 1827 nsize: 7 1828 args: -run_type full -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type newtonls -dm_mat_type is -pc_type bddc -ksp_type cg -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -cells 2,2 1829 1830 test: 1831 TODO: broken 1832 suffix: p4est_fas_q2_conformal_serial 1833 requires: p4est !complex !__float128 1834 args: -run_type full -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type fas -snes_fas_levels 3 -pc_type jacobi -ksp_type gmres -fas_coarse_pc_type svd -fas_coarse_ksp_type gmres -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type svd -fas_levels_ksp_type gmres -fas_levels_snes_monitor_short -simplex 0 -dm_refine_hierarchy 3 -cells 2,2 1835 1836 test: 1837 TODO: broken 1838 suffix: p4est_fas_q2_nonconformal_serial 1839 requires: p4est 1840 args: -run_type full -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type fas -snes_fas_levels 3 -pc_type jacobi -ksp_type gmres -fas_coarse_pc_type jacobi -fas_coarse_ksp_type gmres -fas_coarse_ksp_monitor_true_residual -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_snes_type newtonls -fas_levels_pc_type jacobi -fas_levels_ksp_type gmres -fas_levels_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 1841 1842 test: 1843 suffix: fas_newton_0_p4est 1844 requires: p4est !single !__float128 1845 args: -run_type full -variable_coefficient nonlinear -interpolate 1 -petscspace_degree 1 -snes_type fas -snes_fas_levels 2 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -fas_coarse_snes_linesearch_type basic -snes_converged_reason ::ascii_info_detail -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor_short -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 1846 1847 # Full solve simplicial AMR 1848 test: 1849 suffix: tri_p1_adapt_0 1850 requires: pragmatic 1851 TODO: broken 1852 args: -run_type exact -dim 2 -dm_refine 5 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -variable_coefficient circle -snes_converged_reason ::ascii_info_detail -pc_type lu -adaptor_refinement_factor 1.0 -dm_view -dm_adapt_view -snes_adapt_initial 1 1853 1854 test: 1855 suffix: tri_p1_adapt_1 1856 requires: pragmatic 1857 TODO: broken 1858 args: -run_type exact -dim 2 -dm_refine 5 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -variable_coefficient circle -snes_converged_reason ::ascii_info_detail -pc_type lu -adaptor_refinement_factor 1.0 -dm_view -dm_adapt_iter_view -dm_adapt_view -snes_adapt_sequence 2 1859 1860 test: 1861 suffix: tri_p1_adapt_analytic_0 1862 requires: pragmatic 1863 TODO: broken 1864 args: -run_type exact -dim 2 -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -variable_coefficient cross -snes_adapt_initial 4 -adaptor_target_num 500 -adaptor_monitor -dm_view -dm_adapt_iter_view 1865 1866 # Full solve tensor AMR 1867 test: 1868 suffix: quad_q1_adapt_0 1869 requires: p4est 1870 args: -run_type exact -dim 2 -simplex 0 -dm_plex_convert_type p4est -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -variable_coefficient circle -snes_converged_reason ::ascii_info_detail -pc_type lu -dm_forest_initial_refinement 4 -snes_adapt_initial 1 -dm_view 1871 filter: grep -v DM_ 1872 1873 test: 1874 suffix: amr_0 1875 nsize: 5 1876 args: -run_type test -petscpartitioner_type simple -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dm_refine 1 -cells 2,2 1877 1878 test: 1879 suffix: amr_1 1880 requires: p4est !complex 1881 args: -run_type test -refinement_limit 0.0 -simplex 0 -interpolate -bc_type dirichlet -petscspace_degree 1 -dm_plex_convert_type p4est -dm_p4est_refine_pattern center -dm_forest_maximum_refinement 5 -dm_view vtk:amr.vtu:vtk_vtu -vec_view vtk:amr.vtu:vtk_vtu:append -cells 2,2 1882 1883 test: 1884 suffix: p4est_solve_bddc 1885 requires: p4est !complex 1886 args: -run_type full -variable_coefficient nonlinear -nonzero_initial_guess 1 -interpolate 1 -petscspace_degree 2 -snes_max_it 20 -snes_type newtonls -dm_mat_type is -pc_type bddc -ksp_type cg -snes_monitor_short -ksp_monitor -snes_linesearch_type bt -snes_converged_reason -snes_view -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -petscpartitioner_type simple -pc_bddc_detect_disconnected 1887 nsize: 4 1888 1889 test: 1890 suffix: p4est_solve_fas 1891 requires: p4est 1892 args: -run_type full -variable_coefficient nonlinear -nonzero_initial_guess 1 -interpolate 1 -petscspace_degree 2 -snes_max_it 10 -snes_type fas -snes_linesearch_type bt -snes_fas_levels 3 -fas_coarse_snes_type newtonls -fas_coarse_snes_linesearch_type basic -fas_coarse_ksp_type cg -fas_coarse_pc_type jacobi -fas_coarse_snes_monitor_short -fas_levels_snes_max_it 4 -fas_levels_snes_type newtonls -fas_levels_snes_linesearch_type bt -fas_levels_ksp_type cg -fas_levels_pc_type jacobi -fas_levels_snes_monitor_short -fas_levels_cycle_snes_linesearch_type bt -snes_monitor_short -snes_converged_reason -snes_view -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash 1893 nsize: 4 1894 TODO: identical machine two runs produce slightly different solver trackers 1895 1896 test: 1897 suffix: p4est_convergence_test_1 1898 requires: p4est 1899 args: -quiet -run_type test -interpolate 1 -petscspace_degree 1 -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash 1900 nsize: 4 1901 1902 test: 1903 suffix: p4est_convergence_test_2 1904 requires: p4est 1905 args: -quiet -run_type test -interpolate 1 -petscspace_degree 1 -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 3 -dm_forest_initial_refinement 3 -dm_forest_maximum_refinement 5 -dm_p4est_refine_pattern hash 1906 1907 test: 1908 suffix: p4est_convergence_test_3 1909 requires: p4est 1910 args: -quiet -run_type test -interpolate 1 -petscspace_degree 1 -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 4 -dm_forest_initial_refinement 4 -dm_forest_maximum_refinement 6 -dm_p4est_refine_pattern hash 1911 1912 test: 1913 suffix: p4est_convergence_test_4 1914 requires: p4est 1915 args: -quiet -run_type test -interpolate 1 -petscspace_degree 1 -simplex 0 -petscspace_poly_tensor -dm_plex_convert_type p4est -dm_forest_minimum_refinement 5 -dm_forest_initial_refinement 5 -dm_forest_maximum_refinement 7 -dm_p4est_refine_pattern hash 1916 timeoutfactor: 5 1917 1918 # Serial tests with GLVis visualization 1919 test: 1920 suffix: glvis_2d_tet_p1 1921 args: -quiet -run_type test -interpolate 1 -bc_type dirichlet -petscspace_degree 1 -vec_view glvis: -f ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 1922 test: 1923 suffix: glvis_2d_tet_p2 1924 args: -quiet -run_type test -interpolate 1 -bc_type dirichlet -petscspace_degree 2 -vec_view glvis: -f ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 1925 test: 1926 suffix: glvis_2d_hex_p1 1927 args: -quiet -run_type test -interpolate 1 -bc_type dirichlet -petscspace_degree 1 -vec_view glvis: -simplex 0 -dm_refine 1 1928 test: 1929 suffix: glvis_2d_hex_p2 1930 args: -quiet -run_type test -interpolate 1 -bc_type dirichlet -petscspace_degree 2 -vec_view glvis: -simplex 0 -dm_refine 1 1931 test: 1932 suffix: glvis_2d_hex_p2_p4est 1933 requires: p4est 1934 args: -quiet -run_type test -interpolate 1 -bc_type dirichlet -petscspace_degree 2 -vec_view glvis: -simplex 0 -dm_plex_convert_type p4est -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 1 -dm_forest_maximum_refinement 4 -dm_p4est_refine_pattern hash -cells 2,2 -viewer_glvis_dm_plex_enable_ncmesh 1935 test: 1936 suffix: glvis_2d_tet_p0 1937 args: -run_type exact -interpolate 1 -guess_vec_view glvis: -nonzero_initial_guess 1 -f ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -petscspace_degree 0 1938 test: 1939 suffix: glvis_2d_hex_p0 1940 args: -run_type exact -interpolate 1 -guess_vec_view glvis: -nonzero_initial_guess 1 -cells 5,7 -simplex 0 -petscspace_degree 0 1941 1942 # PCHPDDM tests 1943 testset: 1944 nsize: 4 1945 requires: hpddm slepc !single 1946 args: -run_type test -run_test_check_ksp -quiet -petscspace_degree 1 -interpolate 1 -petscpartitioner_type simple -bc_type none -simplex 0 -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 2 -pc_hpddm_coarse_p 1 -pc_hpddm_coarse_pc_type svd -ksp_rtol 1.e-10 -pc_hpddm_levels_1_st_pc_factor_shift_type INBLOCKS -ksp_converged_reason 1947 test: 1948 suffix: quad_singular_hpddm 1949 args: -cells 6,7 1950 test: 1951 requires: p4est 1952 suffix: p4est_singular_2d_hpddm 1953 args: -dm_plex_convert_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 3 -dm_forest_maximum_refinement 3 1954 test: 1955 requires: p4est 1956 suffix: p4est_nc_singular_2d_hpddm 1957 args: -dm_plex_convert_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 1 -dm_forest_maximum_refinement 3 -dm_p4est_refine_pattern hash 1958 testset: 1959 nsize: 4 1960 requires: hpddm slepc triangle !single 1961 args: -run_type full -petscpartitioner_type simple -dm_refine 2 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -ksp_type gmres -ksp_gmres_restart 100 -pc_type hpddm -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 4 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_rtol 1.e-1 1962 test: 1963 args: -pc_hpddm_coarse_mat_type baij -options_left no 1964 suffix: tri_hpddm_reuse_baij 1965 test: 1966 requires: !complex 1967 suffix: tri_hpddm_reuse 1968 testset: 1969 nsize: 4 1970 requires: hpddm slepc !single 1971 args: -run_type full -petscpartitioner_type simple -cells 7,5 -dm_refine 2 -simplex 0 -bc_type dirichlet -interpolate 1 -petscspace_degree 2 -ksp_type gmres -ksp_gmres_restart 100 -pc_type hpddm -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 4 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_rtol 1.e-1 1972 test: 1973 args: -pc_hpddm_coarse_mat_type baij -options_left no 1974 suffix: quad_hpddm_reuse_baij 1975 test: 1976 requires: !complex 1977 suffix: quad_hpddm_reuse 1978 testset: 1979 nsize: 4 1980 requires: hpddm slepc !single 1981 args: -run_type full -petscpartitioner_type simple -cells 7,5 -dm_refine 2 -simplex 0 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -pc_type hpddm -snes_monitor_short -ksp_monitor_short -snes_converged_reason ::ascii_info_detail -ksp_converged_reason -snes_view -show_solution 0 -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_threshold 0.1 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_rtol 1.e-1 1982 test: 1983 args: -pc_hpddm_coarse_mat_type baij -options_left no 1984 suffix: quad_hpddm_reuse_threshold_baij 1985 test: 1986 requires: !complex 1987 suffix: quad_hpddm_reuse_threshold 1988 testset: 1989 nsize: 4 1990 requires: hpddm slepc parmetis !single 1991 filter: sed -e "s/linear solver iterations=17/linear solver iterations=16/g" 1992 args: -run_type full -petscpartitioner_type parmetis -dm_refine 3 -bc_type dirichlet -interpolate 1 -petscspace_degree 1 -ksp_type gmres -ksp_gmres_restart 100 -pc_type hpddm -snes_monitor_short -snes_converged_reason ::ascii_info_detail -snes_view -show_solution 0 -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type icc -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_rtol 1.e-10 -f ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -pc_hpddm_levels_1_sub_pc_factor_levels 3 -variable_coefficient circle -dm_plex_gmsh_periodic 0 1993 test: 1994 args: -pc_hpddm_coarse_mat_type baij -options_left no 1995 suffix: tri_parmetis_hpddm_baij 1996 test: 1997 requires: !complex 1998 suffix: tri_parmetis_hpddm 1999 2000 # 2D serial P1 tests for adaptive MG 2001 test: 2002 suffix: 2d_p1_adaptmg_0 2003 requires: triangle bamg 2004 args: -dm_refine_hierarchy 3 -cells 4,4 -interpolate -bc_type dirichlet -petscspace_degree 1 \ 2005 -variable_coefficient checkerboard_0 -mat_petscspace_degree 0 -div 16 -k 3 \ 2006 -snes_max_it 1 -ksp_converged_reason \ 2007 -ksp_rtol 1e-8 -pc_type mg 2008 # -ksp_monitor_true_residual -ksp_converged_reason -mg_levels_ksp_monitor_true_residual -pc_mg_mesp_monitor -dm_adapt_interp_view_fine draw -dm_adapt_interp_view_coarse draw -draw_pause 1 2009 test: 2010 suffix: 2d_p1_adaptmg_1 2011 requires: triangle bamg 2012 args: -dm_refine_hierarchy 3 -cells 4,4 -interpolate -bc_type dirichlet -petscspace_degree 1 \ 2013 -variable_coefficient checkerboard_0 -mat_petscspace_degree 0 -div 16 -k 3 \ 2014 -snes_max_it 1 -ksp_converged_reason \ 2015 -ksp_rtol 1e-8 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp -pc_mg_adapt_interp_coarse_space eigenvector -pc_mg_adapt_interp_n 1 \ 2016 -pc_mg_mesp_ksp_type richardson -pc_mg_mesp_ksp_richardson_self_scale -pc_mg_mesp_ksp_max_it 100 -pc_mg_mesp_pc_type none 2017 2018 TEST*/ 2019