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

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

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

309: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
310: {
311:   DMLabel        label;
314:   DMCreateLabel(dm, name);
315:   DMGetLabel(dm, name, &label);
316:   DMPlexMarkBoundaryFaces(dm, 1, label);
317:   DMPlexLabelComplete(dm, label);
318:   return(0);
319: }

321: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
322: {
323:   PetscInt       dim      = ctx->dim;
324:   const char    *filename = ctx->filename;
325:   size_t         len;
326:   PetscMPIInt    numProcs;

330:   MPI_Comm_size(comm, &numProcs);
331:   PetscStrlen(filename, &len);
332:   if (len) {
333:     DMPlexCreateFromFile(comm, filename, PETSC_TRUE, dm);
334:   } else {
335:     PetscInt        d;

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

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

403: static PetscErrorCode Omega_0(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
404: {
405:   u[0] = 0.0;
406:   return 0;
407: }

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

423: static PetscErrorCode initialSolution_n(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
424: {
425:   u[0] = 0.0;
426:   return 0;
427: }

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

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

446: static PetscErrorCode initialSolution_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
447: {
448:   u[0] = 0.0;
449:   return 0;
450: }

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

458: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
459: {
460:   PetscDS        prob;
461:   const PetscInt id = 1;
462:   PetscErrorCode ierr, f;

465:   DMGetDS(dm, &prob);
466:   PetscDSSetResidual(prob, 0, f0_n,     f1_n);
467:   PetscDSSetResidual(prob, 1, f0_Omega, f1_Omega);
468:   PetscDSSetResidual(prob, 2, f0_psi,   f1_psi);
469:   PetscDSSetResidual(prob, 3, f0_phi,   f1_phi);
470:   PetscDSSetResidual(prob, 4, f0_jz,    f1_jz);
471:   ctx->initialFuncs[0] = initialSolution_n;
472:   ctx->initialFuncs[1] = initialSolution_Omega;
473:   ctx->initialFuncs[2] = initialSolution_psi;
474:   ctx->initialFuncs[3] = initialSolution_phi;
475:   ctx->initialFuncs[4] = initialSolution_jz;
476:   for (f = 0; f < 5; ++f) {
477:     PetscDSSetImplicit( prob, f, ctx->implicit);
478:     PetscDSAddBoundary( prob, DM_BC_ESSENTIAL, "wall", "marker", f, 0, NULL, (void (*)(void)) ctx->initialFuncs[f], 1, &id, ctx);
479:   }
480:   PetscDSSetContext(prob, 0, ctx);
481:   return(0);
482: }

484: static PetscErrorCode SetupEquilibriumFields(DM dm, DM dmAux, AppCtx *ctx)
485: {
486:   PetscErrorCode (*eqFuncs[3])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar [], void *) = {log_n_0, Omega_0, psi_0};
487:   Vec            eq;
489:   AppCtx *ctxarr[3];

491:   ctxarr[0] = ctxarr[1] = ctxarr[2] = ctx; /* each variable could have a different context */
493:   DMCreateLocalVector(dmAux, &eq);
494:   DMProjectFunctionLocal(dmAux, 0.0, eqFuncs, (void **)ctxarr, INSERT_ALL_VALUES, eq);
495:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) eq);
496:   if (ctx->plotRef) {  /* plot reference functions */
497:     PetscViewer       viewer = NULL;
498:     PetscBool         isHDF5,isVTK;
499:     char              buf[256];
500:     Vec               global;
501:     DMCreateGlobalVector(dmAux,&global);
502:     VecSet(global,.0); /* BCs! */
503:     DMLocalToGlobalBegin(dmAux,eq,INSERT_VALUES,global);
504:     DMLocalToGlobalEnd(dmAux,eq,INSERT_VALUES,global);
505:     PetscViewerCreate(PetscObjectComm((PetscObject)dmAux),&viewer);
506: #ifdef PETSC_HAVE_HDF5
507:     PetscViewerSetType(viewer,PETSCVIEWERHDF5);
508: #else
509:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
510: #endif
511:     PetscViewerSetFromOptions(viewer);
512:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
513:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
514:     if (isHDF5) {
515:       PetscSNPrintf(buf, 256, "uEquilibrium-%dD.h5", ctx->dim);
516:     } else if (isVTK) {
517:       PetscSNPrintf(buf, 256, "uEquilibrium-%dD.vtu", ctx->dim);
518:       PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
519:     }
520:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
521:     PetscViewerFileSetName(viewer,buf);
522:     if (isHDF5) {DMView(dmAux,viewer);}
523:     /* view equilibrium fields, this will overwrite fine grids with coarse grids! */
524:     PetscObjectSetName((PetscObject) global, "u0");
525:     VecView(global,viewer);
526:     PetscViewerDestroy(&viewer);
527:     VecDestroy(&global);
528:   }
529:   VecDestroy(&eq);
530:   return(0);
531: }

533: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
534: {
535:   DM             dmAux, coordDM;
536:   PetscInt       f;

540:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
541:   DMGetCoordinateDM(dm, &coordDM);
542:   if (!feAux) return(0);
543:   DMClone(dm, &dmAux);
544:   PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
545:   DMSetCoordinateDM(dmAux, coordDM);
546:   for (f = 0; f < NfAux; ++f) {DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);}
547:   DMCreateDS(dmAux);
548:   SetupEquilibriumFields(dm, dmAux, user);
549:   DMDestroy(&dmAux);
550:   return(0);
551: }

