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 int main(int argc, char **argv) 46 { 47 TS ts; 48 Vec x; /*solution vector*/ 49 Mat A; /*Jacobian*/ 50 PetscInt steps,mx,eimex_rowcol[2],two; 51 PetscErrorCode ierr; 52 PetscScalar *x_ptr; 53 PetscReal ftime,dt,norm; 54 Vec ref; 55 struct _User user; /* user-defined work context */ 56 PetscViewer viewer; 57 58 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 59 /* Initialize user application context */ 60 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options",""); 61 user.mu = 1e0; 62 ierr = PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL);CHKERRQ(ierr); 63 ierr = PetscOptionsEnd();CHKERRQ(ierr); 64 65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 Set runtime options 67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 68 /* 69 ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 70 */ 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create necessary matrix and vectors, solve same ODE on every process 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 76 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr); 77 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 78 ierr = MatSetUp(A);CHKERRQ(ierr); 79 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 80 81 ierr = MatCreateVecs(A,&ref,NULL);CHKERRQ(ierr); 82 ierr = VecGetArray(ref,&x_ptr);CHKERRQ(ierr); 83 /* 84 * [0,1], mu=10^-3 85 */ 86 /* 87 x_ptr[0] = -1.8881254106283; 88 x_ptr[1] = 0.7359074233370;*/ 89 90 /* 91 * [0,0.5],mu=10^-3 92 */ 93 /* 94 x_ptr[0] = 1.596980778659137; 95 x_ptr[1] = -1.029103015879544; 96 */ 97 /* 98 * [0,0.5],mu=1 99 */ 100 x_ptr[0] = 1.619084329683235; 101 x_ptr[1] = -0.803530465176385; 102 103 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 104 Create timestepping solver context 105 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 106 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 107 ierr = TSSetType(ts,TSEIMEX);CHKERRQ(ierr); 108 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 109 ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr); 110 ierr = TSSetIJacobian(ts,A,A,IJacobian,&user);CHKERRQ(ierr); 111 112 dt = 0.00001; 113 ftime = 1.1; 114 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 115 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 116 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118 Set initial conditions 119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 120 ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr); 121 x_ptr[0] = 2.; 122 x_ptr[1] = -2./3. + 10./81.*(user.mu) - 292./2187.* (user.mu) * (user.mu) 123 -1814./19683.*(user.mu)*(user.mu)*(user.mu); 124 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 125 ierr = VecGetSize(x,&mx);CHKERRQ(ierr); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Set runtime options 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 131 132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 Solve nonlinear system 134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 135 ierr = TSSolve(ts,x);CHKERRQ(ierr); 136 ierr = TSGetTime(ts,&ftime);CHKERRQ(ierr); 137 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 138 139 ierr = VecAXPY(x,-1.0,ref);CHKERRQ(ierr); 140 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 141 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 142 143 eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2; 144 ierr = PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL);CHKERRQ(ierr); 145 ierr = PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1");CHKERRQ(ierr); 146 ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr); 147 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); 148 ierr = VecRestoreArray(x,&x_ptr);CHKERRQ(ierr); 149 150 /* Write line in convergence log */ 151 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr); 152 ierr = PetscViewerSetType(viewer,PETSCVIEWERASCII);CHKERRQ(ierr); 153 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_APPEND);CHKERRQ(ierr); 154 ierr = PetscViewerFileSetName(viewer,"eimex_nonstiff_vdp.txt");CHKERRQ(ierr); 155 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %10.8f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm);CHKERRQ(ierr); 156 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Free work space. 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 ierr = MatDestroy(&A);CHKERRQ(ierr); 162 ierr = VecDestroy(&x);CHKERRQ(ierr); 163 ierr = VecDestroy(&ref);CHKERRQ(ierr); 164 ierr = TSDestroy(&ts);CHKERRQ(ierr); 165 ierr = PetscFinalize(); 166 return ierr; 167 } 168 169 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 170 { 171 PetscErrorCode ierr; 172 PetscScalar *f; 173 const PetscScalar *x; 174 175 PetscFunctionBegin; 176 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 177 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 178 f[0] = x[1]; 179 f[1] = 0.0; 180 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 181 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr) 186 { 187 User user = (User)ptr; 188 PetscScalar *f; 189 const PetscScalar *x,*xdot; 190 PetscErrorCode ierr; 191 192 PetscFunctionBegin; 193 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 194 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 195 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 196 f[0] = xdot[0]; 197 f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu; 198 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 199 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 200 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr) 205 { 206 PetscErrorCode ierr; 207 User user = (User)ptr; 208 PetscReal mu = user->mu; 209 PetscInt rowcol[] = {0,1}; 210 PetscScalar J[2][2]; 211 const PetscScalar *x; 212 213 PetscFunctionBegin; 214 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 215 J[0][0] = a; 216 J[0][1] = 0; 217 J[1][0] = (2.*x[0]*x[1]+1.)/mu; 218 J[1][1] = a - (1. - x[0]*x[0])/mu; 219 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 220 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 221 222 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 223 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 224 if (A != B) { 225 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 227 } 228 PetscFunctionReturn(0); 229 } 230 231 /*TEST 232 233 test: 234 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 235 requires: x 236 237 test: 238 suffix: adapt 239 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 240 requires: x 241 242 test: 243 suffix: loop 244 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} 245 requires: x 246 247 TEST*/ 248