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