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