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