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";

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$.
\begin{equation}
\begin{aligned}
\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 \\
\partial_t \tilde\Omega &= \left\{ \tilde\Omega, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \mu \nabla^2_\perp \tilde\Omega \\
\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 \\
\nabla^2_\perp\tilde\phi &= \tilde\Omega \\
j_z &= -\nabla^2_\perp \left(\tilde\psi + \psi_0 \right)\\
\end{aligned}
\end{equation}
 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 {DENSITY,OMEGA,PSI,PHI,JZ};

 43: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 44:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 45:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 46:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 47: {
 48:   const PetscScalar *pnDer   = &u_x[uOff_x[DENSITY]];
 49:   const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
 50:   const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
 51:   const PetscScalar *jzDer   = &u_x[uOff_x[JZ]];
 52:   const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
 53:   f0[0] += - poissonBracket(dim,pnDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer) - poissonBracket(dim,logRefDenDer, pphiDer);
 54:   if (u_t) f0[0] += u_t[DENSITY];
 55: }

 57: static void f1_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 58:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 59:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 60:                  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,
 69:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 70:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 71:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 72: {
 73:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 74:   const PetscScalar *ppsiDer   = &u_x[uOff_x[PSI]];
 75:   const PetscScalar *pphiDer   = &u_x[uOff_x[PHI]];
 76:   const PetscScalar *jzDer     = &u_x[uOff_x[JZ]];

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

 82: static void f1_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 83:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 84:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 85:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 86: {
 87:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 88:   PetscInt           d;

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

 93: static void f0_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 94:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 95:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 96:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 97: {
 98:   const PetscScalar *pnDer     = &u_x[uOff_x[DENSITY]];
 99:   const PetscScalar *ppsiDer   = &u_x[uOff_x[PSI]];
100:   const PetscScalar *pphiDer   = &u_x[uOff_x[PHI]];
101:   const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]];
102:   const PetscScalar *logRefDenDer= &a_x[aOff_x[DENSITY]];
103:   PetscScalar       psiDer[3];
104:   PetscScalar       phi_n_Der[3];
105:   PetscInt          d;
106:   if (dim < 2) {MPI_Abort(MPI_COMM_WORLD,1); return;} /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */
107:   for (d = 0; d < dim; ++d) {
108:     psiDer[d]    = refPsiDer[d] + ppsiDer[d];
109:     phi_n_Der[d] = pphiDer[d]   - pnDer[d];
110:   }
111:   f0[0] = - poissonBracket(dim,psiDer, phi_n_Der) + poissonBracket(dim,logRefDenDer, ppsiDer);
112:   if (u_t) f0[0] += u_t[PSI];
113: }

115: static void f1_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
116:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
117:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
118:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
119: {
120:   const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
121:   PetscInt           d;

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

126: static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
127:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
128:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
129:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
130: {
131:   f0[0] = -u[uOff[OMEGA]];
132: }

134: static void f1_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
135:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
136:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
137:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
138: {
139:   const PetscScalar *pphi = &u_x[uOff_x[PHI]];
140:   PetscInt           d;

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

145: static void f0_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
146:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
147:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
148:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
149: {
150:   f0[0] = u[uOff[JZ]];
151: }

153: static void f1_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
154:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
155:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
156:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
157: {
158:   const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
159:   const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */
160:   PetscInt           d;

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

165: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
166: {

170:   options->debug    = 1;
171:   options->plotRef  = PETSC_FALSE;
172:   options->implicit = PETSC_FALSE;
173:   options->mu       = 0;
174:   options->eta      = 0;
175:   options->beta     = 1;
176:   options->a        = 1;
177:   options->b        = PETSC_PI;
178:   options->Jop      = 0;
179:   options->m        = 1;
180:   options->eps      = 1.e-6;

182:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
183:   PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);
184:   PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);
185:   PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);
186:   PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);
187:   PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);
188:   PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);
189:   PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);
190:   PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);
191:   PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);
192:   PetscOptionsEnd();
193:   options->ke = PetscSqrtScalar(options->Jop);
194:   if (options->Jop==0.0) {
195:     options->Jo = 1.0/PetscPowScalar(options->a,2);
196:   } else {
197:     options->Jo = options->Jop*PetscCosReal(options->ke*options->a)/(1.0-PetscCosReal(options->ke*options->a));
198:   }
199:   options->ky = PETSC_PI*options->m/options->b;
200:   if (PetscPowReal(options->ky, 2) < options->Jop) {
201:     options->kx = PetscSqrtScalar(options->Jop-PetscPowScalar(options->ky,2));
202:     options->DeltaPrime = -2.0*options->kx*options->a*PetscCosReal(options->kx*options->a)/PetscSinReal(options->kx*options->a);
203:   } else if (PetscPowReal(options->ky, 2) > options->Jop) {
204:     options->kx = PetscSqrtScalar(PetscPowScalar(options->ky,2)-options->Jop);
205:     options->DeltaPrime = -2.0*options->kx*options->a*PetscCoshReal(options->kx*options->a)/PetscSinhReal(options->kx*options->a);
206:   } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/
207:     options->kx = 0;
208:     options->DeltaPrime = -2.0;
209:   }
210:   PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime);

