1 static char help[] = "Tests implementation of PetscSpace_Sum by solving the Poisson equations using a PetscSpace_Poly and a PetscSpace_Sum and checking that \ 2 solutions agree up to machine precision.\n\n"; 3 4 #include <petscdmplex.h> 5 #include <petscds.h> 6 #include <petscfe.h> 7 #include <petscsnes.h> 8 /* We are solving the system of equations: 9 * \vec{u} = -\grad{p} 10 * \div{u} = f 11 */ 12 13 /* Exact solutions for linear velocity 14 \vec{u} = \vec{x}; 15 p = -0.5*(\vec{x} \cdot \vec{x}); 16 */ 17 static PetscErrorCode linear_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 18 { 19 PetscInt c; 20 21 for (c = 0; c < Nc; ++c) u[c] = x[c]; 22 return PETSC_SUCCESS; 23 } 24 25 static PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 26 { 27 PetscInt d; 28 29 u[0] = 0.; 30 for (d = 0; d < dim; ++d) u[0] += -0.5 * x[d] * x[d]; 31 return PETSC_SUCCESS; 32 } 33 34 static PetscErrorCode linear_divu(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 35 { 36 u[0] = dim; 37 return PETSC_SUCCESS; 38 } 39 40 /* fx_v are the residual functions for the equation \vec{u} = \grad{p}. f0_v is the term <v,u>.*/ 41 static void f0_v(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[]) 42 { 43 PetscInt i; 44 45 for (i = 0; i < dim; ++i) f0[i] = u[uOff[0] + i]; 46 } 47 48 /* f1_v is the term <v,-\grad{p}> but we integrate by parts to get <\grad{v}, -p*I> */ 49 static void f1_v(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[]) 50 { 51 PetscInt c; 52 53 for (c = 0; c < dim; ++c) { 54 PetscInt d; 55 56 for (d = 0; d < dim; ++d) f1[c * dim + d] = (c == d) ? -u[uOff[1]] : 0; 57 } 58 } 59 60 /* Residual function for enforcing \div{u} = f. */ 61 static void f0_q_linear(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[]) 62 { 63 PetscScalar rhs, divu = 0; 64 PetscInt i; 65 66 (void)linear_divu(dim, t, x, dim, &rhs, NULL); 67 for (i = 0; i < dim; ++i) divu += u_x[uOff_x[0] + i * dim + i]; 68 f0[0] = divu - rhs; 69 } 70 71 /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */ 72 static void f0_bd_u_linear(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[]) 73 { 74 PetscScalar pressure; 75 PetscInt d; 76 77 (void)linear_p(dim, t, x, dim, &pressure, NULL); 78 for (d = 0; d < dim; ++d) f0[d] = pressure * n[d]; 79 } 80 81 /* gx_yz are the jacobian functions obtained by taking the derivative of the y residual w.r.t z*/ 82 static void g0_vu(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[]) 83 { 84 PetscInt c; 85 86 for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0; 87 } 88 89 static void g1_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 g1[]) 90 { 91 PetscInt c; 92 93 for (c = 0; c < dim; ++c) g1[c * dim + c] = 1.0; 94 } 95 96 static void g2_vp(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[]) 97 { 98 PetscInt c; 99 100 for (c = 0; c < dim; ++c) g2[c * dim + c] = -1.0; 101 } 102 103 typedef struct { 104 PetscInt dummy; 105 } UserCtx; 106 107 static PetscErrorCode CreateMesh(MPI_Comm comm, UserCtx *user, DM *mesh) 108 { 109 PetscFunctionBegin; 110 PetscCall(DMCreate(comm, mesh)); 111 PetscCall(DMSetType(*mesh, DMPLEX)); 112 PetscCall(DMSetFromOptions(*mesh)); 113 PetscCall(DMSetApplicationContext(*mesh, user)); 114 PetscCall(DMViewFromOptions(*mesh, NULL, "-dm_view")); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 /* Setup the system of equations that we wish to solve */ 119 static PetscErrorCode SetupProblem(DM dm, UserCtx *user) 120 { 121 PetscDS ds; 122 DMLabel label; 123 PetscWeakForm wf; 124 const PetscInt id = 1; 125 PetscInt bd; 126 127 PetscFunctionBegin; 128 PetscCall(DMGetDS(dm, &ds)); 129 /* All of these are independent of the user's choice of solution */ 130 PetscCall(PetscDSSetResidual(ds, 0, f0_v, f1_v)); 131 PetscCall(PetscDSSetResidual(ds, 1, f0_q_linear, NULL)); 132 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_vu, NULL, NULL, NULL)); 133 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_vp, NULL)); 134 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_qu, NULL, NULL)); 135 136 PetscCall(DMGetLabel(dm, "marker", &label)); 137 PetscCall(PetscDSAddBoundary(ds, DM_BC_NATURAL, "Boundary Integral", label, 1, &id, 0, 0, NULL, (PetscFortranCallbackFn *)NULL, NULL, user, &bd)); 138 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 139 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_linear, 0, NULL)); 140 141 PetscCall(PetscDSSetExactSolution(ds, 0, linear_u, NULL)); 142 PetscCall(PetscDSSetExactSolution(ds, 1, linear_p, NULL)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /* Create the finite element spaces we will use for this system */ 147 static PetscErrorCode SetupDiscretization(DM mesh, DM mesh_sum, PetscErrorCode (*setup)(DM, UserCtx *), UserCtx *user) 148 { 149 DM cdm = mesh, cdm_sum = mesh_sum; 150 PetscDS ds; 151 PetscFE u, divu, u_sum, divu_sum; 152 PetscInt dim; 153 PetscBool simplex; 154 155 PetscFunctionBegin; 156 PetscCall(DMGetDimension(mesh, &dim)); 157 PetscCall(DMPlexIsSimplex(mesh, &simplex)); 158 159 { 160 PetscBool force; 161 // Turn off automatic quadrature selection as a test 162 PetscCall(DMGetDS(mesh_sum, &ds)); 163 PetscCall(PetscDSGetForceQuad(ds, &force)); 164 if (force) PetscCall(PetscDSSetForceQuad(ds, PETSC_FALSE)); 165 } 166 167 /* Create FE objects and give them names so that options can be set from 168 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */ 169 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, dim, simplex, "velocity_", -1, &u)); 170 PetscCall(PetscObjectSetName((PetscObject)u, "velocity")); 171 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum), dim, dim, simplex, "velocity_sum_", -1, &u_sum)); 172 PetscCall(PetscObjectSetName((PetscObject)u_sum, "velocity_sum")); 173 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "divu_", -1, &divu)); 174 PetscCall(PetscObjectSetName((PetscObject)divu, "divu")); 175 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum), dim, 1, simplex, "divu_sum_", -1, &divu_sum)); 176 PetscCall(PetscObjectSetName((PetscObject)divu_sum, "divu_sum")); 177 178 PetscCall(PetscFECopyQuadrature(u, divu)); 179 PetscCall(PetscFECopyQuadrature(u_sum, divu_sum)); 180 181 /* Associate the FE objects with the mesh and setup the system */ 182 PetscCall(DMSetField(mesh, 0, NULL, (PetscObject)u)); 183 PetscCall(DMSetField(mesh, 1, NULL, (PetscObject)divu)); 184 PetscCall(DMCreateDS(mesh)); 185 PetscCall((*setup)(mesh, user)); 186 187 PetscCall(DMSetField(mesh_sum, 0, NULL, (PetscObject)u_sum)); 188 PetscCall(DMSetField(mesh_sum, 1, NULL, (PetscObject)divu_sum)); 189 PetscCall(DMCreateDS(mesh_sum)); 190 PetscCall((*setup)(mesh_sum, user)); 191 192 while (cdm) { 193 PetscCall(DMCopyDisc(mesh, cdm)); 194 PetscCall(DMGetCoarseDM(cdm, &cdm)); 195 } 196 197 while (cdm_sum) { 198 PetscCall(DMCopyDisc(mesh_sum, cdm_sum)); 199 PetscCall(DMGetCoarseDM(cdm_sum, &cdm_sum)); 200 } 201 202 /* The Mesh now owns the fields, so we can destroy the FEs created in this 203 * function */ 204 PetscCall(PetscFEDestroy(&u)); 205 PetscCall(PetscFEDestroy(&divu)); 206 PetscCall(PetscFEDestroy(&u_sum)); 207 PetscCall(PetscFEDestroy(&divu_sum)); 208 PetscCall(DMDestroy(&cdm)); 209 PetscCall(DMDestroy(&cdm_sum)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 int main(int argc, char **argv) 214 { 215 UserCtx user; 216 DM dm, dm_sum; 217 SNES snes, snes_sum; 218 Vec u, u_sum; 219 PetscReal errNorm; 220 const PetscReal errTol = PETSC_SMALL; 221 222 PetscFunctionBeginUser; 223 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 224 225 /* Set up a snes for the standard approach, one space with 2 components */ 226 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 227 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 228 PetscCall(SNESSetDM(snes, dm)); 229 230 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */ 231 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_sum)); 232 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm_sum)); 233 PetscCall(SNESSetDM(snes_sum, dm_sum)); 234 PetscCall(SetupDiscretization(dm, dm_sum, SetupProblem, &user)); 235 236 /* Set up and solve the system using standard approach. */ 237 PetscCall(DMCreateGlobalVector(dm, &u)); 238 PetscCall(VecSet(u, 0.0)); 239 PetscCall(PetscObjectSetName((PetscObject)u, "solution")); 240 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user)); 241 PetscCall(SNESSetFromOptions(snes)); 242 PetscCall(DMSNESCheckFromOptions(snes, u)); 243 PetscCall(SNESSolve(snes, NULL, u)); 244 PetscCall(SNESGetSolution(snes, &u)); 245 PetscCall(VecViewFromOptions(u, NULL, "-solution_view")); 246 247 /* Set up and solve the sum space system */ 248 PetscCall(DMCreateGlobalVector(dm_sum, &u_sum)); 249 PetscCall(VecSet(u_sum, 0.0)); 250 PetscCall(PetscObjectSetName((PetscObject)u_sum, "solution_sum")); 251 PetscCall(DMPlexSetSNESLocalFEM(dm_sum, PETSC_FALSE, &user)); 252 PetscCall(SNESSetFromOptions(snes_sum)); 253 PetscCall(DMSNESCheckFromOptions(snes_sum, u_sum)); 254 PetscCall(SNESSolve(snes_sum, NULL, u_sum)); 255 PetscCall(SNESGetSolution(snes_sum, &u_sum)); 256 PetscCall(VecViewFromOptions(u_sum, NULL, "-solution_sum_view")); 257 258 /* Check if standard solution and sum space solution match to machine precision */ 259 PetscCall(VecAXPY(u_sum, -1, u)); 260 PetscCall(VecNorm(u_sum, NORM_2, &errNorm)); 261 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Sum space provides the same solution as a regular space: %s", (errNorm < errTol) ? "true" : "false")); 262 263 /* Cleanup */ 264 PetscCall(VecDestroy(&u_sum)); 265 PetscCall(VecDestroy(&u)); 266 PetscCall(SNESDestroy(&snes_sum)); 267 PetscCall(SNESDestroy(&snes)); 268 PetscCall(DMDestroy(&dm_sum)); 269 PetscCall(DMDestroy(&dm)); 270 PetscCall(PetscFinalize()); 271 return 0; 272 } 273 274 /*TEST 275 test: 276 suffix: 2d_lagrange 277 requires: triangle 278 args: -velocity_petscspace_degree 1 \ 279 -velocity_petscspace_type poly \ 280 -velocity_petscspace_components 2\ 281 -velocity_petscdualspace_type lagrange \ 282 -divu_petscspace_degree 0 \ 283 -divu_petscspace_type poly \ 284 -divu_petscdualspace_lagrange_continuity false \ 285 -velocity_sum_petscfe_default_quadrature_order 1 \ 286 -velocity_sum_petscspace_degree 1 \ 287 -velocity_sum_petscspace_type sum \ 288 -velocity_sum_petscspace_variables 2 \ 289 -velocity_sum_petscspace_components 2 \ 290 -velocity_sum_petscspace_sum_spaces 2 \ 291 -velocity_sum_petscspace_sum_concatenate true \ 292 -velocity_sum_petscdualspace_type lagrange \ 293 -velocity_sum_sumcomp_0_petscspace_type poly \ 294 -velocity_sum_sumcomp_0_petscspace_degree 1 \ 295 -velocity_sum_sumcomp_0_petscspace_variables 2 \ 296 -velocity_sum_sumcomp_0_petscspace_components 1 \ 297 -velocity_sum_sumcomp_1_petscspace_type poly \ 298 -velocity_sum_sumcomp_1_petscspace_degree 1 \ 299 -velocity_sum_sumcomp_1_petscspace_variables 2 \ 300 -velocity_sum_sumcomp_1_petscspace_components 1 \ 301 -divu_sum_petscspace_degree 0 \ 302 -divu_sum_petscspace_type sum \ 303 -divu_sum_petscspace_variables 2 \ 304 -divu_sum_petscspace_components 1 \ 305 -divu_sum_petscspace_sum_spaces 1 \ 306 -divu_sum_petscspace_sum_concatenate true \ 307 -divu_sum_petscdualspace_lagrange_continuity false \ 308 -divu_sum_sumcomp_0_petscspace_type poly \ 309 -divu_sum_sumcomp_0_petscspace_degree 0 \ 310 -divu_sum_sumcomp_0_petscspace_variables 2 \ 311 -divu_sum_sumcomp_0_petscspace_components 1 \ 312 -dm_refine 0 \ 313 -snes_error_if_not_converged \ 314 -ksp_rtol 1e-10 \ 315 -ksp_error_if_not_converged \ 316 -pc_type fieldsplit\ 317 -pc_fieldsplit_type schur\ 318 -divu_sum_petscdualspace_lagrange_continuity false \ 319 -pc_fieldsplit_schur_precondition full 320 TEST*/ 321