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";

 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: }

269: static void f_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
270:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
271:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
272:                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
273: {
274:   const PetscScalar *pn = &u[uOff[DENSITY]];
275:   *f0 = *pn;
276: }

280: static PetscErrorCode PostStep(TS ts)
281: {
282:   PetscErrorCode    ierr;
283:   DM                dm;
284:   AppCtx            *ctx;
285:   PetscInt          stepi,num;
286:   Vec               X;
288:   TSGetApplicationContext(ts, &ctx); assert(ctx);
289:   if (ctx->debug<1) return(0);
290:   TSGetSolution(ts, &X);
291:   VecGetDM(X, &dm);
292:   TSGetStepNumber(ts, &stepi);
293:   DMGetOutputSequenceNumber(dm, &num, NULL);
294:   if (num < 0) {DMSetOutputSequenceNumber(dm, 0, 0.0);}
295:   PetscObjectSetName((PetscObject) X, "u");
296:   VecViewFromOptions(X, NULL, "-vec_view");
297:   /* print integrals */
298:   {
299:     PetscDS          prob;
300:     DM               plex;
301:     PetscScalar den, tt[5];
302:     DMConvert(dm, DMPLEX, &plex);
303:     DMGetDS(plex, &prob);
304:     PetscDSSetObjective(prob, 0, &f_n);
305:     DMPlexComputeIntegralFEM(plex,X,tt,ctx);
306:     den = tt[0];
307:     DMDestroy(&plex);
308:     PetscPrintf(PetscObjectComm((PetscObject)dm), "%D) total perturbed mass = %g\n", stepi, (double) PetscRealPart(den));
309:   }
310:   return(0);
311: }

313: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
314: {
315:   DMLabel        label;
318:   DMCreateLabel(dm, name);
319:   DMGetLabel(dm, name, &label);
320:   DMPlexMarkBoundaryFaces(dm, 1, label);
321:   DMPlexLabelComplete(dm, label);
322:   return(0);
323: }

325: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
326: {
327:   PetscInt       dim      = ctx->dim;
328:   const char    *filename = ctx->filename;
329:   size_t         len;
330:   PetscMPIInt    numProcs;

334:   MPI_Comm_size(comm, &numProcs);
335:   PetscStrlen(filename, &len);
336:   if (len) {
337:     DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);
338:   } else {
339:     PetscInt        d;

341:     /* create DM */
342:     if (ctx->cell_simplex && dim == 3) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh a cylinder with simplices");
343:     if (dim==2) {
344:       PetscInt refineRatio, totCells = 1;
345:       if (ctx->cell_simplex) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot mesh 2D with simplices");
346:       refineRatio = PetscMax((PetscInt) (PetscPowReal(numProcs, 1.0/dim) + 0.1) - 1, 1);
347:       for (d = 0; d < dim; ++d) {
348:         if (ctx->cells[d] < refineRatio) ctx->cells[d] = refineRatio;
349:         if (ctx->periodicity[d]==DM_BOUNDARY_PERIODIC && ctx->cells[d]*refineRatio <= 2) refineRatio = 2;
350:       }
351:       for (d = 0; d < dim; ++d) {
352:         ctx->cells[d] *= refineRatio;
353:         totCells *= ctx->cells[d];
354:       }
355:       if (totCells % numProcs) SETERRQ2(comm,PETSC_ERR_ARG_WRONG,"Total cells %D not divisible by processes %D", totCells, numProcs);
356:       DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, ctx->cells, ctx->domain_lo, ctx->domain_hi, ctx->periodicity, PETSC_TRUE, dm);
357:     } else {
358:       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");
359:       /* we stole dm_refine so clear it */
360:       PetscOptionsClearValue(NULL,"-dm_refine");
361:       DMPlexCreateHexCylinderMesh(comm, ctx->refine, ctx->periodicity[2], dm);
362:     }
363:   }
364:   {
365:     DM distributedMesh = NULL;
366:     /* Distribute mesh over processes */
367:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
368:     if (distributedMesh) {
369:       DMDestroy(dm);
370:       *dm  = distributedMesh;
371:     }
372:   }
373:   {
374:     PetscBool hasLabel;
375:     DMHasLabel(*dm, "marker", &hasLabel);
376:     if (!hasLabel) {CreateBCLabel(*dm, "marker");}
377:   }
378:   {
379:     char      convType[256];
380:     PetscBool flg;
381:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
382:     PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex48",DMList,DMPLEX,convType,256,&flg);
383:     PetscOptionsEnd();
384:     if (flg) {
385:       DM dmConv;
386:       DMConvert(*dm,convType,&dmConv);
387:       if (dmConv) {
388:         DMDestroy(dm);
389:         *dm  = dmConv;
390:       }
391:     }
392:   }
393:   PetscObjectSetName((PetscObject) *dm, "Mesh");
394:   DMSetFromOptions(*dm);
395:   DMLocalizeCoordinates(*dm); /* needed for periodic */
396:   return(0);
397: }

