Actual source code: ex15.c

  1: static char help[] = "Mixed element discretization of the Poisson equation.\n\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscdmswarm.h>
  5: #include <petscds.h>
  6: #include <petscsnes.h>
  7: #include <petscconvest.h>
  8: #include <petscbag.h>

 10: /*
 11: The Poisson equation

 13:   -\Delta\phi = f

 15: can be rewritten in first order form

 17:   q - \nabla\phi  &= 0
 18:   -\nabla \cdot q &= f
 19: */

 21: typedef enum {
 22:   SIGMA,
 23:   NUM_CONSTANTS
 24: } ConstantType;
 25: typedef struct {
 26:   PetscReal sigma; /* Nondimensional charge per length in x */
 27: } Parameter;

 29: typedef enum {
 30:   SOL_CONST,
 31:   SOL_LINEAR,
 32:   SOL_QUADRATIC,
 33:   SOL_TRIG,
 34:   SOL_TRIGX,
 35:   SOL_PARTICLES,
 36:   NUM_SOL_TYPES
 37: } SolType;
 38: static const char *solTypes[] = {"const", "linear", "quadratic", "trig", "trigx", "particles"};

 40: typedef struct {
 41:   SolType   solType; /* MMS solution type */
 42:   PetscBag  bag;     /* Problem parameters */
 43:   PetscBool particleRHS;
 44:   PetscInt  Np;
 45: } AppCtx;

 47: /* SOLUTION CONST: \phi = 1, q = 0, f = 0 */
 48: static PetscErrorCode const_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 49: {
 50:   *u = 1.0;
 51:   return PETSC_SUCCESS;
 52: }

 54: static PetscErrorCode const_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 55: {
 56:   for (PetscInt d = 0; d < dim; ++d) u[d] = 0.0;
 57:   return PETSC_SUCCESS;
 58: }

 60: /* SOLUTION LINEAR: \phi = 2y, q = <0, 2>, f = 0 */
 61: static PetscErrorCode linear_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 62: {
 63:   u[0] = 2. * x[1];
 64:   return PETSC_SUCCESS;
 65: }

 67: static PetscErrorCode linear_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 68: {
 69:   u[0] = 0.;
 70:   u[1] = 2.;
 71:   return PETSC_SUCCESS;
 72: }

 74: /* SOLUTION QUADRATIC: \phi = x (2\pi - x) + (1 + y) (1 - y), q = <2\pi - 2 x, - 2 y> = <2\pi, 0> - 2 x, f = -4 */
 75: static PetscErrorCode quadratic_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 76: {
 77:   u[0] = x[0] * (6.283185307179586 - x[0]) + (1. + x[1]) * (1. - x[1]);
 78:   return PETSC_SUCCESS;
 79: }

 81: static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 82: {
 83:   u[0] = 6.283185307179586 - 2. * x[0];
 84:   u[1] = -2. * x[1];
 85:   return PETSC_SUCCESS;
 86: }

 88: static PetscErrorCode quadratic_q_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 89: {
 90:   u[0] = x[1] > 0. ? -2. * x[1] : 2. * x[1];
 91:   return PETSC_SUCCESS;
 92: }

 94: static void f0_quadratic_phi(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[])
 95: {
 96:   for (PetscInt d = 0; d < dim; ++d) f0[0] -= -2.0;
 97: }

 99: /* SOLUTION TRIG: \phi = sin(x) + (1/3 - y^2), q = <cos(x), -2 y>, f = sin(x) + 2 */
100: static PetscErrorCode trig_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
101: {
102:   u[0] = PetscSinReal(x[0]) + (1. / 3. - x[1] * x[1]);
103:   return PETSC_SUCCESS;
104: }

106: static PetscErrorCode trig_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
107: {
108:   u[0] = PetscCosReal(x[0]);
109:   u[1] = -2. * x[1];
110:   return PETSC_SUCCESS;
111: }

113: static PetscErrorCode trig_q_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
114: {
115:   u[0] = x[1] > 0. ? -2. * x[1] : 2. * x[1];
116:   return PETSC_SUCCESS;
117: }

119: static void f0_trig_phi(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[])
120: {
121:   f0[0] += PetscSinReal(x[0]) + 2.;
122: }

124: /* SOLUTION TRIGX: \phi = sin(x), q = <cos(x), 0>, f = sin(x) */
125: static PetscErrorCode trigx_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
126: {
127:   u[0] = PetscSinReal(x[0]);
128:   return PETSC_SUCCESS;
129: }

