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 **)¶m));
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 **)¶m));
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 **)¶m));
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*/