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, void *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, void *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:   PetscInt d;
 29:   for (d = 0; d < dim; ++d) f0[0] += 2.0;
 30: }

 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: }

 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: }

 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: }

 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: }

 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;

 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];

 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(PETSC_SUCCESS);
 78: }

 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(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
 93: {
 94:   PetscDS        ds;
 95:   PetscWeakForm  wf;
 96:   DMLabel        label;
 97:   const PetscInt id = 1;

 99:   PetscFunctionBeginUser;
100:   PetscCall(DMGetRegionNumDS(dm, 0, &label, NULL, &ds, NULL));
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, NULL));
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(PETSC_SUCCESS);
116: }

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];

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(PETSC_SUCCESS);
152: }

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 */

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: }

185: /*TEST

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

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

197: TEST*/