Actual source code: ex48.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  1: static char help[] = "Evolution of magnetic islands.\n\
  2: The aim of this model is to self-consistently study the interaction between the tearing mode and small scale drift-wave turbulence.\n\n\n";
This is a three field model for the density $\tilde n$, vorticity $\tilde\Omega$, and magnetic flux $\tilde\psi$, using auxiliary variables potential $\tilde\phi$ and current $j_z$.
\begin{equation}
\begin{aligned}
\partial_t \tilde n &= \left\{ \tilde n, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \left\{ \ln n_0, \tilde\phi \right\} + \mu \nabla^2_\perp \tilde n \\
\partial_t \tilde\Omega &= \left\{ \tilde\Omega, \tilde\phi \right\} + \beta \left\{ j_z, \tilde\psi \right\} + \mu \nabla^2_\perp \tilde\Omega \\
\partial_t \tilde\psi &= \left\{ \psi_0 + \tilde\psi, \tilde\phi - \tilde n \right\} - \left\{ \ln n_0, \tilde\psi \right\} + \frac{\eta}{\beta} \nabla^2_\perp \tilde\psi \\
\nabla^2_\perp\tilde\phi &= \tilde\Omega \\
j_z &= -\nabla^2_\perp \left(\tilde\psi + \psi_0 \right)\\
\end{aligned}
\end{equation}
 17: #include <petscdmplex.h>
 18: #include <petscts.h>
 19: #include <petscds.h>
 20: #include <assert.h>

 22: typedef struct {
 23:   PetscInt       debug;             /* The debugging level */
 24:   PetscBool      plotRef;           /* Plot the reference fields */
 25:   /* Domain and mesh definition */
 26:   PetscInt       dim;               /* The topological mesh dimension */
 27:   char           filename[2048];    /* The optional ExodusII file */
 28:   PetscBool      cell_simplex;           /* Simplicial mesh */
 29:   DMBoundaryType boundary_types[3];
 30:   PetscInt       cells[3];
 31:   PetscInt       refine;
 32:   /* geometry  */
 33:   PetscReal      domain_lo[3], domain_hi[3];
 34:   DMBoundaryType periodicity[3];              /* The domain periodicity */
 35:   PetscReal      b0[3]; /* not used */
 36:   /* Problem definition */
 37:   PetscErrorCode (**initialFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 38:   PetscReal      mu, eta, beta;
 39:   PetscReal      a,b,Jo,Jop,m,ke,kx,ky,DeltaPrime,eps;
 40:   /* solver */
 41:   PetscBool      implicit;
 42: } AppCtx;

 44: static AppCtx *s_ctx;

 46: static PetscScalar poissonBracket(PetscInt dim, const PetscScalar df[], const PetscScalar dg[])
 47: {
 48:   PetscScalar ret = df[0]*dg[1] - df[1]*dg[0];
 49:   return ret;
 50: }

 52: enum field_idx {DENSITY,OMEGA,PSI,PHI,JZ};

 54: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 55:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 56:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 57:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 58: {
 59:   const PetscScalar *pnDer   = &u_x[uOff_x[DENSITY]];
 60:   const PetscScalar *ppsiDer = &u_x[uOff_x[PSI]];
 61:   const PetscScalar *pphiDer = &u_x[uOff_x[PHI]];
 62:   const PetscScalar *jzDer   = &u_x[uOff_x[JZ]];
 63:   const PetscScalar *logRefDenDer = &a_x[aOff_x[DENSITY]];
 64:   f0[0] += - poissonBracket(dim,pnDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer) - poissonBracket(dim,logRefDenDer, pphiDer);
 65:   if (u_t) f0[0] += u_t[DENSITY];
 66: }

 68: static void f1_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 69:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 70:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 71:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 72: {
 73:   const PetscScalar *pnDer = &u_x[uOff_x[DENSITY]];
 74:   PetscInt           d;

 76:   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pnDer[d];
 77: }

 79: static void f0_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 80:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 81:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 82:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 83: {
 84:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 85:   const PetscScalar *ppsiDer   = &u_x[uOff_x[PSI]];
 86:   const PetscScalar *pphiDer   = &u_x[uOff_x[PHI]];
 87:   const PetscScalar *jzDer     = &u_x[uOff_x[JZ]];

 89:   f0[0] += - poissonBracket(dim,pOmegaDer, pphiDer) - s_ctx->beta*poissonBracket(dim,jzDer, ppsiDer);
 90:   if (u_t) f0[0] += u_t[OMEGA];
 91: }

 93: static void f1_Omega(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 94:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 95:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 96:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 97: {
 98:   const PetscScalar *pOmegaDer = &u_x[uOff_x[OMEGA]];
 99:   PetscInt           d;

101:   for (d = 0; d < dim-1; ++d) f1[d] = -s_ctx->mu*pOmegaDer[d];
102: }

104: static void f0_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
105:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
106:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
107:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
108: {
109:   const PetscScalar *pnDer     = &u_x[uOff_x[DENSITY]];
110:   const PetscScalar *ppsiDer   = &u_x[uOff_x[PSI]];
111:   const PetscScalar *pphiDer   = &u_x[uOff_x[PHI]];
112:   const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]];
113:   const PetscScalar *logRefDenDer= &a_x[aOff_x[DENSITY]];
114:   PetscScalar       psiDer[3];
115:   PetscScalar       phi_n_Der[3];
116:   PetscInt          d;
117:   if (dim < 2) {MPI_Abort(MPI_COMM_WORLD,1); return;} /* this is needed so that the clang static analyzer does not generate a warning about variables used by not set */
118:   for (d = 0; d < dim; ++d) {
119:     psiDer[d]    = refPsiDer[d] + ppsiDer[d];
120:     phi_n_Der[d] = pphiDer[d]   - pnDer[d];
121:   }
122:   f0[0] = - poissonBracket(dim,psiDer, phi_n_Der) + poissonBracket(dim,logRefDenDer, ppsiDer);
123:   if (u_t) f0[0] += u_t[PSI];
124: }

126: static void f1_psi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
127:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
128:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
129:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
130: {
131:   const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
132:   PetscInt           d;

134:   for (d = 0; d < dim-1; ++d) f1[d] = -(s_ctx->eta/s_ctx->beta)*ppsi[d];
135: }

137: static void f0_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
138:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
139:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
140:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
141: {
142:   f0[0] = -u[uOff[OMEGA]];
143: }

145: static void f1_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux,
146:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
147:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
148:                    PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
149: {
150:   const PetscScalar *pphi = &u_x[uOff_x[PHI]];
151:   PetscInt           d;

153:   for (d = 0; d < dim-1; ++d) f1[d] = pphi[d];
154: }

156: static void f0_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
157:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
158:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
159:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
160: {
161:   f0[0] = u[uOff[JZ]];
162: }

164: static void f1_jz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167:                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
168: {
169:   const PetscScalar *ppsi = &u_x[uOff_x[PSI]];
170:   const PetscScalar *refPsiDer = &a_x[aOff_x[PSI]]; /* aOff_x[PSI] == 2*PSI */
171:   PetscInt           d;

173:   for (d = 0; d < dim-1; ++d) f1[d] = ppsi[d] + refPsiDer[d];
174: }

176: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
177: {
178:   PetscBool      flg;
180:   PetscInt       ii, bd;
182:   options->debug               = 1;
183:   options->plotRef             = PETSC_FALSE;
184:   options->dim                 = 2;
185:   options->filename[0]         = '\0';
186:   options->cell_simplex        = PETSC_FALSE;
187:   options->implicit            = PETSC_FALSE;
188:   options->refine              = 2;
189:   options->domain_lo[0]  = 0.0;
190:   options->domain_lo[1]  = 0.0;
191:   options->domain_lo[2]  = 0.0;
192:   options->domain_hi[0]  = 2.0;
193:   options->domain_hi[1]  = 2.0*PETSC_PI;
194:   options->domain_hi[2]  = 2.0;
195:   options->periodicity[0]    = DM_BOUNDARY_NONE;
196:   options->periodicity[1]    = DM_BOUNDARY_NONE;
197:   options->periodicity[2]    = DM_BOUNDARY_NONE;
198:   options->mu   = 0;
199:   options->eta  = 0;
200:   options->beta = 1;
201:   options->a = 1;
202:   options->b = PETSC_PI;
203:   options->Jop = 0;
204:   options->m = 1;
205:   options->eps = 1.e-6;

207:   for (ii = 0; ii < options->dim; ++ii) options->cells[ii] = 4;
208:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
209:   PetscOptionsInt("-debug", "The debugging level", "ex48.c", options->debug, &options->debug, NULL);
210:   PetscOptionsBool("-plot_ref", "Plot the reference fields", "ex48.c", options->plotRef, &options->plotRef, NULL);
211:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex48.c", options->dim, &options->dim, NULL);
212:   if (options->dim < 2 || options->dim > 3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Dim %D must be 2 or 3",options->dim);
213:   PetscOptionsInt("-dm_refine", "Hack to get refinement level for cylinder", "ex48.c", options->refine, &options->refine, NULL);
214:   PetscOptionsReal("-mu", "mu", "ex48.c", options->mu, &options->mu, NULL);
215:   PetscOptionsReal("-eta", "eta", "ex48.c", options->eta, &options->eta, NULL);
216:   PetscOptionsReal("-beta", "beta", "ex48.c", options->beta, &options->beta, NULL);
217:   PetscOptionsReal("-Jop", "Jop", "ex48.c", options->Jop, &options->Jop, NULL);
218:   PetscOptionsReal("-m", "m", "ex48.c", options->m, &options->m, NULL);
219:   PetscOptionsReal("-eps", "eps", "ex48.c", options->eps, &options->eps, NULL);
220:   PetscOptionsString("-f", "Exodus.II filename to read", "ex48.c", options->filename, options->filename, sizeof(options->filename), &flg);
221:   PetscOptionsBool("-cell_simplex", "Simplicial (true) or tensor (false) mesh", "ex48.c", options->cell_simplex, &options->cell_simplex, NULL);
222:   PetscOptionsBool("-implicit", "Use implicit time integrator", "ex48.c", options->implicit, &options->implicit, NULL);
223:   ii = options->dim;
224:   PetscOptionsRealArray("-domain_hi", "Domain size", "ex48.c", options->domain_hi, &ii, NULL);
225:   ii = options->dim;
226:   PetscOptionsRealArray("-domain_lo", "Domain size", "ex48.c", options->domain_lo, &ii, NULL);
227:   ii = options->dim;
228:   bd = options->periodicity[0];
229:   PetscOptionsEList("-x_periodicity", "The x-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[0]], &bd, NULL);
230:   options->periodicity[0] = (DMBoundaryType) bd;
231:   bd = options->periodicity[1];
232:   PetscOptionsEList("-y_periodicity", "The y-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[1]], &bd, NULL);
233:   options->periodicity[1] = (DMBoundaryType) bd;
234:   bd = options->periodicity[2];
235:   PetscOptionsEList("-z_periodicity", "The z-boundary periodicity", "ex48.c", DMBoundaryTypes, 5, DMBoundaryTypes[options->periodicity[2]], &bd, NULL);
236:   options->periodicity[2] = (DMBoundaryType) bd;
237:   ii = options->dim;
238:   PetscOptionsIntArray("-cells", "Number of cells in each dimension", "ex48.c", options->cells, &ii, NULL);
239:   PetscOptionsEnd();
240:   options->a = (options->domain_hi[0]-options->domain_lo[0])/2.0;
241:   options->b = (options->domain_hi[1]-options->domain_lo[1])/2.0;
242:   for (ii = 0; ii < options->dim; ++ii) {
243:     if (options->domain_hi[ii] <= options->domain_lo[ii]) SETERRQ3(comm,PETSC_ERR_ARG_WRONG,"Domain %D lo=%g hi=%g",ii,options->domain_lo[ii],options->domain_hi[ii]);
244:   }
245:   options->ke = PetscSqrtScalar(options->Jop);
246:   if (options->Jop==0.0) {
247:     options->Jo = 1.0/PetscPowScalar(options->a,2);
248:   } else {
249:     options->Jo = options->Jop*PetscCosReal(options->ke*options->a)/(1.0-PetscCosReal(options->ke*options->a));
250:   }
251:   options->ky = PETSC_PI*options->m/options->b;
252:   if (PetscPowReal(options->ky, 2) < options->Jop) {
253:     options->kx = PetscSqrtScalar(options->Jop-PetscPowScalar(options->ky,2));
254:     options->DeltaPrime = -2.0*options->kx*options->a*PetscCosReal(options->kx*options->a)/PetscSinReal(options->kx*options->a);
255:   } else if (PetscPowReal(options->ky, 2) > options->Jop) {
256:     options->kx = PetscSqrtScalar(PetscPowScalar(options->ky,2)-options->Jop);
257:     options->DeltaPrime = -2.0*options->kx*options->a*PetscCoshReal(options->kx*options->a)/PetscSinhReal(options->kx*options->a);
258:   } else { /*they're equal (or there's a NaN), lim(x*cot(x))_x->0=1*/
259:     options->kx = 0;
260:     options->DeltaPrime = -2.0;
261:   }
262:   PetscPrintf(comm, "DeltaPrime=%g\n",options->DeltaPrime);

264:   return(0);
265: }

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

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

309: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
310: {
311:   DM             plex;
312:   DMLabel        label;

316:   DMCreateLabel(dm, name);
317:   DMGetLabel(dm, name, &label);
318:   DMConvert(dm, DMPLEX, &plex);
319:   DMPlexMarkBoundaryFaces(dm, 1, label);
320:   DMDestroy(&plex);
321:   return(0);
322: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

705: /*TEST

707:   test:
708:     suffix: 0
709:     args: -debug 1 -dim 2 -dm_refine 1 -x_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0
710:   test:
711:     suffix: 1
712:     args: -debug 1 -dim 3 -dm_refine 1 -z_periodicity PERIODIC -ts_max_steps 1 -ts_max_time 10. -ts_dt 1.0 -domain_lo -2,-1,-1 -domain_hi 2,1,1

714: TEST*/