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(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 117 /* determine starting point of each processor */ 118 PetscCall(VecGetOwnershipRange(global, &mybase, &myend)); 119 PetscCall(VecGetLocalSize(global, &locsize)); 120 121 /* Initialize the array */ 122 PetscCall(VecGetArrayWrite(global, &localptr)); 123 for (i = 0; i < locsize; i++) localptr[i] = 1.0; 124 125 if (mybase == 0) localptr[0] = 1.0; 126 127 PetscCall(VecRestoreArrayWrite(global, &localptr)); 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 131 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx) 132 { 133 VecScatter scatter; 134 IS from, to; 135 PetscInt i, n, *idx; 136 Vec tmp_vec; 137 const PetscScalar *tmp; 138 139 PetscFunctionBeginUser; 140 /* Get the size of the vector */ 141 PetscCall(VecGetSize(global, &n)); 142 143 /* Set the index sets */ 144 PetscCall(PetscMalloc1(n, &idx)); 145 for (i = 0; i < n; i++) idx[i] = i; 146 147 /* Create local sequential vectors */ 148 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec)); 149 150 /* Create scatter context */ 151 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 152 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 153 PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter)); 154 PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 155 PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 156 157 PetscCall(VecGetArrayRead(tmp_vec, &tmp)); 158 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]))); 159 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)))); 160 PetscCall(VecRestoreArrayRead(tmp_vec, &tmp)); 161 PetscCall(VecScatterDestroy(&scatter)); 162 PetscCall(ISDestroy(&from)); 163 PetscCall(ISDestroy(&to)); 164 PetscCall(PetscFree(idx)); 165 PetscCall(VecDestroy(&tmp_vec)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx) 170 { 171 PetscScalar *outptr; 172 const PetscScalar *inptr; 173 PetscInt i, n, *idx; 174 IS from, to; 175 VecScatter scatter; 176 Vec tmp_in, tmp_out; 177 178 PetscFunctionBeginUser; 179 /* Get the length of parallel vector */ 180 PetscCall(VecGetSize(globalin, &n)); 181 182 /* Set the index sets */ 183 PetscCall(PetscMalloc1(n, &idx)); 184 for (i = 0; i < n; i++) idx[i] = i; 185 186 /* Create local sequential vectors */ 187 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in)); 188 PetscCall(VecDuplicate(tmp_in, &tmp_out)); 189 190 /* Create scatter context */ 191 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 192 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 193 PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter)); 194 PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 195 PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 196 PetscCall(VecScatterDestroy(&scatter)); 197 198 /*Extract income array */ 199 PetscCall(VecGetArrayRead(tmp_in, &inptr)); 200 201 /* Extract outcome array*/ 202 PetscCall(VecGetArrayWrite(tmp_out, &outptr)); 203 204 outptr[0] = 2.0 * inptr[0] + inptr[1]; 205 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 206 outptr[2] = inptr[1] + 2.0 * inptr[2]; 207 208 PetscCall(VecRestoreArrayRead(tmp_in, &inptr)); 209 PetscCall(VecRestoreArrayWrite(tmp_out, &outptr)); 210 211 PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter)); 212 PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 213 PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 214 215 /* Destroy idx aand scatter */ 216 PetscCall(ISDestroy(&from)); 217 PetscCall(ISDestroy(&to)); 218 PetscCall(VecScatterDestroy(&scatter)); 219 PetscCall(VecDestroy(&tmp_in)); 220 PetscCall(VecDestroy(&tmp_out)); 221 PetscCall(PetscFree(idx)); 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx) 226 { 227 PetscScalar v[3]; 228 const PetscScalar *tmp; 229 PetscInt idx[3], i; 230 231 PetscFunctionBeginUser; 232 idx[0] = 0; 233 idx[1] = 1; 234 idx[2] = 2; 235 PetscCall(VecGetArrayRead(x, &tmp)); 236 237 i = 0; 238 v[0] = 2.0; 239 v[1] = 1.0; 240 v[2] = 0.0; 241 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 242 243 i = 1; 244 v[0] = 1.0; 245 v[1] = 2.0; 246 v[2] = 1.0; 247 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 248 249 i = 2; 250 v[0] = 0.0; 251 v[1] = 1.0; 252 v[2] = 2.0; 253 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 254 255 PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY)); 256 PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY)); 257 258 if (A != BB) { 259 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 260 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 261 } 262 PetscCall(VecRestoreArrayRead(x, &tmp)); 263 264 PetscFunctionReturn(PETSC_SUCCESS); 265 } 266 267 /* 268 The exact solutions 269 */ 270 PetscReal solx(PetscReal t) 271 { 272 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)); 273 } 274 275 PetscReal soly(PetscReal t) 276 { 277 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); 278 } 279 280 PetscReal solz(PetscReal t) 281 { 282 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)); 283 } 284 285 /*TEST 286 287 test: 288 suffix: euler 289 args: -ts_type euler 290 requires: !single 291 292 test: 293 suffix: beuler 294 args: -ts_type beuler 295 requires: !single 296 297 TEST*/ 298