Actual source code: ex77.c
1: static char help[] = "Time-dependent reactive low Mach Flow in 2d and 3d channels with finite elements.\n\
2: We solve the reactive low Mach flow problem in a rectangular domain\n\
3: using a parallel unstructured mesh (DMPLEX) to discretize the flow\n\
4: and particles (DWSWARM) to discretize the chemical species.\n\n\n";
6: /*F
7: This low Mach flow is time-dependent isoviscous Navier-Stokes flow. We discretize using the
8: finite element method on an unstructured mesh. The weak form equations are
10: \begin{align*}
11: < q, \nabla\cdot u > = 0
12: <v, du/dt> + <v, u \cdot \nabla u> + < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > - < v, f > = 0
13: < w, u \cdot \nabla T > + < \nabla w, \alpha \nabla T > - < w, Q > = 0
14: \end{align*}
16: where $\nu$ is the kinematic viscosity and $\alpha$ is thermal diffusivity.
18: For visualization, use
20: -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append
22: The particles can be visualized using
24: -part_dm_view draw -part_dm_view_swarm_radius 0.03
26: F*/
28: #include <petscdmplex.h>
29: #include <petscdmswarm.h>
30: #include <petscts.h>
31: #include <petscds.h>
32: #include <petscbag.h>
34: typedef enum {
35: SOL_TRIG_TRIG,
36: NUM_SOL_TYPES
37: } SolType;
38: const char *solTypes[NUM_SOL_TYPES + 1] = {"trig_trig", "unknown"};
40: typedef enum {
41: PART_LAYOUT_CELL,
42: PART_LAYOUT_BOX,
43: NUM_PART_LAYOUT_TYPES
44: } PartLayoutType;
45: const char *partLayoutTypes[NUM_PART_LAYOUT_TYPES + 1] = {"cell", "box", "unknown"};
47: typedef struct {
48: PetscReal nu; /* Kinematic viscosity */
49: PetscReal alpha; /* Thermal diffusivity */
50: PetscReal T_in; /* Inlet temperature*/
51: PetscReal omega; /* Rotation speed in MMS benchmark */
52: } Parameter;
54: typedef struct {
55: /* Problem definition */
56: PetscBag bag; /* Holds problem parameters */
57: SolType solType; /* MMS solution type */
58: PartLayoutType partLayout; /* Type of particle distribution */
59: PetscInt Npc; /* The initial number of particles per cell */
60: PetscReal partLower[3]; /* Lower left corner of particle box */
61: PetscReal partUpper[3]; /* Upper right corner of particle box */
62: PetscInt Npb; /* The initial number of particles per box dimension */
63: } AppCtx;
65: typedef struct {
66: PetscReal ti; /* The time for ui, at the beginning of the advection solve */
67: PetscReal tf; /* The time for uf, at the end of the advection solve */
68: Vec ui; /* The PDE solution field at ti */
69: Vec uf; /* The PDE solution field at tf */
70: Vec x0; /* The initial particle positions at t = 0 */
71: PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
72: AppCtx *ctx; /* Context for exact solution */
73: } AdvCtx;
75: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
76: {
77: PetscInt d;
78: for (d = 0; d < Nc; ++d) u[d] = 0.0;
79: return PETSC_SUCCESS;
80: }
82: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
83: {
84: PetscInt d;
85: for (d = 0; d < Nc; ++d) u[d] = 1.0;
86: return PETSC_SUCCESS;
87: }
89: /*
90: CASE: trigonometric-trigonometric
91: In 2D we use exact solution:
93: x = r0 cos(w t + theta0) r0 = sqrt(x0^2 + y0^2)
94: y = r0 sin(w t + theta0) theta0 = arctan(y0/x0)
95: u = -w r0 sin(theta0) = -w y
96: v = w r0 cos(theta0) = w x
97: p = x + y - 1
98: T = t + x + y
99: f = <1, 1>
100: Q = 1 + w (x - y)/r
102: so that
104: \nabla \cdot u = 0 + 0 = 0
106: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
107: = <0, 0> + u_i d_i u_j - \nu 0 + <1, 1>
108: = <1, 1> + w^2 <-y, x> . <<0, 1>, <-1, 0>>
109: = <1, 1> + w^2 <-x, -y>
110: = <1, 1> - w^2 <x, y>
112: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
113: = 1 + <u, v> . <1, 1> - \alpha 0
114: = 1 + u + v
115: */
116: static PetscErrorCode trig_trig_x(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *x, PetscCtx ctx)
117: {
118: const PetscReal x0 = X[0];
119: const PetscReal y0 = X[1];
120: const PetscReal R0 = PetscSqrtReal(x0 * x0 + y0 * y0);
121: const PetscReal theta0 = PetscAtan2Real(y0, x0);
122: Parameter *p = (Parameter *)ctx;
124: x[0] = R0 * PetscCosReal(p->omega * time + theta0);
125: x[1] = R0 * PetscSinReal(p->omega * time + theta0);
126: return PETSC_SUCCESS;
127: }
128: static PetscErrorCode trig_trig_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
129: {
130: Parameter *p = (Parameter *)ctx;
132: u[0] = -p->omega * X[1];
133: u[1] = p->omega * X[0];
134: return PETSC_SUCCESS;
135: }
136: static PetscErrorCode trig_trig_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
137: {
138: u[0] = 0.0;
139: u[1] = 0.0;
140: return PETSC_SUCCESS;
141: }
143: static PetscErrorCode trig_trig_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
144: {
145: p[0] = X[0] + X[1] - 1.0;
146: return PETSC_SUCCESS;
147: }
149: static PetscErrorCode trig_trig_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
150: {
151: T[0] = time + X[0] + X[1];
152: return PETSC_SUCCESS;
153: }
154: static PetscErrorCode trig_trig_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
155: {
156: T[0] = 1.0;
157: return PETSC_SUCCESS;
158: }
160: static void f0_trig_trig_v(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[])
161: {
162: const PetscReal omega = PetscRealPart(constants[3]);
163: PetscInt Nc = dim;
164: PetscInt c, d;
166: for (d = 0; d < dim; ++d) f0[d] = u_t[uOff[0] + d];
168: for (c = 0; c < Nc; ++c) {
169: for (d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d];
170: }
171: f0[0] -= 1.0 - omega * omega * X[0];
172: f0[1] -= 1.0 - omega * omega * X[1];
173: }
175: static void f0_trig_trig_w(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[])
176: {
177: const PetscReal omega = PetscRealPart(constants[3]);
178: PetscInt d;
180: for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0] + d] * u_x[uOff_x[2] + d];
181: f0[0] += u_t[uOff[2]] - (1.0 + omega * (X[0] - X[1]));
182: }
184: 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[])
185: {
186: PetscInt d;
187: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
188: }
190: /*f1_v = \nu[grad(u) + grad(u)^T] - pI */
191: static void f1_v(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[])
192: {
193: const PetscReal nu = PetscRealPart(constants[0]);
194: const PetscInt Nc = dim;
195: PetscInt c, d;
197: for (c = 0; c < Nc; ++c) {
198: for (d = 0; d < dim; ++d) f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
199: f1[c * dim + c] -= u[uOff[1]];
200: }
201: }
203: static void f1_w(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[])
204: {
205: const PetscReal alpha = PetscRealPart(constants[1]);
206: PetscInt d;
207: for (d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
208: }
210: /* Jacobians */
211: static void g1_qu(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[])
212: {
213: PetscInt d;
214: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
215: }
217: static void g0_vu(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[])
218: {
219: PetscInt c, d;
220: const PetscInt Nc = dim;
222: for (d = 0; d < dim; ++d) g0[d * dim + d] = u_tShift;
224: for (c = 0; c < Nc; ++c) {
225: for (d = 0; d < dim; ++d) g0[c * Nc + d] += u_x[c * Nc + d];
226: }
227: }
229: static void g1_vu(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[])
230: {
231: PetscInt NcI = dim;
232: PetscInt NcJ = dim;
233: PetscInt c, d, e;
235: for (c = 0; c < NcI; ++c) {
236: for (d = 0; d < NcJ; ++d) {
237: for (e = 0; e < dim; ++e) {
238: if (c == d) g1[(c * NcJ + d) * dim + e] += u[e];
239: }
240: }
241: }
242: }
244: static void g2_vp(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[])
245: {
246: PetscInt d;
247: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
248: }
250: static void g3_vu(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[])
251: {
252: const PetscReal nu = PetscRealPart(constants[0]);
253: const PetscInt Nc = dim;
255: for (PetscInt c = 0; c < Nc; ++c) {
256: for (PetscInt d = 0; d < dim; ++d) {
257: g3[((c * Nc + c) * dim + d) * dim + d] += nu; // gradU
258: g3[((c * Nc + d) * dim + d) * dim + c] += nu; // gradU transpose
259: }
260: }
261: }
263: static void g0_wT(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[])
264: {
265: for (PetscInt d = 0; d < dim; ++d) g0[d] = u_tShift;
266: }
268: static void g0_wu(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[])
269: {
270: for (PetscInt d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
271: }
273: static void g1_wT(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[])
274: {
275: for (PetscInt d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
276: }
278: static void g3_wT(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[])
279: {
280: const PetscReal alpha = PetscRealPart(constants[1]);
282: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
283: }
285: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
286: {
287: PetscInt sol, pl, n;
289: PetscFunctionBeginUser;
290: options->solType = SOL_TRIG_TRIG;
291: options->partLayout = PART_LAYOUT_CELL;
292: options->Npc = 1;
293: options->Npb = 1;
294: options->partLower[0] = options->partLower[1] = options->partLower[2] = 0.;
295: options->partUpper[0] = options->partUpper[1] = options->partUpper[2] = 1.;
296: PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");
297: sol = options->solType;
298: PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex77.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
299: options->solType = (SolType)sol;
300: pl = options->partLayout;
301: PetscCall(PetscOptionsEList("-part_layout_type", "The particle layout type", "ex77.c", partLayoutTypes, NUM_PART_LAYOUT_TYPES, partLayoutTypes[options->partLayout], &pl, NULL));
302: options->partLayout = (PartLayoutType)pl;
303: PetscCall(PetscOptionsInt("-Npc", "The initial number of particles per cell", "ex77.c", options->Npc, &options->Npc, NULL));
304: n = 3;
305: PetscCall(PetscOptionsRealArray("-part_lower", "The lower left corner of the particle box", "ex77.c", options->partLower, &n, NULL));
306: n = 3;
307: PetscCall(PetscOptionsRealArray("-part_upper", "The upper right corner of the particle box", "ex77.c", options->partUpper, &n, NULL));
308: PetscCall(PetscOptionsInt("-Npb", "The initial number of particles per box dimension", "ex77.c", options->Npb, &options->Npb, NULL));
309: PetscOptionsEnd();
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode SetupParameters(AppCtx *user)
314: {
315: PetscBag bag;
316: Parameter *p;
318: PetscFunctionBeginUser;
319: /* setup PETSc parameter bag */
320: PetscCall(PetscBagGetData(user->bag, &p));
321: PetscCall(PetscBagSetName(user->bag, "par", "Low Mach flow parameters"));
322: bag = user->bag;
323: PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
324: PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
325: PetscCall(PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature"));
326: PetscCall(PetscBagRegisterReal(bag, &p->omega, 1.0, "omega", "Rotation speed in MMS benchmark"));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
331: {
332: PetscFunctionBeginUser;
333: PetscCall(DMCreate(comm, dm));
334: PetscCall(DMSetType(*dm, DMPLEX));
335: PetscCall(DMSetFromOptions(*dm));
336: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
341: {
342: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
343: PetscErrorCode (*exactFuncs_t[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
344: PetscDS prob;
345: DMLabel label;
346: Parameter *ctx;
347: PetscInt id;
349: PetscFunctionBeginUser;
350: PetscCall(DMGetLabel(dm, "marker", &label));
351: PetscCall(DMGetDS(dm, &prob));
352: switch (user->solType) {
353: case SOL_TRIG_TRIG:
354: PetscCall(PetscDSSetResidual(prob, 0, f0_trig_trig_v, f1_v));
355: PetscCall(PetscDSSetResidual(prob, 2, f0_trig_trig_w, f1_w));
357: exactFuncs[0] = trig_trig_u;
358: exactFuncs[1] = trig_trig_p;
359: exactFuncs[2] = trig_trig_T;
360: exactFuncs_t[0] = trig_trig_u_t;
361: exactFuncs_t[1] = NULL;
362: exactFuncs_t[2] = trig_trig_T_t;
363: break;
364: default:
365: SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
366: }
368: PetscCall(PetscDSSetResidual(prob, 1, f0_q, NULL));
370: PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu));
371: PetscCall(PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL));
372: PetscCall(PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL));
373: PetscCall(PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL));
374: PetscCall(PetscDSSetJacobian(prob, 2, 2, g0_wT, g1_wT, NULL, g3_wT));
375: /* Setup constants */
376: {
377: Parameter *param;
378: PetscScalar constants[4];
380: PetscCall(PetscBagGetData(user->bag, ¶m));
382: constants[0] = param->nu;
383: constants[1] = param->alpha;
384: constants[2] = param->T_in;
385: constants[3] = param->omega;
386: PetscCall(PetscDSSetConstants(prob, 4, constants));
387: }
388: /* Setup Boundary Conditions */
389: PetscCall(PetscBagGetData(user->bag, &ctx));
390: id = 3;
391: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], (PetscVoidFn *)exactFuncs_t[0], ctx, NULL));
392: id = 1;
393: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], (PetscVoidFn *)exactFuncs_t[0], ctx, NULL));
394: id = 2;
395: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], (PetscVoidFn *)exactFuncs_t[0], ctx, NULL));
396: id = 4;
397: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], (PetscVoidFn *)exactFuncs_t[0], ctx, NULL));
398: id = 3;
399: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], (PetscVoidFn *)exactFuncs_t[2], ctx, NULL));
400: id = 1;
401: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], (PetscVoidFn *)exactFuncs_t[2], ctx, NULL));
402: id = 2;
403: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], (PetscVoidFn *)exactFuncs_t[2], ctx, NULL));
404: id = 4;
405: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], (PetscVoidFn *)exactFuncs_t[2], ctx, NULL));
407: /*setup exact solution.*/
408: PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx));
409: PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx));
410: PetscCall(PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx));
411: PetscCall(PetscDSSetExactSolutionTimeDerivative(prob, 0, exactFuncs_t[0], ctx));
412: PetscCall(PetscDSSetExactSolutionTimeDerivative(prob, 1, exactFuncs_t[1], ctx));
413: PetscCall(PetscDSSetExactSolutionTimeDerivative(prob, 2, exactFuncs_t[2], ctx));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /* x_t = v
419: Note that here we use the velocity field at t_{n+1} to advect the particles from
420: t_n to t_{n+1}. If we use both of these fields, we could use Crank-Nicholson or
421: the method of characteristics.
422: */
423: static PetscErrorCode FreeStreaming(TS ts, PetscReal t, Vec X, Vec F, PetscCtx ctx)
424: {
425: AdvCtx *adv = (AdvCtx *)ctx;
426: Vec u = adv->ui;
427: DM sdm, dm, vdm;
428: Vec vel, locvel, pvel;
429: IS vis;
430: DMInterpolationInfo ictx;
431: const PetscScalar *coords, *v;
432: PetscScalar *f;
433: PetscInt vf[1] = {0};
434: PetscInt dim, Np;
436: PetscFunctionBeginUser;
437: PetscCall(TSGetDM(ts, &sdm));
438: PetscCall(DMSwarmGetCellDM(sdm, &dm));
439: PetscCall(DMGetGlobalVector(sdm, &pvel));
440: PetscCall(DMSwarmGetLocalSize(sdm, &Np));
441: PetscCall(DMGetDimension(dm, &dim));
442: /* Get local velocity */
443: PetscCall(DMCreateSubDM(dm, 1, vf, &vis, &vdm));
444: PetscCall(VecGetSubVector(u, vis, &vel));
445: PetscCall(DMGetLocalVector(vdm, &locvel));
446: PetscCall(DMPlexInsertBoundaryValues(vdm, PETSC_TRUE, locvel, adv->ti, NULL, NULL, NULL));
447: PetscCall(DMGlobalToLocalBegin(vdm, vel, INSERT_VALUES, locvel));
448: PetscCall(DMGlobalToLocalEnd(vdm, vel, INSERT_VALUES, locvel));
449: PetscCall(VecRestoreSubVector(u, vis, &vel));
450: PetscCall(ISDestroy(&vis));
451: /* Interpolate velocity */
452: PetscCall(DMInterpolationCreate(PETSC_COMM_SELF, &ictx));
453: PetscCall(DMInterpolationSetDim(ictx, dim));
454: PetscCall(DMInterpolationSetDof(ictx, dim));
455: PetscCall(VecGetArrayRead(X, &coords));
456: PetscCall(DMInterpolationAddPoints(ictx, Np, (PetscReal *)coords));
457: PetscCall(VecRestoreArrayRead(X, &coords));
458: /* Particles that lie outside the domain should be dropped,
459: whereas particles that move to another partition should trigger a migration */
460: PetscCall(DMInterpolationSetUp(ictx, vdm, PETSC_FALSE, PETSC_TRUE));
461: PetscCall(VecSet(pvel, 0.));
462: PetscCall(DMInterpolationEvaluate(ictx, vdm, locvel, pvel));
463: PetscCall(DMInterpolationDestroy(&ictx));
464: PetscCall(DMRestoreLocalVector(vdm, &locvel));
465: PetscCall(DMDestroy(&vdm));
467: PetscCall(VecGetArray(F, &f));
468: PetscCall(VecGetArrayRead(pvel, &v));
469: PetscCall(PetscArraycpy(f, v, Np * dim));
470: PetscCall(VecRestoreArrayRead(pvel, &v));
471: PetscCall(VecRestoreArray(F, &f));
472: PetscCall(DMRestoreGlobalVector(sdm, &pvel));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: static PetscErrorCode SetInitialParticleConditions(TS ts, Vec u)
477: {
478: AppCtx *user;
479: void *ctx;
480: DM dm;
481: PetscScalar *coords;
482: PetscReal x[3], dx[3];
483: PetscInt n[3];
484: PetscInt dim, d, i, j, k;
486: PetscFunctionBeginUser;
487: PetscCall(TSGetApplicationContext(ts, &ctx));
488: user = ((AdvCtx *)ctx)->ctx;
489: PetscCall(TSGetDM(ts, &dm));
490: PetscCall(DMGetDimension(dm, &dim));
491: switch (user->partLayout) {
492: case PART_LAYOUT_CELL:
493: PetscCall(DMSwarmSetPointCoordinatesRandom(dm, user->Npc));
494: break;
495: case PART_LAYOUT_BOX:
496: for (d = 0; d < dim; ++d) {
497: n[d] = user->Npb;
498: dx[d] = (user->partUpper[d] - user->partLower[d]) / PetscMax(1, n[d] - 1);
499: }
500: PetscCall(VecGetArray(u, &coords));
501: switch (dim) {
502: case 2:
503: x[0] = user->partLower[0];
504: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
505: x[1] = user->partLower[1];
506: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
507: const PetscInt p = j * n[0] + i;
508: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
509: }
510: }
511: break;
512: case 3:
513: x[0] = user->partLower[0];
514: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
515: x[1] = user->partLower[1];
516: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
517: x[2] = user->partLower[2];
518: for (k = 0; k < n[2]; ++k, x[2] += dx[2]) {
519: const PetscInt p = (k * n[1] + j) * n[0] + i;
520: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
521: }
522: }
523: }
524: break;
525: default:
526: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Do not support particle layout in dimension %" PetscInt_FMT, dim);
527: }
528: PetscCall(VecRestoreArray(u, &coords));
529: break;
530: default:
531: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]);
532: }
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: static PetscErrorCode SetupDiscretization(DM dm, DM sdm, AppCtx *user)
537: {
538: DM cdm = dm;
539: DMSwarmCellDM celldm;
540: PetscFE fe[3];
541: Parameter *param;
542: PetscInt *swarm_cellid, n[3];
543: PetscReal x[3], dx[3];
544: PetscScalar *coords;
545: DMPolytopeType ct;
546: PetscInt dim, d, cStart, cEnd, c, Np, p, i, j, k;
547: PetscBool simplex;
548: MPI_Comm comm;
549: const char *cellid;
551: PetscFunctionBeginUser;
552: PetscCall(DMGetDimension(dm, &dim));
553: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
554: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
555: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
556: /* Create finite element */
557: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
558: PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
559: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
561: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
562: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
563: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
565: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
566: PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
567: PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));
569: /* Set discretization and boundary conditions for each mesh */
570: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
571: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
572: PetscCall(DMSetField(dm, 2, NULL, (PetscObject)fe[2]));
573: PetscCall(DMCreateDS(dm));
574: PetscCall(SetupProblem(dm, user));
575: PetscCall(PetscBagGetData(user->bag, ¶m));
576: while (cdm) {
577: PetscCall(DMCopyDisc(dm, cdm));
578: PetscCall(DMGetCoarseDM(cdm, &cdm));
579: }
580: PetscCall(PetscFEDestroy(&fe[0]));
581: PetscCall(PetscFEDestroy(&fe[1]));
582: PetscCall(PetscFEDestroy(&fe[2]));
584: {
585: PetscObject pressure;
586: MatNullSpace nullspacePres;
588: PetscCall(DMGetField(dm, 1, NULL, &pressure));
589: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
590: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
591: PetscCall(MatNullSpaceDestroy(&nullspacePres));
592: }
594: /* Setup particle information */
595: PetscCall(DMSwarmSetType(sdm, DMSWARM_PIC));
596: PetscCall(DMSwarmRegisterPetscDatatypeField(sdm, "mass", 1, PETSC_REAL));
597: PetscCall(DMSwarmFinalizeFieldRegister(sdm));
598: PetscCall(DMSwarmGetCellDMActive(sdm, &celldm));
599: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
600: switch (user->partLayout) {
601: case PART_LAYOUT_CELL:
602: PetscCall(DMSwarmSetLocalSizes(sdm, (cEnd - cStart) * user->Npc, 0));
603: PetscCall(DMSetFromOptions(sdm));
604: PetscCall(DMSwarmGetField(sdm, cellid, NULL, NULL, (void **)&swarm_cellid));
605: for (c = cStart; c < cEnd; ++c) {
606: for (p = 0; p < user->Npc; ++p) {
607: const PetscInt n = c * user->Npc + p;
609: swarm_cellid[n] = c;
610: }
611: }
612: PetscCall(DMSwarmRestoreField(sdm, cellid, NULL, NULL, (void **)&swarm_cellid));
613: PetscCall(DMSwarmSetPointCoordinatesRandom(sdm, user->Npc));
614: break;
615: case PART_LAYOUT_BOX:
616: Np = 1;
617: for (d = 0; d < dim; ++d) {
618: n[d] = user->Npb;
619: dx[d] = (user->partUpper[d] - user->partLower[d]) / PetscMax(1, n[d] - 1);
620: Np *= n[d];
621: }
622: PetscCall(DMSwarmSetLocalSizes(sdm, Np, 0));
623: PetscCall(DMSetFromOptions(sdm));
624: PetscCall(DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
625: switch (dim) {
626: case 2:
627: x[0] = user->partLower[0];
628: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
629: x[1] = user->partLower[1];
630: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
631: const PetscInt p = j * n[0] + i;
632: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
633: }
634: }
635: break;
636: case 3:
637: x[0] = user->partLower[0];
638: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
639: x[1] = user->partLower[1];
640: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
641: x[2] = user->partLower[2];
642: for (k = 0; k < n[2]; ++k, x[2] += dx[2]) {
643: const PetscInt p = (k * n[1] + j) * n[0] + i;
644: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
645: }
646: }
647: }
648: break;
649: default:
650: SETERRQ(comm, PETSC_ERR_SUP, "Do not support particle layout in dimension %" PetscInt_FMT, dim);
651: }
652: PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
653: PetscCall(DMSwarmGetField(sdm, cellid, NULL, NULL, (void **)&swarm_cellid));
654: for (p = 0; p < Np; ++p) swarm_cellid[p] = 0;
655: PetscCall(DMSwarmRestoreField(sdm, cellid, NULL, NULL, (void **)&swarm_cellid));
656: PetscCall(DMSwarmMigrate(sdm, PETSC_TRUE));
657: break;
658: default:
659: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]);
660: }
661: PetscCall(PetscObjectSetName((PetscObject)sdm, "Particles"));
662: PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view"));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
667: {
668: Vec vec;
669: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
671: PetscFunctionBeginUser;
672: PetscCheck(ofield == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %" PetscInt_FMT, ofield);
673: funcs[nfield] = constant;
674: PetscCall(DMCreateGlobalVector(dm, &vec));
675: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
676: PetscCall(VecNormalize(vec, NULL));
677: PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
678: PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
679: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
680: PetscCall(VecDestroy(&vec));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u)
685: {
686: DM dm;
687: MatNullSpace nullsp;
689: PetscFunctionBeginUser;
690: PetscCall(TSGetDM(ts, &dm));
691: PetscCall(CreatePressureNullSpace(dm, 1, 1, &nullsp));
692: PetscCall(MatNullSpaceRemove(nullsp, u));
693: PetscCall(MatNullSpaceDestroy(&nullsp));
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: /* Make the discrete pressure discretely divergence free */
698: static PetscErrorCode RemoveDiscretePressureNullspace(TS ts)
699: {
700: Vec u;
702: PetscFunctionBeginUser;
703: PetscCall(TSGetSolution(ts, &u));
704: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
709: {
710: DM dm;
711: PetscReal t;
713: PetscFunctionBeginUser;
714: PetscCall(TSGetDM(ts, &dm));
715: PetscCall(TSGetTime(ts, &t));
716: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
717: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, PetscCtx ctx)
722: {
723: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
724: void *ctxs[3];
725: DM dm;
726: PetscDS ds;
727: Vec v;
728: PetscReal ferrors[3];
729: PetscInt tl, l, f;
731: PetscFunctionBeginUser;
732: PetscCall(TSGetDM(ts, &dm));
733: PetscCall(DMGetDS(dm, &ds));
735: for (f = 0; f < 3; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]));
736: PetscCall(DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors));
737: PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tl));
738: for (l = 0; l < tl; ++l) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t"));
739: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g]\n", (int)step, (double)crtime, (double)ferrors[0], (double)ferrors[1], (double)ferrors[2]));
741: PetscCall(DMGetGlobalVector(dm, &u));
742: PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
743: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
744: PetscCall(DMRestoreGlobalVector(dm, &u));
746: PetscCall(DMGetGlobalVector(dm, &v));
747: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, v));
748: PetscCall(PetscObjectSetName((PetscObject)v, "Exact Solution"));
749: PetscCall(VecViewFromOptions(v, NULL, "-exact_vec_view"));
750: PetscCall(DMRestoreGlobalVector(dm, &v));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: /* Note that adv->x0 will not be correct after migration */
755: static PetscErrorCode ComputeParticleError(TS ts, Vec u, Vec e)
756: {
757: AdvCtx *adv;
758: DM sdm;
759: Parameter *param;
760: const PetscScalar *xp0, *xp;
761: PetscScalar *ep;
762: PetscReal time;
763: PetscInt dim, Np, p;
764: MPI_Comm comm;
766: PetscFunctionBeginUser;
767: PetscCall(TSGetTime(ts, &time));
768: PetscCall(TSGetApplicationContext(ts, &adv));
769: PetscCall(PetscBagGetData(adv->ctx->bag, ¶m));
770: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
771: PetscCall(TSGetDM(ts, &sdm));
772: PetscCall(DMGetDimension(sdm, &dim));
773: PetscCall(DMSwarmGetLocalSize(sdm, &Np));
774: PetscCall(VecGetArrayRead(adv->x0, &xp0));
775: PetscCall(VecGetArrayRead(u, &xp));
776: PetscCall(VecGetArrayWrite(e, &ep));
777: for (p = 0; p < Np; ++p) {
778: PetscScalar x[3];
779: PetscReal x0[3];
781: for (PetscInt d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p * dim + d]);
782: PetscCall(adv->exact(dim, time, x0, 1, x, param));
783: for (PetscInt d = 0; d < dim; ++d) ep[p * dim + d] += x[d] - xp[p * dim + d];
784: }
785: PetscCall(VecRestoreArrayRead(adv->x0, &xp0));
786: PetscCall(VecRestoreArrayRead(u, &xp));
787: PetscCall(VecRestoreArrayWrite(e, &ep));
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: static PetscErrorCode MonitorParticleError(TS ts, PetscInt step, PetscReal time, Vec u, PetscCtx ctx)
792: {
793: AdvCtx *adv = (AdvCtx *)ctx;
794: DM sdm;
795: Parameter *param;
796: const PetscScalar *xp0, *xp;
797: PetscReal error = 0.0;
798: PetscInt dim, tl, l, Np, p;
799: MPI_Comm comm;
801: PetscFunctionBeginUser;
802: PetscCall(PetscBagGetData(adv->ctx->bag, ¶m));
803: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
804: PetscCall(TSGetDM(ts, &sdm));
805: PetscCall(DMGetDimension(sdm, &dim));
806: PetscCall(DMSwarmGetLocalSize(sdm, &Np));
807: PetscCall(VecGetArrayRead(adv->x0, &xp0));
808: PetscCall(VecGetArrayRead(u, &xp));
809: for (p = 0; p < Np; ++p) {
810: PetscScalar x[3];
811: PetscReal x0[3];
812: PetscReal perror = 0.0;
814: for (PetscInt d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p * dim + d]);
815: PetscCall(adv->exact(dim, time, x0, 1, x, param));
816: for (PetscInt d = 0; d < dim; ++d) perror += PetscSqr(PetscRealPart(x[d] - xp[p * dim + d]));
817: error += perror;
818: }
819: PetscCall(VecRestoreArrayRead(adv->x0, &xp0));
820: PetscCall(VecRestoreArrayRead(u, &xp));
821: PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tl));
822: for (l = 0; l < tl; ++l) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t"));
823: PetscCall(PetscPrintf(comm, "Timestep: %04d time = %-8.4g \t L_2 Particle Error: [%2.3g]\n", (int)step, (double)time, (double)error));
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: static PetscErrorCode AdvectParticles(TS ts)
828: {
829: TS sts;
830: DM sdm;
831: Vec coordinates;
832: AdvCtx *adv;
833: PetscReal time;
834: PetscBool lreset, reset;
835: PetscInt dim, n, N, newn, newN;
837: PetscFunctionBeginUser;
838: PetscCall(PetscObjectQuery((PetscObject)ts, "_SwarmTS", (PetscObject *)&sts));
839: PetscCall(TSGetDM(sts, &sdm));
840: PetscCall(TSGetRHSFunction(sts, NULL, NULL, &adv));
841: PetscCall(DMGetDimension(sdm, &dim));
842: PetscCall(DMSwarmGetSize(sdm, &N));
843: PetscCall(DMSwarmGetLocalSize(sdm, &n));
844: PetscCall(DMSwarmCreateGlobalVectorFromField(sdm, DMSwarmPICField_coor, &coordinates));
845: PetscCall(TSGetTime(ts, &time));
846: PetscCall(TSSetMaxTime(sts, time));
847: adv->tf = time;
848: PetscCall(TSSolve(sts, coordinates));
849: PetscCall(DMSwarmDestroyGlobalVectorFromField(sdm, DMSwarmPICField_coor, &coordinates));
850: PetscCall(VecCopy(adv->uf, adv->ui));
851: adv->ti = adv->tf;
853: PetscCall(DMSwarmMigrate(sdm, PETSC_TRUE));
854: PetscCall(DMSwarmGetSize(sdm, &newN));
855: PetscCall(DMSwarmGetLocalSize(sdm, &newn));
856: lreset = (n != newn || N != newN) ? PETSC_TRUE : PETSC_FALSE;
857: PetscCallMPI(MPIU_Allreduce(&lreset, &reset, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sts)));
858: if (reset) {
859: PetscCall(TSReset(sts));
860: PetscCall(DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor));
861: }
862: PetscCall(DMViewFromOptions(sdm, NULL, "-dm_view"));
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: int main(int argc, char **argv)
867: {
868: DM dm, sdm;
869: TS ts, sts;
870: Vec u, xtmp;
871: AppCtx user;
872: AdvCtx adv;
873: PetscReal t;
874: PetscInt dim;
876: PetscFunctionBeginUser;
877: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
878: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
879: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
880: PetscCall(SetupParameters(&user));
881: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
882: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
883: PetscCall(TSSetDM(ts, dm));
884: PetscCall(DMSetApplicationContext(dm, &user));
885: /* Discretize chemical species */
886: PetscCall(DMCreate(PETSC_COMM_WORLD, &sdm));
887: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sdm, "part_"));
888: PetscCall(DMSetType(sdm, DMSWARM));
889: PetscCall(DMGetDimension(dm, &dim));
890: PetscCall(DMSetDimension(sdm, dim));
891: PetscCall(DMSwarmSetCellDM(sdm, dm));
892: /* Setup problem */
893: PetscCall(SetupDiscretization(dm, sdm, &user));
894: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
896: PetscCall(DMCreateGlobalVector(dm, &u));
897: PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));
899: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
900: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
901: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
902: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
903: PetscCall(TSSetPreStep(ts, RemoveDiscretePressureNullspace));
904: PetscCall(TSMonitorSet(ts, MonitorError, &user, NULL));
905: PetscCall(TSSetFromOptions(ts));
907: PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); /* Must come after SetFromOptions() */
908: PetscCall(SetInitialConditions(ts, u));
909: PetscCall(TSGetTime(ts, &t));
910: PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
911: PetscCall(DMTSCheckFromOptions(ts, u));
913: /* Setup particle position integrator */
914: PetscCall(TSCreate(PETSC_COMM_WORLD, &sts));
915: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)sts, "part_"));
916: PetscCall(PetscObjectIncrementTabLevel((PetscObject)sts, (PetscObject)ts, 1));
917: PetscCall(TSSetDM(sts, sdm));
918: PetscCall(TSSetProblemType(sts, TS_NONLINEAR));
919: PetscCall(TSSetExactFinalTime(sts, TS_EXACTFINALTIME_MATCHSTEP));
920: PetscCall(TSMonitorSet(sts, MonitorParticleError, &adv, NULL));
921: PetscCall(TSSetFromOptions(sts));
922: PetscCall(TSSetApplicationContext(sts, &adv));
923: PetscCall(TSSetComputeExactError(sts, ComputeParticleError));
924: PetscCall(TSSetComputeInitialCondition(sts, SetInitialParticleConditions));
925: adv.ti = t;
926: adv.uf = u;
927: PetscCall(VecDuplicate(adv.uf, &adv.ui));
928: PetscCall(VecCopy(u, adv.ui));
929: PetscCall(TSSetRHSFunction(sts, NULL, FreeStreaming, &adv));
930: PetscCall(TSSetPostStep(ts, AdvectParticles));
931: PetscCall(PetscObjectCompose((PetscObject)ts, "_SwarmTS", (PetscObject)sts));
932: PetscCall(DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor));
933: PetscCall(DMCreateGlobalVector(sdm, &adv.x0));
934: PetscCall(DMSwarmCreateGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp));
935: PetscCall(VecCopy(xtmp, adv.x0));
936: PetscCall(DMSwarmDestroyGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp));
937: switch (user.solType) {
938: case SOL_TRIG_TRIG:
939: adv.exact = trig_trig_x;
940: break;
941: default:
942: SETERRQ(PetscObjectComm((PetscObject)sdm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user.solType, NUM_SOL_TYPES)], user.solType);
943: }
944: adv.ctx = &user;
946: PetscCall(TSSolve(ts, u));
947: PetscCall(DMTSCheckFromOptions(ts, u));
948: PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
950: PetscCall(VecDestroy(&u));
951: PetscCall(VecDestroy(&adv.x0));
952: PetscCall(VecDestroy(&adv.ui));
953: PetscCall(DMDestroy(&dm));
954: PetscCall(DMDestroy(&sdm));
955: PetscCall(TSDestroy(&ts));
956: PetscCall(TSDestroy(&sts));
957: PetscCall(PetscBagDestroy(&user.bag));
958: PetscCall(PetscFinalize());
959: return 0;
960: }
962: /*TEST
964: # Swarm does not work with complex
965: test:
966: suffix: 2d_tri_p2_p1_p1_tconvp
967: requires: triangle !single !complex
968: args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 2 \
969: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
970: -dmts_check .001 -ts_max_steps 4 -ts_time_step 0.1 -ts_monitor_cancel \
971: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
972: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
973: -fieldsplit_0_pc_type lu \
974: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
975: -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \
976: -part_ts_max_steps 2 -part_ts_time_step 0.05 -part_ts_convergence_estimate -convest_num_refine 1 -part_ts_monitor_cancel
977: test:
978: suffix: 2d_tri_p2_p1_p1_exit
979: requires: triangle !single !complex
980: args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 1 \
981: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
982: -dmts_check .001 -ts_max_steps 10 -ts_time_step 0.1 \
983: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
984: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
985: -fieldsplit_0_pc_type lu \
986: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
987: -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \
988: -part_ts_max_steps 20 -part_ts_time_step 0.05
990: TEST*/