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