Actual source code: ex48.c

  1: static char help[] = "Evolution of magnetic islands.\n\
  2: The aim of this model is to self-consistently study the interaction between the tearing mode and small scale drift-wave turbulence.\n\n\n";

  4: /*F
  5: This is a three field model for the density $\tilde n$, vorticity $\tilde\Omega$, and magnetic flux $\tilde\psi$, using auxiliary variables potential $\tilde\phi$ and current $j_z$.
  6: \begin{equation}
  7:   \begin{aligned}
  8:     \partial_t \tilde n       &= \left\{ \tilde n, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \left\{ \ln n_0, \tilde\phi \right\} + \mu \nabla^2_\perp \tilde n \\
  9:   \partial_t \tilde\Omega   &= \left\{ \tilde\Omega, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \mu \nabla^2_\perp \tilde\Omega \\
 10:   \partial_t \tilde\psi     &= \left\{ \psi_0 + \tilde\psi, \tilde\phi - \tilde n \right\} - \left\{ \ln n_0, \tilde\psi \right\} + \frac{\eta}{\beta} \nabla^2_\perp \tilde\psi \\
 11:   \nabla^2_\perp\tilde\phi        &= \tilde\Omega \\
 12:   j_z  &= -\nabla^2_\perp  \left(\tilde\psi + \psi_0  \right)\\
 13:   \end{aligned}
 14: \end{equation}
 15: F*/

 17: #include <petscdmplex.h>
 18: #include <petscts.h>
 19: #include <petscds.h>

 21: typedef struct {
 22:   PetscInt  debug;   /* The debugging level */
 23:   PetscBool plotRef; /* Plot the reference fields */
 24:   PetscReal lower[3], upper[3];
 25:   /* Problem definition */
 26:   PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 27:   PetscReal mu, eta, beta;
 28:   PetscReal a, b, Jo, Jop, m, ke, kx, ky, DeltaPrime, eps;
 29:   /* solver */
 30:   PetscBool implicit;
 31: } AppCtx;

 33: static AppCtx *s_ctx;

 35: static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
 36: {
 37:   PetscScalar ret = df[0] * dg[1] - df[1] * dg[0];
 38:   return ret;
 39: }

 41: enum field_idx {
 42:   DENSITY,
 43:   OMEGA,
 44:   PSI,
 45:   PHI,
 46:   JZ
 47: };

 49: static void f0_n(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[])
 50: {
 51:   const PetscScalar *pnDer        = &u_x[uOff_x[DENSITY]];
 52:   const PetscScalar *ppsiDer      = &u_x[uOff_x[PSI]];
 53:   const PetscScalar *pphiDer      = &u_x[uOff_x[PHI]];
 54:   const PetscScalar *jzDer        = &u_x[uOff_x[JZ]];
 55:   const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
 56:   f0[0] += -poissonBracket(dim, pnDer, pphiDer) - s_ctx->beta * poissonBracket(dim, jzDer, ppsiDer) - poissonBracket(dim, logRefDenDer, pphiDer);
 57:   if (u_t) f0[0] += u_t[DENSITY];
 58: }

 60: static void f1_n(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[])
 61: {
 62:   const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
 63:   PetscInt           d;

 65:   for (d = 0; d < dim - 1; ++d) f1[d] = -s_ctx->mu * pnDer[d];
 66: }

 68: static void f0_Omega(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[])
 69: {
 70:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 71:   const PetscScalar *ppsiDer   = &u_x[uOff_x[PSI]];
 72:   const PetscScalar *pphiDer   = &u_x[uOff_x[PHI]];
 73:   const PetscScalar *jzDer     = &u_x[uOff_x[JZ]];

 75:   f0[0] += -poissonBracket(dim, pOmegaDer, pphiDer) - s_ctx->beta * poissonBracket(dim, jzDer, ppsiDer);
 76:   if (u_t) f0[0] += u_t[OMEGA];
 77: }

 79: static void f1_Omega(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[])
 80: {
 81:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 82:   PetscInt           d;

 84:   for (d = 0; d < dim - 1; ++d) f1[d] = -s_ctx->mu * pOmegaDer[d];
 85: }

 87: static void f0_psi(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[])
 88: {
 89:   const PetscScalar *pnDer        = &u_x[uOff_x[DENSITY]];
 90:   const PetscScalar *ppsiDer      = &u_x[uOff_x[PSI]];
 91:   const PetscScalar *pphiDer      = &u_x[uOff_x[PHI]];
 92:   const PetscScalar *refPsiDer    = &a_x[aOff_x[PSI]];
 93:   const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
 94:   PetscScalar        psiDer[3];
 95:   PetscScalar        phi_n_Der[3];
 96:   PetscInt           d;
 97:   if (dim < 2) {
 98:     MPI_Abort(MPI_COMM_WORLD, 1);
 99:     return;
100:   } /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */
101:   for (d = 0; d < dim; ++d) {
102:     psiDer[d]    = refPsiDer[d] + ppsiDer[d];
103:     phi_n_Der[d] = pphiDer[d] - pnDer[d];
104:   }
105:   f0[0] = -poissonBracket(dim, psiDer, phi_n_Der) + poissonBracket(dim, logRefDenDer, ppsiDer);
106:   if (u_t) f0[0] += u_t[PSI];
107: }