553: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
554: {
555:   DM              cdm = dm;
556:   const PetscInt  dim = ctx->dim;
557:   PetscQuadrature q;
558:   PetscFE         fe[5], feAux[3];
559:   PetscInt        Nf = 5, NfAux = 3, f;
560:   PetscBool       cell_simplex = ctx->cell_simplex;
561:   MPI_Comm        comm;
562:   PetscErrorCode  ierr;

565:   /* Create finite element */
566:   PetscObjectGetComm((PetscObject) dm, &comm);
567:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[0]);
568:   PetscObjectSetName((PetscObject) fe[0], "density");
569:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[1]);
570:   PetscObjectSetName((PetscObject) fe[1], "vorticity");
571:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[2]);
572:   PetscObjectSetName((PetscObject) fe[2], "flux");
573:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[3]);
574:   PetscObjectSetName((PetscObject) fe[3], "potential");
575:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &fe[4]);
576:   PetscObjectSetName((PetscObject) fe[4], "current");

578:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[0]);
579:   PetscFEGetQuadrature(fe[0], &q);
580:   PetscFESetQuadrature(feAux[0], q);
581:   PetscObjectSetName((PetscObject) feAux[0], "n_0");
582:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[1]);
583:   PetscFEGetQuadrature(fe[1], &q);
584:   PetscFESetQuadrature(feAux[1], q);
585:   PetscObjectSetName((PetscObject) feAux[1], "vorticity_0");
586:   PetscFECreateDefault(comm, dim, 1, cell_simplex, NULL, -1, &feAux[2]);
587:   PetscFEGetQuadrature(fe[2], &q);
588:   PetscFESetQuadrature(feAux[2], q);
589:   PetscObjectSetName((PetscObject) feAux[2], "flux_0");
590:   /* Set discretization and boundary conditions for each mesh */
591:   for (f = 0; f < Nf; ++f) {DMSetField(dm, f, NULL, (PetscObject) fe[f]);}
592:   DMCreateDS(dm);
593:   SetupProblem(dm, ctx);
594:   while (cdm) {
595:     DMCopyDisc(dm, cdm);
596:     SetupAuxDM(dm, NfAux, feAux, ctx);
597:     {
598:       PetscBool hasLabel;

600:       DMHasLabel(cdm, "marker", &hasLabel);
601:       if (!hasLabel) {CreateBCLabel(cdm, "marker");}
602:     }
603:     DMGetCoarseDM(cdm, &cdm);
604:   }
605:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&fe[f]);}
606:   for (f = 0; f < NfAux; ++f) {PetscFEDestroy(&feAux[f]);}
607:   return(0);
608: }

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

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

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

702: /*TEST

704:   test:
705:     suffix: 0
706:     args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
707:   test:
708:     suffix: 1
709:     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

711: TEST*/