1 static char help[] = "Poisson Problem in mixed form with 2d and 3d with 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 automatic convergence estimation\n\ 5 and Hdiv elements.\n\n\n"; 6 7 /* 8 9 The mixed form of Poisson's equation, e.g. https://firedrakeproject.org/demos/poisson_mixed.py.html, is given 10 in the strong form by 11 \begin{align} 12 \vb{q} - \nabla u &= 0 \\ 13 \nabla \cdot \vb{q} &= f 14 \end{align} 15 where $u$ is the potential, as in the original problem, but we also discretize the gradient of potential $\vb{q}$, 16 or flux, directly. The weak form is then 17 \begin{align} 18 <t, \vb{q}> + <\nabla \cdot t, u> - <t_n, u>_\Gamma &= 0 \\ 19 <v, \nabla \cdot \vb{q}> &= <v, f> 20 \end{align} 21 For the original Poisson problem, the Dirichlet boundary forces an essential boundary condition on the potential space, 22 and the Neumann boundary gives a natural boundary condition in the weak form. In the mixed formulation, the Neumann 23 boundary gives an essential boundary condition on the flux space, $\vb{q} \cdot \vu{n} = h$, and the Dirichlet condition 24 becomes a natural condition in the weak form, <t_n, g>_\Gamma. 25 */ 26 27 #include <petscdmplex.h> 28 #include <petscsnes.h> 29 #include <petscds.h> 30 #include <petscconvest.h> 31 32 typedef enum { 33 SOL_LINEAR, 34 SOL_QUADRATIC, 35 SOL_QUARTIC, 36 SOL_QUARTIC_NEUMANN, 37 SOL_UNKNOWN, 38 NUM_SOLTYPE 39 } SolType; 40 const char *SolTypeNames[NUM_SOLTYPE + 3] = {"linear", "quadratic", "quartic", "quartic_neumann", "unknown", "SolType", "SOL_", NULL}; 41 42 typedef struct { 43 SolType solType; /* The type of exact solution */ 44 } AppCtx; 45 46 static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 47 PetscInt d; 48 for (d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d]; 49 } 50 51 /* 2D Dirichlet potential example 52 53 u = x 54 q = <1, 0> 55 f = 0 56 57 We will need a boundary integral of u over \Gamma. 58 */ 59 static PetscErrorCode linear_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 60 u[0] = x[0]; 61 return 0; 62 } 63 static PetscErrorCode linear_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 64 PetscInt c; 65 for (c = 0; c < Nc; ++c) u[c] = c ? 0.0 : 1.0; 66 return 0; 67 } 68 69 static void f0_linear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 70 f0[0] = 0.0; 71 f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0); 72 } 73 static void f0_bd_linear_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 74 PetscScalar potential; 75 PetscInt d; 76 77 linear_u(dim, t, x, dim, &potential, NULL); 78 for (d = 0; d < dim; ++d) f0[d] = -potential * n[d]; 79 } 80 81 /* 2D Dirichlet potential example 82 83 u = x^2 + y^2 84 q = <2x, 2y> 85 f = 4 86 87 We will need a boundary integral of u over \Gamma. 88 */ 89 static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 90 PetscInt d; 91 92 u[0] = 0.0; 93 for (d = 0; d < dim; ++d) u[0] += x[d] * x[d]; 94 return 0; 95 } 96 static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 97 PetscInt c; 98 for (c = 0; c < Nc; ++c) u[c] = 2.0 * x[c]; 99 return 0; 100 } 101 102 static void f0_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 103 f0[0] = -4.0; 104 f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0); 105 } 106 static void f0_bd_quadratic_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 107 PetscScalar potential; 108 PetscInt d; 109 110 quadratic_u(dim, t, x, dim, &potential, NULL); 111 for (d = 0; d < dim; ++d) f0[d] = -potential * n[d]; 112 } 113 114 /* 2D Dirichlet potential example 115 116 u = x (1-x) y (1-y) 117 q = <(1-2x) y (1-y), x (1-x) (1-2y)> 118 f = -y (1-y) - x (1-x) 119 120 u|_\Gamma = 0 so that the boundary integral vanishes 121 */ 122 static PetscErrorCode quartic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 123 PetscInt d; 124 125 u[0] = 1.0; 126 for (d = 0; d < dim; ++d) u[0] *= x[d] * (1.0 - x[d]); 127 return 0; 128 } 129 static PetscErrorCode quartic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 130 PetscInt c, d; 131 132 for (c = 0; c < Nc; ++c) { 133 u[c] = 1.0; 134 for (d = 0; d < dim; ++d) { 135 if (c == d) u[c] *= 1 - 2.0 * x[d]; 136 else u[c] *= x[d] * (1.0 - x[d]); 137 } 138 } 139 return 0; 140 } 141 142 static void f0_quartic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 143 PetscInt d; 144 f0[0] = 0.0; 145 for (d = 0; d < dim; ++d) f0[0] += 2.0 * x[d] * (1.0 - x[d]); 146 f0_u(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0); 147 } 148 149 /* 2D Dirichlet potential example 150 151 u = x (1-x) y (1-y) 152 q = <(1-2x) y (1-y), x (1-x) (1-2y)> 153 f = -y (1-y) - x (1-x) 154 155 du/dn_\Gamma = {(1-2x) y (1-y), x (1-x) (1-2y)} . n produces an essential condition on q 156 */ 157 158 static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) { 159 PetscInt c; 160 for (c = 0; c < dim; ++c) f0[c] = u[uOff[0] + c]; 161 } 162 163 /* <\nabla\cdot w, u> = <\nabla w, Iu> */ 164 static void f1_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) { 165 PetscInt c; 166 for (c = 0; c < dim; ++c) f1[c * dim + c] = u[uOff[1]]; 167 } 168 169 /* <t, q> */ 170 static void g0_qq(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) { 171 PetscInt c; 172 for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0; 173 } 174 175 /* <\nabla\cdot t, u> = <\nabla t, Iu> */ 176 static void g2_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) { 177 PetscInt d; 178 for (d = 0; d < dim; ++d) g2[d * dim + d] = 1.0; 179 } 180 181 /* <v, \nabla\cdot q> */ 182 static void g1_uq(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) { 183 PetscInt d; 184 for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; 185 } 186 187 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 188 PetscFunctionBeginUser; 189 options->solType = SOL_LINEAR; 190 PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX"); 191 PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex24.c", SolTypeNames, (PetscEnum)options->solType, (PetscEnum *)&options->solType, NULL)); 192 PetscOptionsEnd(); 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { 197 PetscFunctionBeginUser; 198 PetscCall(DMCreate(comm, dm)); 199 PetscCall(DMSetType(*dm, DMPLEX)); 200 PetscCall(PetscObjectSetName((PetscObject)*dm, "Example Mesh")); 201 PetscCall(DMSetApplicationContext(*dm, user)); 202 PetscCall(DMSetFromOptions(*dm)); 203 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 204 PetscFunctionReturn(0); 205 } 206 207 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) { 208 PetscDS ds; 209 DMLabel label; 210 PetscWeakForm wf; 211 const PetscInt id = 1; 212 PetscInt bd; 213 214 PetscFunctionBeginUser; 215 PetscCall(DMGetLabel(dm, "marker", &label)); 216 PetscCall(DMGetDS(dm, &ds)); 217 PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q)); 218 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL)); 219 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qu, NULL)); 220 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_uq, NULL, NULL)); 221 switch (user->solType) { 222 case SOL_LINEAR: 223 PetscCall(PetscDSSetResidual(ds, 1, f0_linear_u, NULL)); 224 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 225 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 226 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_linear_q, 0, NULL)); 227 PetscCall(PetscDSSetExactSolution(ds, 0, linear_q, user)); 228 PetscCall(PetscDSSetExactSolution(ds, 1, linear_u, user)); 229 break; 230 case SOL_QUADRATIC: 231 PetscCall(PetscDSSetResidual(ds, 1, f0_quadratic_u, NULL)); 232 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 233 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 234 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_quadratic_q, 0, NULL)); 235 PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user)); 236 PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user)); 237 break; 238 case SOL_QUARTIC: 239 PetscCall(PetscDSSetResidual(ds, 1, f0_quartic_u, NULL)); 240 PetscCall(PetscDSSetExactSolution(ds, 0, quartic_q, user)); 241 PetscCall(PetscDSSetExactSolution(ds, 1, quartic_u, user)); 242 break; 243 case SOL_QUARTIC_NEUMANN: 244 PetscCall(PetscDSSetResidual(ds, 1, f0_quartic_u, NULL)); 245 PetscCall(PetscDSSetExactSolution(ds, 0, quartic_q, user)); 246 PetscCall(PetscDSSetExactSolution(ds, 1, quartic_u, user)); 247 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "Flux condition", label, 1, &id, 0, 0, NULL, (void (*)(void))quartic_q, NULL, user, NULL)); 248 break; 249 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]); 250 } 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode SetupDiscretization(DM dm, PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) { 255 DM cdm = dm; 256 PetscFE feq, feu; 257 DMPolytopeType ct; 258 PetscBool simplex; 259 PetscInt dim, cStart; 260 261 PetscFunctionBeginUser; 262 PetscCall(DMGetDimension(dm, &dim)); 263 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 264 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 265 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; 266 /* Create finite element */ 267 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "field_", -1, &feq)); 268 PetscCall(PetscObjectSetName((PetscObject)feq, "field")); 269 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "potential_", -1, &feu)); 270 PetscCall(PetscObjectSetName((PetscObject)feu, "potential")); 271 PetscCall(PetscFECopyQuadrature(feq, feu)); 272 /* Set discretization and boundary conditions for each mesh */ 273 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)feq)); 274 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)feu)); 275 PetscCall(DMCreateDS(dm)); 276 PetscCall((*setup)(dm, user)); 277 while (cdm) { 278 PetscCall(DMCopyDisc(dm, cdm)); 279 PetscCall(DMGetCoarseDM(cdm, &cdm)); 280 } 281 PetscCall(PetscFEDestroy(&feq)); 282 PetscCall(PetscFEDestroy(&feu)); 283 PetscFunctionReturn(0); 284 } 285 286 int main(int argc, char **argv) { 287 DM dm; /* Problem specification */ 288 SNES snes; /* Nonlinear solver */ 289 Vec u; /* Solutions */ 290 AppCtx user; /* User-defined work context */ 291 292 PetscFunctionBeginUser; 293 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 294 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 295 /* Primal system */ 296 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 297 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 298 PetscCall(SNESSetDM(snes, dm)); 299 PetscCall(SetupDiscretization(dm, SetupPrimalProblem, &user)); 300 PetscCall(DMCreateGlobalVector(dm, &u)); 301 PetscCall(VecSet(u, 0.0)); 302 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 303 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 304 PetscCall(SNESSetFromOptions(snes)); 305 PetscCall(DMSNESCheckFromOptions(snes, u)); 306 PetscCall(SNESSolve(snes, NULL, u)); 307 PetscCall(SNESGetSolution(snes, &u)); 308 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 309 /* Cleanup */ 310 PetscCall(VecDestroy(&u)); 311 PetscCall(SNESDestroy(&snes)); 312 PetscCall(DMDestroy(&dm)); 313 PetscCall(PetscFinalize()); 314 return 0; 315 } 316 317 /*TEST 318 319 test: 320 suffix: 2d_bdm1_p0 321 requires: triangle 322 args: -sol_type linear \ 323 -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 1 \ 324 -dmsnes_check .001 -snes_error_if_not_converged \ 325 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 326 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 327 -fieldsplit_field_pc_type lu \ 328 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 329 test: 330 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.0, 1.0] 331 # Using -sol_type quadratic -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: [2.9, 1.0] 332 suffix: 2d_bdm1_p0_conv 333 requires: triangle 334 args: -sol_type quartic \ 335 -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \ 336 -snes_error_if_not_converged \ 337 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 338 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 339 -fieldsplit_field_pc_type lu \ 340 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 341 342 test: 343 # HDF5 output: -dm_view hdf5:${PETSC_DIR}/sol.h5 -potential_view hdf5:${PETSC_DIR}/sol.h5::append -exact_vec_view hdf5:${PETSC_DIR}/sol.h5::append 344 # VTK output: -potential_view vtk: -exact_vec_view vtk: 345 # VTU output: -potential_view vtk:multifield.vtu -exact_vec_view vtk:exact.vtu 346 suffix: 2d_q2_p0 347 requires: triangle 348 args: -sol_type linear -dm_plex_simplex 0 \ 349 -field_petscspace_degree 2 -dm_refine 0 \ 350 -dmsnes_check .001 -snes_error_if_not_converged \ 351 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 352 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 353 -fieldsplit_field_pc_type lu \ 354 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 355 test: 356 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [0.0057, 1.0] 357 suffix: 2d_q2_p0_conv 358 requires: triangle 359 args: -sol_type linear -dm_plex_simplex 0 \ 360 -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \ 361 -snes_error_if_not_converged \ 362 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 363 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 364 -fieldsplit_field_pc_type lu \ 365 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 366 test: 367 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: [-0.014, 0.0066] 368 suffix: 2d_q2_p0_neumann_conv 369 requires: triangle 370 args: -sol_type quartic_neumann -dm_plex_simplex 0 \ 371 -field_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate \ 372 -snes_error_if_not_converged \ 373 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 374 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 375 -fieldsplit_field_pc_type lu \ 376 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type svd 377 378 TEST*/ 379 380 /* These tests will be active once tensor P^- is working 381 382 test: 383 suffix: 2d_bdmq1_p0_0 384 requires: triangle 385 args: -dm_plex_simplex 0 -sol_type linear \ 386 -field_petscspace_poly_type pminus_hdiv -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate \ 387 -dmsnes_check .001 -snes_error_if_not_converged \ 388 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 389 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 390 -fieldsplit_field_pc_type lu \ 391 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 392 test: 393 suffix: 2d_bdmq1_p0_2 394 requires: triangle 395 args: -dm_plex_simplex 0 -sol_type quartic \ 396 -field_petscspace_poly_type_no pminus_hdiv -field_petscspace_degree 1 -field_petscdualspace_type bdm -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate \ 397 -dmsnes_check .001 -snes_error_if_not_converged \ 398 -ksp_rtol 1e-10 -ksp_error_if_not_converged \ 399 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition full \ 400 -fieldsplit_field_pc_type lu \ 401 -fieldsplit_potential_ksp_rtol 1e-10 -fieldsplit_potential_pc_type lu 402 403 */ 404