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 11 When run in parallel, each process solves the same set of equations separately. 12 */ 13 14 static char help[] = "Solves a linear ODE. \n\n"; 15 16 #include <petscts.h> 17 #include <petscpc.h> 18 19 extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 20 extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 21 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); 22 extern PetscErrorCode Initial(Vec, void *); 23 extern PetscErrorCode MyMatMult(Mat, Vec, Vec); 24 25 extern PetscReal solx(PetscReal); 26 extern PetscReal soly(PetscReal); 27 extern PetscReal solz(PetscReal); 28 29 int main(int argc, char **argv) 30 { 31 PetscInt time_steps = 100, steps; 32 Vec global; 33 PetscReal dt, ftime; 34 TS ts; 35 Mat A, S; 36 PetscBool nest = PETSC_FALSE; 37 38 PetscFunctionBeginUser; 39 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 40 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL)); 41 PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL)); 42 43 /* create vector to hold state */ 44 if (nest) { 45 Vec g[3]; 46 47 PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0])); 48 PetscCall(VecSetSizes(g[0], 1, PETSC_DECIDE)); 49 PetscCall(VecSetFromOptions(g[0])); 50 PetscCall(VecDuplicate(g[0], &g[1])); 51 PetscCall(VecDuplicate(g[0], &g[2])); 52 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global)); 53 PetscCall(VecDestroy(&g[0])); 54 PetscCall(VecDestroy(&g[1])); 55 PetscCall(VecDestroy(&g[2])); 56 } else { 57 PetscCall(VecCreate(PETSC_COMM_WORLD, &global)); 58 PetscCall(VecSetSizes(global, 3, PETSC_DECIDE)); 59 PetscCall(VecSetFromOptions(global)); 60 } 61 62 /* set initial conditions */ 63 PetscCall(Initial(global, NULL)); 64 65 /* make timestep context */ 66 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 67 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 68 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL)); 69 dt = 0.001; 70 71 /* 72 The user provides the RHS and Jacobian 73 */ 74 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 75 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 76 PetscCall(MatSetSizes(A, 3, 3, PETSC_DECIDE, PETSC_DECIDE)); 77 PetscCall(MatSetFromOptions(A)); 78 PetscCall(MatSetUp(A)); 79 PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL)); 80 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 81 82 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, NULL, &S)); 83 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MyMatMult)); 84 PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL)); 85 86 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 87 PetscCall(TSSetFromOptions(ts)); 88 89 PetscCall(TSSetTimeStep(ts, dt)); 90 PetscCall(TSSetMaxSteps(ts, time_steps)); 91 PetscCall(TSSetMaxTime(ts, 1)); 92 PetscCall(TSSetSolution(ts, global)); 93 94 PetscCall(TSSolve(ts, global)); 95 PetscCall(TSGetSolveTime(ts, &ftime)); 96 PetscCall(TSGetStepNumber(ts, &steps)); 97 98 /* free the memory */ 99 PetscCall(TSDestroy(&ts)); 100 PetscCall(VecDestroy(&global)); 101 PetscCall(MatDestroy(&A)); 102 PetscCall(MatDestroy(&S)); 103 104 PetscCall(PetscFinalize()); 105 return 0; 106 } 107 108 PetscErrorCode MyMatMult(Mat S, Vec x, Vec y) 109 { 110 const PetscScalar *inptr; 111 PetscScalar *outptr; 112 113 PetscFunctionBeginUser; 114 PetscCall(VecGetArrayRead(x, &inptr)); 115 PetscCall(VecGetArrayWrite(y, &outptr)); 116 117 outptr[0] = 2.0 * inptr[0] + inptr[1]; 118 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 119 outptr[2] = inptr[1] + 2.0 * inptr[2]; 120 121 PetscCall(VecRestoreArrayRead(x, &inptr)); 122 PetscCall(VecRestoreArrayWrite(y, &outptr)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 PetscErrorCode Initial(Vec global, void *ctx) 127 { 128 PetscScalar *localptr; 129 130 PetscFunctionBeginUser; 131 PetscCall(VecGetArrayWrite(global, &localptr)); 132 localptr[0] = solx(0.0); 133 localptr[1] = soly(0.0); 134 localptr[2] = solz(0.0); 135 PetscCall(VecRestoreArrayWrite(global, &localptr)); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx) 140 { 141 const PetscScalar *tmp; 142 PetscScalar exact[] = {solx(time), soly(time), solz(time)}; 143 144 PetscFunctionBeginUser; 145 PetscCall(VecGetArrayRead(global, &tmp)); 146 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]))); 147 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - exact[0]), (double)PetscRealPart(tmp[1] - exact[1]), (double)PetscRealPart(tmp[2] - exact[2]))); 148 PetscCall(VecRestoreArrayRead(global, &tmp)); 149 PetscFunctionReturn(PETSC_SUCCESS); 150 } 151 152 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx) 153 { 154 PetscScalar *outptr; 155 const PetscScalar *inptr; 156 157 PetscFunctionBeginUser; 158 /*Extract income array */ 159 PetscCall(VecGetArrayRead(globalin, &inptr)); 160 161 /* Extract outcome array*/ 162 PetscCall(VecGetArrayWrite(globalout, &outptr)); 163 164 outptr[0] = 2.0 * inptr[0] + inptr[1]; 165 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2]; 166 outptr[2] = inptr[1] + 2.0 * inptr[2]; 167 168 PetscCall(VecRestoreArrayRead(globalin, &inptr)); 169 PetscCall(VecRestoreArrayWrite(globalout, &outptr)); 170 PetscFunctionReturn(PETSC_SUCCESS); 171 } 172 173 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx) 174 { 175 PetscScalar v[3]; 176 PetscInt idx[3], rst; 177 178 PetscFunctionBeginUser; 179 PetscCall(VecGetOwnershipRange(x, &rst, NULL)); 180 idx[0] = 0 + rst; 181 idx[1] = 1 + rst; 182 idx[2] = 2 + rst; 183 184 v[0] = 2.0; 185 v[1] = 1.0; 186 v[2] = 0.0; 187 PetscCall(MatSetValues(BB, 1, idx, 3, idx, v, INSERT_VALUES)); 188 189 v[0] = 1.0; 190 v[1] = 2.0; 191 v[2] = 1.0; 192 PetscCall(MatSetValues(BB, 1, idx + 1, 3, idx, v, INSERT_VALUES)); 193 194 v[0] = 0.0; 195 v[1] = 1.0; 196 v[2] = 2.0; 197 PetscCall(MatSetValues(BB, 1, idx + 2, 3, idx, v, INSERT_VALUES)); 198 199 PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY)); 200 PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY)); 201 202 if (A != BB) { 203 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 204 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 205 } 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208 209 /* 210 The exact solutions 211 */ 212 PetscReal solx(PetscReal t) 213 { 214 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)); 215 } 216 217 PetscReal soly(PetscReal t) 218 { 219 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); 220 } 221 222 PetscReal solz(PetscReal t) 223 { 224 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)); 225 } 226 227 /*TEST 228 229 test: 230 suffix: euler 231 args: -ts_type euler -nest {{0 1}} 232 requires: !single 233 234 test: 235 suffix: beuler 236 args: -ts_type beuler -nest {{0 1}} 237 requires: !single 238 239 test: 240 suffix: rk 241 args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor 242 requires: !single 243 244 test: 245 diff_args: -j 246 requires: double !complex 247 output_file: output/ex2_be_adapt.out 248 suffix: bdf_1_adapt 249 args: -ts_type bdf -ts_bdf_order 1 -ts_adapt_type basic -ts_adapt_clip 0,2 250 251 test: 252 diff_args: -j 253 requires: double !complex 254 suffix: be_adapt 255 args: -ts_type beuler -ts_adapt_type basic -ts_adapt_clip 0,2 256 257 TEST*/ 258