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*/