1 2 static char help[] = "Solves the van der Pol DAE.\n\ 3 Input parameters include:\n"; 4 5 /* 6 Concepts: TS^time-dependent nonlinear problems 7 Concepts: TS^van der Pol DAE 8 Processors: 1 9 */ 10 /* ------------------------------------------------------------------------ 11 12 This program solves the van der Pol DAE 13 y' = -z = f(y,z) (1) 14 0 = y-(z^3/3 - z) = g(y,z) 15 on the domain 0 <= x <= 1, with the boundary conditions 16 y(0) = -2, y'(0) = -2.355301397608119909925287735864250951918 17 This is a nonlinear equation. 18 19 Notes: 20 This code demonstrates the TS solver interface with the Van der Pol DAE, 21 namely it is the case when there is no RHS (meaning the RHS == 0), and the 22 equations are converted to two variants of linear problems, u_t = f(u,t), 23 namely turning (1) into a vector equation in terms of u, 24 25 [ y' + z ] = [ 0 ] 26 [ (z^3/3 - z) - y ] [ 0 ] 27 28 which then we can write as a vector equation 29 30 [ u_1' + u_2 ] = [ 0 ] (2) 31 [ (u_2^3/3 - u_2) - u_1 ] [ 0 ] 32 33 which is now in the desired form of u_t = f(u,t). As this is a DAE, and 34 there is no u_2', there is no need for a split, 35 36 so 37 38 [ G(u',u,t) ] = [ u_1' ] + [ u_2 ] 39 [ 0 ] [ (u_2^3/3 - u_2) - u_1 ] 40 41 Using the definition of the Jacobian of G (from the PETSc user manual), 42 in the equation G(u',u,t) = F(u,t), 43 44 dG dG 45 J(G) = a * -- - -- 46 du' du 47 48 where d is the partial derivative. In this example, 49 50 dG [ 1 ; 0 ] 51 -- = [ ] 52 du' [ 0 ; 0 ] 53 54 dG [ 0 ; 1 ] 55 -- = [ ] 56 du [ -1 ; 1 - u_2^2 ] 57 58 Hence, 59 60 [ a ; -1 ] 61 J(G) = [ ] 62 [ 1 ; u_2^2 - 1 ] 63 64 ------------------------------------------------------------------------- */ 65 66 #include <petscts.h> 67 68 typedef struct _n_User *User; 69 struct _n_User { 70 PetscReal next_output; 71 }; 72 73 /* 74 * User-defined routines 75 */ 76 77 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 78 { 79 PetscErrorCode ierr; 80 PetscScalar *f; 81 const PetscScalar *x,*xdot; 82 83 PetscFunctionBeginUser; 84 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 85 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 86 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 87 f[0] = xdot[0] + x[1]; 88 f[1] = (x[1]*x[1]*x[1]/3.0 - x[1])-x[0]; 89 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 90 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 91 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 96 { 97 PetscErrorCode ierr; 98 PetscInt rowcol[] = {0,1}; 99 PetscScalar J[2][2]; 100 const PetscScalar *x; 101 102 PetscFunctionBeginUser; 103 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 104 J[0][0] = a; J[0][1] = -1.; 105 J[1][0] = 1.; J[1][1] = -1. + x[1]*x[1]; 106 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 107 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 108 109 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 110 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 111 if (A != B) { 112 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 113 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 114 } 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode RegisterMyARK2(void) 119 { 120 PetscErrorCode ierr; 121 122 PetscFunctionBeginUser; 123 { 124 const PetscReal 125 A[3][3] = {{0,0,0}, 126 {0.41421356237309504880,0,0}, 127 {0.75,0.25,0}}, 128 At[3][3] = {{0,0,0}, 129 {0.12132034355964257320,0.29289321881345247560,0}, 130 {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}}, 131 *bembedt = NULL,*bembed = NULL; 132 ierr = TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);CHKERRQ(ierr); 133 } 134 PetscFunctionReturn(0); 135 } 136 137 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 138 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx) 139 { 140 PetscErrorCode ierr; 141 const PetscScalar *x; 142 PetscReal tfinal, dt; 143 User user = (User)ctx; 144 Vec interpolatedX; 145 146 PetscFunctionBeginUser; 147 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 148 ierr = TSGetMaxTime(ts,&tfinal);CHKERRQ(ierr); 149 150 while (user->next_output <= t && user->next_output <= tfinal) { 151 ierr = VecDuplicate(X,&interpolatedX);CHKERRQ(ierr); 152 ierr = TSInterpolate(ts,user->next_output,interpolatedX);CHKERRQ(ierr); 153 ierr = VecGetArrayRead(interpolatedX,&x);CHKERRQ(ierr); 154 ierr = PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %3D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",(double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),(double)PetscRealPart(x[1]));CHKERRQ(ierr); 155 ierr = VecRestoreArrayRead(interpolatedX,&x);CHKERRQ(ierr); 156 ierr = VecDestroy(&interpolatedX);CHKERRQ(ierr); 157 user->next_output += PetscRealConstant(0.1); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 int main(int argc,char **argv) 163 { 164 TS ts; /* nonlinear solver */ 165 Vec x; /* solution, residual vectors */ 166 Mat A; /* Jacobian matrix */ 167 PetscInt steps; 168 PetscReal ftime = 0.5; 169 PetscBool monitor = PETSC_FALSE; 170 PetscScalar *x_ptr; 171 PetscMPIInt size; 172 struct _n_User user; 173 PetscErrorCode ierr; 174 175 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 176 Initialize program 177 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 178 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 179 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 180 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 181 182 ierr = RegisterMyARK2();CHKERRQ(ierr); 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Set runtime options 186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187 188 user.next_output = 0.0; 189 ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 190 191 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192 Create necessary matrix and vectors, solve same ODE on every process 193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 195 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 196 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 197 ierr = MatSetUp(A);CHKERRQ(ierr); 198 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 199 200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 201 Create timestepping solver context 202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 203 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 204 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 205 ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr); 206 ierr = TSSetIJacobian(ts,A,A,IJacobian,&user);CHKERRQ(ierr); 207 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 208 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 209 if (monitor) { 210 ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr); 211 } 212 213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 214 Set initial conditions 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr); 217 x_ptr[0] = -2; x_ptr[1] = -2.355301397608119909925287735864250951918; 218 ierr = VecRestoreArray(x,&x_ptr);CHKERRQ(ierr); 219 ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 220 221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 222 Set runtime options 223 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 224 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 225 226 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227 Solve nonlinear system 228 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 229 ierr = TSSolve(ts,x);CHKERRQ(ierr); 230 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 231 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 232 ierr = PetscPrintf(PETSC_COMM_WORLD,"steps %3D, ftime %g\n",steps,(double)ftime);CHKERRQ(ierr); 233 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 234 235 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 236 Free work space. All PETSc objects should be destroyed when they 237 are no longer needed. 238 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 239 ierr = MatDestroy(&A);CHKERRQ(ierr); 240 ierr = VecDestroy(&x);CHKERRQ(ierr); 241 ierr = TSDestroy(&ts);CHKERRQ(ierr); 242 243 ierr = PetscFinalize(); 244 return ierr; 245 } 246 247 /*TEST 248 249 test: 250 requires: !single 251 suffix: a 252 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp 253 output_file: output/ex19_pi42.out 254 255 test: 256 requires: !single 257 suffix: b 258 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_filter PI42 259 output_file: output/ex19_pi42.out 260 261 test: 262 requires: !single 263 suffix: c 264 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_pid 0.4,0.2 265 output_file: output/ex19_pi42.out 266 267 TEST*/ 268