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