399: static PetscErrorCode log_n_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
400: {
401:   AppCtx *lctx = (AppCtx*)ctx;
402:   assert(ctx);
403:   u[0] = (lctx->domain_hi-lctx->domain_lo)+x[0];
404:   return 0;
405: }

407: static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
408: {
409:   u[0] = 0.0;
410:   return 0;
411: }

413: static PetscErrorCode psi_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
414: {
415:   AppCtx *lctx = (AppCtx*)ctx;
416:   assert(ctx);
417:   /* This sets up a symmetrix By flux aroound the mid point in x, which represents a current density flux along z.  The stability
418:      is analytically known and reported in ProcessOptions. */
419:   if (lctx->ke!=0.0) {
420:     u[0] = (PetscCosReal(lctx->ke*(x[0]-lctx->a))-PetscCosReal(lctx->ke*lctx->a))/(1.0-PetscCosReal(lctx->ke*lctx->a));
421:   } else {
422:     u[0] = 1.0-PetscPowScalar((x[0]-lctx->a)/lctx->a,2);
423:   }
424:   return 0;
425: }

427: static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
428: {
429:   u[0] = 0.0;
430:   return 0;
431: }

433: static PetscErrorCode initialSolution_Omega(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
434: {
435:   u[0] = 0.0;
436:   return 0;
437: }

439: static PetscErrorCode initialSolution_psi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *a_ctx)
440: {
441:   AppCtx *ctx = (AppCtx*)a_ctx;
442:   PetscScalar r = ctx->eps*(PetscScalar) (rand()) / (PetscScalar) (RAND_MAX);
443:   assert(ctx);
444:   if (x[0] == ctx->domain_lo[0] || x[0] == ctx->domain_hi[0]) r = 0;
445:   u[0] = r;
446:   /* PetscPrintf(PETSC_COMM_WORLD, "rand psi %lf\n",u[0]); */
447:   return 0;
448: }

450: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
451: {
452:   u[0] = 0.0;
453:   return 0;
454: }

