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 PetscFunctionBegin; 24 PetscCall(VecGetArrayRead(U,&u)); 25 PetscCall(VecGetArray(F,&f)); 26 27 f[0] = PetscCosReal(t); 28 f[1] = PetscSinReal(u[1]); 29 30 PetscCall(VecRestoreArrayRead(U,&u)); 31 PetscCall(VecRestoreArray(F,&f)); 32 PetscFunctionReturn(0); 33 } 34 35 /* 36 Defines the exact solution. 37 */ 38 static PetscErrorCode ExactSolution(PetscReal t, Vec U) 39 { 40 PetscScalar *u; 41 42 PetscFunctionBegin; 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(0); 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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 68 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 69 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 70 71 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 72 Create timestepping solver context 73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 74 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 75 PetscCall(TSSetType(ts,TSROSW)); 76 77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 78 Set ODE routines 79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 80 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 81 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 82 83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84 Set initial conditions 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 PetscCall(VecCreate(PETSC_COMM_WORLD,&U)); 87 PetscCall(VecSetSizes(U,n,PETSC_DETERMINE)); 88 PetscCall(VecSetUp(U)); 89 PetscCall(VecGetArray(U,&u)); 90 u[0] = 0.0; 91 u[1] = 1.0; 92 PetscCall(VecRestoreArray(U,&u)); 93 PetscCall(TSSetSolution(ts,U)); 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Set solver options 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 PetscCall(TSSetSaveTrajectory(ts)); 99 PetscCall(TSSetMaxTime(ts,final_time)); 100 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 101 PetscCall(TSSetTimeStep(ts,dt)); 102 /* The adaptive time step controller is forced to take constant time steps. */ 103 PetscCall(TSGetAdapt(ts,&adapt)); 104 PetscCall(TSAdaptSetType(adapt,TSADAPTNONE)); 105 106 PetscCall(TSSetFromOptions(ts)); 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Run timestepping solver and compute error 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 PetscCall(TSSolve(ts,U)); 112 PetscCall(TSGetTime(ts,&t)); 113 114 if (PetscAbsReal(t-final_time)>100*PETSC_MACHINE_EPSILON) { 115 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)); 116 } 117 PetscCall(VecDuplicate(U,&Uex)); 118 PetscCall(ExactSolution(t,Uex)); 119 120 PetscCall(VecAYPX(Uex,-1.0,U)); 121 PetscCall(VecNorm(Uex,NORM_2,&error)); 122 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error at final time: %.2E\n",(double)error)); 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Free work space. All PETSc objects should be destroyed when they are no longer needed. 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 PetscCall(VecDestroy(&U)); 128 PetscCall(VecDestroy(&Uex)); 129 PetscCall(TSDestroy(&ts)); 130 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 /*TEST 136 137 test: 138 suffix: 3bs 139 args: -ts_type rk -ts_rk_type 3bs 140 requires: !single 141 142 test: 143 suffix: 5bs 144 args: -ts_type rk -ts_rk_type 5bs 145 requires: !single 146 147 test: 148 suffix: 5dp 149 args: -ts_type rk -ts_rk_type 5dp 150 requires: !single 151 152 test: 153 suffix: 6vr 154 args: -ts_type rk -ts_rk_type 6vr 155 requires: !single 156 157 test: 158 suffix: 7vr 159 args: -ts_type rk -ts_rk_type 7vr 160 requires: !single 161 162 test: 163 suffix: 8vr 164 args: -ts_type rk -ts_rk_type 8vr 165 requires: !single 166 167 TEST*/ 168