131: static PetscErrorCode trigx_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
132: {
133:   u[0] = PetscCosReal(x[0]);
134:   u[1] = 0.;
135:   return PETSC_SUCCESS;
136: }

138: static PetscErrorCode trigx_q_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
139: {
140:   u[0] = x[1] > 0. ? -2. * x[1] : 2. * x[1];
141:   return PETSC_SUCCESS;
142: }

144: static void f0_trigx_phi(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[])
145: {
146:   f0[0] += PetscSinReal(x[0]);
147: }

149: 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[])
150: {
151:   for (PetscInt d = 0; d < dim; ++d) f0[d] += u[uOff[0] + d];
152: }

154: 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[])
155: {
156:   for (PetscInt d = 0; d < dim; ++d) f1[d * dim + d] = u[uOff[1]];
157: }

159: static void f0_phi(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[])
160: {
161:   for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
162: }

164: static void f0_phi_backgroundCharge(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[])
165: {
166:   f0[0] += constants[SIGMA];
167:   for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[uOff_x[0] + d * dim + d];
168: }

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: {
172:   for (PetscInt d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
173: }

175: static void g2_qphi(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[])
176: {
177:   for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = 1.0;
178: }

180: static void g1_phiq(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[])
181: {
182:   for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
183: }

185: /* SOLUTION PARTICLES: \phi = sigma, q = <cos(x), 0>, f = sin(x) */
186: static PetscErrorCode particles_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
187: {
188:   u[0] = 0.0795775;
189:   return PETSC_SUCCESS;
190: }

192: static PetscErrorCode particles_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
193: {
194:   u[0] = 0.;
195:   u[1] = 0.;
196:   return PETSC_SUCCESS;
197: }

199: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
200: {
201:   PetscInt sol;

203:   PetscFunctionBeginUser;
204:   options->solType     = SOL_CONST;
205:   options->particleRHS = PETSC_FALSE;
206:   options->Np          = 100;

208:   PetscOptionsBegin(comm, "", "Mixed Poisson Options", "DMPLEX");
209:   PetscCall(PetscOptionsBool("-particleRHS", "Flag to user particle RHS and background charge", "ex9.c", options->particleRHS, &options->particleRHS, NULL));
210:   sol = options->solType;
211:   PetscCall(PetscOptionsEList("-sol_type", "The MMS solution type", "ex12.c", solTypes, NUM_SOL_TYPES, solTypes[sol], &sol, NULL));
212:   options->solType = (SolType)sol;
213:   PetscOptionsEnd();
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
218: {
219:   PetscFunctionBeginUser;
220:   PetscCall(DMCreate(comm, dm));
221:   PetscCall(DMSetType(*dm, DMPLEX));
222:   PetscCall(DMSetFromOptions(*dm));
223:   PetscCall(DMSetApplicationContext(*dm, user));
224:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
229: {
230:   PetscDS        ds;
231:   PetscWeakForm  wf;
232:   DMLabel        label;
233:   const PetscInt id = 1;

235:   PetscFunctionBeginUser;
236:   PetscCall(DMGetDS(dm, &ds));
237:   PetscCall(PetscDSGetWeakForm(ds, &wf));
238:   PetscCall(DMGetLabel(dm, "marker", &label));
239:   PetscCall(PetscDSSetResidual(ds, 0, f0_q, f1_q));
240:   if (user->particleRHS) {
241:     PetscCall(PetscDSSetResidual(ds, 1, f0_phi_backgroundCharge, NULL));
242:   } else {
243:     PetscCall(PetscDSSetResidual(ds, 1, f0_phi, NULL));
244:   }
245:   PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_qq, NULL, NULL, NULL));
246:   PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_qphi, NULL));
247:   PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_phiq, NULL, NULL));
248:   switch (user->solType) {
249:   case SOL_CONST:
250:     PetscCall(PetscDSSetExactSolution(ds, 0, const_q, user));
251:     PetscCall(PetscDSSetExactSolution(ds, 1, const_phi, user));
252:     break;
253:   case SOL_LINEAR:
254:     PetscCall(PetscDSSetExactSolution(ds, 0, linear_q, user));
255:     PetscCall(PetscDSSetExactSolution(ds, 1, linear_phi, user));
256:     break;
257:   case SOL_QUADRATIC:
258:     PetscCall(PetscWeakFormAddResidual(wf, NULL, 0, 1, 0, f0_quadratic_phi, NULL));
259:     PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user));
260:     PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_phi, user));
261:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_q_bc, NULL, user, NULL));
262:     break;
263:   case SOL_TRIG:
264:     PetscCall(PetscWeakFormAddResidual(wf, NULL, 0, 1, 0, f0_trig_phi, NULL));
265:     PetscCall(PetscDSSetExactSolution(ds, 0, trig_q, user));
266:     PetscCall(PetscDSSetExactSolution(ds, 1, trig_phi, user));
267:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_q_bc, NULL, user, NULL));
268:     break;
269:   case SOL_TRIGX:
270:     PetscCall(PetscWeakFormAddResidual(wf, NULL, 0, 1, 0, f0_trigx_phi, NULL));
271:     PetscCall(PetscDSSetExactSolution(ds, 0, trigx_q, user));
272:     PetscCall(PetscDSSetExactSolution(ds, 1, trigx_phi, user));
273:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trigx_q_bc, NULL, user, NULL));
274:     break;
275:   case SOL_PARTICLES:
276:     PetscCall(PetscDSSetExactSolution(ds, 0, particles_q, user));
277:     PetscCall(PetscDSSetExactSolution(ds, 1, particles_phi, user));
278:     break;
279:   default:
280:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid solution type: %d", user->solType);
281:   }