456: static PetscErrorCode initialSolution_jz(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
457: {
458:   u[0] = 0.0;
459:   return 0;
460: }

462: static PetscErrorCode SetupProblem(PetscDS prob, AppCtx *ctx)
463: {
464:   const PetscInt id = 1;
465:   PetscErrorCode ierr, f;

468:   PetscDSSetResidual(prob, 0, f0_n,     f1_n);
469:   PetscDSSetResidual(prob, 1, f0_Omega, f1_Omega);
470:   PetscDSSetResidual(prob, 2, f0_psi,   f1_psi);
471:   PetscDSSetResidual(prob, 3, f0_phi,   f1_phi);
472:   PetscDSSetResidual(prob, 4, f0_jz,    f1_jz);
473:   ctx->initialFuncs[0] = initialSolution_n;
474:   ctx->initialFuncs[1] = initialSolution_Omega;
475:   ctx->initialFuncs[2] = initialSolution_psi;
476:   ctx->initialFuncs[3] = initialSolution_phi;
477:   ctx->initialFuncs[4] = initialSolution_jz;
478:   for (f = 0; f < 5; ++f) {
479:     PetscDSSetImplicit( prob, f, ctx->implicit);
480:     PetscDSAddBoundary( prob, DM_BC_ESSENTIAL, "wall", "marker", f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], 1, &id, ctx);
481:   }
482:   PetscDSSetContext(prob, 0, ctx);
483:   PetscDSSetFromOptions(prob);
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 SetupDiscretization(DM dm, AppCtx *ctx)
537: {
538:   DM              cdm = dm;
539:   const PetscInt  dim = ctx->dim;
540:   PetscQuadrature q;
541:   PetscFE         fe[5], feAux[3];
542:   PetscDS         prob, probAux;
543:   PetscInt        Nf = 5, NfAux = 3, f;
544:   PetscBool       cell_simplex = ctx->cell_simplex;
545:   MPI_Comm        comm;
546:   PetscErrorCode  ierr;

549:   /* Create finite element */
550:   PetscObjectGetComm((PetscObject) dm, &comm);
551:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);
552:   PetscObjectSetName((PetscObject) fe[0], "density");
553:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);
554:   PetscObjectSetName((PetscObject) fe[1], "vorticity");
555:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);
556:   PetscObjectSetName((PetscObject) fe[2], "flux");
557:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);
558:   PetscObjectSetName((PetscObject) fe[3], "potential");
559:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);
560:   PetscObjectSetName((PetscObject) fe[4], "current");

562:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);
563:   PetscFEGetQuadrature(fe[0], &q);
564:   PetscFESetQuadrature(feAux[0], q);
565:   PetscObjectSetName((PetscObject) feAux[0], "n_0");
566:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);
567:   PetscFEGetQuadrature(fe[1], &q);
568:   PetscFESetQuadrature(feAux[1], q);
569:   PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");
570:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);
571:   PetscFEGetQuadrature(fe[2], &q);
572:   PetscFESetQuadrature(feAux[2], q);
573:   PetscObjectSetName((PetscObject) feAux[2], "flux_0");
574:   /* Set discretization and boundary conditions for each mesh */
575:   DMGetDS(dm, &prob);
576:   for (f = 0; f < Nf; ++f) {PetscDSSetDiscretization(prob, f, (PetscObject) fe[f]);}
577:   PetscDSCreate(PetscObjectComm((PetscObject) dm), &probAux);
578:   for (f = 0; f < NfAux; ++f) {PetscDSSetDiscretization(probAux, f, (PetscObject) feAux[f]);}
579:   SetupProblem(prob, ctx);
580:   while (cdm) {
581:     DM coordDM, dmAux;

583:     DMSetDS(cdm,prob);
584:     DMGetCoordinateDM(cdm,&coordDM);
585:     {
586:       PetscBool hasLabel;

588:       DMHasLabel(cdm, "marker", &hasLabel);
589:       if (!hasLabel) {CreateBCLabel(cdm, "marker");}
590:     }

592:     DMClone(cdm, &dmAux);
593:     DMSetCoordinateDM(dmAux, coordDM);
594:     DMSetDS(dmAux, probAux);
595:     PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
596:     SetupEquilibriumFields(cdm, dmAux, ctx);
597:     DMDestroy(&dmAux);

599:     DMGetCoarseDM(cdm, &cdm);
600:   }
601:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&fe[f]);}
602:   for (f = 0; f < NfAux; ++f) {PetscFEDestroy(&feAux[f]);}
603:   PetscDSDestroy(&probAux);
604:   return(0);
605: }

