Actual source code: extchemfield.c

      1: static const char help[] = "Integrate chemistry using TChem.\n";

      3: #include <petscts.h>
      4: #include <petscdmda.h>

      6: #if defined(PETSC_HAVE_TCHEM)
      7: #if defined(MAX)
      8: #undef MAX
      9: #endif
     10: #if defined(MIN)
     11: #undef MIN
     12: #endif
     13: #  include <TC_params.h>
     14: #  include <TC_interface.h>
     15: #else
     16: #  error TChem is required for this example.  Reconfigure PETSc using --download-tchem.
     17: #endif
     18: /*

     20:     This is an extension of extchem.c to solve the reaction equations independently in each cell of a one dimensional field

     22:     Obtain the three files into this directory

     24:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/grimech30.dat > chem.inp
     25:        curl http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat > therm.dat
     26:        cp $PETSC_DIR/$PETSC_ARCH/externalpackages/tchem/data/periodictable.dat .

     28:     Run with
     29:      ./extchemfield  -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_monitor -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_max_time .005

     31:      Options for visualizing the solution:
     32:         Watch certain variables in each cell evolve with time
     33:         -draw_solution 1 -ts_monitor_lg_solution_variables Temp,H2,O2,H2O,CH4,CO,CO2,C2H2,N2 -lg_use_markers false  -draw_pause -2

     35:         Watch certain variables in all cells evolve with time
     36:         -da_refine 4 -ts_monitor_draw_solution -draw_fields_by_name Temp,H2 -draw_vec_mark_points  -draw_pause -2

     38:         Keep the initial temperature distribution as one monitors the current temperature distribution
     39:         -ts_monitor_draw_solution_initial -draw_bounds .9,1.7 -draw_fields_by_name Temp

     41:         Save the images in a .gif (movie) file
     42:         -draw_save -draw_save_single_file

     44:         Compute the sensitivies of the solution of the first tempature on the initial conditions
     45:         -ts_adjoint_solve  -ts_dt 1.e-5 -ts_type cn -ts_adjoint_view_solution draw

     47:         Turn off diffusion
     48:         -diffusion no

     50:         Turn off reactions
     51:         -reactions no

     54:     The solution for component i = 0 is the temperature.

     56:     The solution, i > 0, is the mass fraction, massf[i], of species i, i.e. mass of species i/ total mass of all species

     58:     The mole fraction molef[i], i > 0, is the number of moles of a species/ total number of moles of all species
     59:         Define M[i] = mass per mole of species i then
     60:         molef[i] = massf[i]/(M[i]*(sum_j massf[j]/M[j]))

     62:     FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species.

     64: */
     65: typedef struct _User *User;
     66: struct _User {
     67:   PetscReal pressure;
     68:   int       Nspec;
     69:   int       Nreac;
     70:   PetscReal Tini,dx;
     71:   PetscReal diffus;
     72:   DM        dm;
     73:   PetscBool diffusion,reactions;
     74:   double    *tchemwork;
     75:   double    *Jdense;        /* Dense array workspace where Tchem computes the Jacobian */
     76:   PetscInt  *rows;
     77: };

     79: static PetscErrorCode MonitorCell(TS,User,PetscInt);
     80: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
     81: static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
     82: static PetscErrorCode FormInitialSolution(TS,Vec,void*);

     84: #define TCCHKERRQ(ierr) do {if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} while (0)

     86: int main(int argc,char **argv)
     87: {
     88:   TS                ts;         /* time integrator */
     89:   TSAdapt           adapt;
     90:   Vec               X;          /* solution vector */
     91:   Mat               J;          /* Jacobian matrix */
     92:   PetscInt          steps,ncells,xs,xm,i;
     93:   PetscErrorCode    ierr;
     94:   PetscReal         ftime,dt;
     95:   char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp",thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";
     96:   struct _User      user;
     97:   TSConvergedReason reason;
     98:   PetscBool         showsolutions = PETSC_FALSE;
     99:   char              **snames,*names;
    100:   Vec               lambda;     /* used with TSAdjoint for sensitivities */

    102:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
    103:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");
    104:   PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);
    105:   PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);
    106:   user.pressure = 1.01325e5;    /* Pascal */
    107:   PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);
    108:   user.Tini   = 1550;
    109:   PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);
    110:   user.diffus = 100;
    111:   PetscOptionsReal("-diffus","Diffusion constant","",user.diffus,&user.diffus,NULL);
    112:   PetscOptionsBool("-draw_solution","Plot the solution for each cell","",showsolutions,&showsolutions,NULL);
    113:   user.diffusion = PETSC_TRUE;
    114:   PetscOptionsBool("-diffusion","Have diffusion","",user.diffusion,&user.diffusion,NULL);
    115:   user.reactions = PETSC_TRUE;
    116:   PetscOptionsBool("-reactions","Have reactions","",user.reactions,&user.reactions,NULL);
    117:   PetscOptionsEnd();

    119:   TC_initChem(chemfile, thermofile, 0, 1.0);TC
    120:   user.Nspec = TC_getNspec();
    121:   user.Nreac = TC_getNreac();

    123:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,10,user.Nspec+1,1,NULL,&user.dm);
    124:   DMSetFromOptions(user.dm);
    125:   DMSetUp(user.dm);
    126:   DMDAGetInfo(user.dm,NULL,&ncells,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    127:   user.dx = 1.0/ncells;  /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */
    128:   DMDASetUniformCoordinates(user.dm,0.0,1.0,0.0,1.0,0.0,1.0);

    130:   /* set the names of each field in the DMDA based on the species name */
    131:   PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);
    132:   PetscStrcpy(names,"Temp");
    133:   TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);
    134:   PetscMalloc1((user.Nspec+2),&snames);
    135:   for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;
    136:   snames[user.Nspec+1] = NULL;
    137:   DMDASetFieldNames(user.dm,(const char * const *)snames);
    138:   PetscFree(snames);
    139:   PetscFree(names);

    142:   DMCreateMatrix(user.dm,&J);
    143:   DMCreateGlobalVector(user.dm,&X);

    145:   PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);

    147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:      Create timestepping solver context
    149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    150:   TSCreate(PETSC_COMM_WORLD,&ts);
    151:   TSSetDM(ts,user.dm);
    152:   TSSetType(ts,TSARKIMEX);
    153:   TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
    154:   TSARKIMEXSetType(ts,TSARKIMEX4);
    155:   TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
    156:   TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);

    158:   ftime = 1.0;
    159:   TSSetMaxTime(ts,ftime);
    160:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

    162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    163:      Set initial conditions
    164:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    165:   FormInitialSolution(ts,X,&user);
    166:   TSSetSolution(ts,X);
    167:   dt   = 1e-10;                 /* Initial time step */
    168:   TSSetTimeStep(ts,dt);
    169:   TSGetAdapt(ts,&adapt);
    170:   TSAdaptSetStepLimits(adapt,1e-12,1e-4); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
    171:   TSSetMaxSNESFailures(ts,-1);            /* Retry step an unlimited number of times */

    174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175:      Pass information to graphical monitoring routine
    176:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    177:   if (showsolutions) {
    178:     DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL);
    179:     for (i=xs;i<xs+xm;i++) {
    180:       MonitorCell(ts,&user,i);
    181:     }
    182:   }

    184:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    185:      Set runtime options
    186:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    187:   TSSetFromOptions(ts);

    189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    190:      Set final conditions for sensitivities
    191:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    192:   DMCreateGlobalVector(user.dm,&lambda);
    193:   TSSetCostGradients(ts,1,&lambda,NULL);
    194:   VecSetValue(lambda,0,1.0,INSERT_VALUES);
    195:   VecAssemblyBegin(lambda);
    196:   VecAssemblyEnd(lambda);

    198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    199:      Solve ODE
    200:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    201:   TSSolve(ts,X);
    202:   TSGetSolveTime(ts,&ftime);
    203:   TSGetStepNumber(ts,&steps);
    204:   TSGetConvergedReason(ts,&reason);
    205:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);

    207:   {
    208:     Vec                max;
    209:     const char * const *names;
    210:     PetscInt           i;
    211:     const PetscReal    *bmax;

    213:     TSMonitorEnvelopeGetBounds(ts,&max,NULL);
    214:     if (max) {
    215:       TSMonitorLGGetVariableNames(ts,&names);
    216:       if (names) {
    217:         VecGetArrayRead(max,&bmax);
    218:         PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");
    219:         for (i=1; i<user.Nspec; i++) {
    220:           if (bmax[i] > .01) {PetscPrintf(PETSC_COMM_SELF,"%s %g\n",names[i],bmax[i]);}
    221:         }
    222:         VecRestoreArrayRead(max,&bmax);
    223:       }
    224:     }
    225:   }

    227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    228:      Free work space.
    229:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    230:   TC_reset();
    231:   DMDestroy(&user.dm);
    232:   MatDestroy(&J);
    233:   VecDestroy(&X);
    234:   VecDestroy(&lambda);
    235:   TSDestroy(&ts);
    236:   PetscFree3(user.tchemwork,user.Jdense,user.rows);
    237:   PetscFinalize();
    238:   return ierr;
    239: }

    241: /*
    242:    Applies the second order centered difference diffusion operator on a one dimensional periodic domain
    243: */
    244: static PetscErrorCode FormDiffusionFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
    245: {
    246:   User              user = (User)ptr;
    247:   PetscErrorCode    ierr;
    248:   PetscScalar       **f;
    249:   const PetscScalar **x;
    250:   DM                dm;
    251:   PetscInt          i,xs,xm,j,dof;
    252:   Vec               Xlocal;
    253:   PetscReal         idx;

    256:   TSGetDM(ts,&dm);
    257:   DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
    258:   DMGetLocalVector(dm,&Xlocal);
    259:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xlocal);
    260:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xlocal);
    261:   DMDAVecGetArrayDOFRead(dm,Xlocal,&x);
    262:   DMDAVecGetArrayDOF(dm,F,&f);
    263:   DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);

    265:   idx = 1.0*user->diffus/user->dx;
    266:   for (i=xs; i<xs+xm; i++) {
    267:     for (j=0; j<dof; j++) {
    268:       f[i][j] += idx*(x[i+1][j] - 2.0*x[i][j] + x[i-1][j]);
    269:     }
    270:   }
    271:   DMDAVecRestoreArrayDOFRead(dm,Xlocal,&x);
    272:   DMDAVecRestoreArrayDOF(dm,F,&f);
    273:   DMRestoreLocalVector(dm,&Xlocal);
    274:   return(0);
    275: }

    277: /*
    278:    Produces the second order centered difference diffusion operator on a one dimensional periodic domain
    279: */
    280: static PetscErrorCode FormDiffusionJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
    281: {
    282:   User              user = (User)ptr;
    283:   PetscErrorCode    ierr;
    284:   DM                dm;
    285:   PetscInt          i,xs,xm,j,dof;
    286:   PetscReal         idx,values[3];
    287:   MatStencil        row,col[3];

    290:   TSGetDM(ts,&dm);
    291:   DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
    292:   DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);

    294:   idx = 1.0*user->diffus/user->dx;
    295:   values[0] = idx;
    296:   values[1] = -2.0*idx;
    297:   values[2] = idx;
    298:   for (i=xs; i<xs+xm; i++) {
    299:     for (j=0; j<dof; j++) {
    300:       row.i = i;      row.c = j;
    301:       col[0].i = i-1; col[0].c = j;
    302:       col[1].i = i;   col[1].c = j;
    303:       col[2].i = i+1; col[2].c = j;
    304:       MatSetValuesStencil(Pmat,1,&row,3,col,values,ADD_VALUES);
    305:     }
    306:   }
    307:   MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
    308:   MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
    309:   return(0);
    310: }

    312: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
    313: {
    314:   User              user = (User)ptr;
    315:   PetscErrorCode    ierr;
    316:   PetscScalar       **f;
    317:   const PetscScalar **x;
    318:   DM                dm;
    319:   PetscInt          i,xs,xm;

    322:   if (user->reactions) {
    323:     TSGetDM(ts,&dm);
    324:     DMDAVecGetArrayDOFRead(dm,X,&x);
    325:     DMDAVecGetArrayDOF(dm,F,&f);
    326:     DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);

    328:     for (i=xs; i<xs+xm; i++) {
    329:       PetscArraycpy(user->tchemwork,x[i],user->Nspec+1);
    330:       user->tchemwork[0] *= user->Tini; /* Dimensionalize */
    331:       TC_getSrc(user->tchemwork,user->Nspec+1,f[i]);TC
    332:       f[i][0] /= user->Tini;           /* Non-dimensionalize */
    333:     }

    335:     DMDAVecRestoreArrayDOFRead(dm,X,&x);
    336:     DMDAVecRestoreArrayDOF(dm,F,&f);
    337:   } else {
    338:     VecZeroEntries(F);
    339:   }
    340:   if (user->diffusion) {
    341:     FormDiffusionFunction(ts,t,X,F,ptr);
    342:   }
    343:   return(0);
    344: }

    346: static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
    347: {
    348:   User              user = (User)ptr;
    349:   PetscErrorCode    ierr;
    350:   const PetscScalar **x;
    351:   PetscInt          M = user->Nspec+1,i,j,xs,xm;
    352:   DM                dm;

    355:   if (user->reactions) {
    356:     TSGetDM(ts,&dm);
    357:     MatZeroEntries(Pmat);
    358:     MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);
    359:     MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
    360:     DMDAVecGetArrayDOFRead(dm,X,&x);
    361:     DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);

    363:     for (i=xs; i<xs+xm; i++) {
    364:       PetscArraycpy(user->tchemwork,x[i],user->Nspec+1);
    365:       user->tchemwork[0] *= user->Tini;  /* Dimensionalize temperature (first row) because that is what Tchem wants */
    366:       TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);

    368:       for (j=0; j<M; j++) user->Jdense[j + 0*M] /= user->Tini; /* Non-dimensionalize first column */
    369:       for (j=0; j<M; j++) user->Jdense[0 + j*M] /= user->Tini; /* Non-dimensionalize first row */
    370:       for (j=0; j<M; j++) user->rows[j] = i*M+j;
    371:       MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);
    372:     }
    373:     DMDAVecRestoreArrayDOFRead(dm,X,&x);
    374:     MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
    375:     MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
    376:   } else {
    377:     MatZeroEntries(Pmat);
    378:   }
    379:   if (user->diffusion) {
    380:     FormDiffusionJacobian(ts,t,X,Amat,Pmat,ptr);
    381:   }
    382:   if (Amat != Pmat) {
    383:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
    384:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
    385:   }
    386:   return(0);
    387: }

    389: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
    390: {
    391:   PetscScalar    **x,*xc;
    393:   struct {const char *name; PetscReal massfrac;} initial[] = {
    394:     {"CH4", 0.0948178320887},
    395:     {"O2", 0.189635664177},
    396:     {"N2", 0.706766236705},
    397:     {"AR", 0.00878026702874}
    398:   };
    399:   PetscInt       i,j,xs,xm;
    400:   DM             dm;

    403:   VecZeroEntries(X);
    404:   TSGetDM(ts,&dm);
    405:   DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);

    407:   DMDAGetCoordinateArray(dm,&xc);
    408:   DMDAVecGetArrayDOF(dm,X,&x);
    409:   for (i=xs; i<xs+xm; i++) {
    410:     x[i][0] = 1.0 + .05*PetscSinScalar(2.*PETSC_PI*xc[i]);  /* Non-dimensionalized by user->Tini */
    411:     for (j=0; j<sizeof(initial)/sizeof(initial[0]); j++) {
    412:       int ispec = TC_getSpos(initial[j].name, strlen(initial[j].name));
    413:       if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",initial[j].name);
    414:       PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",j,initial[j].name,initial[j].massfrac);
    415:       x[i][1+ispec] = initial[j].massfrac;
    416:     }
    417:   }
    418:   DMDAVecRestoreArrayDOF(dm,X,&x);
    419:   DMDARestoreCoordinateArray(dm,&xc);
    420:   return(0);
    421: }

    423: /*
    424:     Routines for displaying the solutions
    425: */
    426: typedef struct {
    427:   PetscInt cell;
    428:   User     user;
    429: } UserLGCtx;

    431: static PetscErrorCode FormMoleFraction(UserLGCtx *ctx,Vec massf,Vec *molef)
    432: {
    433:   User              user = ctx->user;
    434:   PetscErrorCode    ierr;
    435:   PetscReal         *M,tM=0;
    436:   PetscInt          i,n = user->Nspec+1;
    437:   PetscScalar       *mof;
    438:   const PetscScalar **maf;

    441:   VecCreateSeq(PETSC_COMM_SELF,n,molef);
    442:   PetscMalloc1(user->Nspec,&M);
    443:   TC_getSmass(user->Nspec, M);
    444:   DMDAVecGetArrayDOFRead(user->dm,massf,&maf);
    445:   VecGetArray(*molef,&mof);
    446:   mof[0] = maf[ctx->cell][0]; /* copy over temperature */
    447:   for (i=1; i<n; i++) tM += maf[ctx->cell][i]/M[i-1];
    448:   for (i=1; i<n; i++) {
    449:     mof[i] = maf[ctx->cell][i]/(M[i-1]*tM);
    450:   }
    451:   DMDAVecRestoreArrayDOFRead(user->dm,massf,&maf);
    452:   VecRestoreArray(*molef,&mof);
    453:   PetscFree(M);
    454:   return(0);
    455: }

    457: static PetscErrorCode MonitorCellDestroy(UserLGCtx *uctx)
    458: {

    462:   PetscFree(uctx);
    463:   return(0);
    464: }

    466: /*
    467:    Use TSMonitorLG to monitor the reactions in a particular cell
    468: */
    469: static PetscErrorCode MonitorCell(TS ts,User user,PetscInt cell)
    470: {
    472:   TSMonitorLGCtx ctx;
    473:   char           **snames;
    474:   UserLGCtx      *uctx;
    475:   char           label[128];
    476:   PetscReal      temp,*xc;
    477:   PetscMPIInt    rank;

    480:   DMDAGetCoordinateArray(user->dm,&xc);
    481:   temp = 1.0 + .05*PetscSinScalar(2.*PETSC_PI*xc[cell]);  /* Non-dimensionalized by user->Tini */
    482:   DMDARestoreCoordinateArray(user->dm,&xc);
    483:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    484:   PetscSNPrintf(label,sizeof(label),"Initial Temperature %g Cell %d Rank %d",(double)user->Tini*temp,(int)cell,rank);
    485:   TSMonitorLGCtxCreate(PETSC_COMM_SELF,NULL,label,PETSC_DECIDE,PETSC_DECIDE,600,400,1,&ctx);
    486:   DMDAGetFieldNames(user->dm,(const char * const **)&snames);
    487:   TSMonitorLGCtxSetVariableNames(ctx,(const char * const *)snames);
    488:   PetscNew(&uctx);
    489:   uctx->cell = cell;
    490:   uctx->user = user;
    491:   TSMonitorLGCtxSetTransform(ctx,(PetscErrorCode (*)(void*,Vec,Vec*))FormMoleFraction,(PetscErrorCode (*)(void*))MonitorCellDestroy,uctx);
    492:   TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
    493:   return(0);
    494: }