static char help[] = "Nonlinear Reaction Problem from Chemistry.\n";
/*F
This directory contains examples based on the PDES/ODES given in the book
Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
W. Hundsdorf and J.G. Verwer
Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
\begin{eqnarray}
{U_1}_t - k U_1 U_2 & = & 0 \\
{U_2}_t - k U_1 U_2 & = & 0 \\
{U_3}_t - k U_1 U_2 & = & 0
\end{eqnarray}
Helpful runtime monitoring options:
-ts_view - prints information about the solver being used
-ts_monitor - prints the progress of the solver
-ts_adapt_monitor - prints the progress of the time-step adaptor
-ts_monitor_lg_timestep - plots the size of each timestep (at each time-step)
-ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep)
-ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep)
-draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process
-ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process)
-ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process)
-ts_monitor_lg_error -1 - plots each component of the error in the solution as a function of time (at the end of the solution process)
-lg_use_markers false - do NOT show the data points on the plots
-draw_save - save the timestep and solution plot as a .Gif image file
F*/
/*
Project: Generate a nicely formatted HTML page using
1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html
2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and
3) the text output (output.txt) generated by running the following commands.
4)
rm -rf *.Gif
./ex1 -ts_monitor_lg_error -1 -ts_monitor_lg_solution -1 -draw_pause -2 -ts_max_steps 100 -ts_monitor_lg_timestep -1 -draw_save -draw_virtual -ts_monitor -ts_adapt_monitor -ts_view > output.txt
For example something like
PETSc Example -- Nonlinear Reaction Problem from Chemistry
*/
/*
Include "petscts.h" so that we can use TS solvers. Note that this
file automatically includes:
petscsys.h - base PETSc routines petscvec.h - vectors
petscmat.h - matrices
petscis.h - index sets petscksp.h - Krylov subspace methods
petscviewer.h - viewers petscpc.h - preconditioners
petscksp.h - linear solvers
*/
#include
typedef struct {
PetscScalar k;
Vec initialsolution;
} AppCtx;
PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v)
{
PetscFunctionBegin;
PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
PetscFunctionReturn(PETSC_SUCCESS);
}
PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v)
{
PetscFunctionBegin;
PetscCall(PetscNew(ctx));
PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR));
PetscFunctionReturn(PETSC_SUCCESS);
}
/*
Defines the ODE passed to the ODE solver
*/
PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
{
PetscScalar *f;
const PetscScalar *u, *udot;
PetscFunctionBegin;
/* The next three lines allow us to access the entries of the vectors directly */
PetscCall(VecGetArrayRead(U, &u));
PetscCall(VecGetArrayRead(Udot, &udot));
PetscCall(VecGetArrayWrite(F, &f));
f[0] = udot[0] + ctx->k * u[0] * u[1];
f[1] = udot[1] + ctx->k * u[0] * u[1];
f[2] = udot[2] - ctx->k * u[0] * u[1];
PetscCall(VecRestoreArrayRead(U, &u));
PetscCall(VecRestoreArrayRead(Udot, &udot));
PetscCall(VecRestoreArrayWrite(F, &f));
PetscFunctionReturn(PETSC_SUCCESS);
}
/*
Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
*/
PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
{
PetscInt rowcol[] = {0, 1, 2};
PetscScalar J[3][3];
const PetscScalar *u, *udot;
PetscFunctionBegin;
PetscCall(VecGetArrayRead(U, &u));
PetscCall(VecGetArrayRead(Udot, &udot));
J[0][0] = a + ctx->k * u[1];
J[0][1] = ctx->k * u[0];
J[0][2] = 0.0;
J[1][0] = ctx->k * u[1];
J[1][1] = a + ctx->k * u[0];
J[1][2] = 0.0;
J[2][0] = -ctx->k * u[1];
J[2][1] = -ctx->k * u[0];
J[2][2] = a;
PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
PetscCall(VecRestoreArrayRead(U, &u));
PetscCall(VecRestoreArrayRead(Udot, &udot));
PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
if (A != B) {
PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
}
PetscFunctionReturn(PETSC_SUCCESS);
}
/*
Defines the exact (analytic) solution to the ODE
*/
static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
{
const PetscScalar *uinit;
PetscScalar *u, d0, q;
PetscFunctionBegin;
PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit));
PetscCall(VecGetArrayWrite(U, &u));
d0 = uinit[0] - uinit[1];
if (d0 == 0.0) q = ctx->k * t;
else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0;
u[0] = uinit[0] / (1.0 + uinit[1] * q);
u[1] = u[0] - d0;
u[2] = uinit[1] + uinit[2] - u[1];
PetscCall(VecRestoreArrayWrite(U, &u));
PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit));
PetscFunctionReturn(PETSC_SUCCESS);
}
int main(int argc, char **argv)
{
TS ts; /* ODE integrator */
Vec U; /* solution will be stored here */
Mat A; /* Jacobian matrix */
PetscMPIInt size;
PetscInt n = 3;
AppCtx ctx;
PetscScalar *u;
const char *const names[] = {"U1", "U2", "U3", NULL};
PetscBool mf = PETSC_FALSE;
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Initialize program
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &argv, NULL, help));
PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &mf, NULL));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create necessary matrix and vectors
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
PetscCall(MatSetFromOptions(A));
PetscCall(MatSetUp(A));
PetscCall(MatCreateVecs(A, &U, NULL));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ctx.k = .9;
PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL));
PetscCall(VecDuplicate(U, &ctx.initialsolution));
PetscCall(VecGetArrayWrite(ctx.initialsolution, &u));
u[0] = 1;
u[1] = .7;
u[2] = 0;
PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u));
PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create timestepping solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
PetscCall(TSSetType(ts, TSROSW));
PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx));
if (!mf) PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));
else PetscCall(TSSetIJacobian(ts, NULL, NULL, (TSIJacobianFn *)IJacobian, &ctx));
PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx));
{
DM dm;
void *ptr;
PetscCall(TSGetDM(ts, &dm));
PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr));
PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr));
PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set initial conditions
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(Solution(ts, 0, U, &ctx));
PetscCall(TSSetSolution(ts, U));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set solver options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(TSSetTimeStep(ts, .001));
PetscCall(TSSetMaxSteps(ts, 1000));
PetscCall(TSSetMaxTime(ts, 20.0));
PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
PetscCall(TSSetFromOptions(ts));
PetscCall(TSMonitorLGSetVariableNames(ts, names));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(TSSolve(ts, U));
PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free work space. All PETSc objects should be destroyed when they are no longer needed.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscCall(VecDestroy(&ctx.initialsolution));
PetscCall(MatDestroy(&A));
PetscCall(VecDestroy(&U));
PetscCall(TSDestroy(&ts));
PetscCall(PetscFinalize());
return 0;
}
/*TEST
test:
args: -ts_view
requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
test:
suffix: 2
args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view
requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
output_file: output/ex1_1.out
test:
requires: !single
suffix: 3
args: -ts_view -snes_mf_operator
requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
TEST*/