Actual source code: ex15.c

  1: static char help[] = "Test conservation properties for 2-variable system\n\n";

  3: /*F
  4: We consider a linear reaction system with two concentrations

  6: \begin{align}
  7:   \frac{\partial c_0}{\partial t} &= -c_0 \\
  8:   \frac{\partial c_1}{\partial t} &= c_0,
  9: \end{align}

 11: wherethe sum $c_0 + c_1$ is conserved, as can be seen by adding the two equations.

 13: We now consider a different set of variables, defined implicitly by $c(u) = e^u$.  This type of transformation is
 14: sometimes used to ensure positivity, and related transformations are sometimes used to develop a well-conditioned
 15: formulation in limits such as zero Mach number.  In this instance, the relation is explicitly invertible, but that is
 16: not always the case.  We can rewrite the differential equation in terms of non-conservative variables u,

 18: \begin{align}
 19:   \frac{\partial c_0}{\partial u_0} \frac{\partial u_0}{\partial t} &= -c_0(u_0) \\
 20:   \frac{\partial c_1}{\partial u_1} \frac{\partial u_1}{\partial t} &= c_0(u_0).
 21: \end{align}

 23: We'll consider this three ways, each using an IFunction

 25: 1. CONSERVATIVE: standard integration in conservative variables: F(C, Cdot) = 0
 26: 2. NONCONSERVATIVE: chain rule formulation entirely in primitive variables: F(U, Udot) = 0
 27: 3. TRANSIENTVAR: Provide function C(U) and solve F(U, Cdot) = 0, where the time integrators handles the transformation

 29: We will see that 1 and 3 are conservative (up to machine precision/solver tolerance, independent of temporal
 30: discretization error) while 2 is not conservative (i.e., scales with temporal discretization error).

 32: F*/

 34: #include <petscts.h>

 36: typedef enum {
 37:   VAR_CONSERVATIVE,
 38:   VAR_NONCONSERVATIVE,
 39:   VAR_TRANSIENTVAR
 40: } VarMode;
 41: static const char *const VarModes[] = {"CONSERVATIVE", "NONCONSERVATIVE", "TRANSIENTVAR", "VarMode", "VAR_", NULL};

 43: static PetscErrorCode IFunction_Conservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
 44: {
 45:   const PetscScalar *u, *udot;
 46:   PetscScalar       *f;

 49:   VecGetArrayRead(U, &u);
 50:   VecGetArrayRead(Udot, &udot);
 51:   VecGetArray(F, &f);

 53:   f[0] = udot[0] + u[0];
 54:   f[1] = udot[1] - u[0];

 56:   VecRestoreArrayRead(U, &u);
 57:   VecRestoreArrayRead(Udot, &udot);
 58:   VecRestoreArray(F, &f);
 59:   return 0;
 60: }

 62: static PetscErrorCode IFunction_Nonconservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
 63: {
 64:   const PetscScalar *u, *udot;
 65:   PetscScalar       *f;

 68:   VecGetArrayRead(U, &u);
 69:   VecGetArrayRead(Udot, &udot);
 70:   VecGetArray(F, &f);

 72:   f[0] = PetscExpScalar(u[0]) * udot[0] + PetscExpScalar(u[0]);
 73:   f[1] = PetscExpScalar(u[1]) * udot[1] - PetscExpScalar(u[0]);

 75:   VecRestoreArrayRead(U, &u);
 76:   VecRestoreArrayRead(Udot, &udot);
 77:   VecRestoreArray(F, &f);
 78:   return 0;
 79: }

 81: static PetscErrorCode IFunction_TransientVar(TS ts, PetscReal t, Vec U, Vec Cdot, Vec F, void *ctx)
 82: {
 83:   const PetscScalar *u, *cdot;
 84:   PetscScalar       *f;

 87:   VecGetArrayRead(U, &u);
 88:   VecGetArrayRead(Cdot, &cdot);
 89:   VecGetArray(F, &f);

 91:   f[0] = cdot[0] + PetscExpScalar(u[0]);
 92:   f[1] = cdot[1] - PetscExpScalar(u[0]);

 94:   VecRestoreArrayRead(U, &u);
 95:   VecRestoreArrayRead(Cdot, &cdot);
 96:   VecRestoreArray(F, &f);
 97:   return 0;
 98: }

100: static PetscErrorCode TransientVar(TS ts, Vec U, Vec C, void *ctx)
101: {
103:   VecCopy(U, C);
104:   VecExp(C);
105:   return 0;
106: }

108: int main(int argc, char *argv[])
109: {
110:   TS          ts;
111:   DM          dm;
112:   Vec         U;
113:   VarMode     var = VAR_CONSERVATIVE;
114:   PetscScalar sum;

117:   PetscInitialize(&argc, &argv, NULL, help);
118:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "TS conservation example", "");
119:   PetscOptionsEnum("-var", "Variable formulation", NULL, VarModes, (PetscEnum)var, (PetscEnum *)&var, NULL);
120:   PetscOptionsEnd();

122:   TSCreate(PETSC_COMM_WORLD, &ts);
123:   TSSetType(ts, TSBDF);
124:   TSGetDM(ts, &dm);
125:   VecCreateSeq(PETSC_COMM_SELF, 2, &U);
126:   VecSetValue(U, 0, 2., INSERT_VALUES);
127:   VecSetValue(U, 1, 1., INSERT_VALUES);
128:   switch (var) {
129:   case VAR_CONSERVATIVE:
130:     DMTSSetIFunction(dm, IFunction_Conservative, NULL);
131:     break;
132:   case VAR_NONCONSERVATIVE:
133:     VecLog(U);
134:     DMTSSetIFunction(dm, IFunction_Nonconservative, NULL);
135:     break;
136:   case VAR_TRANSIENTVAR:
137:     VecLog(U);
138:     DMTSSetIFunction(dm, IFunction_TransientVar, NULL);
139:     DMTSSetTransientVariable(dm, TransientVar, NULL);
140:   }
141:   TSSetMaxTime(ts, 1.);
142:   TSSetFromOptions(ts);

144:   TSSolve(ts, U);
145:   switch (var) {
146:   case VAR_CONSERVATIVE:
147:     break;
148:   case VAR_NONCONSERVATIVE:
149:   case VAR_TRANSIENTVAR:
150:     VecExp(U);
151:     break;
152:   }
153:   VecView(U, PETSC_VIEWER_STDOUT_WORLD);
154:   VecSum(U, &sum);
155:   PetscPrintf(PETSC_COMM_WORLD, "Conservation error %g\n", (double)PetscRealPart(sum - 3.));

157:   VecDestroy(&U);
158:   TSDestroy(&ts);
159:   PetscFinalize();
160:   return 0;
161: }

163: /*TEST

165:   test:
166:     suffix: conservative
167:     args: -snes_fd -var conservative
168:   test:
169:     suffix: nonconservative
170:     args: -snes_fd -var nonconservative
171:   test:
172:     suffix: transientvar
173:     args: -snes_fd -var transientvar

175: TEST*/