109: static void f1_psi(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[])
110: {
111:   const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
112:   PetscInt           d;

114:   for (d = 0; d < dim - 1; ++d) f1[d] = -(s_ctx->eta / s_ctx->beta) * ppsi[d];
115: }

117: static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
118: {
119:   f0[0] = -u[uOff[OMEGA]];
120: }

122: static void f1_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
123: {
124:   const PetscScalar *pphi = &u_x[uOff_x[PHI]];
125:   PetscInt           d;

127:   for (d = 0; d < dim - 1; ++d) f1[d] = pphi[d];
128: }

130: static void f0_jz(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[])
131: {
132:   f0[0] = u[uOff[JZ]];
133: }

135: static void f1_jz(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[])
136: {
137:   const PetscScalar *ppsi      = &u_x[uOff_x[PSI]];
138:   const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */
139:   PetscInt           d;

141:   for (d = 0; d < dim - 1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
142: }

144: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
145: {
147:   options->debug    = 1;
148:   options->plotRef  = PETSC_FALSE;
149:   options->implicit = PETSC_FALSE;
150:   options->mu       = 0;
151:   options->eta      = 0;
152:   options->beta     = 1;
153:   options->a        = 1;
154:   options->b        = PETSC_PI;
155:   options->Jop      = 0;
156:   options->m        = 1;
157:   options->eps      = 1.e-6;

159:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
160:   PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);
161:   PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);
162:   PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);
163:   PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);
164:   PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);
165:   PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);
166:   PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);
167:   PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);
168:   PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);
169:   PetscOptionsEnd();
170:   options->ke = PetscSqrtScalar(options->Jop);
171:   if (options->Jop == 0.0) {
172:     options->Jo = 1.0 / PetscPowScalar(options->a, 2);
173:   } else {
174:     options->Jo = options->Jop * PetscCosReal(options->ke * options->a) / (1.0 - PetscCosReal(options->ke * options->a));
175:   }
176:   options->ky = PETSC_PI * options->m / options->b;
177:   if (PetscPowReal(options->ky, 2) < options->Jop) {
178:     options->kx         = PetscSqrtScalar(options->Jop - PetscPowScalar(options->ky, 2));
179:     options->DeltaPrime = -2.0 * options->kx * options->a * PetscCosReal(options->kx * options->a) / PetscSinReal(options->kx * options->a);
180:   } else if (PetscPowReal(options->ky, 2) > options->Jop) {
181:     options->kx         = PetscSqrtScalar(PetscPowScalar(options->ky, 2) - options->Jop);
182:     options->DeltaPrime = -2.0 * options->kx * options->a * PetscCoshReal(options->kx * options->a) / PetscSinhReal(options->kx * options->a);
183:   } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/
184:     options->kx         = 0;
185:     options->DeltaPrime = -2.0;
186:   }
187:   PetscPrintf(comm, "DeltaPrime=%g\n", (double)options->DeltaPrime);

189:   return 0;
190: }

192: static void f_n(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)
193: {
194:   const PetscScalar *pn = &u[uOff[DENSITY]];
195:   *f0                   = *pn;
196: }

