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