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

#include <petscts.h>
#include <petscdmda.h>

#if defined(PETSC_HAVE_TCHEM)
  #if defined(MAX)
    #undef MAX
  #endif
  #if defined(MIN)
    #undef MIN
  #endif
  #include <TC_params.h>
  #include <TC_interface.h>
#else
  #error TChem is required for this example.  Reconfigure PETSc using --download-tchem.
#endif
/*

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

    Obtain the three files into this directory

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

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

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

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

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

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

        Compute the sensitivities of the solution of the first temperature on the initial conditions
        -ts_adjoint_solve  -ts_time_step 1.e-5 -ts_type cn -ts_adjoint_view_solution draw

        Turn off diffusion
        -diffusion no

        Turn off reactions
        -reactions no

    The solution for component i = 0 is the temperature.

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

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

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

*/
typedef struct _User *User;
struct _User {
  PetscReal pressure;
  int       Nspec;
  int       Nreac;
  PetscReal Tini, dx;
  PetscReal diffus;
  DM        dm;
  PetscBool diffusion, reactions;
  double   *tchemwork;
  double   *Jdense; /* Dense array workspace where Tchem computes the Jacobian */
  PetscInt *rows;
};

static PetscErrorCode MonitorCell(TS, User, PetscInt);
static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
static PetscErrorCode FormRHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
static PetscErrorCode FormInitialSolution(TS, Vec, void *);

