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, void *ctx)
76: {
77: PetscInt d;
78: for (d = 0; d < Nc; ++d) u[d] = 0.0;
79: return 0;
80: }
82: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
83: {
84: PetscInt d;
85: for (d = 0; d < Nc; ++d) u[d] = 1.0;
86: return 0;
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, void *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 0;
127: }
128: static PetscErrorCode trig_trig_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
129: {
130: Parameter *p = (Parameter *)ctx;
132: u[0] = -p->omega * X[1];
133: u[1] = p->omega * X[0];
134: return 0;
135: }
136: static PetscErrorCode trig_trig_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
137: {
138: u[0] = 0.0;
139: u[1] = 0.0;
140: return 0;
141: }
143: static PetscErrorCode trig_trig_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
144: {
145: p[0] = X[0] + X[1] - 1.0;
146: return 0;
147: }
149: static PetscErrorCode trig_trig_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
150: {
151: T[0] = time + X[0] + X[1];
152: return 0;
153: }
154: static PetscErrorCode trig_trig_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
155: {
156: T[0] = 1.0;
157: return 0;
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;
254: PetscInt c, d;
256: for (c = 0; c < Nc; ++c) {
257: for (d = 0; d < dim; ++d) {
258: g3[((c * Nc + c) * dim + d) * dim + d] += nu; // gradU
259: g3[((c * Nc + d) * dim + d) * dim + c] += nu; // gradU transpose
260: }
261: }
262: }
264: 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[])
265: {
266: PetscInt d;
267: for (d = 0; d < dim; ++d) g0[d] = u_tShift;
268: }
270: 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[])
271: {
272: PetscInt d;
273: for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
274: }
276: 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[])
277: {
278: PetscInt d;
279: for (d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
280: }
282: 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[])
283: {
284: const PetscReal alpha = PetscRealPart(constants[1]);
285: PetscInt d;
287: for (d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
288: }
290: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
291: {
292: PetscInt sol, pl, n;
295: options->solType = SOL_TRIG_TRIG;
296: options->partLayout = PART_LAYOUT_CELL;
297: options->Npc = 1;
298: options->Npb = 1;
299: options->partLower[0] = options->partLower[1] = options->partLower[2] = 0.;
300: options->partUpper[0] = options->partUpper[1] = options->partUpper[2] = 1.;
301: PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");
302: sol = options->solType;
303: PetscOptionsEList("-sol_type", "The solution type", "ex77.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL);
304: options->solType = (SolType)sol;
305: pl = options->partLayout;
306: PetscOptionsEList("-part_layout_type", "The particle layout type", "ex77.c", partLayoutTypes, NUM_PART_LAYOUT_TYPES, partLayoutTypes[options->partLayout], &pl, NULL);
307: options->partLayout = (PartLayoutType)pl;
308: PetscOptionsInt("-Npc", "The initial number of particles per cell", "ex77.c", options->Npc, &options->Npc, NULL);
309: n = 3;
310: PetscOptionsRealArray("-part_lower", "The lower left corner of the particle box", "ex77.c", options->partLower, &n, NULL);
311: n = 3;
312: PetscOptionsRealArray("-part_upper", "The upper right corner of the particle box", "ex77.c", options->partUpper, &n, NULL);
313: PetscOptionsInt("-Npb", "The initial number of particles per box dimension", "ex77.c", options->Npb, &options->Npb, NULL);
314: PetscOptionsEnd();
315: return 0;
316: }
318: static PetscErrorCode SetupParameters(AppCtx *user)
319: {
320: PetscBag bag;
321: Parameter *p;
324: /* setup PETSc parameter bag */
325: PetscBagGetData(user->bag, (void **)&p);
326: PetscBagSetName(user->bag, "par", "Low Mach flow parameters");
327: bag = user->bag;
328: PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");
329: PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity");
330: PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature");
331: PetscBagRegisterReal(bag, &p->omega, 1.0, "omega", "Rotation speed in MMS benchmark");
332: return 0;
333: }
335: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
336: {
338: DMCreate(comm, dm);
339: DMSetType(*dm, DMPLEX);
340: DMSetFromOptions(*dm);
341: DMViewFromOptions(*dm, NULL, "-dm_view");
342: return 0;
343: }
345: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
346: {
347: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
348: PetscErrorCode (*exactFuncs_t[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
349: PetscDS prob;
350: DMLabel label;
351: Parameter *ctx;
352: PetscInt id;
355: DMGetLabel(dm, "marker", &label);
356: DMGetDS(dm, &prob);
357: switch (user->solType) {
358: case SOL_TRIG_TRIG:
359: PetscDSSetResidual(prob, 0, f0_trig_trig_v, f1_v);
360: PetscDSSetResidual(prob, 2, f0_trig_trig_w, f1_w);
362: exactFuncs[0] = trig_trig_u;
363: exactFuncs[1] = trig_trig_p;
364: exactFuncs[2] = trig_trig_T;
365: exactFuncs_t[0] = trig_trig_u_t;
366: exactFuncs_t[1] = NULL;
367: exactFuncs_t[2] = trig_trig_T_t;
368: break;
369: default:
370: SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
371: }
373: PetscDSSetResidual(prob, 1, f0_q, NULL);
375: PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu);
376: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL);
377: PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);
378: PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL);
379: PetscDSSetJacobian(prob, 2, 2, g0_wT, g1_wT, NULL, g3_wT);
380: /* Setup constants */
381: {
382: Parameter *param;
383: PetscScalar constants[4];
385: PetscBagGetData(user->bag, (void **)¶m);
387: constants[0] = param->nu;
388: constants[1] = param->alpha;
389: constants[2] = param->T_in;
390: constants[3] = param->omega;
391: PetscDSSetConstants(prob, 4, constants);
392: }
393: /* Setup Boundary Conditions */
394: PetscBagGetData(user->bag, (void **)&ctx);
395: id = 3;
396: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], (void (*)(void))exactFuncs_t[0], ctx, NULL);
397: id = 1;
398: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], (void (*)(void))exactFuncs_t[0], ctx, NULL);
399: id = 2;
400: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], (void (*)(void))exactFuncs_t[0], ctx, NULL);
401: id = 4;
402: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], (void (*)(void))exactFuncs_t[0], ctx, NULL);
403: id = 3;
404: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], (void (*)(void))exactFuncs_t[2], ctx, NULL);
405: id = 1;
406: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], (void (*)(void))exactFuncs_t[2], ctx, NULL);
407: id = 2;
408: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], (void (*)(void))exactFuncs_t[2], ctx, NULL);
409: id = 4;
410: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], (void (*)(void))exactFuncs_t[2], ctx, NULL);
412: /*setup exact solution.*/
413: PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx);
414: PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx);
415: PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx);
416: PetscDSSetExactSolutionTimeDerivative(prob, 0, exactFuncs_t[0], ctx);
417: PetscDSSetExactSolutionTimeDerivative(prob, 1, exactFuncs_t[1], ctx);
418: PetscDSSetExactSolutionTimeDerivative(prob, 2, exactFuncs_t[2], ctx);
419: return 0;
420: }
422: /* x_t = v
424: Note that here we use the velocity field at t_{n+1} to advect the particles from
425: t_n to t_{n+1}. If we use both of these fields, we could use Crank-Nicholson or
426: the method of characteristics.
427: */
428: static PetscErrorCode FreeStreaming(TS ts, PetscReal t, Vec X, Vec F, void *ctx)
429: {
430: AdvCtx *adv = (AdvCtx *)ctx;
431: Vec u = adv->ui;
432: DM sdm, dm, vdm;
433: Vec vel, locvel, pvel;
434: IS vis;
435: DMInterpolationInfo ictx;
436: const PetscScalar *coords, *v;
437: PetscScalar *f;
438: PetscInt vf[1] = {0};
439: PetscInt dim, Np;
442: TSGetDM(ts, &sdm);
443: DMSwarmGetCellDM(sdm, &dm);
444: DMGetGlobalVector(sdm, &pvel);
445: DMSwarmGetLocalSize(sdm, &Np);
446: DMGetDimension(dm, &dim);
447: /* Get local velocity */
448: DMCreateSubDM(dm, 1, vf, &vis, &vdm);
449: VecGetSubVector(u, vis, &vel);
450: DMGetLocalVector(vdm, &locvel);
451: DMPlexInsertBoundaryValues(vdm, PETSC_TRUE, locvel, adv->ti, NULL, NULL, NULL);
452: DMGlobalToLocalBegin(vdm, vel, INSERT_VALUES, locvel);
453: DMGlobalToLocalEnd(vdm, vel, INSERT_VALUES, locvel);
454: VecRestoreSubVector(u, vis, &vel);
455: ISDestroy(&vis);
456: /* Interpolate velocity */
457: DMInterpolationCreate(PETSC_COMM_SELF, &ictx);
458: DMInterpolationSetDim(ictx, dim);
459: DMInterpolationSetDof(ictx, dim);
460: VecGetArrayRead(X, &coords);
461: DMInterpolationAddPoints(ictx, Np, (PetscReal *)coords);
462: VecRestoreArrayRead(X, &coords);
463: /* Particles that lie outside the domain should be dropped,
464: whereas particles that move to another partition should trigger a migration */
465: DMInterpolationSetUp(ictx, vdm, PETSC_FALSE, PETSC_TRUE);
466: VecSet(pvel, 0.);
467: DMInterpolationEvaluate(ictx, vdm, locvel, pvel);
468: DMInterpolationDestroy(&ictx);
469: DMRestoreLocalVector(vdm, &locvel);
470: DMDestroy(&vdm);
472: VecGetArray(F, &f);
473: VecGetArrayRead(pvel, &v);
474: PetscArraycpy(f, v, Np * dim);
475: VecRestoreArrayRead(pvel, &v);
476: VecRestoreArray(F, &f);
477: DMRestoreGlobalVector(sdm, &pvel);
478: return 0;
479: }
481: static PetscErrorCode SetInitialParticleConditions(TS ts, Vec u)
482: {
483: AppCtx *user;
484: void *ctx;
485: DM dm;
486: PetscScalar *coords;
487: PetscReal x[3], dx[3];
488: PetscInt n[3];
489: PetscInt dim, d, i, j, k;
492: TSGetApplicationContext(ts, &ctx);
493: user = ((AdvCtx *)ctx)->ctx;
494: TSGetDM(ts, &dm);
495: DMGetDimension(dm, &dim);
496: switch (user->partLayout) {
497: case PART_LAYOUT_CELL:
498: DMSwarmSetPointCoordinatesRandom(dm, user->Npc);
499: break;
500: case PART_LAYOUT_BOX:
501: for (d = 0; d < dim; ++d) {
502: n[d] = user->Npb;
503: dx[d] = (user->partUpper[d] - user->partLower[d]) / PetscMax(1, n[d] - 1);
504: }
505: VecGetArray(u, &coords);
506: switch (dim) {
507: case 2:
508: x[0] = user->partLower[0];
509: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
510: x[1] = user->partLower[1];
511: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
512: const PetscInt p = j * n[0] + i;
513: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
514: }
515: }
516: break;
517: case 3:
518: x[0] = user->partLower[0];
519: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
520: x[1] = user->partLower[1];
521: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
522: x[2] = user->partLower[2];
523: for (k = 0; k < n[2]; ++k, x[2] += dx[2]) {
524: const PetscInt p = (k * n[1] + j) * n[0] + i;
525: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
526: }
527: }
528: }
529: break;
530: default:
531: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Do not support particle layout in dimension %" PetscInt_FMT, dim);
532: }
533: VecRestoreArray(u, &coords);
534: break;
535: default:
536: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]);
537: }
538: return 0;
539: }
541: static PetscErrorCode SetupDiscretization(DM dm, DM sdm, AppCtx *user)
542: {
543: DM cdm = dm;
544: PetscFE fe[3];
545: Parameter *param;
546: PetscInt *cellid, n[3];
547: PetscReal x[3], dx[3];
548: PetscScalar *coords;
549: DMPolytopeType ct;
550: PetscInt dim, d, cStart, cEnd, c, Np, p, i, j, k;
551: PetscBool simplex;
552: MPI_Comm comm;
555: DMGetDimension(dm, &dim);
556: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
557: DMPlexGetCellType(dm, cStart, &ct);
558: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
559: /* Create finite element */
560: PetscObjectGetComm((PetscObject)dm, &comm);
561: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);
562: PetscObjectSetName((PetscObject)fe[0], "velocity");
564: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
565: PetscFECopyQuadrature(fe[0], fe[1]);
566: PetscObjectSetName((PetscObject)fe[1], "pressure");
568: PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]);
569: PetscFECopyQuadrature(fe[0], fe[2]);
570: PetscObjectSetName((PetscObject)fe[2], "temperature");
572: /* Set discretization and boundary conditions for each mesh */
573: DMSetField(dm, 0, NULL, (PetscObject)fe[0]);
574: DMSetField(dm, 1, NULL, (PetscObject)fe[1]);
575: DMSetField(dm, 2, NULL, (PetscObject)fe[2]);
576: DMCreateDS(dm);
577: SetupProblem(dm, user);
578: PetscBagGetData(user->bag, (void **)¶m);
579: while (cdm) {
580: DMCopyDisc(dm, cdm);
581: DMGetCoarseDM(cdm, &cdm);
582: }
583: PetscFEDestroy(&fe[0]);
584: PetscFEDestroy(&fe[1]);
585: PetscFEDestroy(&fe[2]);
587: {
588: PetscObject pressure;
589: MatNullSpace nullspacePres;
591: DMGetField(dm, 1, NULL, &pressure);
592: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres);
593: PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres);
594: MatNullSpaceDestroy(&nullspacePres);
595: }
597: /* Setup particle information */
598: DMSwarmSetType(sdm, DMSWARM_PIC);
599: DMSwarmRegisterPetscDatatypeField(sdm, "mass", 1, PETSC_REAL);
600: DMSwarmFinalizeFieldRegister(sdm);
601: switch (user->partLayout) {
602: case PART_LAYOUT_CELL:
603: DMSwarmSetLocalSizes(sdm, (cEnd - cStart) * user->Npc, 0);
604: DMSetFromOptions(sdm);
605: DMSwarmGetField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
606: for (c = cStart; c < cEnd; ++c) {
607: for (p = 0; p < user->Npc; ++p) {
608: const PetscInt n = c * user->Npc + p;
610: cellid[n] = c;
611: }
612: }
613: DMSwarmRestoreField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
614: DMSwarmSetPointCoordinatesRandom(sdm, user->Npc);
615: break;
616: case PART_LAYOUT_BOX:
617: Np = 1;
618: for (d = 0; d < dim; ++d) {
619: n[d] = user->Npb;
620: dx[d] = (user->partUpper[d] - user->partLower[d]) / PetscMax(1, n[d] - 1);
621: Np *= n[d];
622: }
623: DMSwarmSetLocalSizes(sdm, Np, 0);
624: DMSetFromOptions(sdm);
625: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
626: switch (dim) {
627: case 2:
628: x[0] = user->partLower[0];
629: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
630: x[1] = user->partLower[1];
631: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
632: const PetscInt p = j * n[0] + i;
633: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
634: }
635: }
636: break;
637: case 3:
638: x[0] = user->partLower[0];
639: for (i = 0; i < n[0]; ++i, x[0] += dx[0]) {
640: x[1] = user->partLower[1];
641: for (j = 0; j < n[1]; ++j, x[1] += dx[1]) {
642: x[2] = user->partLower[2];
643: for (k = 0; k < n[2]; ++k, x[2] += dx[2]) {
644: const PetscInt p = (k * n[1] + j) * n[0] + i;
645: for (d = 0; d < dim; ++d) coords[p * dim + d] = x[d];
646: }
647: }
648: }
649: break;
650: default:
651: SETERRQ(comm, PETSC_ERR_SUP, "Do not support particle layout in dimension %" PetscInt_FMT, dim);
652: }
653: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
654: DMSwarmGetField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
655: for (p = 0; p < Np; ++p) cellid[p] = 0;
656: DMSwarmRestoreField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid);
657: DMSwarmMigrate(sdm, PETSC_TRUE);
658: break;
659: default:
660: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]);
661: }
662: PetscObjectSetName((PetscObject)sdm, "Particles");
663: DMViewFromOptions(sdm, NULL, "-dm_view");
664: return 0;
665: }
667: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
668: {
669: Vec vec;
670: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
674: funcs[nfield] = constant;
675: DMCreateGlobalVector(dm, &vec);
676: DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);
677: VecNormalize(vec, NULL);
678: PetscObjectSetName((PetscObject)vec, "Pressure Null Space");
679: VecViewFromOptions(vec, NULL, "-pressure_nullspace_view");
680: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
681: VecDestroy(&vec);
682: return 0;
683: }
685: static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u)
686: {
687: DM dm;
688: MatNullSpace nullsp;
691: TSGetDM(ts, &dm);
692: CreatePressureNullSpace(dm, 1, 1, &nullsp);
693: MatNullSpaceRemove(nullsp, u);
694: MatNullSpaceDestroy(&nullsp);
695: return 0;
696: }
698: /* Make the discrete pressure discretely divergence free */
699: static PetscErrorCode RemoveDiscretePressureNullspace(TS ts)
700: {
701: Vec u;
704: TSGetSolution(ts, &u);
705: RemoveDiscretePressureNullspace_Private(ts, u);
706: return 0;
707: }
709: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
710: {
711: DM dm;
712: PetscReal t;
715: TSGetDM(ts, &dm);
716: TSGetTime(ts, &t);
717: DMComputeExactSolution(dm, t, u, NULL);
718: RemoveDiscretePressureNullspace_Private(ts, u);
719: return 0;
720: }
722: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
723: {
724: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
725: void *ctxs[3];
726: DM dm;
727: PetscDS ds;
728: Vec v;
729: PetscReal ferrors[3];
730: PetscInt tl, l, f;
733: TSGetDM(ts, &dm);
734: DMGetDS(dm, &ds);
736: for (f = 0; f < 3; ++f) PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]);
737: DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors);
738: PetscObjectGetTabLevel((PetscObject)ts, &tl);
739: for (l = 0; l < tl; ++l) PetscPrintf(PETSC_COMM_WORLD, "\t");
740: 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]);
742: DMGetGlobalVector(dm, &u);
743: PetscObjectSetName((PetscObject)u, "Numerical Solution");
744: VecViewFromOptions(u, NULL, "-sol_vec_view");
745: DMRestoreGlobalVector(dm, &u);
747: DMGetGlobalVector(dm, &v);
748: DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, v);
749: PetscObjectSetName((PetscObject)v, "Exact Solution");
750: VecViewFromOptions(v, NULL, "-exact_vec_view");
751: DMRestoreGlobalVector(dm, &v);
753: return 0;
754: }
756: /* Note that adv->x0 will not be correct after migration */
757: static PetscErrorCode ComputeParticleError(TS ts, Vec u, Vec e)
758: {
759: AdvCtx *adv;
760: DM sdm;
761: Parameter *param;
762: const PetscScalar *xp0, *xp;
763: PetscScalar *ep;
764: PetscReal time;
765: PetscInt dim, Np, p;
766: MPI_Comm comm;
769: TSGetTime(ts, &time);
770: TSGetApplicationContext(ts, &adv);
771: PetscBagGetData(adv->ctx->bag, (void **)¶m);
772: PetscObjectGetComm((PetscObject)ts, &comm);
773: TSGetDM(ts, &sdm);
774: DMGetDimension(sdm, &dim);
775: DMSwarmGetLocalSize(sdm, &Np);
776: VecGetArrayRead(adv->x0, &xp0);
777: VecGetArrayRead(u, &xp);
778: VecGetArrayWrite(e, &ep);
779: for (p = 0; p < Np; ++p) {
780: PetscScalar x[3];
781: PetscReal x0[3];
782: PetscInt d;
784: for (d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p * dim + d]);
785: adv->exact(dim, time, x0, 1, x, param);
786: for (d = 0; d < dim; ++d) ep[p * dim + d] += x[d] - xp[p * dim + d];
787: }
788: VecRestoreArrayRead(adv->x0, &xp0);
789: VecRestoreArrayRead(u, &xp);
790: VecRestoreArrayWrite(e, &ep);
791: return 0;
792: }
794: static PetscErrorCode MonitorParticleError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
795: {
796: AdvCtx *adv = (AdvCtx *)ctx;
797: DM sdm;
798: Parameter *param;
799: const PetscScalar *xp0, *xp;
800: PetscReal error = 0.0;
801: PetscInt dim, tl, l, Np, p;
802: MPI_Comm comm;
805: PetscBagGetData(adv->ctx->bag, (void **)¶m);
806: PetscObjectGetComm((PetscObject)ts, &comm);
807: TSGetDM(ts, &sdm);
808: DMGetDimension(sdm, &dim);
809: DMSwarmGetLocalSize(sdm, &Np);
810: VecGetArrayRead(adv->x0, &xp0);
811: VecGetArrayRead(u, &xp);
812: for (p = 0; p < Np; ++p) {
813: PetscScalar x[3];
814: PetscReal x0[3];
815: PetscReal perror = 0.0;
816: PetscInt d;
818: for (d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p * dim + d]);
819: adv->exact(dim, time, x0, 1, x, param);
820: for (d = 0; d < dim; ++d) perror += PetscSqr(PetscRealPart(x[d] - xp[p * dim + d]));
821: error += perror;
822: }
823: VecRestoreArrayRead(adv->x0, &xp0);
824: VecRestoreArrayRead(u, &xp);
825: PetscObjectGetTabLevel((PetscObject)ts, &tl);
826: for (l = 0; l < tl; ++l) PetscPrintf(PETSC_COMM_WORLD, "\t");
827: PetscPrintf(comm, "Timestep: %04d time = %-8.4g \t L_2 Particle Error: [%2.3g]\n", (int)step, (double)time, (double)error);
828: return 0;
829: }
831: static PetscErrorCode AdvectParticles(TS ts)
832: {
833: TS sts;
834: DM sdm;
835: Vec coordinates;
836: AdvCtx *adv;
837: PetscReal time;
838: PetscBool lreset, reset;
839: PetscInt dim, n, N, newn, newN;
842: PetscObjectQuery((PetscObject)ts, "_SwarmTS", (PetscObject *)&sts);
843: TSGetDM(sts, &sdm);
844: TSGetRHSFunction(sts, NULL, NULL, (void **)&adv);
845: DMGetDimension(sdm, &dim);
846: DMSwarmGetSize(sdm, &N);
847: DMSwarmGetLocalSize(sdm, &n);
848: DMSwarmCreateGlobalVectorFromField(sdm, DMSwarmPICField_coor, &coordinates);
849: TSGetTime(ts, &time);
850: TSSetMaxTime(sts, time);
851: adv->tf = time;
852: TSSolve(sts, coordinates);
853: DMSwarmDestroyGlobalVectorFromField(sdm, DMSwarmPICField_coor, &coordinates);
854: VecCopy(adv->uf, adv->ui);
855: adv->ti = adv->tf;
857: DMSwarmMigrate(sdm, PETSC_TRUE);
858: DMSwarmGetSize(sdm, &newN);
859: DMSwarmGetLocalSize(sdm, &newn);
860: lreset = (n != newn || N != newN) ? PETSC_TRUE : PETSC_FALSE;
861: MPI_Allreduce(&lreset, &reset, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sts));
862: if (reset) {
863: TSReset(sts);
864: DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);
865: }
866: DMViewFromOptions(sdm, NULL, "-dm_view");
867: return 0;
868: }
870: int main(int argc, char **argv)
871: {
872: DM dm, sdm;
873: TS ts, sts;
874: Vec u, xtmp;
875: AppCtx user;
876: AdvCtx adv;
877: PetscReal t;
878: PetscInt dim;
881: PetscInitialize(&argc, &argv, NULL, help);
882: ProcessOptions(PETSC_COMM_WORLD, &user);
883: PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);
884: SetupParameters(&user);
885: TSCreate(PETSC_COMM_WORLD, &ts);
886: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
887: TSSetDM(ts, dm);
888: DMSetApplicationContext(dm, &user);
889: /* Discretize chemical species */
890: DMCreate(PETSC_COMM_WORLD, &sdm);
891: PetscObjectSetOptionsPrefix((PetscObject)sdm, "part_");
892: DMSetType(sdm, DMSWARM);
893: DMGetDimension(dm, &dim);
894: DMSetDimension(sdm, dim);
895: DMSwarmSetCellDM(sdm, dm);
896: /* Setup problem */
897: SetupDiscretization(dm, sdm, &user);
898: DMPlexCreateClosureIndex(dm, NULL);
900: DMCreateGlobalVector(dm, &u);
901: DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace);
903: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user);
904: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user);
905: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user);
906: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
907: TSSetPreStep(ts, RemoveDiscretePressureNullspace);
908: TSMonitorSet(ts, MonitorError, &user, NULL);
909: TSSetFromOptions(ts);
911: TSSetComputeInitialCondition(ts, SetInitialConditions); /* Must come after SetFromOptions() */
912: SetInitialConditions(ts, u);
913: TSGetTime(ts, &t);
914: DMSetOutputSequenceNumber(dm, 0, t);
915: DMTSCheckFromOptions(ts, u);
917: /* Setup particle position integrator */
918: TSCreate(PETSC_COMM_WORLD, &sts);
919: PetscObjectSetOptionsPrefix((PetscObject)sts, "part_");
920: PetscObjectIncrementTabLevel((PetscObject)sts, (PetscObject)ts, 1);
921: TSSetDM(sts, sdm);
922: TSSetProblemType(sts, TS_NONLINEAR);
923: TSSetExactFinalTime(sts, TS_EXACTFINALTIME_MATCHSTEP);
924: TSMonitorSet(sts, MonitorParticleError, &adv, NULL);
925: TSSetFromOptions(sts);
926: TSSetApplicationContext(sts, &adv);
927: TSSetComputeExactError(sts, ComputeParticleError);
928: TSSetComputeInitialCondition(sts, SetInitialParticleConditions);
929: adv.ti = t;
930: adv.uf = u;
931: VecDuplicate(adv.uf, &adv.ui);
932: VecCopy(u, adv.ui);
933: TSSetRHSFunction(sts, NULL, FreeStreaming, &adv);
934: TSSetPostStep(ts, AdvectParticles);
935: PetscObjectCompose((PetscObject)ts, "_SwarmTS", (PetscObject)sts);
936: DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);
937: DMCreateGlobalVector(sdm, &adv.x0);
938: DMSwarmCreateGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp);
939: VecCopy(xtmp, adv.x0);
940: DMSwarmDestroyGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp);
941: switch (user.solType) {
942: case SOL_TRIG_TRIG:
943: adv.exact = trig_trig_x;
944: break;
945: default:
946: SETERRQ(PetscObjectComm((PetscObject)sdm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user.solType, NUM_SOL_TYPES)], user.solType);
947: }
948: adv.ctx = &user;
950: TSSolve(ts, u);
951: DMTSCheckFromOptions(ts, u);
952: PetscObjectSetName((PetscObject)u, "Numerical Solution");
954: VecDestroy(&u);
955: VecDestroy(&adv.x0);
956: VecDestroy(&adv.ui);
957: DMDestroy(&dm);
958: DMDestroy(&sdm);
959: TSDestroy(&ts);
960: TSDestroy(&sts);
961: PetscBagDestroy(&user.bag);
962: PetscFinalize();
963: return 0;
964: }
966: /*TEST
968: # Swarm does not work with complex
969: test:
970: suffix: 2d_tri_p2_p1_p1_tconvp
971: requires: triangle !single !complex
972: args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 2 \
973: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
974: -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 -ts_monitor_cancel \
975: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
976: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
977: -fieldsplit_0_pc_type lu \
978: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
979: -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \
980: -part_ts_max_steps 2 -part_ts_dt 0.05 -part_ts_convergence_estimate -convest_num_refine 1 -part_ts_monitor_cancel
981: test:
982: suffix: 2d_tri_p2_p1_p1_exit
983: requires: triangle !single !complex
984: args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 1 \
985: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
986: -dmts_check .001 -ts_max_steps 10 -ts_dt 0.1 \
987: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
988: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
989: -fieldsplit_0_pc_type lu \
990: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \
991: -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \
992: -part_ts_max_steps 20 -part_ts_dt 0.05
994: TEST*/