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: PetscErrorCode ierr;
227: DM dm;
228: AppCtx *ctx;
229: PetscInt stepi,num;
230: Vec X;
232: TSGetApplicationContext(ts, &ctx);
233: if (ctx->debug<1) return(0);
234: TSGetSolution(ts, &X);
235: VecGetDM(X, &dm);
236: TSGetStepNumber(ts, &stepi);
237: DMGetOutputSequenceNumber(dm, &num, NULL);
238: if (num < 0) {DMSetOutputSequenceNumber(dm, 0, 0.0);}
239: PetscObjectSetName((PetscObject) X, "u");
240: VecViewFromOptions(X, NULL, "-vec_view");
241: /* print integrals */
242: {
243: PetscDS prob;
244: DM plex;
245: PetscScalar den, tt[5];
246: DMConvert(dm, DMPLEX, &plex);
247: DMGetDS(plex, &prob);
248: PetscDSSetObjective(prob, 0, &f_n);
249: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
250: den = tt[0];
251: DMDestroy(&plex);
252: PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));
253: }
254: return(0);
255: }
257: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
258: {
262: DMCreate(comm, dm);
263: DMSetType(*dm, DMPLEX);
264: DMSetFromOptions(*dm);
265: DMViewFromOptions(*dm, NULL, "-dm_view");
267: DMGetBoundingBox(*dm, ctx->lower, ctx->upper);
268: ctx->a = (ctx->upper[0] - ctx->lower[0])/2.0;
269: ctx->b = (ctx->upper[1] - ctx->lower[1])/2.0;
270: return(0);
271: }
273: static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
274: {
275: AppCtx *lctx = (AppCtx*)ctx;
276: u[0] = 2.*lctx->a + x[0];
277: return 0;
278: }
280: static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
281: {
282: u[0] = 0.0;
283: return 0;
284: }
286: static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
287: {
288: AppCtx *lctx = (AppCtx*)ctx;
289: /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z. The stability
290: is analytically known and reported in ProcessOptions. */
291: if (lctx->ke!=0.0) {
292: u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
293: } else {
294: u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
295: }
296: return 0;
297: }
299: static PetscErrorCode initialSolution_n(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 initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
306: {
307: u[0] = 0.0;
308: return 0;
309: }
311: static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
312: {
313: AppCtx *ctx = (AppCtx*)a_ctx;
314: PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
315: if (x[0] == ctx->lower[0] || x[0] == ctx->upper[0]) r = 0;
316: u[0] = r;
317: /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
318: return 0;
319: }
321: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
322: {
323: u[0] = 0.0;
324: return 0;
325: }
327: static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
328: {
329: u[0] = 0.0;
330: return 0;
331: }
333: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
334: {
335: PetscDS ds;
336: DMLabel label;
337: const PetscInt id = 1;
338: PetscErrorCode ierr, f;
341: DMGetLabel(dm, "marker", &label);
342: DMGetDS(dm, &ds);
343: PetscDSSetResidual(ds, 0, f0_n, f1_n);
344: PetscDSSetResidual(ds, 1, f0_Omega, f1_Omega);
345: PetscDSSetResidual(ds, 2, f0_psi, f1_psi);
346: PetscDSSetResidual(ds, 3, f0_phi, f1_phi);
347: PetscDSSetResidual(ds, 4, f0_jz, f1_jz);
348: ctx->initialFuncs[0] = initialSolution_n;
349: ctx->initialFuncs[1] = initialSolution_Omega;
350: ctx->initialFuncs[2] = initialSolution_psi;
351: ctx->initialFuncs[3] = initialSolution_phi;
352: ctx->initialFuncs[4] = initialSolution_jz;
353: for (f = 0; f < 5; ++f) {
354: PetscDSSetImplicit(ds, f, ctx->implicit);
355: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, ctx, NULL);
356: }
357: PetscDSSetContext(ds, 0, ctx);
358: return(0);
359: }
361: static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
362: {
363: PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
364: Vec eq;
366: AppCtx *ctxarr[3];
368: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
370: DMCreateLocalVector(dmAux, &eq);
371: DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);
372: DMSetAuxiliaryVec(dm, NULL, 0, eq);
373: if (ctx->plotRef) { /* plot reference functions */
374: PetscViewer viewer = NULL;
375: PetscBool isHDF5,isVTK;
376: char buf[256];
377: Vec global;
378: PetscInt dim;
380: DMGetDimension(dm, &dim);
381: DMCreateGlobalVector(dmAux,&global);
382: VecSet(global,.0); /* BCs! */
383: DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);
384: DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);
385: PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);
386: #ifdef PETSC_HAVE_HDF5
387: PetscViewerSetType(viewer,PETSCVIEWERHDF5);
388: #else
389: PetscViewerSetType(viewer,PETSCVIEWERVTK);
390: #endif
391: PetscViewerSetFromOptions(viewer);
392: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
393: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
394: if (isHDF5) {
395: PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", dim);
396: } else if (isVTK) {
397: PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", dim);
398: PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
399: }
400: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
401: PetscViewerFileSetName(viewer,buf);
402: if (isHDF5) {DMView(dmAux,viewer);}
403: /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
404: PetscObjectSetName((PetscObject) global, "u0");
405: VecView(global,viewer);
406: PetscViewerDestroy(&viewer);
407: VecDestroy(&global);
408: }
409: VecDestroy(&eq);
410: return(0);
411: }
413: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
414: {
415: DM dmAux, coordDM;
416: PetscInt f;
420: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
421: DMGetCoordinateDM(dm, &coordDM);
422: if (!feAux) return(0);
423: DMClone(dm, &dmAux);
424: DMSetCoordinateDM(dmAux, coordDM);
425: for (f = 0; f < NfAux; ++f) {DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);}
426: DMCreateDS(dmAux);
427: SetupEquilibriumFields(dm, dmAux, user);
428: DMDestroy(&dmAux);
429: return(0);
430: }
432: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
433: {
434: DM cdm = dm;
435: PetscFE fe[5], feAux[3];
436: PetscInt dim, Nf = 5, NfAux = 3, f;
437: PetscBool simplex;
438: MPI_Comm comm;
439: PetscErrorCode ierr;
442: /* Create finite element */
443: PetscObjectGetComm((PetscObject) dm, &comm);
444: DMGetDimension(dm, &dim);
445: DMPlexIsSimplex(dm, &simplex);
446: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[0]);
447: PetscObjectSetName((PetscObject) fe[0], "density");
448: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[1]);
449: PetscObjectSetName((PetscObject) fe[1], "vorticity");
450: PetscFECopyQuadrature(fe[0], fe[1]);
451: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[2]);
452: PetscObjectSetName((PetscObject) fe[2], "flux");
453: PetscFECopyQuadrature(fe[0], fe[2]);
454: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[3]);
455: PetscObjectSetName((PetscObject) fe[3], "potential");
456: PetscFECopyQuadrature(fe[0], fe[3]);
457: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[4]);
458: PetscObjectSetName((PetscObject) fe[4], "current");
459: PetscFECopyQuadrature(fe[0], fe[4]);
461: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[0]);
462: PetscObjectSetName((PetscObject) feAux[0], "n_0");
463: PetscFECopyQuadrature(fe[0], feAux[0]);
464: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[1]);
465: PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");
466: PetscFECopyQuadrature(fe[0], feAux[1]);
467: PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[2]);
468: PetscObjectSetName((PetscObject) feAux[2], "flux_0");
469: PetscFECopyQuadrature(fe[0], feAux[2]);
470: /* Set discretization and boundary conditions for each mesh */
471: for (f = 0; f < Nf; ++f) {DMSetField(dm, f, NULL, (PetscObject) fe[f]);}
472: DMCreateDS(dm);
473: SetupProblem(dm, ctx);
474: while (cdm) {
475: SetupAuxDM(dm, NfAux, feAux, ctx);
476: DMCopyDisc(dm, cdm);
477: DMGetCoarseDM(cdm, &cdm);
478: }
479: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&fe[f]);}
480: for (f = 0; f < NfAux; ++f) {PetscFEDestroy(&feAux[f]);}
481: return(0);
482: }
484: int main(int argc, char **argv)
485: {
486: DM dm;
487: TS ts;
488: Vec u, r;
489: AppCtx ctx;
490: PetscReal t = 0.0;
491: PetscReal L2error = 0.0;
493: AppCtx *ctxarr[5];
495: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
496: s_ctx = &ctx;
497: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
498: ProcessOptions(PETSC_COMM_WORLD, &ctx);
499: /* create mesh and problem */
500: CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
501: DMSetApplicationContext(dm, &ctx);
502: PetscMalloc1(5, &ctx.initialFuncs);
503: SetupDiscretization(dm, &ctx);
504: DMCreateGlobalVector(dm, &u);
505: PetscObjectSetName((PetscObject) u, "u");
506: VecDuplicate(u, &r);
507: PetscObjectSetName((PetscObject) r, "r");
508: /* create TS */
509: TSCreate(PETSC_COMM_WORLD, &ts);
510: TSSetDM(ts, dm);
511: TSSetApplicationContext(ts, &ctx);
512: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
513: if (ctx.implicit) {
514: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
515: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
516: } else {
517: DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);
518: }
519: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
520: TSSetFromOptions(ts);
521: TSSetPostStep(ts, PostStep);
522: /* make solution & solve */
523: DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);
524: TSSetSolution(ts,u);
525: DMViewFromOptions(dm, NULL, "-dm_view");
526: PostStep(ts); /* print the initial state */
527: TSSolve(ts, u);
528: TSGetTime(ts, &t);
529: DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);
530: if (L2error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
531: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);}
532: VecDestroy(&u);
533: VecDestroy(&r);
534: TSDestroy(&ts);
535: DMDestroy(&dm);
536: PetscFree(ctx.initialFuncs);
537: PetscFinalize();
538: return ierr;
539: }
541: /*TEST
543: test:
544: suffix: 0
545: 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 \
546: -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
547: test:
548: # Remapping with periodicity is broken
549: suffix: 1
550: 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 \
551: -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
553: TEST*/