607: int main(int argc, char **argv)
608: {
609:   DM             dm;
610:   TS             ts;
611:   Vec            u, r;
612:   AppCtx         ctx;
613:   PetscReal      t       = 0.0;
614:   PetscReal      L2error = 0.0;
616:   AppCtx        *ctxarr[5];

618:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctxarr[3] = ctxarr[4] = &ctx; /* each variable could have a different context */
619:   s_ctx = &ctx;
620:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
621:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
622:   /* create mesh and problem */
623:   CreateMesh(PETSC_COMM_WORLD, &ctx, &dm);
624:   DMSetApplicationContext(dm, &ctx);
625:   PetscMalloc1(5, &ctx.initialFuncs);
626:   SetupDiscretization(dm, &ctx);
627:   DMCreateGlobalVector(dm, &u);
628:   PetscObjectSetName((PetscObject) u, "u");
629:   VecDuplicate(u, &r);
630:   PetscObjectSetName((PetscObject) r, "r");
631:   /* create TS */
632:   TSCreate(PETSC_COMM_WORLD, &ts);
633:   TSSetDM(ts, dm);
634:   TSSetApplicationContext(ts, &ctx);
635:   DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
636:   if (ctx.implicit) {
637:     DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
638:     DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
639:   } else {
640:     DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &ctx);
641:   }
642:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
643:   TSSetFromOptions(ts);
644:   TSSetPostStep(ts, PostStep);
645:   /* make solution & solve */
646:   DMProjectFunction(dm, t, ctx.initialFuncs, (void **)ctxarr, INSERT_ALL_VALUES, u);
647:   TSSetSolution(ts,u);
648:   DMViewFromOptions(dm, NULL, "-dm_view");
649:   PostStep(ts); /* print the initial state */
650:   TSSolve(ts, u);
651:   TSGetTime(ts, &t);
652:   DMComputeL2Diff(dm, t, ctx.initialFuncs, (void **)ctxarr, u, &L2error);
653:   if (L2error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
654:   else                   {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", L2error);}
655: #if 0
656:   {
657:     PetscReal res = 0.0;
658:     /* Check discretization error */
659:     VecViewFromOptions(u, NULL, "-initial_guess_view");
660:     DMComputeL2Diff(dm, 0.0, ctx.exactFuncs, NULL, u, &error);
661:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
662:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
663:     /* Check residual */
664:     SNESComputeFunction(snes, u, r);
665:     VecChop(r, 1.0e-10);
666:     VecViewFromOptions(r, NULL, "-initial_residual_view");
667:     VecNorm(r, NORM_2, &res);
668:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
669:     /* Check Jacobian */
670:     {
671:       Mat A;
672:       Vec b;

674:       SNESGetJacobian(snes, &A, NULL, NULL, NULL);
675:       SNESComputeJacobian(snes, u, A, A);
676:       VecDuplicate(u, &b);
677:       VecSet(r, 0.0);
678:       SNESComputeFunction(snes, r, b);
679:       MatMult(A, u, r);
680:       VecAXPY(r, 1.0, b);
681:       VecDestroy(&b);
682:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
683:       VecChop(r, 1.0e-10);
684:       VecViewFromOptions(r, NULL, "-linear_residual_view");
685:       VecNorm(r, NORM_2, &res);
686:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
687:     }
688:   }
689: #endif
690:   VecDestroy(&u);
691:   VecDestroy(&r);
692:   TSDestroy(&ts);
693:   DMDestroy(&dm);
694:   PetscFree(ctx.initialFuncs);
695:   PetscFinalize();
696:   return ierr;
697: }

699: /*TEST

701:   test:
702:     suffix: 0
703:     args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_final_time 10. -ts_dt 1.0
704:   test:
705:     suffix: 1
706:     args: -debug 1 -dim 3 -dm_refine 1 -z_periodicity PERIODIC -ts_max_steps 1 -ts_final_time 10. -ts_dt 1.0 -domain_lo -2,-1,-1 -domain_hi 2,1,1

708: TEST*/