1 static char help[] = "Small ODE to test TS accuracy.\n"; 2 3 /* 4 The ODE 5 u1_t = cos(t), 6 u2_t = sin(u2) 7 with analytical solution 8 u1(t) = sin(t), 9 u2(t) = 2 * atan(exp(t) * tan(0.5)) 10 is used to test the accuracy of TS schemes. 11 */ 12 13 #include <petscts.h> 14 15 /* 16 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U) 17 */ 18 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s) 19 { 20 PetscScalar *f; 21 const PetscScalar *u; 22 23 PetscFunctionBeginUser; 24 PetscCall(VecGetArrayRead(U, &u)); 25 PetscCall(VecGetArray(F, &f)); 26 27 f[0] = PetscCosReal(t); 28 f[1] = PetscSinScalar(u[1]); 29 30 PetscCall(VecRestoreArrayRead(U, &u)); 31 PetscCall(VecRestoreArray(F, &f)); 32 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 /* 36 Defines the exact solution. 37 */ 38 static PetscErrorCode ExactSolution(PetscReal t, Vec U) 39 { 40 PetscScalar *u; 41 42 PetscFunctionBeginUser; 43 PetscCall(VecGetArray(U, &u)); 44 45 u[0] = PetscSinReal(t); 46 u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5)); 47 48 PetscCall(VecRestoreArray(U, &u)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 int main(int argc, char **argv) 53 { 54 TS ts; /* ODE integrator */ 55 Vec U; /* numerical solution will be stored here */ 56 Vec Uex; /* analytical (exact) solution will be stored here */ 57 PetscMPIInt size; 58 PetscInt n = 2; 59 PetscScalar *u; 60 PetscReal t, final_time = 1.0, dt = 0.25; 61 PetscReal error; 62 TSAdapt adapt; 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Initialize program 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 67 PetscFunctionBeginUser; 68 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 69 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 70 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create timestepping solver context 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 76 PetscCall(TSSetType(ts, TSROSW)); 77 78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 Set ODE routines 80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 81 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 82 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Set initial conditions 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 88 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 89 PetscCall(VecSetUp(U)); 90 PetscCall(VecGetArray(U, &u)); 91 u[0] = 0.0; 92 u[1] = 1.0; 93 PetscCall(VecRestoreArray(U, &u)); 94 PetscCall(TSSetSolution(ts, U)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Set solver options 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 PetscCall(TSSetSaveTrajectory(ts)); 100 PetscCall(TSSetMaxTime(ts, final_time)); 101 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 102 PetscCall(TSSetTimeStep(ts, dt)); 103 /* The adaptive time step controller is forced to take constant time steps. */ 104 PetscCall(TSGetAdapt(ts, &adapt)); 105 PetscCall(TSAdaptSetType(adapt, TSADAPTNONE)); 106 107 PetscCall(TSSetFromOptions(ts)); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Run timestepping solver and compute error 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 PetscCall(TSSolve(ts, U)); 113 PetscCall(TSGetTime(ts, &t)); 114 115 if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) { 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference of %g between the prescribed final time %g and the actual final time.\n", (double)(final_time - t), (double)final_time)); 117 } 118 PetscCall(VecDuplicate(U, &Uex)); 119 PetscCall(ExactSolution(t, Uex)); 120 121 PetscCall(VecAYPX(Uex, -1.0, U)); 122 PetscCall(VecNorm(Uex, NORM_2, &error)); 123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error)); 124 125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126 Free work space. All PETSc objects should be destroyed when they are no longer needed. 127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 128 PetscCall(VecDestroy(&U)); 129 PetscCall(VecDestroy(&Uex)); 130 PetscCall(TSDestroy(&ts)); 131 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 test: 139 suffix: 3bs 140 args: -ts_type rk -ts_rk_type 3bs 141 requires: !single 142 143 test: 144 suffix: 5bs 145 args: -ts_type rk -ts_rk_type 5bs 146 requires: !single 147 148 test: 149 suffix: 5dp 150 args: -ts_type rk -ts_rk_type 5dp 151 requires: !single 152 153 test: 154 suffix: 6vr 155 args: -ts_type rk -ts_rk_type 6vr 156 requires: !single 157 158 test: 159 suffix: 7vr 160 args: -ts_type rk -ts_rk_type 7vr 161 requires: !single 162 163 test: 164 suffix: 8vr 165 args: -ts_type rk -ts_rk_type 8vr 166 requires: !single 167 168 TEST*/ 169