212:   return 0;
213: }

215: static void f_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
216:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
217:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
218:                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
219: {
220:   const PetscScalar *pn = &u[uOff[DENSITY]];
221:   *f0 = *pn;
222: }

224: static PetscErrorCode PostStep(TS ts)
225: {
226:   DM        dm;
227:   AppCtx   *ctx;
228:   PetscInt  stepi,num;
229:   Vec       X;

231:   TSGetApplicationContext(ts, &ctx);
232:   if (ctx->debug<1) return 0;
233:   TSGetSolution(ts, &X);
234:   VecGetDM(X, &dm);
235:   TSGetStepNumber(ts, &stepi);
236:   DMGetOutputSequenceNumber(dm, &num, NULL);
237:   if (num < 0) DMSetOutputSequenceNumber(dm, 0, 0.0);
238:   PetscObjectSetName((PetscObject) X, "u");
239:   VecViewFromOptions(X, NULL, "-vec_view");
240:   /* print integrals */
241:   {
242:     PetscDS          prob;
243:     DM               plex;
244:     PetscScalar den, tt[5];
245:     DMConvert(dm, DMPLEX, &plex);
246:     DMGetDS(plex, &prob);
247:     PetscDSSetObjective(prob, 0, &f_n);
248:     DMPlexComputeIntegralFEM(plex,X,tt,ctx);
249:     den = tt[0];
250:     DMDestroy(&plex);
251:     PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));
252:   }
253:   return 0;
254: }

256: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
257: {
259:   DMCreate(comm, dm);
260:   DMSetType(*dm, DMPLEX);
261:   DMSetFromOptions(*dm);
262:   DMViewFromOptions(*dm, NULL, "-dm_view");

264:   DMGetBoundingBox(*dm, ctx->lower, ctx->upper);
265:   ctx->a = (ctx->upper[0] - ctx->lower[0])/2.0;
266:   ctx->b = (ctx->upper[1] - ctx->lower[1])/2.0;
267:   return 0;
268: }

270: static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
271: {
272:   AppCtx *lctx = (AppCtx*)ctx;
273:   u[0] = 2.*lctx->a + x[0];
274:   return 0;
275: }

