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 CHKERRQ(VecGetArrayRead(U,&u)); 25 CHKERRQ(VecGetArray(F,&f)); 26 27 f[0] = PetscCosReal(t); 28 f[1] = PetscSinReal(u[1]); 29 30 CHKERRQ(VecRestoreArrayRead(U,&u)); 31 CHKERRQ(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 CHKERRQ(VecGetArray(U,&u)); 44 45 u[0] = PetscSinReal(t); 46 u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5)); 47 48 CHKERRQ(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 PetscErrorCode ierr; 58 PetscMPIInt size; 59 PetscInt n = 2; 60 PetscScalar *u; 61 PetscReal t, final_time = 1.0, dt = 0.25; 62 PetscReal error; 63 TSAdapt adapt; 64 65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 Initialize program 67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 68 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 69 CHKERRMPI(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 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 76 CHKERRQ(TSSetType(ts,TSROSW)); 77 78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 Set ODE routines 80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 81 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 82 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Set initial conditions 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&U)); 88 CHKERRQ(VecSetSizes(U,n,PETSC_DETERMINE)); 89 CHKERRQ(VecSetUp(U)); 90 CHKERRQ(VecGetArray(U,&u)); 91 u[0] = 0.0; 92 u[1] = 1.0; 93 CHKERRQ(VecRestoreArray(U,&u)); 94 CHKERRQ(TSSetSolution(ts,U)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Set solver options 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 CHKERRQ(TSSetSaveTrajectory(ts)); 100 CHKERRQ(TSSetMaxTime(ts,final_time)); 101 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 102 CHKERRQ(TSSetTimeStep(ts,dt)); 103 /* The adaptive time step controller is forced to take constant time steps. */ 104 CHKERRQ(TSGetAdapt(ts,&adapt)); 105 CHKERRQ(TSAdaptSetType(adapt,TSADAPTNONE)); 106 107 CHKERRQ(TSSetFromOptions(ts)); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Run timestepping solver and compute error 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 CHKERRQ(TSSolve(ts,U)); 113 CHKERRQ(TSGetTime(ts,&t)); 114 115 if (PetscAbsReal(t-final_time)>100*PETSC_MACHINE_EPSILON) { 116 CHKERRQ(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 CHKERRQ(VecDuplicate(U,&Uex)); 119 CHKERRQ(ExactSolution(t,Uex)); 120 121 CHKERRQ(VecAYPX(Uex,-1.0,U)); 122 CHKERRQ(VecNorm(Uex,NORM_2,&error)); 123 CHKERRQ(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 CHKERRQ(VecDestroy(&U)); 129 CHKERRQ(VecDestroy(&Uex)); 130 CHKERRQ(TSDestroy(&ts)); 131 132 ierr = PetscFinalize(); 133 return ierr; 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