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 PetscBool nest; 38 39 PetscFunctionBeginUser; 40 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 41 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL)); 42 PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL)); 43 44 /* create vector to hold state */ 45 if (nest) { 46 Vec g[3]; 47 48 PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0])); 49 PetscCall(VecSetSizes(g[0], PETSC_DECIDE, 1)); 50 PetscCall(VecSetFromOptions(g[0])); 51 PetscCall(VecDuplicate(g[0], &g[1])); 52 PetscCall(VecDuplicate(g[0], &g[2])); 53 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global)); 54 PetscCall(VecDestroy(&g[0])); 55 PetscCall(VecDestroy(&g[1])); 56 PetscCall(VecDestroy(&g[2])); 57 } else { 58 PetscCall(VecCreate(PETSC_COMM_WORLD, &global)); 59 PetscCall(VecSetSizes(global, PETSC_DECIDE, 3)); 60 PetscCall(VecSetFromOptions(global)); 61 } 62 63 /* set initial conditions */ 64 PetscCall(Initial(global, NULL)); 65 66 /* make timestep context */ 67 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 68 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 69 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL)); 70 dt = 0.001; 71 72 /* 73 The user provides the RHS and Jacobian 74 */ 75 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 76 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 77 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3)); 78 PetscCall(MatSetFromOptions(A)); 79 PetscCall(MatSetUp(A)); 80 PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL)); 81 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 82 83 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, 3, 3, NULL, &S)); 84 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult)); 85 PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL)); 86 87 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 88 PetscCall(TSSetFromOptions(ts)); 89 90 PetscCall(TSSetTimeStep(ts, dt)); 91 PetscCall(TSSetMaxSteps(ts, time_steps)); 92 PetscCall(TSSetMaxTime(ts, 1)); 93 PetscCall(TSSetSolution(ts, global)); 94 95 PetscCall(TSSolve(ts, global)); 96 PetscCall(TSGetSolveTime(ts, &ftime)); 97 PetscCall(TSGetStepNumber(ts, &steps)); 98 99 /* free the memories */ 100 101 PetscCall(TSDestroy(&ts)); 102 PetscCall(VecDestroy(&global)); 103 PetscCall(MatDestroy(&A)); 104 PetscCall(MatDestroy(&S)); 105 106 PetscCall(PetscFinalize()); 107 return 0; 108 } 109 110 PetscErrorCode MyMatMult(Mat S, Vec x, Vec y) 111 { 112 const PetscScalar *inptr; 113 PetscScalar *outptr; 114 115 PetscFunctionBeginUser; 116 PetscCall(VecGetArrayRead(x, &inptr)); 117 PetscCall(VecGetArrayWrite(y, &outptr)); 118 119 outptr[0] = 2.0 * inptr[0] + inptr[1]; 120 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 121 outptr[2] = inptr[1] + 2.0 * inptr[2]; 122 123 PetscCall(VecRestoreArrayRead(x, &inptr)); 124 PetscCall(VecRestoreArrayWrite(y, &outptr)); 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127 128 /* this test problem has initial values (1,1,1). */ 129 PetscErrorCode Initial(Vec global, void *ctx) 130 { 131 PetscScalar *localptr; 132 PetscInt i, mybase, myend, locsize; 133 134 PetscFunctionBeginUser; 135 /* determine starting point of each processor */ 136 PetscCall(VecGetOwnershipRange(global, &mybase, &myend)); 137 PetscCall(VecGetLocalSize(global, &locsize)); 138 139 /* Initialize the array */ 140 PetscCall(VecGetArrayWrite(global, &localptr)); 141 for (i = 0; i < locsize; i++) localptr[i] = 1.0; 142 143 if (mybase == 0) localptr[0] = 1.0; 144 145 PetscCall(VecRestoreArrayWrite(global, &localptr)); 146 PetscFunctionReturn(PETSC_SUCCESS); 147 } 148 149 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx) 150 { 151 VecScatter scatter; 152 IS from, to; 153 PetscInt i, n, *idx; 154 Vec tmp_vec; 155 const PetscScalar *tmp; 156 157 PetscFunctionBeginUser; 158 /* Get the size of the vector */ 159 PetscCall(VecGetSize(global, &n)); 160 161 /* Set the index sets */ 162 PetscCall(PetscMalloc1(n, &idx)); 163 for (i = 0; i < n; i++) idx[i] = i; 164 165 /* Create local sequential vectors */ 166 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec)); 167 168 /* Create scatter context */ 169 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 170 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 171 PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter)); 172 PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 173 PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD)); 174 175 PetscCall(VecGetArrayRead(tmp_vec, &tmp)); 176 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]))); 177 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)))); 178 PetscCall(VecRestoreArrayRead(tmp_vec, &tmp)); 179 PetscCall(VecScatterDestroy(&scatter)); 180 PetscCall(ISDestroy(&from)); 181 PetscCall(ISDestroy(&to)); 182 PetscCall(PetscFree(idx)); 183 PetscCall(VecDestroy(&tmp_vec)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx) 188 { 189 PetscScalar *outptr; 190 const PetscScalar *inptr; 191 PetscInt i, n, *idx; 192 IS from, to; 193 VecScatter scatter; 194 Vec tmp_in, tmp_out; 195 196 PetscFunctionBeginUser; 197 /* Get the length of parallel vector */ 198 PetscCall(VecGetSize(globalin, &n)); 199 200 /* Set the index sets */ 201 PetscCall(PetscMalloc1(n, &idx)); 202 for (i = 0; i < n; i++) idx[i] = i; 203 204 /* Create local sequential vectors */ 205 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in)); 206 PetscCall(VecDuplicate(tmp_in, &tmp_out)); 207 208 /* Create scatter context */ 209 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from)); 210 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to)); 211 PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter)); 212 PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 213 PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD)); 214 PetscCall(VecScatterDestroy(&scatter)); 215 216 /*Extract income array */ 217 PetscCall(VecGetArrayRead(tmp_in, &inptr)); 218 219 /* Extract outcome array*/ 220 PetscCall(VecGetArrayWrite(tmp_out, &outptr)); 221 222 outptr[0] = 2.0 * inptr[0] + inptr[1]; 223 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 224 outptr[2] = inptr[1] + 2.0 * inptr[2]; 225 226 PetscCall(VecRestoreArrayRead(tmp_in, &inptr)); 227 PetscCall(VecRestoreArrayWrite(tmp_out, &outptr)); 228 229 PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter)); 230 PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 231 PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD)); 232 233 /* Destroy idx aand scatter */ 234 PetscCall(ISDestroy(&from)); 235 PetscCall(ISDestroy(&to)); 236 PetscCall(VecScatterDestroy(&scatter)); 237 PetscCall(VecDestroy(&tmp_in)); 238 PetscCall(VecDestroy(&tmp_out)); 239 PetscCall(PetscFree(idx)); 240 PetscFunctionReturn(PETSC_SUCCESS); 241 } 242 243 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx) 244 { 245 PetscScalar v[3]; 246 const PetscScalar *tmp; 247 PetscInt idx[3], i; 248 249 PetscFunctionBeginUser; 250 idx[0] = 0; 251 idx[1] = 1; 252 idx[2] = 2; 253 PetscCall(VecGetArrayRead(x, &tmp)); 254 255 i = 0; 256 v[0] = 2.0; 257 v[1] = 1.0; 258 v[2] = 0.0; 259 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 260 261 i = 1; 262 v[0] = 1.0; 263 v[1] = 2.0; 264 v[2] = 1.0; 265 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 266 267 i = 2; 268 v[0] = 0.0; 269 v[1] = 1.0; 270 v[2] = 2.0; 271 PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES)); 272 273 PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY)); 274 PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY)); 275 276 if (A != BB) { 277 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 278 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 279 } 280 PetscCall(VecRestoreArrayRead(x, &tmp)); 281 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 /* 286 The exact solutions 287 */ 288 PetscReal solx(PetscReal t) 289 { 290 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)); 291 } 292 293 PetscReal soly(PetscReal t) 294 { 295 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); 296 } 297 298 PetscReal solz(PetscReal t) 299 { 300 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)); 301 } 302 303 /*TEST 304 305 test: 306 suffix: euler 307 args: -ts_type euler -nest {{0 1}} 308 requires: !single 309 310 test: 311 suffix: beuler 312 args: -ts_type beuler -nest {{0 1}} 313 requires: !single 314 315 test: 316 suffix: rk 317 args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor 318 requires: !single 319 320 TEST*/ 321