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