xref: /petsc/src/ts/tutorials/ex48.c (revision bef158480efac06de457f7a665168877ab3c2fd7) 
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 #include <assert.h>
21 
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;
43 
44 static AppCtx *s_ctx;
45 
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 }
51 
52 enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};
53 
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 }
67 
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;
75 
76   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
77 }
78 
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]];
88 
89   f0[0] += - poissonBracket(dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer);
90   if (u_t) f0[0] += u_t[OMEGA];
91 }
92 
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;
100 
101   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
102 }
103 
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 }
125 
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;
133 
134   for (d = 0; d < dim-1; ++d) f1[d] = -(s_ctx->eta/s_ctx->beta)*ppsi[d];
135 }
136 
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 }
144 
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;
152 
153   for (d = 0; d < dim-1; ++d) f1[d] = pphi[d];
154 }
155 
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 }
163 
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;
172 
173   for (d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
174 }
175 
176 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
177 {
178   PetscBool      flg;
179   PetscErrorCode ierr;
180   PetscInt       ii, bd;
181   PetscFunctionBeginUser;
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;
206 
207   for (ii = 0; ii < options->dim; ++ii) options->cells[ii] = 4;
208   ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr);
209   ierr = PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);CHKERRQ(ierr);
210   ierr = PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);CHKERRQ(ierr);
211   ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex48.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
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);CHKERRQ(ierr);
213   ierr = PetscOptionsInt("-dm_refine", "Hack to get refinement level for cylinder", "ex48.c", options->refine, &options->refine, NULL);CHKERRQ(ierr);
214   ierr = PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);CHKERRQ(ierr);
215   ierr = PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);CHKERRQ(ierr);
216   ierr = PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);CHKERRQ(ierr);
217   ierr = PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);CHKERRQ(ierr);
218   ierr = PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);CHKERRQ(ierr);
219   ierr = PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);CHKERRQ(ierr);
220   ierr = PetscOptionsString("-f", "Exodus.II filename to read", "ex48.c", options->filename, options->filename, sizeof(options->filename), &flg);CHKERRQ(ierr);
221   ierr = PetscOptionsBool("-cell_simplex", "Simplicial (true) or tensor (false) mesh", "ex48.c", options->cell_simplex, &options->cell_simplex, NULL);CHKERRQ(ierr);
222   ierr = PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);CHKERRQ(ierr);
223   ii = options->dim;
224   ierr = PetscOptionsRealArray("-domain_hi", "Domain size", "ex48.c", options->domain_hi, &ii, NULL);CHKERRQ(ierr);
225   ii = options->dim;
226   ierr = PetscOptionsRealArray("-domain_lo", "Domain size", "ex48.c", options->domain_lo, &ii, NULL);CHKERRQ(ierr);
227   ii = options->dim;
228   bd = options->periodicity[0];
229   ierr = PetscOptionsEList("-x_periodicity", "The x-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);CHKERRQ(ierr);
230   options->periodicity[0] = (DMBoundaryType) bd;
231   bd = options->periodicity[1];
232   ierr = PetscOptionsEList("-y_periodicity", "The y-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);CHKERRQ(ierr);
233   options->periodicity[1] = (DMBoundaryType) bd;
234   bd = options->periodicity[2];
235   ierr = PetscOptionsEList("-z_periodicity", "The z-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);CHKERRQ(ierr);
236   options->periodicity[2] = (DMBoundaryType) bd;
237   ii = options->dim;
238   ierr = PetscOptionsIntArray("-cells", "Number of cells in each dimension", "ex48.c", options->cells, &ii, NULL);CHKERRQ(ierr);
239   ierr = PetscOptionsEnd();CHKERRQ(ierr);
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   ierr = PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime);CHKERRQ(ierr);
263 
264   PetscFunctionReturn(0);
265 }
266 
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 }
275 
276 static PetscErrorCode PostStep(TS ts)
277 {
278   PetscErrorCode    ierr;
279   DM                dm;
280   AppCtx            *ctx;
281   PetscInt          stepi,num;
282   Vec               X;
283   PetscFunctionBegin;
284   ierr = TSGetApplicationContext(ts, &ctx);CHKERRQ(ierr); assert(ctx);
285   if (ctx->debug<1) PetscFunctionReturn(0);
286   ierr = TSGetSolution(ts, &X);CHKERRQ(ierr);
287   ierr = VecGetDM(X, &dm);CHKERRQ(ierr);
288   ierr = TSGetStepNumber(ts, &stepi);CHKERRQ(ierr);
289   ierr = DMGetOutputSequenceNumber(dm, &num, NULL);CHKERRQ(ierr);
290   if (num < 0) {ierr = DMSetOutputSequenceNumber(dm, 0, 0.0);CHKERRQ(ierr);}
291   ierr = PetscObjectSetName((PetscObject) X, "u");CHKERRQ(ierr);
292   ierr = VecViewFromOptions(X, NULL, "-vec_view");CHKERRQ(ierr);
293   /* print integrals */
294   {
295     PetscDS          prob;
296     DM               plex;
297     PetscScalar den, tt[5];
298     ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
299     ierr = DMGetDS(plex, &prob);CHKERRQ(ierr);
300     ierr = PetscDSSetObjective(prob, 0, &f_n);CHKERRQ(ierr);
301     ierr = DMPlexComputeIntegralFEM(plex,X,tt,ctx);CHKERRQ(ierr);
302     den = tt[0];
303     ierr = DMDestroy(&plex);CHKERRQ(ierr);
304     PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));CHKERRQ(ierr);
305   }
306   PetscFunctionReturn(0);
307 }
308 
309 static PetscErrorCode CreateBCLabel(DM dm, const char name[])
310 {
311   DM             plex;
312   DMLabel        label;
313   PetscErrorCode ierr;
314 
315   PetscFunctionBeginUser;
316   ierr = DMCreateLabel(dm, name);CHKERRQ(ierr);
317   ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr);
318   ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
319   ierr = DMPlexMarkBoundaryFaces(dm, 1, label);CHKERRQ(ierr);
320   ierr = DMDestroy(&plex);CHKERRQ(ierr);
321   PetscFunctionReturn(0);
322 }
323 
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;
330   PetscErrorCode ierr;
331 
332   PetscFunctionBeginUser;
333   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
334   ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
335   if (len) {
336     ierr = DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);CHKERRQ(ierr);
337   } else {
338     PetscInt        d;
339 
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       ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, ctx->cells, ctx->domain_lo, ctx->domain_hi, ctx->periodicity, PETSC_TRUE, dm);CHKERRQ(ierr);
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       ierr = PetscOptionsClearValue(NULL,"-dm_refine");CHKERRQ(ierr);
360       ierr = DMPlexCreateHexCylinderMesh(comm, ctx->refine, ctx->periodicity[2], dm);CHKERRQ(ierr);
361     }
362   }
363   {
364     DM distributedMesh = NULL;
365     /* Distribute mesh over processes */
366     ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
367     if (distributedMesh) {
368       ierr = DMDestroy(dm);CHKERRQ(ierr);
369       *dm  = distributedMesh;
370     }
371   }
372   {
373     PetscBool hasLabel;
374     ierr = DMHasLabel(*dm, "marker", &hasLabel);CHKERRQ(ierr);
375     if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);}
376   }
377   {
378     char      convType[256];
379     PetscBool flg;
380     ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr);
381     ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex48",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr);
382     ierr = PetscOptionsEnd();CHKERRQ(ierr);
383     if (flg) {
384       DM dmConv;
385       ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr);
386       if (dmConv) {
387         ierr = DMDestroy(dm);CHKERRQ(ierr);
388         *dm  = dmConv;
389       }
390     }
391   }
392   ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
393   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
394   ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */
395   PetscFunctionReturn(0);
396 }
397 
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 }
405 
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 }
411 
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 }
425 
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 }
431 
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 }
437 
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 }
448 
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 }
454 
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 }
460 
461 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
462 {
463   PetscDS        prob;
464   const PetscInt id = 1;
465   PetscErrorCode ierr, f;
466 
467   PetscFunctionBeginUser;
468   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
469   ierr = PetscDSSetResidual(prob, 0, f0_n,     f1_n);CHKERRQ(ierr);
470   ierr = PetscDSSetResidual(prob, 1, f0_Omega, f1_Omega);CHKERRQ(ierr);
471   ierr = PetscDSSetResidual(prob, 2, f0_psi,   f1_psi);CHKERRQ(ierr);
472   ierr = PetscDSSetResidual(prob, 3, f0_phi,   f1_phi);CHKERRQ(ierr);
473   ierr = PetscDSSetResidual(prob, 4, f0_jz,    f1_jz);CHKERRQ(ierr);
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     ierr = PetscDSSetImplicit(prob, f, ctx->implicit);CHKERRQ(ierr);
481     ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], NULL, 1, &id, ctx);CHKERRQ(ierr);
482   }
483   ierr = PetscDSSetContext(prob, 0, ctx);CHKERRQ(ierr);
484   PetscFunctionReturn(0);
485 }
486 
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;
491   PetscErrorCode ierr;
492   AppCtx *ctxarr[3];
493 
494   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
495   PetscFunctionBegin;
496   ierr = DMCreateLocalVector(dmAux, &eq);CHKERRQ(ierr);
497   ierr = DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);CHKERRQ(ierr);
498   ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) eq);CHKERRQ(ierr);
499   if (ctx->plotRef) {  /* plot reference functions */
500     PetscViewer       viewer = NULL;
501     PetscBool         isHDF5,isVTK;
502     char              buf[256];
503     Vec               global;
504     ierr = DMCreateGlobalVector(dmAux,&global);CHKERRQ(ierr);
505     ierr = VecSet(global,.0);CHKERRQ(ierr); /* BCs! */
506     ierr = DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);CHKERRQ(ierr);
507     ierr = DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);CHKERRQ(ierr);
508     ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);CHKERRQ(ierr);
509 #ifdef PETSC_HAVE_HDF5
510     ierr = PetscViewerSetType(viewer,PETSCVIEWERHDF5);CHKERRQ(ierr);
511 #else
512     ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
513 #endif
514     ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr);
515     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);CHKERRQ(ierr);
516     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);CHKERRQ(ierr);
517     if (isHDF5) {
518       ierr = PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", ctx->dim);CHKERRQ(ierr);
519     } else if (isVTK) {
520       ierr = PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", ctx->dim);CHKERRQ(ierr);
521       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
522     }
523     ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
524     ierr = PetscViewerFileSetName(viewer,buf);CHKERRQ(ierr);
525     if (isHDF5) {ierr = DMView(dmAux,viewer);CHKERRQ(ierr);}
526     /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
527     ierr = PetscObjectSetName((PetscObject) global, "u0");CHKERRQ(ierr);
528     ierr = VecView(global,viewer);CHKERRQ(ierr);
529     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
530     ierr = VecDestroy(&global);CHKERRQ(ierr);
531   }
532   ierr = VecDestroy(&eq);CHKERRQ(ierr);
533   PetscFunctionReturn(0);
534 }
535 
536 static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
537 {
538   DM             dmAux, coordDM;
539   PetscInt       f;
540   PetscErrorCode ierr;
541 
542   PetscFunctionBegin;
543   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
544   ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr);
545   if (!feAux) PetscFunctionReturn(0);
546   ierr = DMClone(dm, &dmAux);CHKERRQ(ierr);
547   ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr);
548   ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr);
549   for (f = 0; f < NfAux; ++f) {ierr = DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);CHKERRQ(ierr);}
550   ierr = DMCreateDS(dmAux);CHKERRQ(ierr);
551   ierr = SetupEquilibriumFields(dm, dmAux, user);CHKERRQ(ierr);
552   ierr = DMDestroy(&dmAux);CHKERRQ(ierr);
553   PetscFunctionReturn(0);
554 }
555 
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;
565 
566   PetscFunctionBeginUser;
567   /* Create finite element */
568   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
569   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);CHKERRQ(ierr);
570   ierr = PetscObjectSetName((PetscObject) fe[0], "density");CHKERRQ(ierr);
571   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);CHKERRQ(ierr);
572   ierr = PetscObjectSetName((PetscObject) fe[1], "vorticity");CHKERRQ(ierr);
573   ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr);
574   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);CHKERRQ(ierr);
575   ierr = PetscObjectSetName((PetscObject) fe[2], "flux");CHKERRQ(ierr);
576   ierr = PetscFECopyQuadrature(fe[0], fe[2]);CHKERRQ(ierr);
577   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);CHKERRQ(ierr);
578   ierr = PetscObjectSetName((PetscObject) fe[3], "potential");CHKERRQ(ierr);
579   ierr = PetscFECopyQuadrature(fe[0], fe[3]);CHKERRQ(ierr);
580   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);CHKERRQ(ierr);
581   ierr = PetscObjectSetName((PetscObject) fe[4], "current");CHKERRQ(ierr);
582   ierr = PetscFECopyQuadrature(fe[0], fe[4]);CHKERRQ(ierr);
583 
584   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);CHKERRQ(ierr);
585   ierr = PetscObjectSetName((PetscObject) feAux[0], "n_0");CHKERRQ(ierr);
586   ierr = PetscFECopyQuadrature(fe[0], feAux[0]);CHKERRQ(ierr);
587   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);CHKERRQ(ierr);
588   ierr = PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");CHKERRQ(ierr);
589   ierr = PetscFECopyQuadrature(fe[0], feAux[1]);CHKERRQ(ierr);
590   ierr = PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);CHKERRQ(ierr);
591   ierr = PetscObjectSetName((PetscObject) feAux[2], "flux_0");CHKERRQ(ierr);
592   ierr = PetscFECopyQuadrature(fe[0], feAux[2]);CHKERRQ(ierr);
593   /* Set discretization and boundary conditions for each mesh */
594   for (f = 0; f < Nf; ++f) {ierr = DMSetField(dm, f, NULL, (PetscObject) fe[f]);CHKERRQ(ierr);}
595   ierr = DMCreateDS(dm);CHKERRQ(ierr);
596   ierr = SetupProblem(dm, ctx);CHKERRQ(ierr);
597   while (cdm) {
598     ierr = SetupAuxDM(dm, NfAux, feAux, ctx);CHKERRQ(ierr);
599     {
600       PetscBool hasLabel;
601 
602       ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr);
603       if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);}
604     }
605     ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
606     ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
607   }
608   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&fe[f]);CHKERRQ(ierr);}
609   for (f = 0; f < NfAux; ++f) {ierr = PetscFEDestroy(&feAux[f]);CHKERRQ(ierr);}
610   PetscFunctionReturn(0);
611 }
612 
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;
621   PetscErrorCode ierr;
622   AppCtx        *ctxarr[5];
623 
624   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
625   s_ctx = &ctx;
626   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
627   ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
628   /* create mesh and problem */
629   ierr = CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);CHKERRQ(ierr);
630   ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);
631   ierr = PetscMalloc1(5, &ctx.initialFuncs);CHKERRQ(ierr);
632   ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr);
633   ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
634   ierr = PetscObjectSetName((PetscObject) u, "u");CHKERRQ(ierr);
635   ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
636   ierr = PetscObjectSetName((PetscObject) r, "r");CHKERRQ(ierr);
637   /* create TS */
638   ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
639   ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
640   ierr = TSSetApplicationContext(ts, &ctx);CHKERRQ(ierr);
641   ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr);
642   if (ctx.implicit) {
643     ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr);
644     ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr);
645   } else {
646     ierr = DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);CHKERRQ(ierr);
647   }
648   ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
649   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
650   ierr = TSSetPostStep(ts, PostStep);CHKERRQ(ierr);
651   /* make solution & solve */
652   ierr = DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
653   ierr = TSSetSolution(ts,u);CHKERRQ(ierr);
654   ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
655   ierr = PostStep(ts);CHKERRQ(ierr); /* print the initial state */
656   ierr = TSSolve(ts, u);CHKERRQ(ierr);
657   ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
658   ierr = DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);CHKERRQ(ierr);
659   if (L2error < 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);}
660   else                   {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);CHKERRQ(ierr);}
661 #if 0
662   {
663     PetscReal res = 0.0;
664     /* Check discretization error */
665     ierr = VecViewFromOptions(u, NULL, "-initial_guess_view");CHKERRQ(ierr);
666     ierr = DMComputeL2Diff(dm, 0.0, ctx.exactFuncs, NULL, u, &error);CHKERRQ(ierr);
667     if (error < 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);}
668     else                 {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);CHKERRQ(ierr);}
669     /* Check residual */
670     ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr);
671     ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
672     ierr = VecViewFromOptions(r, NULL, "-initial_residual_view");CHKERRQ(ierr);
673     ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
674     ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);CHKERRQ(ierr);
675     /* Check Jacobian */
676     {
677       Mat A;
678       Vec b;
679 
680       ierr = SNESGetJacobian(snes, &A, NULL, NULL, NULL);CHKERRQ(ierr);
681       ierr = SNESComputeJacobian(snes, u, A, A);CHKERRQ(ierr);
682       ierr = VecDuplicate(u, &b);CHKERRQ(ierr);
683       ierr = VecSet(r, 0.0);CHKERRQ(ierr);
684       ierr = SNESComputeFunction(snes, r, b);CHKERRQ(ierr);
685       ierr = MatMult(A, u, r);CHKERRQ(ierr);
686       ierr = VecAXPY(r, 1.0, b);CHKERRQ(ierr);
687       ierr = VecDestroy(&b);CHKERRQ(ierr);
688       ierr = PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");CHKERRQ(ierr);
689       ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
690       ierr = VecViewFromOptions(r, NULL, "-linear_residual_view");CHKERRQ(ierr);
691       ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
692       ierr = PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);CHKERRQ(ierr);
693     }
694   }
695 #endif
696   ierr = VecDestroy(&u);CHKERRQ(ierr);
697   ierr = VecDestroy(&r);CHKERRQ(ierr);
698   ierr = TSDestroy(&ts);CHKERRQ(ierr);
699   ierr = DMDestroy(&dm);CHKERRQ(ierr);
700   ierr = PetscFree(ctx.initialFuncs);CHKERRQ(ierr);
701   ierr = PetscFinalize();
702   return ierr;
703 }
704 
705 /*TEST
706 
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
713 
714 TEST*/
715