283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: static PetscErrorCode SetupDiscretization(DM dm, PetscInt Nf, const char *names[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
287: {
288:   DM             cdm = dm;
289:   PetscFE        fe;
290:   DMPolytopeType ct;
291:   PetscInt       dim, cStart;
292:   char           prefix[PETSC_MAX_PATH_LEN];

294:   PetscFunctionBeginUser;
295:   PetscCall(DMGetDimension(dm, &dim));
296:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
297:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
298:   for (PetscInt f = 0; f < Nf; ++f) {
299:     PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", names[f]));
300:     PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, -1, &fe));
301:     PetscCall(PetscObjectSetName((PetscObject)fe, names[f]));
302:     if (f > 0) {
303:       PetscFE fe0;

305:       PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe0));
306:       PetscCall(PetscFECopyQuadrature(fe0, fe));
307:     }
308:     PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
309:     PetscCall(PetscFEDestroy(&fe));
310:   }
311:   PetscCall(DMCreateDS(dm));
312:   PetscCall((*setup)(dm, user));
313:   while (cdm) {
314:     PetscCall(DMCopyDisc(dm, cdm));
315:     PetscCall(DMGetCoarseDM(cdm, &cdm));
316:   }
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: static PetscErrorCode InitializeParticlesAndWeights(DM sw, AppCtx *user)
321: {
322:   DM           dm;
323:   PetscScalar *weight;
324:   PetscInt     Np, Npc, p, dim, c, cStart, cEnd, q, *cellid;
325:   PetscReal    weightsum = 0.0;
326:   PetscMPIInt  size, rank;
327:   Parameter   *param;

329:   PetscFunctionBegin;
330:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
331:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
332:   PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
333:   PetscCall(PetscOptionsInt("-dm_swarm_num_particles", "The target number of particles", "", user->Np, &user->Np, NULL));
334:   PetscOptionsEnd();

336:   Np = user->Np;
337:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Np = %" PetscInt_FMT "\n", Np));
338:   PetscCall(DMGetDimension(sw, &dim));
339:   PetscCall(DMSwarmGetCellDM(sw, &dm));
340:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));

342:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
343:   PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));

345:   Npc = Np / (cEnd - cStart);
346:   PetscCall(DMSwarmGetField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
347:   for (c = 0, p = 0; c < cEnd - cStart; ++c) {
348:     for (q = 0; q < Npc; ++q, ++p) cellid[p] = c;
349:   }
350:   PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));

352:   PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
353:   PetscCall(DMSwarmSortGetAccess(sw));
354:   for (p = 0; p < Np; ++p) {
355:     weight[p] = 1.0 / Np;
356:     weightsum += PetscRealPart(weight[p]);
357:   }

359:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "weightsum = %1.10f\n", (double)weightsum));
360:   PetscCall(DMSwarmSortRestoreAccess(sw));
361:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }

365: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
366: {
367:   PetscInt dim;

369:   PetscFunctionBeginUser;
370:   PetscCall(DMGetDimension(dm, &dim));
371:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
372:   PetscCall(DMSetType(*sw, DMSWARM));
373:   PetscCall(DMSetDimension(*sw, dim));
374:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
375:   PetscCall(DMSwarmSetCellDM(*sw, dm));
376:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));

378:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));

380:   PetscCall(InitializeParticlesAndWeights(*sw, user));

382:   PetscCall(DMSetFromOptions(*sw));
383:   PetscCall(DMSetApplicationContext(*sw, user));
384:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
385:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));