277: static PetscErrorCode Omega_0(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 psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
284: {
285:   AppCtx *lctx = (AppCtx*)ctx;
286:   /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z.  The stability
287:      is analytically known and reported in ProcessOptions. */
288:   if (lctx->ke!=0.0) {
289:     u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
290:   } else {
291:     u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
292:   }
293:   return 0;
294: }

296: static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
297: {
298:   u[0] = 0.0;
299:   return 0;
300: }

302: static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
303: {
304:   u[0] = 0.0;
305:   return 0;
306: }

308: static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
309: {
310:   AppCtx *ctx = (AppCtx*)a_ctx;
311:   PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
312:   if (x[0] == ctx->lower[0] || x[0] == ctx->upper[0]) r = 0;
313:   u[0] = r;
314:   /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
315:   return 0;
316: }

318: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
319: {
320:   u[0] = 0.0;
321:   return 0;
322: }

324: static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
325: {
326:   u[0] = 0.0;
327:   return 0;
328: }

330: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
331: {
332:   PetscDS        ds;
333:   DMLabel        label;
334:   const PetscInt id = 1;

337:   DMGetLabel(dm, "marker", &label);
338:   DMGetDS(dm, &ds);
339:   PetscDSSetResidual(ds, 0, f0_n,     f1_n);
340:   PetscDSSetResidual(ds, 1, f0_Omega, f1_Omega);
341:   PetscDSSetResidual(ds, 2, f0_psi,   f1_psi);
342:   PetscDSSetResidual(ds, 3, f0_phi,   f1_phi);
343:   PetscDSSetResidual(ds, 4, f0_jz,    f1_jz);
344:   ctx->initialFuncs[0] = initialSolution_n;
345:   ctx->initialFuncs[1] = initialSolution_Omega;
346:   ctx->initialFuncs[2] = initialSolution_psi;
347:   ctx->initialFuncs[3] = initialSolution_phi;
348:   ctx->initialFuncs[4] = initialSolution_jz;
349:   for (PetscInt f = 0; f < 5; ++f) {
350:     PetscDSSetImplicit(ds, f, ctx->implicit);
351:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, ctx, NULL);
352:   }
353:   PetscDSSetContext(ds, 0, ctx);
354:   return 0;
355: }

357: static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
358: {
359:   PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
360:   Vec            eq;
361:   AppCtx *ctxarr[3];

363:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
364:   DMCreateLocalVector(dmAux, &eq);
365:   DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);
366:   DMSetAuxiliaryVec(dm, NULL, 0, 0, eq);
367:   if (ctx->plotRef) {  /* plot reference functions */
368:     PetscViewer       viewer = NULL;
369:     PetscBool         isHDF5,isVTK;
370:     char              buf[256];
371:     Vec               global;
372:     PetscInt          dim;

374:     DMGetDimension(dm, &dim);
375:     DMCreateGlobalVector(dmAux,&global);
376:     VecSet(global,.0); /* BCs! */
377:     DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);
378:     DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);
379:     PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);
380: #ifdef PETSC_HAVE_HDF5
381:     PetscViewerSetType(viewer,PETSCVIEWERHDF5);
382: #else
383:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
384: #endif
385:     PetscViewerSetFromOptions(viewer);
386:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
387:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
388:     if (isHDF5) {
389:       PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", dim);
390:     } else if (isVTK) {
391:       PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", dim);
392:       PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
393:     }
394:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
395:     PetscViewerFileSetName(viewer,buf);
396:     if (isHDF5) DMView(dmAux,viewer);
397:     /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
398:     PetscObjectSetName((PetscObject) global, "u0");
399:     VecView(global,viewer);
400:     PetscViewerDestroy(&viewer);
401:     VecDestroy(&global);
402:   }
403:   VecDestroy(&eq);
404:   return 0;
405: }

407: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
408: {
409:   DM             dmAux, coordDM;
410:   PetscInt       f;

412:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
413:   DMGetCoordinateDM(dm, &coordDM);
414:   if (!feAux) return 0;
415:   DMClone(dm, &dmAux);
416:   DMSetCoordinateDM(dmAux, coordDM);
417:   for (f = 0; f < NfAux; ++f) DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);
418:   DMCreateDS(dmAux);
419:   SetupEquilibriumFields(dm, dmAux, user);
420:   DMDestroy(&dmAux);
421:   return 0;
422: }