#define PetscCallTC(ierr) \
  do { \
    PetscCheck(!ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in TChem library, return code %d", ierr); \
  } while (0)

int main(int argc, char **argv)
{
  TS                ts; /* time integrator */
  TSAdapt           adapt;
  Vec               X; /* solution vector */
  Mat               J; /* Jacobian matrix */
  PetscInt          steps, ncells, xs, xm, i;
  PetscReal         ftime, dt;
  char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp", thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";
  struct _User      user;
  TSConvergedReason reason;
  PetscBool         showsolutions = PETSC_FALSE;
  char            **snames, *names;
  Vec               lambda; /* used with TSAdjoint for sensitivities */

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Chemistry solver options", "");
  PetscCall(PetscOptionsString("-chem", "CHEMKIN input file", "", chemfile, chemfile, sizeof(chemfile), NULL));
  PetscCall(PetscOptionsString("-thermo", "NASA thermo input file", "", thermofile, thermofile, sizeof(thermofile), NULL));
  user.pressure = 1.01325e5; /* Pascal */
  PetscCall(PetscOptionsReal("-pressure", "Pressure of reaction [Pa]", "", user.pressure, &user.pressure, NULL));
  user.Tini = 1550;
  PetscCall(PetscOptionsReal("-Tini", "Initial temperature [K]", "", user.Tini, &user.Tini, NULL));
  user.diffus = 100;
  PetscCall(PetscOptionsReal("-diffus", "Diffusion constant", "", user.diffus, &user.diffus, NULL));
  PetscCall(PetscOptionsBool("-draw_solution", "Plot the solution for each cell", "", showsolutions, &showsolutions, NULL));
  user.diffusion = PETSC_TRUE;
  PetscCall(PetscOptionsBool("-diffusion", "Have diffusion", "", user.diffusion, &user.diffusion, NULL));
  user.reactions = PETSC_TRUE;
  PetscCall(PetscOptionsBool("-reactions", "Have reactions", "", user.reactions, &user.reactions, NULL));
  PetscOptionsEnd();

  PetscCallTC(TC_initChem(chemfile, thermofile, 0, 1.0));
  user.Nspec = TC_getNspec();
  user.Nreac = TC_getNreac();

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

  /* set the names of each field in the DMDA based on the species name */
  PetscCall(PetscMalloc1((user.Nspec + 1) * LENGTHOFSPECNAME, &names));
  PetscCall(PetscStrncpy(names, "Temp", (user.Nspec + 1) * LENGTHOFSPECNAME));
  TC_getSnames(user.Nspec, names + LENGTHOFSPECNAME);
  PetscCall(PetscMalloc1(user.Nspec + 2, &snames));
  for (i = 0; i < user.Nspec + 1; i++) snames[i] = names + i * LENGTHOFSPECNAME;
  snames[user.Nspec + 1] = NULL;
  PetscCall(DMDASetFieldNames(user.dm, (const char *const *)snames));
  PetscCall(PetscFree(snames));
  PetscCall(PetscFree(names));

  PetscCall(DMCreateMatrix(user.dm, &J));
  PetscCall(DMCreateGlobalVector(user.dm, &X));

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

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create timestepping solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
  PetscCall(TSSetDM(ts, user.dm));
  PetscCall(TSSetType(ts, TSARKIMEX));
  PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE));
  PetscCall(TSARKIMEXSetType(ts, TSARKIMEX4));
  PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
  PetscCall(TSSetRHSJacobian(ts, J, J, FormRHSJacobian, &user));

  ftime = 1.0;
  PetscCall(TSSetMaxTime(ts, ftime));
  PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set initial conditions
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(FormInitialSolution(ts, X, &user));
  PetscCall(TSSetSolution(ts, X));
  dt = 1e-10; /* Initial time step */
  PetscCall(TSSetTimeStep(ts, dt));
  PetscCall(TSGetAdapt(ts, &adapt));
  PetscCall(TSAdaptSetStepLimits(adapt, 1e-12, 1e-4)); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
  PetscCall(TSSetMaxSNESFailures(ts, -1));             /* Retry step an unlimited number of times */

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Pass information to graphical monitoring routine
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  if (showsolutions) {
    PetscCall(DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL));
    for (i = xs; i < xs + xm; i++) PetscCall(MonitorCell(ts, &user, i));
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSSetFromOptions(ts));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set final conditions for sensitivities
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(DMCreateGlobalVector(user.dm, &lambda));
  PetscCall(TSSetCostGradients(ts, 1, &lambda, NULL));
  PetscCall(VecSetValue(lambda, 0, 1.0, INSERT_VALUES));
  PetscCall(VecAssemblyBegin(lambda));
  PetscCall(VecAssemblyEnd(lambda));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Solve ODE
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSSolve(ts, X));
  PetscCall(TSGetSolveTime(ts, &ftime));
  PetscCall(TSGetStepNumber(ts, &steps));
  PetscCall(TSGetConvergedReason(ts, &reason));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

  {
    Vec                max;
    const char *const *names;
    PetscInt           i;
    const PetscReal   *bmax;

    PetscCall(TSMonitorEnvelopeGetBounds(ts, &max, NULL));
    if (max) {
      PetscCall(TSMonitorLGGetVariableNames(ts, &names));
      if (names) {
        PetscCall(VecGetArrayRead(max, &bmax));
        PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species - maximum mass fraction\n"));
        for (i = 1; i < user.Nspec; i++) {
          if (bmax[i] > .01) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s %g\n", names[i], (double)bmax[i]));
        }
        PetscCall(VecRestoreArrayRead(max, &bmax));
      }
    }
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Free work space.
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  TC_reset();
  PetscCall(DMDestroy(&user.dm));
  PetscCall(MatDestroy(&J));
  PetscCall(VecDestroy(&X));
  PetscCall(VecDestroy(&lambda));
  PetscCall(TSDestroy(&ts));
  PetscCall(PetscFree3(user.tchemwork, user.Jdense, user.rows));
  PetscCall(PetscFinalize());
  return 0;
}

