1 /* 2 * ex_vdp.c 3 * 4 * Created on: Jun 1, 2012 5 * Author: Hong Zhang 6 */ 7 static char help[] = "Solves the van der Pol equation. \n Input parameters include:\n"; 8 9 /* 10 * Processors:1 11 */ 12 13 /* 14 * This program solves the van der Pol equation 15 * y' = z (1) 16 * z' = (((1-y^2)*z-y)/eps (2) 17 * on the domain 0<=x<=0.5, with the initial conditions 18 * y(0) = 2, 19 * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3 20 * IMEX schemes are applied to the splitted equation 21 * [y'] = [z] + [0 ] 22 * [z'] [0] [(((1-y^2)*z-y)/eps] 23 * 24 * F(x)= [z] 25 * [0] 26 * 27 * G(x)= [y'] - [0 ] 28 * [z'] [(((1-y^2)*z-y)/eps] 29 * 30 * JG(x) = G_x + a G_xdot 31 */ 32 33 #include <petscdmda.h> 34 #include <petscts.h> 35 36 typedef struct _User *User; 37 struct _User { 38 PetscReal mu; /*stiffness control coefficient: epsilon*/ 39 }; 40 41 static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 42 static PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*); 43 static PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 44 45 46 int main(int argc, char **argv) 47 { 48 TS ts; 49 Vec x; /*solution vector*/ 50 Mat A; /*Jacobian*/ 51 PetscInt steps,mx,eimex_rowcol[2],two; 52 PetscErrorCode ierr; 53 PetscScalar *x_ptr; 54 PetscReal ftime,dt,norm; 55 Vec ref; 56 struct _User user; /* user-defined work context */ 57 PetscViewer viewer; 58 59 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 60 /* Initialize user application context */ 61 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options",""); 62 user.mu = 1e0; 63 ierr = PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL);CHKERRQ(ierr); 64 ierr = PetscOptionsEnd();CHKERRQ(ierr); 65 66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67 Set runtime options 68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 69 /* 70 ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 71 */ 72 73 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 Create necessary matrix and vectors, solve same ODE on every process 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 77 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 78 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 79 ierr = MatSetUp(A);CHKERRQ(ierr); 80 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 81 82 ierr = MatCreateVecs(A,&ref,NULL);CHKERRQ(ierr); 83 ierr = VecGetArray(ref,&x_ptr);CHKERRQ(ierr); 84 /* 85 * [0,1], mu=10^-3 86 */ 87 /* 88 x_ptr[0] = -1.8881254106283; 89 x_ptr[1] = 0.7359074233370;*/ 90 91 /* 92 * [0,0.5],mu=10^-3 93 */ 94 /* 95 x_ptr[0] = 1.596980778659137; 96 x_ptr[1] = -1.029103015879544; 97 */ 98 /* 99 * [0,0.5],mu=1 100 */ 101 x_ptr[0] = 1.619084329683235; 102 x_ptr[1] = -0.803530465176385; 103 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Create timestepping solver context 106 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 107 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 108 ierr = TSSetType(ts,TSEIMEX);CHKERRQ(ierr); 109 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 110 ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr); 111 ierr = TSSetIJacobian(ts,A,A,IJacobian,&user);CHKERRQ(ierr); 112 113 dt = 0.00001; 114 ftime = 1.1; 115 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 116 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 117 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Set initial conditions 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr); 122 x_ptr[0] = 2.; 123 x_ptr[1] = -2./3. + 10./81.*(user.mu) - 292./2187.* (user.mu) * (user.mu) 124 -1814./19683.*(user.mu)*(user.mu)*(user.mu); 125 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 126 ierr = VecGetSize(x,&mx);CHKERRQ(ierr); 127 128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 Set runtime options 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 132 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Solve nonlinear system 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 ierr = TSSolve(ts,x);CHKERRQ(ierr); 137 ierr = TSGetTime(ts,&ftime);CHKERRQ(ierr); 138 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 139 140 ierr = VecAXPY(x,-1.0,ref);CHKERRQ(ierr); 141 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 142 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 143 144 eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2; 145 ierr = PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL);CHKERRQ(ierr); 146 ierr = PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1");CHKERRQ(ierr); 147 ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr); 148 ierr = PetscPrintf(PETSC_COMM_WORLD,"(%D,%D) %10.8f %18.15f %18.15f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm,(double)PetscRealPart(x_ptr[0]),(double)PetscRealPart(x_ptr[1]));CHKERRQ(ierr); 149 ierr = VecRestoreArray(x,&x_ptr);CHKERRQ(ierr); 150 151 /* Write line in convergence log */ 152 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr); 153 ierr = PetscViewerSetType(viewer,PETSCVIEWERASCII);CHKERRQ(ierr); 154 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_APPEND);CHKERRQ(ierr); 155 ierr = PetscViewerFileSetName(viewer,"eimex_nonstiff_vdp.txt");CHKERRQ(ierr); 156 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %10.8f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm);CHKERRQ(ierr); 157 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 158 159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 Free work space. 161 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 162 ierr = MatDestroy(&A);CHKERRQ(ierr); 163 ierr = VecDestroy(&x);CHKERRQ(ierr); 164 ierr = VecDestroy(&ref);CHKERRQ(ierr); 165 ierr = TSDestroy(&ts);CHKERRQ(ierr); 166 ierr = PetscFinalize(); 167 return ierr; 168 } 169 170 171 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 172 { 173 PetscErrorCode ierr; 174 PetscScalar *f; 175 const PetscScalar *x; 176 177 PetscFunctionBegin; 178 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 179 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 180 f[0] = x[1]; 181 f[1] = 0.0; 182 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 183 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr) 188 { 189 User user = (User)ptr; 190 PetscScalar *f; 191 const PetscScalar *x,*xdot; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 196 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 197 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 198 f[0] = xdot[0]; 199 f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu; 200 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 201 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 202 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr) 207 { 208 PetscErrorCode ierr; 209 User user = (User)ptr; 210 PetscReal mu = user->mu; 211 PetscInt rowcol[] = {0,1}; 212 PetscScalar J[2][2]; 213 const PetscScalar *x; 214 215 PetscFunctionBegin; 216 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 217 J[0][0] = a; 218 J[0][1] = 0; 219 J[1][0] = (2.*x[0]*x[1]+1.)/mu; 220 J[1][1] = a - (1. - x[0]*x[0])/mu; 221 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 222 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 223 224 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 if (A != B) { 227 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 228 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 229 } 230 PetscFunctionReturn(0); 231 } 232 233 /*TEST 234 235 test: 236 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution 237 requires: x 238 239 test: 240 suffix: adapt 241 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution 242 requires: x 243 244 test: 245 suffix: loop 246 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt {{0.005 0.001 0.0005}separate output} -ts_max_steps {{100 500 1000}separate output} -ts_eimex_row_col {{1,1 2,1 3,1 2,2 3,2 3,3}separate output} 247 requires: x 248 249 TEST*/ 250 251