424: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
425: {
426:   DM              cdm = dm;
427:   PetscFE         fe[5], feAux[3];
428:   PetscInt        dim, Nf = 5, NfAux = 3, f;
429:   PetscBool       simplex;
430:   MPI_Comm        comm;

433:   /* Create finite element */
434:   PetscObjectGetComm((PetscObject) dm, &comm);
435:   DMGetDimension(dm, &dim);
436:   DMPlexIsSimplex(dm, &simplex);
437:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[0]);
438:   PetscObjectSetName((PetscObject) fe[0], "density");
439:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[1]);
440:   PetscObjectSetName((PetscObject) fe[1], "vorticity");
441:   PetscFECopyQuadrature(fe[0], fe[1]);
442:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[2]);
443:   PetscObjectSetName((PetscObject) fe[2], "flux");
444:   PetscFECopyQuadrature(fe[0], fe[2]);
445:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[3]);
446:   PetscObjectSetName((PetscObject) fe[3], "potential");
447:   PetscFECopyQuadrature(fe[0], fe[3]);
448:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[4]);
449:   PetscObjectSetName((PetscObject) fe[4], "current");
450:   PetscFECopyQuadrature(fe[0], fe[4]);

452:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[0]);
453:   PetscObjectSetName((PetscObject) feAux[0], "n_0");
454:   PetscFECopyQuadrature(fe[0], feAux[0]);
455:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[1]);
456:   PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");
457:   PetscFECopyQuadrature(fe[0], feAux[1]);
458:   PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[2]);
459:   PetscObjectSetName((PetscObject) feAux[2], "flux_0");
460:   PetscFECopyQuadrature(fe[0], feAux[2]);
461:   /* Set discretization and boundary conditions for each mesh */
462:   for (f = 0; f < Nf; ++f) DMSetField(dm, f, NULL, (PetscObject) fe[f]);
463:   DMCreateDS(dm);
464:   SetupProblem(dm, ctx);
465:   while (cdm) {
466:     SetupAuxDM(dm, NfAux, feAux, ctx);
467:     DMCopyDisc(dm, cdm);
468:     DMGetCoarseDM(cdm, &cdm);
469:   }
470:   for (f = 0; f < Nf; ++f) PetscFEDestroy(&fe[f]);
471:   for (f = 0; f < NfAux; ++f) PetscFEDestroy(&feAux[f]);
472:   return 0;
473: }

475: int main(int argc, char **argv)
476: {
477:   DM             dm;
478:   TS             ts;
479:   Vec            u, r;
480:   AppCtx         ctx;
481:   PetscReal      t       = 0.0;
482:   PetscReal      L2error = 0.0;
483:   AppCtx        *ctxarr[5];

485:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
486:   s_ctx = &ctx;
487:   PetscInitialize(&argc, &argv, NULL,help);
488:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
489:   /* create mesh and problem */
490:   CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
491:   DMSetApplicationContext(dm, &ctx);
492:   PetscMalloc1(5, &ctx.initialFuncs);
493:   SetupDiscretization(dm, &ctx);
494:   DMCreateGlobalVector(dm, &u);
495:   PetscObjectSetName((PetscObject) u, "u");
496:   VecDuplicate(u, &r);
497:   PetscObjectSetName((PetscObject) r, "r");
498:   /* create TS */
499:   TSCreate(PETSC_COMM_WORLD, &ts);
500:   TSSetDM(ts, dm);
501:   TSSetApplicationContext(ts, &ctx);
502:   DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
503:   if (ctx.implicit) {
504:     DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
505:     DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
506:   } else {
507:     DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);
508:   }
509:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
510:   TSSetFromOptions(ts);
511:   TSSetPostStep(ts, PostStep);
512:   /* make solution & solve */
513:   DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);
514:   TSSetSolution(ts,u);
515:   DMViewFromOptions(dm, NULL, "-dm_view");
516:   PostStep(ts); /* print the initial state */
517:   TSSolve(ts, u);
518:   TSGetTime(ts, &t);
519:   DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);
520:   if (L2error < 1.0e-11) PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");
521:   else                   PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);
522:   VecDestroy(&u);
523:   VecDestroy(&r);
524:   TSDestroy(&ts);
525:   DMDestroy(&dm);
526:   PetscFree(ctx.initialFuncs);
527:   PetscFinalize();
528:   return 0;
529: }

531: /*TEST

533:   test:
534:     suffix: 0
535:     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 \
536:           -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
537:   test:
538:     # Remapping with periodicity is broken
539:     suffix: 1
540:     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 \
541:            -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0

543: TEST*/