1 static char help[] = "Test conservation properties for 2-variable system\n\n";
2
3 /*F
4 We consider a linear reaction system with two concentrations
5
6 \begin{align}
7 \frac{\partial c_0}{\partial t} &= -c_0 \\
8 \frac{\partial c_1}{\partial t} &= c_0,
9 \end{align}
10
11 wherethe sum $c_0 + c_1$ is conserved, as can be seen by adding the two equations.
12
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,
17
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}
22
23 We'll consider this three ways, each using an IFunction
24
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
28
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).
31
32 F*/
33
34 #include <petscts.h>
35
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};
42
IFunction_Conservative(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)43 static PetscErrorCode IFunction_Conservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
44 {
45 const PetscScalar *u, *udot;
46 PetscScalar *f;
47
48 PetscFunctionBeginUser;
49 PetscCall(VecGetArrayRead(U, &u));
50 PetscCall(VecGetArrayRead(Udot, &udot));
51 PetscCall(VecGetArray(F, &f));
52
53 f[0] = udot[0] + u[0];
54 f[1] = udot[1] - u[0];
55
56 PetscCall(VecRestoreArrayRead(U, &u));
57 PetscCall(VecRestoreArrayRead(Udot, &udot));
58 PetscCall(VecRestoreArray(F, &f));
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
IFunction_Nonconservative(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)62 static PetscErrorCode IFunction_Nonconservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
63 {
64 const PetscScalar *u, *udot;
65 PetscScalar *f;
66
67 PetscFunctionBeginUser;
68 PetscCall(VecGetArrayRead(U, &u));
69 PetscCall(VecGetArrayRead(Udot, &udot));
70 PetscCall(VecGetArray(F, &f));
71
72 f[0] = PetscExpScalar(u[0]) * udot[0] + PetscExpScalar(u[0]);
73 f[1] = PetscExpScalar(u[1]) * udot[1] - PetscExpScalar(u[0]);
74
75 PetscCall(VecRestoreArrayRead(U, &u));
76 PetscCall(VecRestoreArrayRead(Udot, &udot));
77 PetscCall(VecRestoreArray(F, &f));
78 PetscFunctionReturn(PETSC_SUCCESS);
79 }
80
IFunction_TransientVar(TS ts,PetscReal t,Vec U,Vec Cdot,Vec F,PetscCtx ctx)81 static PetscErrorCode IFunction_TransientVar(TS ts, PetscReal t, Vec U, Vec Cdot, Vec F, PetscCtx ctx)
82 {
83 const PetscScalar *u, *cdot;
84 PetscScalar *f;
85
86 PetscFunctionBeginUser;
87 PetscCall(VecGetArrayRead(U, &u));
88 PetscCall(VecGetArrayRead(Cdot, &cdot));
89 PetscCall(VecGetArray(F, &f));
90
91 f[0] = cdot[0] + PetscExpScalar(u[0]);
92 f[1] = cdot[1] - PetscExpScalar(u[0]);
93
94 PetscCall(VecRestoreArrayRead(U, &u));
95 PetscCall(VecRestoreArrayRead(Cdot, &cdot));
96 PetscCall(VecRestoreArray(F, &f));
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
TransientVar(TS ts,Vec U,Vec C,PetscCtx ctx)100 static PetscErrorCode TransientVar(TS ts, Vec U, Vec C, PetscCtx ctx)
101 {
102 PetscFunctionBeginUser;
103 PetscCall(VecCopy(U, C));
104 PetscCall(VecExp(C));
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
main(int argc,char * argv[])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;
115
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();
121
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));
143
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)));
156
157 PetscCall(VecDestroy(&U));
158 PetscCall(TSDestroy(&ts));
159 PetscCall(PetscFinalize());
160 return 0;
161 }
162
163 /*TEST
164
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
174
175 TEST*/
176