xref: /petsc/src/ts/tutorials/ex48.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254) 
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";
3 
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*/
16 
17 #include <petscdmplex.h>
18 #include <petscts.h>
19 #include <petscds.h>
20 
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;
32 
33 static AppCtx *s_ctx;
34 
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 }
40 
41 enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};
42 
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 }
56 
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;
64 
65   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
66 }
67 
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]];
77 
78   f0[0] += - poissonBracket(dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer);
79   if (u_t) f0[0] += u_t[OMEGA];
80 }
81 
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;
89 
90   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
91 }
92 
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 }
114 
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;
122 
123   for (d = 0; d < dim-1; ++d) f1[d] = -(s_ctx->eta/s_ctx->beta)*ppsi[d];
124 }
125 
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 }
133 
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;
141 
142   for (d = 0; d < dim-1; ++d) f1[d] = pphi[d];
143 }
144 
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 }
152 
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;
161 
162   for (d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
163 }
164 
165 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
166 {
167   PetscErrorCode ierr;
168 
169   PetscFunctionBeginUser;
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;
181 
182   ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr);
183   CHKERRQ(PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL));
184   CHKERRQ(PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL));
185   CHKERRQ(PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL));
186   CHKERRQ(PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL));
187   CHKERRQ(PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL));
188   CHKERRQ(PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL));
189   CHKERRQ(PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL));
190   CHKERRQ(PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL));
191   CHKERRQ(PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL));
192   ierr = PetscOptionsEnd();CHKERRQ(ierr);
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   CHKERRQ(PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime));
211 
212   PetscFunctionReturn(0);
213 }
214 
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 }
223 
224 static PetscErrorCode PostStep(TS ts)
225 {
226   DM        dm;
227   AppCtx   *ctx;
228   PetscInt  stepi,num;
229   Vec       X;
230 
231   PetscFunctionBegin;
232   CHKERRQ(TSGetApplicationContext(ts, &ctx));
233   if (ctx->debug<1) PetscFunctionReturn(0);
234   CHKERRQ(TSGetSolution(ts, &X));
235   CHKERRQ(VecGetDM(X, &dm));
236   CHKERRQ(TSGetStepNumber(ts, &stepi));
237   CHKERRQ(DMGetOutputSequenceNumber(dm, &num, NULL));
238   if (num < 0) CHKERRQ(DMSetOutputSequenceNumber(dm, 0, 0.0));
239   CHKERRQ(PetscObjectSetName((PetscObject) X, "u"));
240   CHKERRQ(VecViewFromOptions(X, NULL, "-vec_view"));
241   /* print integrals */
242   {
243     PetscDS          prob;
244     DM               plex;
245     PetscScalar den, tt[5];
246     CHKERRQ(DMConvert(dm, DMPLEX, &plex));
247     CHKERRQ(DMGetDS(plex, &prob));
248     CHKERRQ(PetscDSSetObjective(prob, 0, &f_n));
249     CHKERRQ(DMPlexComputeIntegralFEM(plex,X,tt,ctx));
250     den = tt[0];
251     CHKERRQ(DMDestroy(&plex));
252     CHKERRQ(PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den)));
253   }
254   PetscFunctionReturn(0);
255 }
256 
257 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
258 {
259   PetscFunctionBeginUser;
260   CHKERRQ(DMCreate(comm, dm));
261   CHKERRQ(DMSetType(*dm, DMPLEX));
262   CHKERRQ(DMSetFromOptions(*dm));
263   CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view"));
264 
265   CHKERRQ(DMGetBoundingBox(*dm, ctx->lower, ctx->upper));
266   ctx->a = (ctx->upper[0] - ctx->lower[0])/2.0;
267   ctx->b = (ctx->upper[1] - ctx->lower[1])/2.0;
268   PetscFunctionReturn(0);
269 }
270 
271 static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
272 {
273   AppCtx *lctx = (AppCtx*)ctx;
274   u[0] = 2.*lctx->a + x[0];
275   return 0;
276 }
277 
278 static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
279 {
280   u[0] = 0.0;
281   return 0;
282 }
283 
284 static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
285 {
286   AppCtx *lctx = (AppCtx*)ctx;
287   /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z.  The stability
288      is analytically known and reported in ProcessOptions. */
289   if (lctx->ke!=0.0) {
290     u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
291   } else {
292     u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
293   }
294   return 0;
295 }
296 
297 static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
298 {
299   u[0] = 0.0;
300   return 0;
301 }
302 
303 static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
304 {
305   u[0] = 0.0;
306   return 0;
307 }
308 
309 static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
310 {
311   AppCtx *ctx = (AppCtx*)a_ctx;
312   PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
313   if (x[0] == ctx->lower[0] || x[0] == ctx->upper[0]) r = 0;
314   u[0] = r;
315   /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
316   return 0;
317 }
318 
319 static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
320 {
321   u[0] = 0.0;
322   return 0;
323 }
324 
325 static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
326 {
327   u[0] = 0.0;
328   return 0;
329 }
330 
331 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
332 {
333   PetscDS        ds;
334   DMLabel        label;
335   const PetscInt id = 1;
336 
337   PetscFunctionBeginUser;
338   CHKERRQ(DMGetLabel(dm, "marker", &label));
339   CHKERRQ(DMGetDS(dm, &ds));
340   CHKERRQ(PetscDSSetResidual(ds, 0, f0_n,     f1_n));
341   CHKERRQ(PetscDSSetResidual(ds, 1, f0_Omega, f1_Omega));
342   CHKERRQ(PetscDSSetResidual(ds, 2, f0_psi,   f1_psi));
343   CHKERRQ(PetscDSSetResidual(ds, 3, f0_phi,   f1_phi));
344   CHKERRQ(PetscDSSetResidual(ds, 4, f0_jz,    f1_jz));
345   ctx->initialFuncs[0] = initialSolution_n;
346   ctx->initialFuncs[1] = initialSolution_Omega;
347   ctx->initialFuncs[2] = initialSolution_psi;
348   ctx->initialFuncs[3] = initialSolution_phi;
349   ctx->initialFuncs[4] = initialSolution_jz;
350   for (PetscInt f = 0; f < 5; ++f) {
351     CHKERRQ(PetscDSSetImplicit(ds, f, ctx->implicit));
352     CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, ctx, NULL));
353   }
354   CHKERRQ(PetscDSSetContext(ds, 0, ctx));
355   PetscFunctionReturn(0);
356 }
357 
358 static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
359 {
360   PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
361   Vec            eq;
362   AppCtx *ctxarr[3];
363 
364   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
365   PetscFunctionBegin;
366   CHKERRQ(DMCreateLocalVector(dmAux, &eq));
367   CHKERRQ(DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq));
368   CHKERRQ(DMSetAuxiliaryVec(dm, NULL, 0, 0, eq));
369   if (ctx->plotRef) {  /* plot reference functions */
370     PetscViewer       viewer = NULL;
371     PetscBool         isHDF5,isVTK;
372     char              buf[256];
373     Vec               global;
374     PetscInt          dim;
375 
376     CHKERRQ(DMGetDimension(dm, &dim));
377     CHKERRQ(DMCreateGlobalVector(dmAux,&global));
378     CHKERRQ(VecSet(global,.0)); /* BCs! */
379     CHKERRQ(DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global));
380     CHKERRQ(DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global));
381     CHKERRQ(PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer));
382 #ifdef PETSC_HAVE_HDF5
383     CHKERRQ(PetscViewerSetType(viewer,PETSCVIEWERHDF5));
384 #else
385     CHKERRQ(PetscViewerSetType(viewer,PETSCVIEWERVTK));
386 #endif
387     CHKERRQ(PetscViewerSetFromOptions(viewer));
388     CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5));
389     CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK));
390     if (isHDF5) {
391       CHKERRQ(PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", dim));
392     } else if (isVTK) {
393       CHKERRQ(PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", dim));
394       CHKERRQ(PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU));
395     }
396     CHKERRQ(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE));
397     CHKERRQ(PetscViewerFileSetName(viewer,buf));
398     if (isHDF5) CHKERRQ(DMView(dmAux,viewer));
399     /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
400     CHKERRQ(PetscObjectSetName((PetscObject) global, "u0"));
401     CHKERRQ(VecView(global,viewer));
402     CHKERRQ(PetscViewerDestroy(&viewer));
403     CHKERRQ(VecDestroy(&global));
404   }
405   CHKERRQ(VecDestroy(&eq));
406   PetscFunctionReturn(0);
407 }
408 
409 static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
410 {
411   DM             dmAux, coordDM;
412   PetscInt       f;
413 
414   PetscFunctionBegin;
415   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
416   CHKERRQ(DMGetCoordinateDM(dm, &coordDM));
417   if (!feAux) PetscFunctionReturn(0);
418   CHKERRQ(DMClone(dm, &dmAux));
419   CHKERRQ(DMSetCoordinateDM(dmAux, coordDM));
420   for (f = 0; f < NfAux; ++f) CHKERRQ(DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]));
421   CHKERRQ(DMCreateDS(dmAux));
422   CHKERRQ(SetupEquilibriumFields(dm, dmAux, user));
423   CHKERRQ(DMDestroy(&dmAux));
424   PetscFunctionReturn(0);
425 }
426 
427 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
428 {
429   DM              cdm = dm;
430   PetscFE         fe[5], feAux[3];
431   PetscInt        dim, Nf = 5, NfAux = 3, f;
432   PetscBool       simplex;
433   MPI_Comm        comm;
434 
435   PetscFunctionBeginUser;
436   /* Create finite element */
437   CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm));
438   CHKERRQ(DMGetDimension(dm, &dim));
439   CHKERRQ(DMPlexIsSimplex(dm, &simplex));
440   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[0]));
441   CHKERRQ(PetscObjectSetName((PetscObject) fe[0], "density"));
442   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[1]));
443   CHKERRQ(PetscObjectSetName((PetscObject) fe[1], "vorticity"));
444   CHKERRQ(PetscFECopyQuadrature(fe[0], fe[1]));
445   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[2]));
446   CHKERRQ(PetscObjectSetName((PetscObject) fe[2], "flux"));
447   CHKERRQ(PetscFECopyQuadrature(fe[0], fe[2]));
448   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[3]));
449   CHKERRQ(PetscObjectSetName((PetscObject) fe[3], "potential"));
450   CHKERRQ(PetscFECopyQuadrature(fe[0], fe[3]));
451   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &fe[4]));
452   CHKERRQ(PetscObjectSetName((PetscObject) fe[4], "current"));
453   CHKERRQ(PetscFECopyQuadrature(fe[0], fe[4]));
454 
455   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[0]));
456   CHKERRQ(PetscObjectSetName((PetscObject) feAux[0], "n_0"));
457   CHKERRQ(PetscFECopyQuadrature(fe[0], feAux[0]));
458   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[1]));
459   CHKERRQ(PetscObjectSetName((PetscObject) feAux[1], "vorticity_0"));
460   CHKERRQ(PetscFECopyQuadrature(fe[0], feAux[1]));
461   CHKERRQ(PetscFECreateDefault(comm, dim, 1, simplex, NULL, -1, &feAux[2]));
462   CHKERRQ(PetscObjectSetName((PetscObject) feAux[2], "flux_0"));
463   CHKERRQ(PetscFECopyQuadrature(fe[0], feAux[2]));
464   /* Set discretization and boundary conditions for each mesh */
465   for (f = 0; f < Nf; ++f) CHKERRQ(DMSetField(dm, f, NULL, (PetscObject) fe[f]));
466   CHKERRQ(DMCreateDS(dm));
467   CHKERRQ(SetupProblem(dm, ctx));
468   while (cdm) {
469     CHKERRQ(SetupAuxDM(dm, NfAux, feAux, ctx));
470     CHKERRQ(DMCopyDisc(dm, cdm));
471     CHKERRQ(DMGetCoarseDM(cdm, &cdm));
472   }
473   for (f = 0; f < Nf; ++f) CHKERRQ(PetscFEDestroy(&fe[f]));
474   for (f = 0; f < NfAux; ++f) CHKERRQ(PetscFEDestroy(&feAux[f]));
475   PetscFunctionReturn(0);
476 }
477 
478 int main(int argc, char **argv)
479 {
480   DM             dm;
481   TS             ts;
482   Vec            u, r;
483   AppCtx         ctx;
484   PetscReal      t       = 0.0;
485   PetscReal      L2error = 0.0;
486   AppCtx        *ctxarr[5];
487 
488   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
489   s_ctx = &ctx;
490   CHKERRQ(PetscInitialize(&argc, &argv, NULL,help));
491   CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &ctx));
492   /* create mesh and problem */
493   CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
494   CHKERRQ(DMSetApplicationContext(dm, &ctx));
495   CHKERRQ(PetscMalloc1(5, &ctx.initialFuncs));
496   CHKERRQ(SetupDiscretization(dm, &ctx));
497   CHKERRQ(DMCreateGlobalVector(dm, &u));
498   CHKERRQ(PetscObjectSetName((PetscObject) u, "u"));
499   CHKERRQ(VecDuplicate(u, &r));
500   CHKERRQ(PetscObjectSetName((PetscObject) r, "r"));
501   /* create TS */
502   CHKERRQ(TSCreate(PETSC_COMM_WORLD, &ts));
503   CHKERRQ(TSSetDM(ts, dm));
504   CHKERRQ(TSSetApplicationContext(ts, &ctx));
505   CHKERRQ(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
506   if (ctx.implicit) {
507     CHKERRQ(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
508     CHKERRQ(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
509   } else {
510     CHKERRQ(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx));
511   }
512   CHKERRQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
513   CHKERRQ(TSSetFromOptions(ts));
514   CHKERRQ(TSSetPostStep(ts, PostStep));
515   /* make solution & solve */
516   CHKERRQ(DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u));
517   CHKERRQ(TSSetSolution(ts,u));
518   CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view"));
519   CHKERRQ(PostStep(ts)); /* print the initial state */
520   CHKERRQ(TSSolve(ts, u));
521   CHKERRQ(TSGetTime(ts, &t));
522   CHKERRQ(DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error));
523   if (L2error < 1.0e-11) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n"));
524   else                   CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error));
525   CHKERRQ(VecDestroy(&u));
526   CHKERRQ(VecDestroy(&r));
527   CHKERRQ(TSDestroy(&ts));
528   CHKERRQ(DMDestroy(&dm));
529   CHKERRQ(PetscFree(ctx.initialFuncs));
530   CHKERRQ(PetscFinalize());
531   return 0;
532 }
533 
534 /*TEST
535 
536   test:
537     suffix: 0
538     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 \
539           -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
540   test:
541     # Remapping with periodicity is broken
542     suffix: 1
543     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 \
544            -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
545 
546 TEST*/
547