1 2 static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n"; 3 /* 4 p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy 5 xmin < x < xmax, ymin < y < ymax; 6 7 Boundary conditions Neumman using mirror values 8 9 Note that x_t and y_t in the above are given functions of x and y; they are not derivatives of x and y. 10 x_t = (y - ws) y_t = (ws/2H)*(Pm - Pmax*sin(x)) 11 12 */ 13 14 #include <petscdm.h> 15 #include <petscdmda.h> 16 #include <petscts.h> 17 18 /* 19 User-defined data structures and routines 20 */ 21 typedef struct { 22 PetscScalar ws; /* Synchronous speed */ 23 PetscScalar H; /* Inertia constant */ 24 PetscScalar D; /* Damping constant */ 25 PetscScalar Pmax; /* Maximum power output of generator */ 26 PetscScalar PM_min; /* Mean mechanical power input */ 27 PetscScalar lambda; /* correlation time */ 28 PetscScalar q; /* noise strength */ 29 PetscScalar mux; /* Initial average angle */ 30 PetscScalar sigmax; /* Standard deviation of initial angle */ 31 PetscScalar muy; /* Average speed */ 32 PetscScalar sigmay; /* standard deviation of initial speed */ 33 PetscScalar rho; /* Cross-correlation coefficient */ 34 PetscScalar xmin; /* left boundary of angle */ 35 PetscScalar xmax; /* right boundary of angle */ 36 PetscScalar ymin; /* bottom boundary of speed */ 37 PetscScalar ymax; /* top boundary of speed */ 38 PetscScalar dx; /* x step size */ 39 PetscScalar dy; /* y step size */ 40 PetscScalar disper_coe; /* Dispersion coefficient */ 41 DM da; 42 PetscInt st_width; /* Stencil width */ 43 DMBoundaryType bx; /* x boundary type */ 44 DMBoundaryType by; /* y boundary type */ 45 PetscBool nonoiseinitial; 46 } AppCtx; 47 48 PetscErrorCode Parameter_settings(AppCtx*); 49 PetscErrorCode ini_bou(Vec,AppCtx*); 50 PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*); 51 PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 52 PetscErrorCode PostStep(TS); 53 54 int main(int argc, char **argv) 55 { 56 Vec x; /* Solution vector */ 57 TS ts; /* Time-stepping context */ 58 AppCtx user; /* Application context */ 59 PetscViewer viewer; 60 61 PetscFunctionBeginUser; 62 PetscCall(PetscInitialize(&argc,&argv,"petscopt_ex7", help)); 63 64 /* Get physics and time parameters */ 65 PetscCall(Parameter_settings(&user)); 66 /* Create a 2D DA with dof = 1 */ 67 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,user.bx,user.by,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,user.st_width,NULL,NULL,&user.da)); 68 PetscCall(DMSetFromOptions(user.da)); 69 PetscCall(DMSetUp(user.da)); 70 /* Set x and y coordinates */ 71 PetscCall(DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0,0)); 72 PetscCall(DMDASetCoordinateName(user.da,0,"X - the angle")); 73 PetscCall(DMDASetCoordinateName(user.da,1,"Y - the speed")); 74 75 /* Get global vector x from DM */ 76 PetscCall(DMCreateGlobalVector(user.da,&x)); 77 78 PetscCall(ini_bou(x,&user)); 79 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer)); 80 PetscCall(VecView(x,viewer)); 81 PetscCall(PetscViewerDestroy(&viewer)); 82 83 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 84 PetscCall(TSSetDM(ts,user.da)); 85 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 86 PetscCall(TSSetType(ts,TSARKIMEX)); 87 PetscCall(TSSetIFunction(ts,NULL,IFunction,&user)); 88 /* PetscCall(TSSetIJacobian(ts,NULL,NULL,IJacobian,&user)); */ 89 PetscCall(TSSetApplicationContext(ts,&user)); 90 PetscCall(TSSetTimeStep(ts,.005)); 91 PetscCall(TSSetFromOptions(ts)); 92 PetscCall(TSSetPostStep(ts,PostStep)); 93 PetscCall(TSSolve(ts,x)); 94 95 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer)); 96 PetscCall(VecView(x,viewer)); 97 PetscCall(PetscViewerDestroy(&viewer)); 98 99 PetscCall(VecDestroy(&x)); 100 PetscCall(DMDestroy(&user.da)); 101 PetscCall(TSDestroy(&ts)); 102 PetscCall(PetscFinalize()); 103 return 0; 104 } 105 106 PetscErrorCode PostStep(TS ts) 107 { 108 Vec X,gc; 109 AppCtx *user; 110 PetscScalar sum = 0,asum; 111 PetscReal t,**p; 112 DMDACoor2d **coors; 113 DM cda; 114 PetscInt i,j,xs,ys,xm,ym; 115 116 PetscFunctionBegin; 117 PetscCall(TSGetApplicationContext(ts,&user)); 118 PetscCall(TSGetTime(ts,&t)); 119 PetscCall(TSGetSolution(ts,&X)); 120 121 PetscCall(DMGetCoordinateDM(user->da,&cda)); 122 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 123 PetscCall(DMGetCoordinates(user->da,&gc)); 124 PetscCall(DMDAVecGetArrayRead(cda,gc,&coors)); 125 PetscCall(DMDAVecGetArrayRead(user->da,X,&p)); 126 for (i=xs; i < xs+xm; i++) { 127 for (j=ys; j < ys+ym; j++) { 128 if (coors[j][i].y < 5) sum += p[j][i]; 129 } 130 } 131 PetscCall(DMDAVecRestoreArrayRead(cda,gc,&coors)); 132 PetscCall(DMDAVecRestoreArrayRead(user->da,X,&p)); 133 PetscCallMPI(MPI_Allreduce(&sum,&asum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ts))); 134 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %f = %f\n",(double)t,(double)(asum))); 135 if (sum < 1.0e-2) { 136 PetscCall(TSSetConvergedReason(ts,TS_CONVERGED_USER)); 137 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Exiting TS as the integral of PDF is almost zero\n")); 138 } 139 PetscFunctionReturn(0); 140 } 141 142 PetscErrorCode ini_bou(Vec X,AppCtx* user) 143 { 144 DM cda; 145 DMDACoor2d **coors; 146 PetscScalar **p; 147 Vec gc; 148 PetscInt i,j; 149 PetscInt xs,ys,xm,ym,M,N; 150 PetscScalar xi,yi; 151 PetscScalar sigmax=user->sigmax,sigmay=user->sigmay; 152 PetscScalar rho =user->rho; 153 PetscScalar muy=user->muy,mux; 154 PetscMPIInt rank; 155 PetscScalar sum; 156 157 PetscFunctionBeginUser; 158 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 159 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 160 user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1); 161 PetscCall(DMGetCoordinateDM(user->da,&cda)); 162 PetscCall(DMGetCoordinates(user->da,&gc)); 163 PetscCall(DMDAVecGetArray(cda,gc,&coors)); 164 PetscCall(DMDAVecGetArray(user->da,X,&p)); 165 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 166 167 /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable 168 muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point 169 in the y-direction. We only modify mux here 170 */ 171 mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */ 172 if (user->nonoiseinitial) { 173 for (i=xs; i < xs+xm; i++) { 174 for (j=ys; j < ys+ym; j++) { 175 xi = coors[j][i].x; yi = coors[j][i].y; 176 if ((xi == mux) && (yi == muy)) { 177 p[j][i] = 1.0; 178 } 179 } 180 } 181 } else { 182 /* Change PM_min accordingly */ 183 user->PM_min = user->Pmax*PetscSinScalar(mux); 184 for (i=xs; i < xs+xm; i++) { 185 for (j=ys; j < ys+ym; j++) { 186 xi = coors[j][i].x; yi = coors[j][i].y; 187 p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay))); 188 } 189 } 190 } 191 PetscCall(DMDAVecRestoreArray(cda,gc,&coors)); 192 PetscCall(DMDAVecRestoreArray(user->da,X,&p)); 193 PetscCall(VecSum(X,&sum)); 194 PetscCall(VecScale(X,1.0/sum)); 195 PetscFunctionReturn(0); 196 } 197 198 /* First advection term */ 199 PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user) 200 { 201 PetscScalar f,fpos,fneg; 202 PetscFunctionBegin; 203 f = (y - user->ws); 204 fpos = PetscMax(f,0); 205 fneg = PetscMin(f,0); 206 if (user->st_width == 1) { 207 *p1 = fpos*(p[j][i] - p[j][i-1])/user->dx + fneg*(p[j][i+1] - p[j][i])/user->dx; 208 } else if (user->st_width == 2) { 209 *p1 = fpos*(3*p[j][i] - 4*p[j][i-1] + p[j][i-2])/(2*user->dx) + fneg*(-p[j][i+2] + 4*p[j][i+1] - 3*p[j][i])/(2*user->dx); 210 } else if (user->st_width == 3) { 211 *p1 = fpos*(2*p[j][i+1] + 3*p[j][i] - 6*p[j][i-1] + p[j][i-2])/(6*user->dx) + fneg*(-p[j][i+2] + 6*p[j][i+1] - 3*p[j][i] - 2*p[j][i-1])/(6*user->dx); 212 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils"); 213 PetscFunctionReturn(0); 214 } 215 216 /* Second advection term */ 217 PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user) 218 { 219 PetscScalar f,fpos,fneg; 220 PetscFunctionBegin; 221 f = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x)); 222 fpos = PetscMax(f,0); 223 fneg = PetscMin(f,0); 224 if (user->st_width == 1) { 225 *p2 = fpos*(p[j][i] - p[j-1][i])/user->dy + fneg*(p[j+1][i] - p[j][i])/user->dy; 226 } else if (user->st_width ==2) { 227 *p2 = fpos*(3*p[j][i] - 4*p[j-1][i] + p[j-2][i])/(2*user->dy) + fneg*(-p[j+2][i] + 4*p[j+1][i] - 3*p[j][i])/(2*user->dy); 228 } else if (user->st_width == 3) { 229 *p2 = fpos*(2*p[j+1][i] + 3*p[j][i] - 6*p[j-1][i] + p[j-2][i])/(6*user->dy) + fneg*(-p[j+2][i] + 6*p[j+1][i] - 3*p[j][i] - 2*p[j-1][i])/(6*user->dy); 230 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils"); 231 PetscFunctionReturn(0); 232 } 233 234 /* Diffusion term */ 235 PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user) 236 { 237 PetscFunctionBeginUser; 238 if (user->st_width == 1) { 239 *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy)); 240 } else if (user->st_width == 2) { 241 *p_diff = user->disper_coe*((-p[j-2][i] + 16*p[j-1][i] - 30*p[j][i] + 16*p[j+1][i] - p[j+2][i])/(12.0*user->dy*user->dy)); 242 } else if (user->st_width == 3) { 243 *p_diff = user->disper_coe*((2*p[j-3][i] - 27*p[j-2][i] + 270*p[j-1][i] - 490*p[j][i] + 270*p[j+1][i] - 27*p[j+2][i] + 2*p[j+3][i])/(180.0*user->dy*user->dy)); 244 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils"); 245 PetscFunctionReturn(0); 246 } 247 248 PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 249 { 250 AppCtx *user=(AppCtx*)ctx; 251 DM cda; 252 DMDACoor2d **coors; 253 PetscScalar **p,**f,**pdot; 254 PetscInt i,j; 255 PetscInt xs,ys,xm,ym,M,N; 256 Vec localX,gc,localXdot; 257 PetscScalar p_adv1 = 0.0,p_adv2 = 0.0,p_diff; /* initialize to prevent incorrect compiler warnings */ 258 259 PetscFunctionBeginUser; 260 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 261 PetscCall(DMGetCoordinateDM(user->da,&cda)); 262 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 263 264 PetscCall(DMGetLocalVector(user->da,&localX)); 265 PetscCall(DMGetLocalVector(user->da,&localXdot)); 266 267 PetscCall(DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX)); 268 PetscCall(DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX)); 269 PetscCall(DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot)); 270 PetscCall(DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot)); 271 272 PetscCall(DMGetCoordinatesLocal(user->da,&gc)); 273 274 PetscCall(DMDAVecGetArrayRead(cda,gc,&coors)); 275 PetscCall(DMDAVecGetArrayRead(user->da,localX,&p)); 276 PetscCall(DMDAVecGetArrayRead(user->da,localXdot,&pdot)); 277 PetscCall(DMDAVecGetArray(user->da,F,&f)); 278 279 user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)); 280 for (i=xs; i < xs+xm; i++) { 281 for (j=ys; j < ys+ym; j++) { 282 PetscCall(adv1(p,coors[j][i].y,i,j,M,&p_adv1,user)); 283 PetscCall(adv2(p,coors[j][i].x,i,j,N,&p_adv2,user)); 284 PetscCall(diffuse(p,i,j,t,&p_diff,user)); 285 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i]; 286 } 287 } 288 PetscCall(DMDAVecRestoreArrayRead(user->da,localX,&p)); 289 PetscCall(DMDAVecRestoreArrayRead(user->da,localX,&pdot)); 290 PetscCall(DMRestoreLocalVector(user->da,&localX)); 291 PetscCall(DMRestoreLocalVector(user->da,&localXdot)); 292 PetscCall(DMDAVecRestoreArray(user->da,F,&f)); 293 PetscCall(DMDAVecRestoreArrayRead(cda,gc,&coors)); 294 295 PetscFunctionReturn(0); 296 } 297 298 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx) 299 { 300 AppCtx *user=(AppCtx*)ctx; 301 DM cda; 302 DMDACoor2d **coors; 303 PetscInt i,j; 304 PetscInt xs,ys,xm,ym,M,N; 305 Vec gc; 306 PetscScalar val[5],xi,yi; 307 MatStencil row,col[5]; 308 PetscScalar c1,c3,c5,c1pos,c1neg,c3pos,c3neg; 309 310 PetscFunctionBeginUser; 311 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 312 PetscCall(DMGetCoordinateDM(user->da,&cda)); 313 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 314 315 PetscCall(DMGetCoordinatesLocal(user->da,&gc)); 316 PetscCall(DMDAVecGetArrayRead(cda,gc,&coors)); 317 for (i=xs; i < xs+xm; i++) { 318 for (j=ys; j < ys+ym; j++) { 319 PetscInt nc = 0; 320 xi = coors[j][i].x; yi = coors[j][i].y; 321 row.i = i; row.j = j; 322 c1 = (yi-user->ws)/user->dx; 323 c1pos = PetscMax(c1,0); 324 c1neg = PetscMin(c1,0); 325 c3 = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/user->dy; 326 c3pos = PetscMax(c3,0); 327 c3neg = PetscMin(c3,0); 328 c5 = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy); 329 col[nc].i = i-1; col[nc].j = j; val[nc++] = c1pos; 330 col[nc].i = i+1; col[nc].j = j; val[nc++] = -c1neg; 331 col[nc].i = i; col[nc].j = j-1; val[nc++] = c3pos + c5; 332 col[nc].i = i; col[nc].j = j+1; val[nc++] = -c3neg + c5; 333 col[nc].i = i; col[nc].j = j; val[nc++] = -c1pos + c1neg -c3pos + c3neg -2*c5 -a; 334 PetscCall(MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES)); 335 } 336 } 337 PetscCall(DMDAVecRestoreArrayRead(cda,gc,&coors)); 338 339 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY)); 340 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY)); 341 if (J != Jpre) { 342 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 343 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode Parameter_settings(AppCtx *user) 349 { 350 PetscBool flg; 351 352 PetscFunctionBeginUser; 353 354 /* Set default parameters */ 355 user->ws = 1.0; 356 user->H = 5.0; 357 user->Pmax = 2.1; 358 user->PM_min = 1.0; 359 user->lambda = 0.1; 360 user->q = 1.0; 361 user->mux = PetscAsinScalar(user->PM_min/user->Pmax); 362 user->sigmax = 0.1; 363 user->sigmay = 0.1; 364 user->rho = 0.0; 365 user->xmin = -PETSC_PI; 366 user->xmax = PETSC_PI; 367 user->bx = DM_BOUNDARY_PERIODIC; 368 user->by = DM_BOUNDARY_MIRROR; 369 user->nonoiseinitial = PETSC_FALSE; 370 371 /* 372 ymin of -3 seems to let the unstable solution move up and leave a zero in its wake 373 with an ymin of -1 the wake is never exactly zero 374 */ 375 user->ymin = -3.0; 376 user->ymax = 10.0; 377 user->st_width = 1; 378 379 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg)); 380 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg)); 381 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg)); 382 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg)); 383 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg)); 384 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg)); 385 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg)); 386 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg)); 387 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg)); 388 if (flg == 0) { 389 user->muy = user->ws; 390 } 391 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg)); 392 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg)); 393 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg)); 394 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg)); 395 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg)); 396 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg)); 397 PetscCall(PetscOptionsGetInt(NULL,NULL,"-stencil_width",&user->st_width,&flg)); 398 PetscCall(PetscOptionsGetEnum(NULL,NULL,"-bx",DMBoundaryTypes,(PetscEnum*)&user->bx,&flg)); 399 PetscCall(PetscOptionsGetEnum(NULL,NULL,"-by",DMBoundaryTypes,(PetscEnum*)&user->by,&flg)); 400 PetscCall(PetscOptionsGetBool(NULL,NULL,"-nonoiseinitial",&user->nonoiseinitial,&flg)); 401 402 PetscFunctionReturn(0); 403 } 404 405 /*TEST 406 407 build: 408 requires: !complex !single 409 410 test: 411 args: -ts_max_steps 2 412 localrunfiles: petscopt_ex7 413 414 test: 415 suffix: 2 416 args: -ts_max_steps 2 -snes_mf_operator 417 output_file: output/ex7_1.out 418 localrunfiles: petscopt_ex7 419 timeoutfactor: 2 420 421 test: 422 suffix: 3 423 args: -ts_max_steps 2 -snes_mf -pc_type none 424 output_file: output/ex7_1.out 425 localrunfiles: petscopt_ex7 426 timeoutfactor: 2 427 428 TEST*/ 429