387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
391: {
392:   PetscBag   bag;
393:   Parameter *p;

395:   PetscFunctionBeginUser;
396:   /* setup PETSc parameter bag */
397:   PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
398:   PetscCall(PetscBagSetName(ctx->bag, "par", "Parameters"));
399:   bag = ctx->bag;
400:   PetscCall(PetscBagRegisterScalar(bag, &p->sigma, 1.0, "sigma", "Charge per unit area, C/m^3"));
401:   PetscCall(PetscBagSetFromOptions(bag));
402:   {
403:     PetscViewer       viewer;
404:     PetscViewerFormat format;
405:     PetscBool         flg;

407:     PetscCall(PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
408:     if (flg) {
409:       PetscCall(PetscViewerPushFormat(viewer, format));
410:       PetscCall(PetscBagView(bag, viewer));
411:       PetscCall(PetscViewerFlush(viewer));
412:       PetscCall(PetscViewerPopFormat(viewer));
413:       PetscCall(PetscViewerDestroy(&viewer));
414:     }
415:   }
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: static PetscErrorCode InitializeConstants(DM sw, AppCtx *user)
420: {
421:   DM         dm;
422:   PetscReal *weight, totalCharge, totalWeight = 0., gmin[3], gmax[3];
423:   PetscInt   Np, p, dim;

425:   PetscFunctionBegin;
426:   PetscCall(DMSwarmGetCellDM(sw, &dm));
427:   PetscCall(DMGetDimension(sw, &dim));
428:   PetscCall(DMSwarmGetLocalSize(sw, &Np));
429:   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
430:   PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&weight));
431:   for (p = 0; p < Np; ++p) totalWeight += weight[p];
432:   totalCharge = -1.0 * totalWeight;
433:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
434:   {
435:     Parameter *param;
436:     PetscReal  Area;

438:     PetscCall(PetscBagGetData(user->bag, (void **)&param));
439:     switch (dim) {
440:     case 1:
441:       Area = (gmax[0] - gmin[0]);
442:       break;
443:     case 2:
444:       Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]);
445:       break;
446:     case 3:
447:       Area = (gmax[0] - gmin[0]) * (gmax[1] - gmin[1]) * (gmax[2] - gmin[2]);
448:       break;
449:     default:
450:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", dim);
451:     }
452:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim = %" PetscInt_FMT "\ttotalWeight = %f\ttotalCharge = %f, Total Area = %f\n", dim, (double)totalWeight, (double)totalCharge, (double)Area));
453:     param->sigma = PetscAbsReal(totalCharge / (Area));

455:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "sigma: %g\n", (double)param->sigma));
456:   }
457:   /* Setup Constants */
458:   {
459:     PetscDS    ds;
460:     Parameter *param;
461:     PetscCall(PetscBagGetData(user->bag, (void **)&param));
462:     PetscScalar constants[NUM_CONSTANTS];
463:     constants[SIGMA] = param->sigma;
464:     PetscCall(DMGetDS(dm, &ds));
465:     PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants));
466:   }
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: int main(int argc, char **argv)
471: {
472:   DM          dm, sw;
473:   SNES        snes;
474:   Vec         u;
475:   AppCtx      user;
476:   const char *names[] = {"q", "phi"};

478:   PetscFunctionBeginUser;
479:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
480:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
481:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
482:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
483:   PetscCall(SNESSetDM(snes, dm));
484:   PetscCall(SetupDiscretization(dm, 2, names, SetupPrimalProblem, &user));
485:   if (user.particleRHS) {
486:     PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag));
487:     PetscCall(CreateSwarm(dm, &user, &sw));
488:     PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
489:     PetscCall(InitializeConstants(sw, &user));
490:   }
491:   PetscCall(DMCreateGlobalVector(dm, &u));
492:   PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
493:   PetscCall(SNESSetFromOptions(snes));
494:   PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user));
495:   PetscCall(DMSNESCheckFromOptions(snes, u));
496:   if (user.particleRHS) {
497:     DM       potential_dm;
498:     IS       potential_IS;
499:     Mat      M_p;
500:     Vec      rho, f, temp_rho;
501:     PetscInt fields = 1;

503:     PetscCall(DMGetGlobalVector(dm, &rho));
504:     PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
505:     PetscCall(DMCreateSubDM(dm, 1, &fields, &potential_IS, &potential_dm));
506:     PetscCall(DMCreateMassMatrix(sw, potential_dm, &M_p));
507:     PetscCall(MatViewFromOptions(M_p, NULL, "-mp_view"));
508:     PetscCall(DMGetGlobalVector(potential_dm, &temp_rho));
509:     PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
510:     PetscCall(PetscObjectSetName((PetscObject)f, "particle weight"));
511:     PetscCall(VecViewFromOptions(f, NULL, "-weights_view"));
512:     PetscCall(MatMultTranspose(M_p, f, temp_rho));
513:     PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
514:     PetscCall(MatDestroy(&M_p));
515:     PetscCall(PetscObjectSetName((PetscObject)rho, "rho"));
516:     PetscCall(VecViewFromOptions(rho, NULL, "-poisson_rho_view"));
517:     PetscCall(VecISCopy(rho, potential_IS, SCATTER_FORWARD, temp_rho));
518:     PetscCall(VecViewFromOptions(temp_rho, NULL, "-rho_view"));
519:     PetscCall(DMRestoreGlobalVector(potential_dm, &temp_rho));
520:     PetscCall(DMDestroy(&potential_dm));
521:     PetscCall(ISDestroy(&potential_IS));

523:     PetscCall(SNESSolve(snes, rho, u));
524:     PetscCall(DMRestoreGlobalVector(dm, &rho));
525:   } else {
526:     PetscCall(SNESSolve(snes, NULL, u));
527:   }
528:   PetscCall(VecDestroy(&u));
529:   PetscCall(SNESDestroy(&snes));
530:   PetscCall(DMDestroy(&dm));
531:   if (user.particleRHS) {
532:     PetscCall(DMDestroy(&sw));
533:     PetscCall(PetscBagDestroy(&user.bag));
534:   }
535:   PetscCall(PetscFinalize());
536:   return PETSC_SUCCESS;
537: }

