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