1 static char help[] = "Poisson Problem with a split domain.\n\ 2 We solve the Poisson problem in two halves of a domain.\n\ 3 In one half, we include an additional field.\n\n\n"; 4 5 #include <petscdmplex.h> 6 #include <petscsnes.h> 7 #include <petscds.h> 8 #include <petscconvest.h> 9 10 typedef struct { 11 PetscInt dummy; 12 } AppCtx; 13 14 static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 15 { 16 u[0] = x[0] * x[0] + x[1] * x[1]; 17 return 0; 18 } 19 20 static PetscErrorCode quad_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 21 { 22 u[0] = 2.0; 23 return 0; 24 } 25 26 static void f0_quad_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[]) 27 { 28 PetscInt d; 29 for (d = 0; d < dim; ++d) f0[0] += 2.0; 30 } 31 32 static void f0_quad_p(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[]) 33 { 34 f0[0] = u[uOff[1]] - 2.0; 35 } 36 37 static void f1_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 f1[]) 38 { 39 PetscInt d; 40 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 41 } 42 43 static void g3_uu(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 g3[]) 44 { 45 PetscInt d; 46 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 47 } 48 49 static void g0_pp(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[]) 50 { 51 g0[0] = 1.0; 52 } 53 54 static PetscErrorCode DivideDomain(DM dm, AppCtx *user) 55 { 56 DMLabel top, bottom; 57 PetscReal low[3], high[3], midy; 58 PetscInt cStart, cEnd, c; 59 60 PetscFunctionBeginUser; 61 PetscCall(DMCreateLabel(dm, "top")); 62 PetscCall(DMCreateLabel(dm, "bottom")); 63 PetscCall(DMGetLabel(dm, "top", &top)); 64 PetscCall(DMGetLabel(dm, "bottom", &bottom)); 65 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 66 PetscCall(DMGetBoundingBox(dm, low, high)); 67 midy = 0.5 * (high[1] - low[1]); 68 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 69 for (c = cStart; c < cEnd; ++c) { 70 PetscReal centroid[3]; 71 72 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 73 if (centroid[1] > midy) PetscCall(DMLabelSetValue(top, c, 1)); 74 else PetscCall(DMLabelSetValue(bottom, c, 1)); 75 } 76 PetscCall(DMPlexLabelComplete(dm, top)); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 81 { 82 PetscFunctionBeginUser; 83 PetscCall(DMCreate(comm, dm)); 84 PetscCall(DMSetType(*dm, DMPLEX)); 85 PetscCall(DMSetFromOptions(*dm)); 86 PetscCall(DivideDomain(*dm, user)); 87 PetscCall(DMSetApplicationContext(*dm, user)); 88 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 93 { 94 PetscDS ds; 95 PetscWeakForm wf; 96 DMLabel label; 97 const PetscInt id = 1; 98 99 PetscFunctionBeginUser; 100 PetscCall(DMGetRegionNumDS(dm, 0, &label, NULL, &ds)); 101 PetscCall(PetscDSGetWeakForm(ds, &wf)); 102 PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u)); 103 PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu)); 104 PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user)); 105 PetscCall(DMGetRegionNumDS(dm, 1, &label, NULL, &ds)); 106 PetscCall(PetscDSGetWeakForm(ds, &wf)); 107 PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u)); 108 PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu)); 109 PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, 0, f0_quad_p, 0, NULL)); 110 PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL)); 111 PetscCall(PetscDSSetExactSolution(ds, 0, quad_u, user)); 112 PetscCall(PetscDSSetExactSolution(ds, 1, quad_p, user)); 113 PetscCall(DMGetLabel(dm, "marker", &label)); 114 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL)); 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 119 { 120 DM cdm = dm; 121 DMLabel top; 122 PetscFE fe, feTop; 123 PetscQuadrature q; 124 PetscInt dim; 125 PetscBool simplex; 126 const char *nameTop = "pressure"; 127 char prefix[PETSC_MAX_PATH_LEN]; 128 129 PetscFunctionBeginUser; 130 PetscCall(DMGetDimension(dm, &dim)); 131 PetscCall(DMPlexIsSimplex(dm, &simplex)); 132 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 133 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 134 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 135 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 136 PetscCall(PetscFEGetQuadrature(fe, &q)); 137 PetscCall(PetscFEDestroy(&fe)); 138 PetscCall(DMGetLabel(dm, "top", &top)); 139 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop)); 140 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &feTop)); 141 PetscCall(PetscObjectSetName((PetscObject)feTop, nameTop)); 142 PetscCall(PetscFESetQuadrature(feTop, q)); 143 PetscCall(DMSetField(dm, 1, top, (PetscObject)feTop)); 144 PetscCall(PetscFEDestroy(&feTop)); 145 PetscCall(DMCreateDS(dm)); 146 PetscCall((*setup)(dm, user)); 147 while (cdm) { 148 PetscCall(DMCopyDisc(dm, cdm)); 149 PetscCall(DMGetCoarseDM(cdm, &cdm)); 150 } 151 PetscFunctionReturn(0); 152 } 153 154 int main(int argc, char **argv) 155 { 156 DM dm; /* Problem specification */ 157 SNES snes; /* Nonlinear solver */ 158 Vec u; /* Solutions */ 159 AppCtx user; /* User-defined work context */ 160 161 PetscFunctionBeginUser; 162 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 163 /* Primal system */ 164 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 165 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 166 PetscCall(SNESSetDM(snes, dm)); 167 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 168 PetscCall(DMCreateGlobalVector(dm, &u)); 169 PetscCall(VecSet(u, 0.0)); 170 PetscCall(PetscObjectSetName((PetscObject)u, "solution")); 171 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 172 PetscCall(SNESSetFromOptions(snes)); 173 PetscCall(DMSNESCheckFromOptions(snes, u)); 174 PetscCall(SNESSolve(snes, NULL, u)); 175 PetscCall(SNESGetSolution(snes, &u)); 176 PetscCall(VecViewFromOptions(u, NULL, "-sol_view")); 177 /* Cleanup */ 178 PetscCall(VecDestroy(&u)); 179 PetscCall(SNESDestroy(&snes)); 180 PetscCall(DMDestroy(&dm)); 181 PetscCall(PetscFinalize()); 182 return 0; 183 } 184 185 /*TEST 186 187 test: 188 suffix: 2d_p1_0 189 requires: triangle 190 args: -potential_petscspace_degree 2 -pressure_top_petscspace_degree 0 -dm_refine 0 -dmsnes_check 191 192 test: 193 suffix: 2d_p1_1 194 requires: triangle 195 args: -potential_petscspace_degree 1 -pressure_top_petscspace_degree 0 -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate 196 197 TEST*/ 198