198: static PetscErrorCode PostStep(TS ts)
199: {
200:   DM       dm;
201:   AppCtx  *ctx;
202:   PetscInt stepi, num;
203:   Vec      X;

206:   TSGetApplicationContext(ts, &ctx);
207:   if (ctx->debug < 1) return 0;
208:   TSGetSolution(ts, &X);
209:   VecGetDM(X, &dm);
210:   TSGetStepNumber(ts, &stepi);
211:   DMGetOutputSequenceNumber(dm, &num, NULL);
212:   if (num < 0) DMSetOutputSequenceNumber(dm, 0, 0.0);
213:   PetscObjectSetName((PetscObject)X, "u");
214:   VecViewFromOptions(X, NULL, "-vec_view");
215:   /* print integrals */
216:   {
217:     PetscDS     prob;
218:     DM          plex;
219:     PetscScalar den, tt[5];
220:     DMConvert(dm, DMPLEX, &plex);
221:     DMGetDS(plex, &prob);
222:     PetscDSSetObjective(prob, 0, &f_n);
223:     DMPlexComputeIntegralFEM(plex, X, tt, ctx);
224:     den = tt[0];
225:     DMDestroy(&plex);
226:     PetscPrintf(PetscObjectComm((PetscObject)dm), "%" PetscInt_FMT ") total perturbed mass = %g\n", stepi, (double)PetscRealPart(den));
227:   }
228:   return 0;
229: }

231: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
232: {
234:   DMCreate(comm, dm);
235:   DMSetType(*dm, DMPLEX);
236:   DMSetFromOptions(*dm);
237:   DMViewFromOptions(*dm, NULL, "-dm_view");

239:   DMGetBoundingBox(*dm, ctx->lower, ctx->upper);
240:   ctx->a = (ctx->upper[0] - ctx->lower[0]) / 2.0;
241:   ctx->b = (ctx->upper[1] - ctx->lower[1]) / 2.0;
242:   return 0;
243: }

245: static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
246: {
247:   AppCtx *lctx = (AppCtx *)ctx;
248:   u[0]         = 2. * lctx->a + x[0];
249:   return 0;
250: }

252: static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
253: {
254:   u[0] = 0.0;
255:   return 0;
256: }

258: static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
259: {
260:   AppCtx *lctx = (AppCtx *)ctx;
261:   /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z.  The stability
262:      is analytically known and reported in ProcessOptions. */
263:   if (lctx->ke != 0.0) {
264:     u[0] = (PetscCosReal(lctx->ke * (x[0] - lctx->a)) - PetscCosReal(lctx->ke * lctx->a)) / (1.0 - PetscCosReal(lctx->ke * lctx->a));
265:   } else {
266:     u[0] = 1.0 - PetscPowScalar((x[0] - lctx->a) / lctx->a, 2);
267:   }
268:   return 0;
269: }

271: static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
272: {
273:   u[0] = 0.0;
274:   return 0;
275: }

277: static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
278: {
279:   u[0] = 0.0;
280:   return 0;
281: }

283: static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
284: {
285:   AppCtx     *ctx = (AppCtx *)a_ctx;
286:   PetscScalar r   = ctx->eps * (PetscScalar)(rand()) / (PetscScalar)(RAND_MAX);
287:   if (x[0] == ctx->lower[0] || x[0] == ctx->upper[0]) r = 0;
288:   u[0] = r;
289:   /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
290:   return 0;
291: }

293: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
294: {
295:   u[0] = 0.0;
296:   return 0;
297: }

299: static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
300: {
301:   u[0] = 0.0;
302:   return 0;
303: }

305: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
306: {
307:   PetscDS        ds;
308:   DMLabel        label;
309:   const PetscInt id = 1;

312:   DMGetLabel(dm, "marker", &label);
313:   DMGetDS(dm, &ds);
314:   PetscDSSetResidual(ds, 0, f0_n, f1_n);
315:   PetscDSSetResidual(ds, 1, f0_Omega, f1_Omega);
316:   PetscDSSetResidual(ds, 2, f0_psi, f1_psi);
317:   PetscDSSetResidual(ds, 3, f0_phi, f1_phi);
318:   PetscDSSetResidual(ds, 4, f0_jz, f1_jz);
319:   ctx->initialFuncs[0] = initialSolution_n;
320:   ctx->initialFuncs[1] = initialSolution_Omega;
321:   ctx->initialFuncs[2] = initialSolution_psi;
322:   ctx->initialFuncs[3] = initialSolution_phi;
323:   ctx->initialFuncs[4] = initialSolution_jz;
324:   for (PetscInt f = 0; f < 5; ++f) {
325:     PetscDSSetImplicit(ds, f, ctx->implicit);
326:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, f, 0, NULL, (void (*)(void))ctx->initialFuncs[f], NULL, ctx, NULL);
327:   }
328:   PetscDSSetContext(ds, 0, ctx);
329:   return 0;
330: }

