1 /* 2 Formatted test for TS routines. 3 4 Solves U_t=F(t,u) 5 Where: 6 7 [2*u1+u2 8 F(t,u)= [u1+2*u2+u3 9 [ u2+2*u3 10 We can compare the solutions from euler, beuler and SUNDIALS to 11 see what is the difference. 12 13 */ 14 15 static char help[] = "Solves a linear ODE. \n\n"; 16 17 #include <petscts.h> 18 #include <petscpc.h> 19 20 extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 21 extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 22 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); 23 extern PetscErrorCode Initial(Vec, void *); 24 extern PetscErrorCode MyMatMult(Mat, Vec, Vec); 25 26 extern PetscReal solx(PetscReal); 27 extern PetscReal soly(PetscReal); 28 extern PetscReal solz(PetscReal); 29 30 int main(int argc, char **argv) 31 { 32 PetscInt time_steps = 100, steps; 33 Vec global; 34 PetscReal dt, ftime; 35 TS ts; 36 Mat A = 0, S; 37 38 PetscFunctionBeginUser; 39 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 40 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL)); 41 42 /* set initial conditions */ 43 PetscCall(VecCreate(PETSC_COMM_WORLD, &global)); 44 PetscCall(VecSetSizes(global, PETSC_DECIDE, 3)); 45 PetscCall(VecSetFromOptions(global)); 46 PetscCall(Initial(global, NULL)); 47 48 /* make timestep context */ 49 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 50 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 51 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL)); 52 dt = 0.001; 53 54 /* 55 The user provides the RHS and Jacobian 56 */ 57 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 58 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 59 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3)); 60 PetscCall(MatSetFromOptions(A)); 61 PetscCall(MatSetUp(A)); 62 PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL)); 63 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 64 65 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, 3, 3, NULL, &S)); 66 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult)); 67 PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL)); 68 69 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 70 PetscCall(TSSetFromOptions(ts)); 71 72 PetscCall(TSSetTimeStep(ts, dt)); 73 PetscCall(TSSetMaxSteps(ts, time_steps)); 74 PetscCall(TSSetMaxTime(ts, 1)); 75 PetscCall(TSSetSolution(ts, global)); 76 77 PetscCall(TSSolve(ts, global)); 78 PetscCall(TSGetSolveTime(ts, &ftime)); 79 PetscCall(TSGetStepNumber(ts, &steps)); 80 81 /* free the memories */ 82 83 PetscCall(TSDestroy(&ts)); 84 PetscCall(VecDestroy(&global)); 85 PetscCall(MatDestroy(&A)); 86 PetscCall(MatDestroy(&S)); 87 88 PetscCall(PetscFinalize()); 89 return 0; 90 } 91 92 PetscErrorCode MyMatMult(Mat S, Vec x, Vec y) 93 { 94 const PetscScalar *inptr; 95 PetscScalar *outptr; 96 97 PetscFunctionBeginUser; 98 PetscCall(VecGetArrayRead(x, &inptr)); 99 PetscCall(VecGetArrayWrite(y, &outptr)); 100 101 outptr[0] = 2.0 * inptr[0] + inptr[1]; 102 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 103 outptr[2] = inptr[1] + 2.0 * inptr[2]; 104 105 PetscCall(VecRestoreArrayRead(x, &inptr)); 106 PetscCall(VecRestoreArrayWrite(y, &outptr)); 107 PetscFunctionReturn(0); 108 } 109 110 /* this test problem has initial values (1,1,1). */ 111 PetscErrorCode Initial(Vec global, void *ctx) 112 { 113 PetscScalar *localptr; 114 PetscInt i, mybase, myend, locsize; 115 116 /* determine starting point of each processor */ 117 PetscCall(VecGetOwnershipRange(global, &mybase, &myend)); 118 PetscCall(VecGetLocalSize(global, &locsize)); 119 120 /* Initialize the array */ 121 PetscCall(VecGetArrayWrite(global, &localptr)); 122 for (i = 0; i < locsize; i++) localptr[i] = 1.0; 123 124 if (mybase == 0) localptr[0] = 1.0; 125 126 PetscCall(VecRestoreArrayWrite(global, &localptr)); 127 return 0; 128 } 129 130 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx) 131 { 132 VecScatter scatter; 133 IS from, to; 134 PetscInt i, n, *idx; 135 Vec tmp_vec; 136 const PetscScalar *tmp; 137 138 /* Get the size of the vector */ 139 PetscCall(VecGetSize(global, &n)); 140 141 /* Set the index sets */ 142 PetscCall(PetscMalloc1(n, &idx)); 143 for (i = 0; i < n; i++) idx[i] = i; 144 145 /* Create local sequential vectors */ 146 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec)); 147 148 /* Create scatter context */ 149 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 150 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 151 PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter)); 152 PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 153 PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 154 155 PetscCall(VecGetArrayRead(tmp_vec, &tmp)); 156 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e u = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0]), (double)PetscRealPart(tmp[1]), (double)PetscRealPart(tmp[2]))); 157 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - solx(time)), (double)PetscRealPart(tmp[1] - soly(time)), (double)PetscRealPart(tmp[2] - solz(time)))); 158 PetscCall(VecRestoreArrayRead(tmp_vec, &tmp)); 159 PetscCall(VecScatterDestroy(&scatter)); 160 PetscCall(ISDestroy(&from)); 161 PetscCall(ISDestroy(&to)); 162 PetscCall(PetscFree(idx)); 163 PetscCall(VecDestroy(&tmp_vec)); 164 return 0; 165 } 166 167 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx) 168 { 169 PetscScalar *outptr; 170 const PetscScalar *inptr; 171 PetscInt i, n, *idx; 172 IS from, to; 173 VecScatter scatter; 174 Vec tmp_in, tmp_out; 175 176 /* Get the length of parallel vector */ 177 PetscCall(VecGetSize(globalin, &n)); 178 179 /* Set the index sets */ 180 PetscCall(PetscMalloc1(n, &idx)); 181 for (i = 0; i < n; i++) idx[i] = i; 182 183 /* Create local sequential vectors */ 184 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in)); 185 PetscCall(VecDuplicate(tmp_in, &tmp_out)); 186 187 /* Create scatter context */ 188 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 189 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 190 PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter)); 191 PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 192 PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 193 PetscCall(VecScatterDestroy(&scatter)); 194 195 /*Extract income array */ 196 PetscCall(VecGetArrayRead(tmp_in, &inptr)); 197 198 /* Extract outcome array*/ 199 PetscCall(VecGetArrayWrite(tmp_out, &outptr)); 200 201 outptr[0] = 2.0 * inptr[0] + inptr[1]; 202 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 203 outptr[2] = inptr[1] + 2.0 * inptr[2]; 204 205 PetscCall(VecRestoreArrayRead(tmp_in, &inptr)); 206 PetscCall(VecRestoreArrayWrite(tmp_out, &outptr)); 207 208 PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter)); 209 PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 210 PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 211 212 /* Destroy idx aand scatter */ 213 PetscCall(ISDestroy(&from)); 214 PetscCall(ISDestroy(&to)); 215 PetscCall(VecScatterDestroy(&scatter)); 216 PetscCall(VecDestroy(&tmp_in)); 217 PetscCall(VecDestroy(&tmp_out)); 218 PetscCall(PetscFree(idx)); 219 return 0; 220 } 221 222 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx) 223 { 224 PetscScalar v[3]; 225 const PetscScalar *tmp; 226 PetscInt idx[3], i; 227 228 idx[0] = 0; 229 idx[1] = 1; 230 idx[2] = 2; 231 PetscCall(VecGetArrayRead(x, &tmp)); 232 233 i = 0; 234 v[0] = 2.0; 235 v[1] = 1.0; 236 v[2] = 0.0; 237 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 238 239 i = 1; 240 v[0] = 1.0; 241 v[1] = 2.0; 242 v[2] = 1.0; 243 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 244 245 i = 2; 246 v[0] = 0.0; 247 v[1] = 1.0; 248 v[2] = 2.0; 249 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 250 251 PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY)); 252 PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY)); 253 254 if (A != BB) { 255 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 256 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 257 } 258 PetscCall(VecRestoreArrayRead(x, &tmp)); 259 260 return 0; 261 } 262 263 /* 264 The exact solutions 265 */ 266 PetscReal solx(PetscReal t) 267 { 268 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)); 269 } 270 271 PetscReal soly(PetscReal t) 272 { 273 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0); 274 } 275 276 PetscReal solz(PetscReal t) 277 { 278 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)); 279 } 280 281 /*TEST 282 283 test: 284 suffix: euler 285 args: -ts_type euler 286 requires: !single 287 288 test: 289 suffix: beuler 290 args: -ts_type beuler 291 requires: !single 292 293 TEST*/ 294