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;
48: PetscFunctionBeginUser;
49: PetscCall(VecGetArrayRead(U, &u));
50: PetscCall(VecGetArrayRead(Udot, &udot));
51: PetscCall(VecGetArray(F, &f));
53: f[0] = udot[0] + u[0];
54: f[1] = udot[1] - u[0];
56: PetscCall(VecRestoreArrayRead(U, &u));
57: PetscCall(VecRestoreArrayRead(Udot, &udot));
58: PetscCall(VecRestoreArray(F, &f));
59: PetscFunctionReturn(PETSC_SUCCESS);
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;
67: PetscFunctionBeginUser;
68: PetscCall(VecGetArrayRead(U, &u));
69: PetscCall(VecGetArrayRead(Udot, &udot));
70: PetscCall(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: PetscCall(VecRestoreArrayRead(U, &u));
76: PetscCall(VecRestoreArrayRead(Udot, &udot));
77: PetscCall(VecRestoreArray(F, &f));
78: PetscFunctionReturn(PETSC_SUCCESS);
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;
86: PetscFunctionBeginUser;
87: PetscCall(VecGetArrayRead(U, &u));
88: PetscCall(VecGetArrayRead(Cdot, &cdot));
89: PetscCall(VecGetArray(F, &f));
91: f[0] = cdot[0] + PetscExpScalar(u[0]);
92: f[1] = cdot[1] - PetscExpScalar(u[0]);
94: PetscCall(VecRestoreArrayRead(U, &u));
95: PetscCall(VecRestoreArrayRead(Cdot, &cdot));
96: PetscCall(VecRestoreArray(F, &f));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode TransientVar(TS ts, Vec U, Vec C, void *ctx)
101: {
102: PetscFunctionBeginUser;
103: PetscCall(VecCopy(U, C));
104: PetscCall(VecExp(C));
105: PetscFunctionReturn(PETSC_SUCCESS);
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;
116: PetscFunctionBeginUser;
117: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
118: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "TS conservation example", "");
119: PetscCall(PetscOptionsEnum("-var", "Variable formulation", NULL, VarModes, (PetscEnum)var, (PetscEnum *)&var, NULL));
120: PetscOptionsEnd();
122: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
123: PetscCall(TSSetType(ts, TSBDF));
124: PetscCall(TSGetDM(ts, &dm));
125: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &U));
126: PetscCall(VecSetValue(U, 0, 2., INSERT_VALUES));
127: PetscCall(VecSetValue(U, 1, 1., INSERT_VALUES));
128: switch (var) {
129: case VAR_CONSERVATIVE:
130: PetscCall(DMTSSetIFunction(dm, IFunction_Conservative, NULL));
131: break;
132: case VAR_NONCONSERVATIVE:
133: PetscCall(VecLog(U));
134: PetscCall(DMTSSetIFunction(dm, IFunction_Nonconservative, NULL));
135: break;
136: case VAR_TRANSIENTVAR:
137: PetscCall(VecLog(U));
138: PetscCall(DMTSSetIFunction(dm, IFunction_TransientVar, NULL));
139: PetscCall(DMTSSetTransientVariable(dm, TransientVar, NULL));
140: }
141: PetscCall(TSSetMaxTime(ts, 1.));
142: PetscCall(TSSetFromOptions(ts));
144: PetscCall(TSSolve(ts, U));
145: switch (var) {
146: case VAR_CONSERVATIVE:
147: break;
148: case VAR_NONCONSERVATIVE:
149: case VAR_TRANSIENTVAR:
150: PetscCall(VecExp(U));
151: break;
152: }
153: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
154: PetscCall(VecSum(U, &sum));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Conservation error %g\n", (double)PetscRealPart(sum - 3.)));
157: PetscCall(VecDestroy(&U));
158: PetscCall(TSDestroy(&ts));
159: PetscCall(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*/