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 { 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(0); 60 } 61 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; 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(0); 79 } 80 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; 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(0); 98 } 99 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(0); 106 } 107 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