332: static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
333: {
334:   PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {log_n_0, Omega_0, psi_0};
335:   Vec     eq;
336:   AppCtx *ctxarr[3];

338:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
340:   DMCreateLocalVector(dmAux, &eq);
341:   DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);
342:   DMSetAuxiliaryVec(dm, NULL, 0, 0, eq);
343:   if (ctx->plotRef) { /* plot reference functions */
344:     PetscViewer viewer = NULL;
345:     PetscBool   isHDF5, isVTK;
346:     char        buf[256];
347:     Vec         global;
348:     PetscInt    dim;

350:     DMGetDimension(dm, &dim);
351:     DMCreateGlobalVector(dmAux, &global);
352:     VecSet(global, .0); /* BCs! */
353:     DMLocalToGlobalBegin(dmAux, eq, INSERT_VALUES, global);
354:     DMLocalToGlobalEnd(dmAux, eq, INSERT_VALUES, global);
355:     PetscViewerCreate(PetscObjectComm((PetscObject)dmAux), &viewer);
356: #ifdef PETSC_HAVE_HDF5
357:     PetscViewerSetType(viewer, PETSCVIEWERHDF5);
358: #else
359:     PetscViewerSetType(viewer, PETSCVIEWERVTK);
360: #endif
361:     PetscViewerSetFromOptions(viewer);
362:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &isHDF5);
363:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isVTK);
364:     if (isHDF5) {
365:       PetscSNPrintf(buf, 256, "uEquilibrium-%" PetscInt_FMT "D.h5", dim);
366:     } else if (isVTK) {
367:       PetscSNPrintf(buf, 256, "uEquilibrium-%" PetscInt_FMT "D.vtu", dim);
368:       PetscViewerPushFormat(viewer, PETSC_VIEWER_VTK_VTU);
369:     }
370:     PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
371:     PetscViewerFileSetName(viewer, buf);
372:     if (isHDF5) DMView(dmAux, viewer);
373:     /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
374:     PetscObjectSetName((PetscObject)global, "u0");
375:     VecView(global, viewer);
376:     PetscViewerDestroy(&viewer);
377:     VecDestroy(&global);
378:   }
379:   VecDestroy(&eq);
380:   return 0;
381: }

383: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
384: {
385:   DM       dmAux, coordDM;
386:   PetscInt f;

389:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
390:   DMGetCoordinateDM(dm, &coordDM);
391:   if (!feAux) return 0;
392:   DMClone(dm, &dmAux);
393:   DMSetCoordinateDM(dmAux, coordDM);
394:   for (f = 0; f < NfAux; ++f) DMSetField(dmAux, f, NULL, (PetscObject)feAux[f]);
395:   DMCreateDS(dmAux);
396:   SetupEquilibriumFields(dm, dmAux, user);
397:   DMDestroy(&dmAux);
398:   return 0;
399: }

401: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
402: {
403:   DM        cdm = dm;
404:   PetscFE   fe[5], feAux[3];
405:   PetscInt  dim, Nf = 5, NfAux = 3, f;
406:   PetscBool simplex;
407:   MPI_Comm  comm;

410:   /* Create finite element */
411:   PetscObjectGetComm((PetscObject)dm, &comm);
412:   DMGetDimension(dm, &dim);
413:   DMPlexIsSimplex(dm, &simplex);
414:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[0]);
415:   PetscObjectSetName((PetscObject)fe[0], "density");
416:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[1]);
417:   PetscObjectSetName((PetscObject)fe[1], "vorticity");
418:   PetscFECopyQuadrature(fe[0], fe[1]);
419:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[2]);
420:   PetscObjectSetName((PetscObject)fe[2], "flux");
421:   PetscFECopyQuadrature(fe[0], fe[2]);
422:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[3]);
423:   PetscObjectSetName((PetscObject)fe[3], "potential");
424:   PetscFECopyQuadrature(fe[0], fe[3]);
425:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[4]);
426:   PetscObjectSetName((PetscObject)fe[4], "current");
427:   PetscFECopyQuadrature(fe[0], fe[4]);

429:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[0]);
430:   PetscObjectSetName((PetscObject)feAux[0], "n_0");
431:   PetscFECopyQuadrature(fe[0], feAux[0]);
432:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[1]);
433:   PetscObjectSetName((PetscObject)feAux[1], "vorticity_0");
434:   PetscFECopyQuadrature(fe[0], feAux[1]);
435:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[2]);
436:   PetscObjectSetName((PetscObject)feAux[2], "flux_0");
437:   PetscFECopyQuadrature(fe[0], feAux[2]);
438:   /* Set discretization and boundary conditions for each mesh */
439:   for (f = 0; f < Nf; ++f) DMSetField(dm, f, NULL, (PetscObject)fe[f]);
440:   DMCreateDS(dm);
441:   SetupProblem(dm, ctx);
442:   while (cdm) {
443:     SetupAuxDM(dm, NfAux, feAux, ctx);
444:     DMCopyDisc(dm, cdm);
445:     DMGetCoarseDM(cdm, &cdm);
446:   }
447:   for (f = 0; f < Nf; ++f) PetscFEDestroy(&fe[f]);
448:   for (f = 0; f < NfAux; ++f) PetscFEDestroy(&feAux[f]);
449:   return 0;
450: }

452: int main(int argc, char **argv)
453: {
454:   DM        dm;
455:   TS        ts;
456:   Vec       u, r;
457:   AppCtx    ctx;
458:   PetscReal t       = 0.0;
459:   PetscReal L2error = 0.0;
460:   AppCtx   *ctxarr[5];

462:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
463:   s_ctx                                                     = &ctx;
465:   PetscInitialize(&argc, &argv, NULL, help);
466:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
467:   /* create mesh and problem */
468:   CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
469:   DMSetApplicationContext(dm, &ctx);
470:   PetscMalloc1(5, &ctx.initialFuncs);
471:   SetupDiscretization(dm, &ctx);
472:   DMCreateGlobalVector(dm, &u);
473:   PetscObjectSetName((PetscObject)u, "u");
474:   VecDuplicate(u, &r);
475:   PetscObjectSetName((PetscObject)r, "r");
476:   /* create TS */
477:   TSCreate(PETSC_COMM_WORLD, &ts);
478:   TSSetDM(ts, dm);
479:   TSSetApplicationContext(ts, &ctx);
480:   DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
481:   if (ctx.implicit) {
482:     DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
483:     DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
484:   } else {
485:     DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);
486:   }
487:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
488:   TSSetFromOptions(ts);
489:   TSSetPostStep(ts, PostStep);
490:   /* make solution & solve */
491:   DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);
492:   TSSetSolution(ts, u);
493:   DMViewFromOptions(dm, NULL, "-dm_view");
494:   PostStep(ts); /* print the initial state */
495:   TSSolve(ts, u);
496:   TSGetTime(ts, &t);
497:   DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);
498:   if (L2error < 1.0e-11) PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");
499:   else PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)L2error);
500:   VecDestroy(&u);
501:   VecDestroy(&r);
502:   TSDestroy(&ts);
503:   DMDestroy(&dm);
504:   PetscFree(ctx.initialFuncs);
505:   PetscFinalize();
506:   return 0;
507: }

509: /*TEST

511:   test:
512:     suffix: 0
513:     args: -debug 1 -dm_refine 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_box_bd periodic,none -dm_plex_box_upper 2.0,6.283185307179586 \
514:           -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
515:   test:
516:     # Remapping with periodicity is broken
517:     suffix: 1
518:     args: -debug 1 -dm_plex_shape cylinder -dm_plex_dim 3 -dm_refine 1 -dm_refine_remap 0 -dm_plex_cylinder_bd periodic -dm_plex_boundary_label marker \
519:            -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0

521: TEST*/