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