xref: /petsc/src/snes/tests/ex15.c (revision 21e3ffae2f3b73c0bd738cf6d0a809700fc04bb0) 
1 static char help[] = "Mixed element discretization of the Poisson equation.\n\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscdmswarm.h>
5 #include <petscds.h>
6 #include <petscsnes.h>
7 #include <petscconvest.h>
8 #include <petscbag.h>
9 
10 /*
11 The Poisson equation
12 
13   -\Delta\phi = f
14 
15 can be rewritten in first order form
16 
17   q - \nabla\phi  &= 0
18   -\nabla \cdot q &= f
19 */
20 
21 typedef enum {
22   SIGMA,
23   NUM_CONSTANTS
24 } ConstantType;
25 typedef struct {
26   PetscReal sigma; /* Nondimensional charge per length in x */
27 } Parameter;
28 
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"};
39 
40 typedef struct {
41   SolType   solType; /* MMS solution type */
42   PetscBag  bag;     /* Problem parameters */
43   PetscBool particleRHS;
44   PetscInt  Np;
45 } AppCtx;
46 
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 }
53 
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 }
59 
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 }
66 
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 }
73 
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 }
80 
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 }
87 
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 }
93 
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 }
98 
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 }
105 
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 }
112 
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 }
118 
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 }
123 
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 }
130 
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 }
137 
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 }
143 
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 }
148 
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 }
153 
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 }
158 
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 }
163 
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 }
169 
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 }
174 
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 }
179 
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 }
184 
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 }
191 
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 }
198 
199 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
200 {
201   PetscInt sol;
202 
203   PetscFunctionBeginUser;
204   options->solType     = SOL_CONST;
205   options->particleRHS = PETSC_FALSE;
206   options->Np          = 100;
207 
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 }
216 
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 }
227 
228 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
229 {
230   PetscDS        ds;
231   PetscWeakForm  wf;
232   DMLabel        label;
233   const PetscInt id = 1;
234 
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   }
282 
283   PetscFunctionReturn(PETSC_SUCCESS);
284 }
285 
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];
293 
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;
304 
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 }
319 
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;
328 
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();
335 
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));
341 
342   PetscCall(PetscBagGetData(user->bag, (void **)&param));
343   PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
344 
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));
351 
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   }
358 
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 }
364 
365 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
366 {
367   PetscInt dim;
368 
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));
377 
378   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
379 
380   PetscCall(InitializeParticlesAndWeights(*sw, user));
381 
382   PetscCall(DMSetFromOptions(*sw));
383   PetscCall(DMSetApplicationContext(*sw, user));
384   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
385   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
386 
387   PetscFunctionReturn(PETSC_SUCCESS);
388 }
389 
390 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
391 {
392   PetscBag   bag;
393   Parameter *p;
394 
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;
406 
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 }
418 
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;
424 
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;
437 
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));
454 
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 }
469 
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"};
477 
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;
502 
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));
522 
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 }
538 
539 /*TEST
540 
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
572 
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.*''//"
578 
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.*''//"
584 
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
589 
590     test:
591           suffix: quad_hdiv_3
592           args: -sol_type trig \
593                 -fieldsplit_q_pc_type lu -fieldsplit_phi_pc_type svd
594 
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
600 
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
606 
607 TEST*/
608