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>
20: #include <assert.h>
22: typedef struct {
23: PetscInt debug; /* The debugging level */
24: PetscBool plotRef; /* Plot the reference fields */
25: /* Domain and mesh definition */
26: PetscInt dim; /* The topological mesh dimension */
27: char filename[2048]; /* The optional ExodusII file */
28: PetscBool cell_simplex; /* Simplicial mesh */
29: DMBoundaryType boundary_types[3];
30: PetscInt cells[3];
31: PetscInt refine;
32: /* geometry */
33: PetscReal domain_lo[3], domain_hi[3];
34: DMBoundaryType periodicity[3]; /* The domain periodicity */
35: PetscReal b0[3]; /* not used */
36: /* Problem definition */
37: PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
38: PetscReal mu, eta, beta;
39: PetscReal a,b,Jo,Jop,m,ke,kx,ky,DeltaPrime,eps;
40: /* solver */
41: PetscBool implicit;
42: } AppCtx;
44: static AppCtx *s_ctx;
46: static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
47: {
48: PetscScalar ret = df[0]*dg[1] - df[1]*dg[0];
49: return ret;
50: }
52: enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};
54: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
55: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
56: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
57: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
58: {
59: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
60: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
61: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
62: const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
63: const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
64: f0[0] += - poissonBracket(dim,pnDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer) - poissonBracket(dim,logRefDenDer, pphiDer);
65: if (u_t) f0[0] += u_t[DENSITY];
66: }
68: static void f1_n(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 f1[])
72: {
73: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
74: PetscInt d;
76: for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
77: }
79: static void f0_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
80: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
81: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
82: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
83: {
84: const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
85: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
86: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
87: const PetscScalar *jzDer = &u_x[uOff_x[JZ]];
89: f0[0] += - poissonBracket(dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer);
90: if (u_t) f0[0] += u_t[OMEGA];
91: }
93: static void f1_Omega(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 f1[])
97: {
98: const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
99: PetscInt d;
101: for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
102: }
104: static void f0_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
105: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
106: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
107: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
108: {
109: const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
110: const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
111: const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
112: const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]];
113: const PetscScalar *logRefDenDer= &a_x[aOff_x[DENSITY]];
114: PetscScalar psiDer[3];
115: PetscScalar phi_n_Der[3];
116: PetscInt d;
117: 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 */
118: for (d = 0; d < dim; ++d) {
119: psiDer[d] = refPsiDer[d] + ppsiDer[d];
120: phi_n_Der[d] = pphiDer[d] - pnDer[d];
121: }
122: f0[0] = - poissonBracket(dim,psiDer, phi_n_Der) + poissonBracket(dim,logRefDenDer, ppsiDer);
123: if (u_t) f0[0] += u_t[PSI];
124: }
126: static void f1_psi(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 f1[])
130: {
131: const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
132: PetscInt d;
134: for (d = 0; d < dim-1; ++d) f1[d] = -(s_ctx->eta/s_ctx->beta)*ppsi[d];
135: }
137: static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
138: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
139: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
140: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
141: {
142: f0[0] = -u[uOff[OMEGA]];
143: }
145: static void f1_phi(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 f1[])
149: {
150: const PetscScalar *pphi = &u_x[uOff_x[PHI]];
151: PetscInt d;
153: for (d = 0; d < dim-1; ++d) f1[d] = pphi[d];
154: }
156: static void f0_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
157: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
158: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
159: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
160: {
161: f0[0] = u[uOff[JZ]];
162: }
164: static void f1_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
168: {
169: const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
170: const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */
171: PetscInt d;
173: for (d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
174: }
176: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
177: {
178: PetscBool flg;
180: PetscInt ii, bd;
182: options->debug = 1;
183: options->plotRef = PETSC_FALSE;
184: options->dim = 2;
185: options->filename[0] = '\0';
186: options->cell_simplex = PETSC_FALSE;
187: options->implicit = PETSC_FALSE;
188: options->refine = 2;
189: options->domain_lo[0] = 0.0;
190: options->domain_lo[1] = 0.0;
191: options->domain_lo[2] = 0.0;
192: options->domain_hi[0] = 2.0;
193: options->domain_hi[1] = 2.0*PETSC_PI;
194: options->domain_hi[2] = 2.0;
195: options->periodicity[0] = DM_BOUNDARY_NONE;
196: options->periodicity[1] = DM_BOUNDARY_NONE;
197: options->periodicity[2] = DM_BOUNDARY_NONE;
198: options->mu = 0;
199: options->eta = 0;
200: options->beta = 1;
201: options->a = 1;
202: options->b = PETSC_PI;
203: options->Jop = 0;
204: options->m = 1;
205: options->eps = 1.e-6;
207: for (ii = 0; ii < options->dim; ++ii) options->cells[ii] = 4;
208: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
209: PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);
210: PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);
211: PetscOptionsInt("-dim", "The topological mesh dimension", "ex48.c", options->dim, &options->dim, NULL);
212: if (options->dim < 2 || options->dim > 3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Dim %D must be 2 or 3",options->dim);
213: PetscOptionsInt("-dm_refine", "Hack to get refinement level for cylinder", "ex48.c", options->refine, &options->refine, NULL);
214: PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);
215: PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);
216: PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);
217: PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);
218: PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);
219: PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);
220: PetscOptionsString("-f", "Exodus.II filename to read", "ex48.c", options->filename, options->filename, sizeof(options->filename), &flg);
221: PetscOptionsBool("-cell_simplex", "Simplicial (true) or tensor (false) mesh", "ex48.c", options->cell_simplex, &options->cell_simplex, NULL);
222: PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);
223: ii = options->dim;
224: PetscOptionsRealArray("-domain_hi", "Domain size", "ex48.c", options->domain_hi, &ii, NULL);
225: ii = options->dim;
226: PetscOptionsRealArray("-domain_lo", "Domain size", "ex48.c", options->domain_lo, &ii, NULL);
227: ii = options->dim;
228: bd = options->periodicity[0];
229: PetscOptionsEList("-x_periodicity", "The x-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);
230: options->periodicity[0] = (DMBoundaryType) bd;
231: bd = options->periodicity[1];
232: PetscOptionsEList("-y_periodicity", "The y-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);
233: options->periodicity[1] = (DMBoundaryType) bd;
234: bd = options->periodicity[2];
235: PetscOptionsEList("-z_periodicity", "The z-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);
236: options->periodicity[2] = (DMBoundaryType) bd;
237: ii = options->dim;
238: PetscOptionsIntArray("-cells", "Number of cells in each dimension", "ex48.c", options->cells, &ii, NULL);
239: PetscOptionsEnd();
240: options->a = (options->domain_hi[0]-options->domain_lo[0])/2.0;
241: options->b = (options->domain_hi[1]-options->domain_lo[1])/2.0;
242: for (ii = 0; ii < options->dim; ++ii) {
243: if (options->domain_hi[ii] <= options->domain_lo[ii]) SETERRQ3(comm,PETSC_ERR_ARG_WRONG,"Domain %D lo=%g hi=%g",ii,options->domain_lo[ii],options->domain_hi[ii]);
244: }
245: options->ke = PetscSqrtScalar(options->Jop);
246: if (options->Jop==0.0) {
247: options->Jo = 1.0/PetscPowScalar(options->a,2);
248: } else {
249: options->Jo = options->Jop*PetscCosReal(options->ke*options->a)/(1.0-PetscCosReal(options->ke*options->a));
250: }
251: options->ky = PETSC_PI*options->m/options->b;
252: if (PetscPowReal(options->ky, 2) < options->Jop) {
253: options->kx = PetscSqrtScalar(options->Jop-PetscPowScalar(options->ky,2));
254: options->DeltaPrime = -2.0*options->kx*options->a*PetscCosReal(options->kx*options->a)/PetscSinReal(options->kx*options->a);
255: } else if (PetscPowReal(options->ky, 2) > options->Jop) {
256: options->kx = PetscSqrtScalar(PetscPowScalar(options->ky,2)-options->Jop);
257: options->DeltaPrime = -2.0*options->kx*options->a*PetscCoshReal(options->kx*options->a)/PetscSinhReal(options->kx*options->a);
258: } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/
259: options->kx = 0;
260: options->DeltaPrime = -2.0;
261: }
262: PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime);
264: return(0);
265: }
267: static void f_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
268: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
269: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
270: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
271: {
272: const PetscScalar *pn = &u[uOff[DENSITY]];
273: *f0 = *pn;
274: }
276: static PetscErrorCode PostStep(TS ts)
277: {
278: PetscErrorCode ierr;
279: DM dm;
280: AppCtx *ctx;
281: PetscInt stepi,num;
282: Vec X;
284: TSGetApplicationContext(ts, &ctx); assert(ctx);
285: if (ctx->debug<1) return(0);
286: TSGetSolution(ts, &X);
287: VecGetDM(X, &dm);
288: TSGetStepNumber(ts, &stepi);
289: DMGetOutputSequenceNumber(dm, &num, NULL);
290: if (num < 0) {DMSetOutputSequenceNumber(dm, 0, 0.0);}
291: PetscObjectSetName((PetscObject) X, "u");
292: VecViewFromOptions(X, NULL, "-vec_view");
293: /* print integrals */
294: {
295: PetscDS prob;
296: DM plex;
297: PetscScalar den, tt[5];
298: DMConvert(dm, DMPLEX, &plex);
299: DMGetDS(plex, &prob);
300: PetscDSSetObjective(prob, 0, &f_n);
301: DMPlexComputeIntegralFEM(plex,X,tt,ctx);
302: den = tt[0];
303: DMDestroy(&plex);
304: PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));
305: }
306: return(0);
307: }
309: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
310: {
311: DM plex;
312: DMLabel label;
316: DMCreateLabel(dm, name);
317: DMGetLabel(dm, name, &label);
318: DMConvert(dm, DMPLEX, &plex);
319: DMPlexMarkBoundaryFaces(dm, 1, label);
320: DMDestroy(&plex);
321: return(0);
322: }
324: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
325: {
326: PetscInt dim = ctx->dim;
327: const char *filename = ctx->filename;
328: size_t len;
329: PetscMPIInt numProcs;
333: MPI_Comm_size(comm, &numProcs);
334: PetscStrlen(filename, &len);
335: if (len) {
336: DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);
337: } else {
338: PetscInt d;
340: /* create DM */
341: if (ctx->cell_simplex && dim == 3) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh a cylinder with simplices");
342: if (dim==2) {
343: PetscInt refineRatio, totCells = 1;
344: if (ctx->cell_simplex) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh 2D with simplices");
345: refineRatio = PetscMax((PetscInt) (PetscPowReal(numProcs, 1.0/dim) + 0.1) - 1, 1);
346: for (d = 0; d < dim; ++d) {
347: if (ctx->cells[d] < refineRatio) ctx->cells[d] = refineRatio;
348: if (ctx->periodicity[d]==DM_BOUNDARY_PERIODIC && ctx->cells[d]*refineRatio <= 2) refineRatio = 2;
349: }
350: for (d = 0; d < dim; ++d) {
351: ctx->cells[d] *= refineRatio;
352: totCells *= ctx->cells[d];
353: }
354: if (totCells % numProcs) SETERRQ2(comm,PETSC_ERR_ARG_WRONG,"Total cells %D not divisible by processes %D", totCells, numProcs);
355: DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, ctx->cells, ctx->domain_lo, ctx->domain_hi, ctx->periodicity, PETSC_TRUE, dm);
356: } else {
357: if (ctx->periodicity[0]==DM_BOUNDARY_PERIODIC || ctx->periodicity[1]==DM_BOUNDARY_PERIODIC) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot do periodic in x or y in a cylinder");
358: /* we stole dm_refine so clear it */
359: PetscOptionsClearValue(NULL,"-dm_refine");
360: DMPlexCreateHexCylinderMesh(comm, ctx->refine, ctx->periodicity[2], dm);
361: }
362: }
363: {
364: DM distributedMesh = NULL;
365: /* Distribute mesh over processes */
366: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
367: if (distributedMesh) {
368: DMDestroy(dm);
369: *dm = distributedMesh;
370: }
371: }
372: {
373: PetscBool hasLabel;
374: DMHasLabel(*dm, "marker", &hasLabel);
375: if (!hasLabel) {CreateBCLabel(*dm, "marker");}
376: }
377: {
378: char convType[256];
379: PetscBool flg;
380: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
381: PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex48",DMList,DMPLEX,convType,256,&flg);
382: PetscOptionsEnd();
383: if (flg) {
384: DM dmConv;
385: DMConvert(*dm,convType,&dmConv);
386: if (dmConv) {
387: DMDestroy(dm);
388: *dm = dmConv;
389: }
390: }
391: }
392: PetscObjectSetName((PetscObject) *dm, "Mesh");
393: DMSetFromOptions(*dm);
394: DMLocalizeCoordinates(*dm); /* needed for periodic */
395: return(0);
396: }
398: static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
399: {
400: AppCtx *lctx = (AppCtx*)ctx;
401: assert(ctx);
402: u[0] = (lctx->domain_hi-lctx->domain_lo)+x[0];
403: return 0;
404: }
406: static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
407: {
408: u[0] = 0.0;
409: return 0;
410: }
412: static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
413: {
414: AppCtx *lctx = (AppCtx*)ctx;
415: assert(ctx);
416: /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z. The stability
417: is analytically known and reported in ProcessOptions. */
418: if (lctx->ke!=0.0) {
419: u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
420: } else {
421: u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
422: }
423: return 0;
424: }
426: static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
427: {
428: u[0] = 0.0;
429: return 0;
430: }
432: static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
433: {
434: u[0] = 0.0;
435: return 0;
436: }
438: static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
439: {
440: AppCtx *ctx = (AppCtx*)a_ctx;
441: PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
442: assert(ctx);
443: if (x[0] == ctx->domain_lo[0] || x[0] == ctx->domain_hi[0]) r = 0;
444: u[0] = r;
445: /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
446: return 0;
447: }
449: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
450: {
451: u[0] = 0.0;
452: return 0;
453: }
455: static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
456: {
457: u[0] = 0.0;
458: return 0;
459: }
461: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
462: {
463: PetscDS prob;
464: const PetscInt id = 1;
465: PetscErrorCode ierr, f;
468: DMGetDS(dm, &prob);
469: PetscDSSetResidual(prob, 0, f0_n, f1_n);
470: PetscDSSetResidual(prob, 1, f0_Omega, f1_Omega);
471: PetscDSSetResidual(prob, 2, f0_psi, f1_psi);
472: PetscDSSetResidual(prob, 3, f0_phi, f1_phi);
473: PetscDSSetResidual(prob, 4, f0_jz, f1_jz);
474: ctx->initialFuncs[0] = initialSolution_n;
475: ctx->initialFuncs[1] = initialSolution_Omega;
476: ctx->initialFuncs[2] = initialSolution_psi;
477: ctx->initialFuncs[3] = initialSolution_phi;
478: ctx->initialFuncs[4] = initialSolution_jz;
479: for (f = 0; f < 5; ++f) {
480: PetscDSSetImplicit(prob, f, ctx->implicit);
481: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, 1, &id, ctx);
482: }
483: PetscDSSetContext(prob, 0, ctx);
484: return(0);
485: }
487: static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
488: {
489: PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
490: Vec eq;
492: AppCtx *ctxarr[3];
494: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
496: DMCreateLocalVector(dmAux, &eq);
497: DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);
498: PetscObjectCompose((PetscObject) dm, "A", (PetscObject) eq);
499: if (ctx->plotRef) { /* plot reference functions */
500: PetscViewer viewer = NULL;
501: PetscBool isHDF5,isVTK;
502: char buf[256];
503: Vec global;
504: DMCreateGlobalVector(dmAux,&global);
505: VecSet(global,.0); /* BCs! */
506: DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);
507: DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);
508: PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);
509: #ifdef PETSC_HAVE_HDF5
510: PetscViewerSetType(viewer,PETSCVIEWERHDF5);
511: #else
512: PetscViewerSetType(viewer,PETSCVIEWERVTK);
513: #endif
514: PetscViewerSetFromOptions(viewer);
515: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
516: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
517: if (isHDF5) {
518: PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", ctx->dim);
519: } else if (isVTK) {
520: PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", ctx->dim);
521: PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
522: }
523: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
524: PetscViewerFileSetName(viewer,buf);
525: if (isHDF5) {DMView(dmAux,viewer);}
526: /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
527: PetscObjectSetName((PetscObject) global, "u0");
528: VecView(global,viewer);
529: PetscViewerDestroy(&viewer);
530: VecDestroy(&global);
531: }
532: VecDestroy(&eq);
533: return(0);
534: }
536: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
537: {
538: DM dmAux, coordDM;
539: PetscInt f;
543: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
544: DMGetCoordinateDM(dm, &coordDM);
545: if (!feAux) return(0);
546: DMClone(dm, &dmAux);
547: PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
548: DMSetCoordinateDM(dmAux, coordDM);
549: for (f = 0; f < NfAux; ++f) {DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);}
550: DMCreateDS(dmAux);
551: SetupEquilibriumFields(dm, dmAux, user);
552: DMDestroy(&dmAux);
553: return(0);
554: }
556: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
557: {
558: DM cdm = dm;
559: const PetscInt dim = ctx->dim;
560: PetscFE fe[5], feAux[3];
561: PetscInt Nf = 5, NfAux = 3, f;
562: PetscBool cell_simplex = ctx->cell_simplex;
563: MPI_Comm comm;
564: PetscErrorCode ierr;
567: /* Create finite element */
568: PetscObjectGetComm((PetscObject) dm, &comm);
569: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);
570: PetscObjectSetName((PetscObject) fe[0], "density");
571: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);
572: PetscObjectSetName((PetscObject) fe[1], "vorticity");
573: PetscFECopyQuadrature(fe[0], fe[1]);
574: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);
575: PetscObjectSetName((PetscObject) fe[2], "flux");
576: PetscFECopyQuadrature(fe[0], fe[2]);
577: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);
578: PetscObjectSetName((PetscObject) fe[3], "potential");
579: PetscFECopyQuadrature(fe[0], fe[3]);
580: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);
581: PetscObjectSetName((PetscObject) fe[4], "current");
582: PetscFECopyQuadrature(fe[0], fe[4]);
584: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);
585: PetscObjectSetName((PetscObject) feAux[0], "n_0");
586: PetscFECopyQuadrature(fe[0], feAux[0]);
587: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);
588: PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");
589: PetscFECopyQuadrature(fe[0], feAux[1]);
590: PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);
591: PetscObjectSetName((PetscObject) feAux[2], "flux_0");
592: PetscFECopyQuadrature(fe[0], feAux[2]);
593: /* Set discretization and boundary conditions for each mesh */
594: for (f = 0; f < Nf; ++f) {DMSetField(dm, f, NULL, (PetscObject) fe[f]);}
595: DMCreateDS(dm);
596: SetupProblem(dm, ctx);
597: while (cdm) {
598: SetupAuxDM(dm, NfAux, feAux, ctx);
599: {
600: PetscBool hasLabel;
602: DMHasLabel(cdm, "marker", &hasLabel);
603: if (!hasLabel) {CreateBCLabel(cdm, "marker");}
604: }
605: DMCopyDisc(dm, cdm);
606: DMGetCoarseDM(cdm, &cdm);
607: }
608: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&fe[f]);}
609: for (f = 0; f < NfAux; ++f) {PetscFEDestroy(&feAux[f]);}
610: return(0);
611: }
613: int main(int argc, char **argv)
614: {
615: DM dm;
616: TS ts;
617: Vec u, r;
618: AppCtx ctx;
619: PetscReal t = 0.0;
620: PetscReal L2error = 0.0;
622: AppCtx *ctxarr[5];
624: ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
625: s_ctx = &ctx;
626: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
627: ProcessOptions(PETSC_COMM_WORLD, &ctx);
628: /* create mesh and problem */
629: CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
630: DMSetApplicationContext(dm, &ctx);
631: PetscMalloc1(5, &ctx.initialFuncs);
632: SetupDiscretization(dm, &ctx);
633: DMCreateGlobalVector(dm, &u);
634: PetscObjectSetName((PetscObject) u, "u");
635: VecDuplicate(u, &r);
636: PetscObjectSetName((PetscObject) r, "r");
637: /* create TS */
638: TSCreate(PETSC_COMM_WORLD, &ts);
639: TSSetDM(ts, dm);
640: TSSetApplicationContext(ts, &ctx);
641: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
642: if (ctx.implicit) {
643: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
644: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
645: } else {
646: DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);
647: }
648: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
649: TSSetFromOptions(ts);
650: TSSetPostStep(ts, PostStep);
651: /* make solution & solve */
652: DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);
653: TSSetSolution(ts,u);
654: DMViewFromOptions(dm, NULL, "-dm_view");
655: PostStep(ts); /* print the initial state */
656: TSSolve(ts, u);
657: TSGetTime(ts, &t);
658: DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);
659: if (L2error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
660: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);}
661: #if 0
662: {
663: PetscReal res = 0.0;
664: /* Check discretization error */
665: VecViewFromOptions(u, NULL, "-initial_guess_view");
666: DMComputeL2Diff(dm, 0.0, ctx.exactFuncs, NULL, u, &error);
667: if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
668: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
669: /* Check residual */
670: SNESComputeFunction(snes, u, r);
671: VecChop(r, 1.0e-10);
672: VecViewFromOptions(r, NULL, "-initial_residual_view");
673: VecNorm(r, NORM_2, &res);
674: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
675: /* Check Jacobian */
676: {
677: Mat A;
678: Vec b;
680: SNESGetJacobian(snes, &A, NULL, NULL, NULL);
681: SNESComputeJacobian(snes, u, A, A);
682: VecDuplicate(u, &b);
683: VecSet(r, 0.0);
684: SNESComputeFunction(snes, r, b);
685: MatMult(A, u, r);
686: VecAXPY(r, 1.0, b);
687: VecDestroy(&b);
688: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
689: VecChop(r, 1.0e-10);
690: VecViewFromOptions(r, NULL, "-linear_residual_view");
691: VecNorm(r, NORM_2, &res);
692: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
693: }
694: }
695: #endif
696: VecDestroy(&u);
697: VecDestroy(&r);
698: TSDestroy(&ts);
699: DMDestroy(&dm);
700: PetscFree(ctx.initialFuncs);
701: PetscFinalize();
702: return ierr;
703: }
705: /*TEST
707: test:
708: suffix: 0
709: args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
710: test:
711: suffix: 1
712: args: -debug 1 -dim 3 -dm_refine 1 -z_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 -domain_lo -2,-1,-1 -domain_hi 2,1,1
714: TEST*/