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 PetscFunctionBeginUser; 68 PetscCall(PetscInitialize(&argc,&argv,(char*)0,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