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