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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 59 /* Initialize user application context */ 60 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options","");PetscCall(ierr); 61 user.mu = 1e0; 62 PetscCall(PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL)); 63 ierr = PetscOptionsEnd();PetscCall(ierr); 64 65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 Set runtime options 67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 68 /* 69 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 70 */ 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create necessary matrix and vectors, solve same ODE on every process 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 76 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 77 PetscCall(MatSetFromOptions(A)); 78 PetscCall(MatSetUp(A)); 79 PetscCall(MatCreateVecs(A,&x,NULL)); 80 81 PetscCall(MatCreateVecs(A,&ref,NULL)); 82 PetscCall(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 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 107 PetscCall(TSSetType(ts,TSEIMEX)); 108 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,&user)); 109 PetscCall(TSSetIFunction(ts,NULL,IFunction,&user)); 110 PetscCall(TSSetIJacobian(ts,A,A,IJacobian,&user)); 111 112 dt = 0.00001; 113 ftime = 1.1; 114 PetscCall(TSSetTimeStep(ts,dt)); 115 PetscCall(TSSetMaxTime(ts,ftime)); 116 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118 Set initial conditions 119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 120 PetscCall(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 PetscCall(TSSetSolution(ts,x)); 125 PetscCall(VecGetSize(x,&mx)); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Set runtime options 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 PetscCall(TSSetFromOptions(ts)); 131 132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 Solve nonlinear system 134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 135 PetscCall(TSSolve(ts,x)); 136 PetscCall(TSGetTime(ts,&ftime)); 137 PetscCall(TSGetStepNumber(ts,&steps)); 138 139 PetscCall(VecAXPY(x,-1.0,ref)); 140 PetscCall(VecNorm(x,NORM_2,&norm)); 141 PetscCall(TSGetTimeStep(ts,&dt)); 142 143 eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2; 144 PetscCall(PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL)); 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1")); 146 PetscCall(VecGetArray(x,&x_ptr)); 147 PetscCall(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 PetscCall(VecRestoreArray(x,&x_ptr)); 149 150 /* Write line in convergence log */ 151 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&viewer)); 152 PetscCall(PetscViewerSetType(viewer,PETSCVIEWERASCII)); 153 PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_APPEND)); 154 PetscCall(PetscViewerFileSetName(viewer,"eimex_nonstiff_vdp.txt")); 155 PetscCall(PetscViewerASCIIPrintf(viewer,"%D %D %10.8f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm)); 156 PetscCall(PetscViewerDestroy(&viewer)); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Free work space. 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 PetscCall(MatDestroy(&A)); 162 PetscCall(VecDestroy(&x)); 163 PetscCall(VecDestroy(&ref)); 164 PetscCall(TSDestroy(&ts)); 165 PetscCall(PetscFinalize()); 166 return 0; 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 PetscCall(VecGetArrayRead(X,&x)); 176 PetscCall(VecGetArray(F,&f)); 177 f[0] = x[1]; 178 f[1] = 0.0; 179 PetscCall(VecRestoreArrayRead(X,&x)); 180 PetscCall(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 PetscCall(VecGetArrayRead(X,&x)); 192 PetscCall(VecGetArrayRead(Xdot,&xdot)); 193 PetscCall(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 PetscCall(VecRestoreArrayRead(X,&x)); 197 PetscCall(VecRestoreArrayRead(Xdot,&xdot)); 198 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 217 PetscCall(VecRestoreArrayRead(X,&x)); 218 219 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 220 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 221 if (A != B) { 222 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 223 PetscCall(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