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 **)&param);

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 **)&param);
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 **)&param);
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 **)&param);
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*/