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