539: /*TEST

541:   # RT1-P0 on quads
542:   testset:
543:     args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,1 \
544:           -dm_plex_box_lower 0,-1 -dm_plex_box_upper 6.283185307179586,1\
545:           -phi_petscspace_degree 0 \
546:           -phi_petscdualspace_lagrange_use_moments \
547:           -phi_petscdualspace_lagrange_moment_order 2 \
548:           -q_petscfe_default_quadrature_order 1 \
549:           -q_petscspace_type sum \
550:           -q_petscspace_variables 2 \
551:           -q_petscspace_components 2 \
552:           -q_petscspace_sum_spaces 2 \
553:           -q_petscspace_sum_concatenate true \
554:           -q_sumcomp_0_petscspace_variables 2 \
555:           -q_sumcomp_0_petscspace_type tensor \
556:           -q_sumcomp_0_petscspace_tensor_spaces 2 \
557:           -q_sumcomp_0_petscspace_tensor_uniform false \
558:           -q_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
559:           -q_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
560:           -q_sumcomp_1_petscspace_variables 2 \
561:           -q_sumcomp_1_petscspace_type tensor \
562:           -q_sumcomp_1_petscspace_tensor_spaces 2 \
563:           -q_sumcomp_1_petscspace_tensor_uniform false \
564:           -q_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
565:           -q_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
566:           -q_petscdualspace_form_degree -1 \
567:           -q_petscdualspace_order 1 \
568:           -q_petscdualspace_lagrange_trimmed true \
569:           -ksp_error_if_not_converged \
570:           -pc_type fieldsplit -pc_fieldsplit_type schur \
571:           -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full

573:     # The Jacobian test is meaningless here
574:     test:
575:           suffix: quad_hdiv_0
576:           args: -dmsnes_check
577:           filter: sed -e "s/Taylor approximation converging at order.*''//"

579:     # The Jacobian test is meaningless here
580:     test:
581:           suffix: quad_hdiv_1
582:           args: -sol_type linear -dmsnes_check
583:           filter: sed -e "s/Taylor approximation converging at order.*''//"

585:     test:
586:           suffix: quad_hdiv_2
587:           args: -sol_type quadratic -dmsnes_check \
588:                 -fieldsplit_q_pc_type lu -fieldsplit_phi_pc_type svd

590:     test:
591:           suffix: quad_hdiv_3
592:           args: -sol_type trig \
593:                 -fieldsplit_q_pc_type lu -fieldsplit_phi_pc_type svd

595:     test:
596:           suffix: quad_hdiv_4
597:           requires: !single
598:           args: -sol_type trigx \
599:                 -fieldsplit_q_pc_type lu -fieldsplit_phi_pc_type svd

601:     test:
602:           suffix: particle_hdiv_5
603:           requires: !complex
604:           args: -dm_swarm_num_particles 100 -particleRHS -sol_type particles \
605:                 -fieldsplit_q_pc_type lu -fieldsplit_phi_pc_type svd

607: TEST*/