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 CHKERRQ(VecGetArrayRead(U,&u)); 46 CHKERRQ(VecGetArrayRead(Udot,&udot)); 47 CHKERRQ(VecGetArray(F,&f)); 48 49 f[0] = udot[0] + u[0]; 50 f[1] = udot[1] - u[0]; 51 52 CHKERRQ(VecRestoreArrayRead(U,&u)); 53 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 54 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 65 CHKERRQ(VecGetArrayRead(Udot,&udot)); 66 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 72 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 73 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 84 CHKERRQ(VecGetArrayRead(Cdot,&cdot)); 85 CHKERRQ(VecGetArray(F,&f)); 86 87 f[0] = cdot[0] + PetscExpScalar(u[0]); 88 f[1] = cdot[1] - PetscExpScalar(u[0]); 89 90 CHKERRQ(VecRestoreArrayRead(U,&u)); 91 CHKERRQ(VecRestoreArrayRead(Cdot,&cdot)); 92 CHKERRQ(VecRestoreArray(F,&f)); 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode TransientVar(TS ts,Vec U,Vec C,void *ctx) 97 { 98 PetscFunctionBegin; 99 CHKERRQ(VecCopy(U,C)); 100 CHKERRQ(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 PetscErrorCode ierr; 112 113 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 114 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS conservation example","");CHKERRQ(ierr); 115 CHKERRQ(PetscOptionsEnum("-var","Variable formulation",NULL,VarModes,(PetscEnum)var,(PetscEnum*)&var,NULL)); 116 ierr = PetscOptionsEnd();CHKERRQ(ierr); 117 118 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 119 CHKERRQ(TSSetType(ts,TSBDF)); 120 CHKERRQ(TSGetDM(ts,&dm)); 121 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,2,&U)); 122 CHKERRQ(VecSetValue(U,0,2.,INSERT_VALUES)); 123 CHKERRQ(VecSetValue(U,1,1.,INSERT_VALUES)); 124 switch (var) { 125 case VAR_CONSERVATIVE: 126 CHKERRQ(DMTSSetIFunction(dm,IFunction_Conservative,NULL)); 127 break; 128 case VAR_NONCONSERVATIVE: 129 CHKERRQ(VecLog(U)); 130 CHKERRQ(DMTSSetIFunction(dm,IFunction_Nonconservative,NULL)); 131 break; 132 case VAR_TRANSIENTVAR: 133 CHKERRQ(VecLog(U)); 134 CHKERRQ(DMTSSetIFunction(dm,IFunction_TransientVar,NULL)); 135 CHKERRQ(DMTSSetTransientVariable(dm,TransientVar,NULL)); 136 } 137 CHKERRQ(TSSetMaxTime(ts,1.)); 138 CHKERRQ(TSSetFromOptions(ts)); 139 140 CHKERRQ(TSSolve(ts,U)); 141 switch (var) { 142 case VAR_CONSERVATIVE: 143 break; 144 case VAR_NONCONSERVATIVE: 145 case VAR_TRANSIENTVAR: 146 CHKERRQ(VecExp(U)); 147 break; 148 } 149 CHKERRQ(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 150 CHKERRQ(VecSum(U,&sum)); 151 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Conservation error %g\n", PetscRealPart(sum - 3.))); 152 153 CHKERRQ(VecDestroy(&U)); 154 CHKERRQ(TSDestroy(&ts)); 155 CHKERRQ(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