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