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 43 static PetscErrorCode IFunction_Conservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) { 44 const PetscScalar *u, *udot; 45 PetscScalar *f; 46 47 PetscFunctionBeginUser; 48 PetscCall(VecGetArrayRead(U, &u)); 49 PetscCall(VecGetArrayRead(Udot, &udot)); 50 PetscCall(VecGetArray(F, &f)); 51 52 f[0] = udot[0] + u[0]; 53 f[1] = udot[1] - u[0]; 54 55 PetscCall(VecRestoreArrayRead(U, &u)); 56 PetscCall(VecRestoreArrayRead(Udot, &udot)); 57 PetscCall(VecRestoreArray(F, &f)); 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode IFunction_Nonconservative(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) { 62 const PetscScalar *u, *udot; 63 PetscScalar *f; 64 65 PetscFunctionBeginUser; 66 PetscCall(VecGetArrayRead(U, &u)); 67 PetscCall(VecGetArrayRead(Udot, &udot)); 68 PetscCall(VecGetArray(F, &f)); 69 70 f[0] = PetscExpScalar(u[0]) * udot[0] + PetscExpScalar(u[0]); 71 f[1] = PetscExpScalar(u[1]) * udot[1] - PetscExpScalar(u[0]); 72 73 PetscCall(VecRestoreArrayRead(U, &u)); 74 PetscCall(VecRestoreArrayRead(Udot, &udot)); 75 PetscCall(VecRestoreArray(F, &f)); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode IFunction_TransientVar(TS ts, PetscReal t, Vec U, Vec Cdot, Vec F, void *ctx) { 80 const PetscScalar *u, *cdot; 81 PetscScalar *f; 82 83 PetscFunctionBeginUser; 84 PetscCall(VecGetArrayRead(U, &u)); 85 PetscCall(VecGetArrayRead(Cdot, &cdot)); 86 PetscCall(VecGetArray(F, &f)); 87 88 f[0] = cdot[0] + PetscExpScalar(u[0]); 89 f[1] = cdot[1] - PetscExpScalar(u[0]); 90 91 PetscCall(VecRestoreArrayRead(U, &u)); 92 PetscCall(VecRestoreArrayRead(Cdot, &cdot)); 93 PetscCall(VecRestoreArray(F, &f)); 94 PetscFunctionReturn(0); 95 } 96 97 static PetscErrorCode TransientVar(TS ts, Vec U, Vec C, void *ctx) { 98 PetscFunctionBeginUser; 99 PetscCall(VecCopy(U, C)); 100 PetscCall(VecExp(C)); 101 PetscFunctionReturn(0); 102 } 103 104 int main(int argc, char *argv[]) { 105 TS ts; 106 DM dm; 107 Vec U; 108 VarMode var = VAR_CONSERVATIVE; 109 PetscScalar sum; 110 111 PetscFunctionBeginUser; 112 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 113 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "TS conservation example", ""); 114 PetscCall(PetscOptionsEnum("-var", "Variable formulation", NULL, VarModes, (PetscEnum)var, (PetscEnum *)&var, NULL)); 115 PetscOptionsEnd(); 116 117 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 118 PetscCall(TSSetType(ts, TSBDF)); 119 PetscCall(TSGetDM(ts, &dm)); 120 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &U)); 121 PetscCall(VecSetValue(U, 0, 2., INSERT_VALUES)); 122 PetscCall(VecSetValue(U, 1, 1., INSERT_VALUES)); 123 switch (var) { 124 case VAR_CONSERVATIVE: PetscCall(DMTSSetIFunction(dm, IFunction_Conservative, NULL)); break; 125 case VAR_NONCONSERVATIVE: 126 PetscCall(VecLog(U)); 127 PetscCall(DMTSSetIFunction(dm, IFunction_Nonconservative, NULL)); 128 break; 129 case VAR_TRANSIENTVAR: 130 PetscCall(VecLog(U)); 131 PetscCall(DMTSSetIFunction(dm, IFunction_TransientVar, NULL)); 132 PetscCall(DMTSSetTransientVariable(dm, TransientVar, NULL)); 133 } 134 PetscCall(TSSetMaxTime(ts, 1.)); 135 PetscCall(TSSetFromOptions(ts)); 136 137 PetscCall(TSSolve(ts, U)); 138 switch (var) { 139 case VAR_CONSERVATIVE: break; 140 case VAR_NONCONSERVATIVE: 141 case VAR_TRANSIENTVAR: PetscCall(VecExp(U)); break; 142 } 143 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD)); 144 PetscCall(VecSum(U, &sum)); 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Conservation error %g\n", (double)PetscRealPart(sum - 3.))); 146 147 PetscCall(VecDestroy(&U)); 148 PetscCall(TSDestroy(&ts)); 149 PetscCall(PetscFinalize()); 150 return 0; 151 } 152 153 /*TEST 154 155 test: 156 suffix: conservative 157 args: -snes_fd -var conservative 158 test: 159 suffix: nonconservative 160 args: -snes_fd -var nonconservative 161 test: 162 suffix: transientvar 163 args: -snes_fd -var transientvar 164 165 TEST*/ 166