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