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 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: 125 PetscCall(DMTSSetIFunction(dm,IFunction_Conservative,NULL)); 126 break; 127 case VAR_NONCONSERVATIVE: 128 PetscCall(VecLog(U)); 129 PetscCall(DMTSSetIFunction(dm,IFunction_Nonconservative,NULL)); 130 break; 131 case VAR_TRANSIENTVAR: 132 PetscCall(VecLog(U)); 133 PetscCall(DMTSSetIFunction(dm,IFunction_TransientVar,NULL)); 134 PetscCall(DMTSSetTransientVariable(dm,TransientVar,NULL)); 135 } 136 PetscCall(TSSetMaxTime(ts,1.)); 137 PetscCall(TSSetFromOptions(ts)); 138 139 PetscCall(TSSolve(ts,U)); 140 switch (var) { 141 case VAR_CONSERVATIVE: 142 break; 143 case VAR_NONCONSERVATIVE: 144 case VAR_TRANSIENTVAR: 145 PetscCall(VecExp(U)); 146 break; 147 } 148 PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 149 PetscCall(VecSum(U,&sum)); 150 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Conservation error %g\n", (double)PetscRealPart(sum - 3.))); 151 152 PetscCall(VecDestroy(&U)); 153 PetscCall(TSDestroy(&ts)); 154 PetscCall(PetscFinalize()); 155 return 0; 156 } 157 158 /*TEST 159 160 test: 161 suffix: conservative 162 args: -snes_fd -var conservative 163 test: 164 suffix: nonconservative 165 args: -snes_fd -var nonconservative 166 test: 167 suffix: transientvar 168 args: -snes_fd -var transientvar 169 170 TEST*/ 171