Actual source code: ex48.c

      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*/