/*
   Applies the second order centered difference diffusion operator on a one dimensional periodic domain
*/
static PetscErrorCode FormDiffusionFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
{
  User                user = (User)ptr;
  PetscScalar       **f;
  const PetscScalar **x;
  DM                  dm;
  PetscInt            i, xs, xm, j, dof;
  Vec                 Xlocal;
  PetscReal           idx;

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &dm));
  PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
  PetscCall(DMGetLocalVector(dm, &Xlocal));
  PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xlocal));
  PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xlocal));
  PetscCall(DMDAVecGetArrayDOFRead(dm, Xlocal, &x));
  PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
  PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

  idx = 1.0 * user->diffus / user->dx;
  for (i = xs; i < xs + xm; i++) {
    for (j = 0; j < dof; j++) f[i][j] += idx * (x[i + 1][j] - 2.0 * x[i][j] + x[i - 1][j]);
  }
  PetscCall(DMDAVecRestoreArrayDOFRead(dm, Xlocal, &x));
  PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
  PetscCall(DMRestoreLocalVector(dm, &Xlocal));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  PetscCall(TSGetDM(ts, &dm));
  PetscCall(DMDAGetInfo(dm, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
  PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

  idx       = 1.0 * user->diffus / user->dx;
  values[0] = idx;
  values[1] = -2.0 * idx;
  values[2] = idx;
  for (i = xs; i < xs + xm; i++) {
    for (j = 0; j < dof; j++) {
      row.i    = i;
      row.c    = j;
      col[0].i = i - 1;
      col[0].c = j;
      col[1].i = i;
      col[1].c = j;
      col[2].i = i + 1;
      col[2].c = j;
      PetscCall(MatSetValuesStencil(Pmat, 1, &row, 3, col, values, ADD_VALUES));
    }
  }
  PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
{
  User                user = (User)ptr;
  PetscScalar       **f;
  const PetscScalar **x;
  DM                  dm;
  PetscInt            i, xs, xm;

  PetscFunctionBeginUser;
  if (user->reactions) {
    PetscCall(TSGetDM(ts, &dm));
    PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
    PetscCall(DMDAVecGetArrayDOF(dm, F, &f));
    PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

    for (i = xs; i < xs + xm; i++) {
      PetscCall(PetscArraycpy(user->tchemwork, x[i], user->Nspec + 1));
      user->tchemwork[0] *= user->Tini; /* Dimensionalize */
      PetscCallTC(TC_getSrc(user->tchemwork, user->Nspec + 1, f[i]));
      f[i][0] /= user->Tini; /* Non-dimensionalize */
    }

    PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
    PetscCall(DMDAVecRestoreArrayDOF(dm, F, &f));
  } else {
    PetscCall(VecZeroEntries(F));
  }
  if (user->diffusion) PetscCall(FormDiffusionFunction(ts, t, X, F, ptr));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBeginUser;
  if (user->reactions) {
    PetscCall(TSGetDM(ts, &dm));
    PetscCall(MatZeroEntries(Pmat));
    PetscCall(MatSetOption(Pmat, MAT_ROW_ORIENTED, PETSC_FALSE));
    PetscCall(MatSetOption(Pmat, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
    PetscCall(DMDAVecGetArrayDOFRead(dm, X, &x));
    PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

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

      for (j = 0; j < M; j++) user->Jdense[j + 0 * M] /= user->Tini; /* Non-dimensionalize first column */
      for (j = 0; j < M; j++) user->Jdense[0 + j * M] /= user->Tini; /* Non-dimensionalize first row */
      for (j = 0; j < M; j++) user->rows[j] = i * M + j;
      PetscCall(MatSetValues(Pmat, M, user->rows, M, user->rows, user->Jdense, INSERT_VALUES));
    }
    PetscCall(DMDAVecRestoreArrayDOFRead(dm, X, &x));
    PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
  } else {
    PetscCall(MatZeroEntries(Pmat));
  }
  if (user->diffusion) PetscCall(FormDiffusionJacobian(ts, t, X, Amat, Pmat, ptr));
  if (Amat != Pmat) {
    PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
{
  PetscScalar **x, *xc;
  struct {
    const char *name;
    PetscReal   massfrac;
  } initial[] = {
    {"CH4", 0.0948178320887 },
    {"O2",  0.189635664177  },
    {"N2",  0.706766236705  },
    {"AR",  0.00878026702874}
  };
  PetscInt i, j, xs, xm;
  DM       dm;

  PetscFunctionBeginUser;
  PetscCall(VecZeroEntries(X));
  PetscCall(TSGetDM(ts, &dm));
  PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xm, NULL, NULL));

  PetscCall(DMDAGetCoordinateArray(dm, &xc));
  PetscCall(DMDAVecGetArrayDOF(dm, X, &x));
  for (i = xs; i < xs + xm; i++) {
    x[i][0] = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[i]); /* Non-dimensionalized by user->Tini */
    for (j = 0; j < PETSC_STATIC_ARRAY_LENGTH(initial); j++) {
      int ispec = TC_getSpos(initial[j].name, (int)strlen(initial[j].name));
      PetscCheck(ispec >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Could not find species %s", initial[j].name);
      PetscCall(PetscPrintf(PETSC_COMM_SELF, "Species %d: %s %g\n", j, initial[j].name, (double)initial[j].massfrac));
      x[i][1 + ispec] = initial[j].massfrac;
    }
  }
  PetscCall(DMDAVecRestoreArrayDOF(dm, X, &x));
  PetscCall(DMDARestoreCoordinateArray(dm, &xc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
    Routines for displaying the solutions
*/
typedef struct {
  PetscInt cell;
  User     user;
} UserLGCtx;

static PetscErrorCode FormMoleFraction(UserLGCtx *ctx, Vec massf, Vec *molef)
{
  User                user = ctx->user;
  PetscReal          *M, tM = 0;
  PetscInt            i, n  = user->Nspec + 1;
  PetscScalar        *mof;
  const PetscScalar **maf;

  PetscFunctionBeginUser;
  PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, molef));
  PetscCall(PetscMalloc1(user->Nspec, &M));
  TC_getSmass(user->Nspec, M);
  PetscCall(DMDAVecGetArrayDOFRead(user->dm, massf, &maf));
  PetscCall(VecGetArray(*molef, &mof));
  mof[0] = maf[ctx->cell][0]; /* copy over temperature */
  for (i = 1; i < n; i++) tM += maf[ctx->cell][i] / M[i - 1];
  for (i = 1; i < n; i++) mof[i] = maf[ctx->cell][i] / (M[i - 1] * tM);
  PetscCall(DMDAVecRestoreArrayDOFRead(user->dm, massf, &maf));
  PetscCall(VecRestoreArray(*molef, &mof));
  PetscCall(PetscFree(M));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MonitorCellDestroy(UserLGCtx **uctx)
{
  PetscFunctionBeginUser;
  PetscCall(PetscFree(*uctx));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
   Use TSMonitorLG to monitor the reactions in a particular cell
*/
static PetscErrorCode MonitorCell(TS ts, User user, PetscInt cell)
{
  TSMonitorLGCtx ctx;
  char         **snames;
  UserLGCtx     *uctx;
  char           label[128];
  PetscReal      temp, *xc;
  PetscMPIInt    rank;

  PetscFunctionBeginUser;
  PetscCall(DMDAGetCoordinateArray(user->dm, &xc));
  temp = 1.0 + .05 * PetscSinScalar(2. * PETSC_PI * xc[cell]); /* Non-dimensionalized by user->Tini */
  PetscCall(DMDARestoreCoordinateArray(user->dm, &xc));
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
  PetscCall(PetscSNPrintf(label, sizeof(label), "Initial Temperature %g Cell %d Rank %d", (double)user->Tini * temp, (int)cell, rank));
  PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, label, PETSC_DECIDE, PETSC_DECIDE, 600, 400, 1, &ctx));
  PetscCall(DMDAGetFieldNames(user->dm, (const char *const **)&snames));
  PetscCall(TSMonitorLGCtxSetVariableNames(ctx, (const char *const *)snames));
  PetscCall(PetscNew(&uctx));
  uctx->cell = cell;
  uctx->user = user;
  PetscCall(TSMonitorLGCtxSetTransform(ctx, (PetscErrorCode (*)(void *, Vec, Vec *))FormMoleFraction, (PetscCtxDestroyFn *)MonitorCellDestroy, uctx));
  PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
  PetscFunctionReturn(PETSC_SUCCESS);
}