Actual source code: ex7.c

      1: static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n";
      2: /*
      3:    p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy
      4:    xmin < x < xmax, ymin < y < ymax;

      6:    Boundary conditions Neumman using mirror values

      8:    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.
      9:    x_t = (y - ws)  y_t = (ws/2H)*(Pm - Pmax*sin(x))

     11: */

     13: #include <petscdm.h>
     14: #include <petscdmda.h>
     15: #include <petscts.h>

     17: /*
     18:    User-defined data structures and routines
     19: */
     20: typedef struct {
     21:   PetscScalar    ws;         /* Synchronous speed */
     22:   PetscScalar    H;          /* Inertia constant */
     23:   PetscScalar    D;          /* Damping constant */
     24:   PetscScalar    Pmax;       /* Maximum power output of generator */
     25:   PetscScalar    PM_min;     /* Mean mechanical power input */
     26:   PetscScalar    lambda;     /* correlation time */
     27:   PetscScalar    q;          /* noise strength */
     28:   PetscScalar    mux;        /* Initial average angle */
     29:   PetscScalar    sigmax;     /* Standard deviation of initial angle */
     30:   PetscScalar    muy;        /* Average speed */
     31:   PetscScalar    sigmay;     /* standard deviation of initial speed */
     32:   PetscScalar    rho;        /* Cross-correlation coefficient */
     33:   PetscScalar    xmin;       /* left boundary of angle */
     34:   PetscScalar    xmax;       /* right boundary of angle */
     35:   PetscScalar    ymin;       /* bottom boundary of speed */
     36:   PetscScalar    ymax;       /* top boundary of speed */
     37:   PetscScalar    dx;         /* x step size */
     38:   PetscScalar    dy;         /* y step size */
     39:   PetscScalar    disper_coe; /* Dispersion coefficient */
     40:   DM             da;
     41:   PetscInt       st_width; /* Stencil width */
     42:   DMBoundaryType bx;       /* x boundary type */
     43:   DMBoundaryType by;       /* y boundary type */
     44:   PetscBool      nonoiseinitial;
     45: } AppCtx;

     47: PetscErrorCode Parameter_settings(AppCtx *);
     48: PetscErrorCode ini_bou(Vec, AppCtx *);
     49: PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     50: PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
     51: PetscErrorCode PostStep(TS);

     53: int main(int argc, char **argv)
     54: {
     55:   Vec         x;    /* Solution vector */
     56:   TS          ts;   /* Time-stepping context */
     57:   AppCtx      user; /* Application context */
     58:   PetscViewer viewer;

     60:   PetscFunctionBeginUser;
     61:   PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex7", help));

     63:   /* Get physics and time parameters */
     64:   PetscCall(Parameter_settings(&user));
     65:   /* Create a 2D DA with dof = 1 */
     66:   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));
     67:   PetscCall(DMSetFromOptions(user.da));
     68:   PetscCall(DMSetUp(user.da));
     69:   /* Set x and y coordinates */
     70:   PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0));
     71:   PetscCall(DMDASetCoordinateName(user.da, 0, "X - the angle"));
     72:   PetscCall(DMDASetCoordinateName(user.da, 1, "Y - the speed"));

     74:   /* Get global vector x from DM  */
     75:   PetscCall(DMCreateGlobalVector(user.da, &x));

     77:   PetscCall(ini_bou(x, &user));
     78:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer));
     79:   PetscCall(VecView(x, viewer));
     80:   PetscCall(PetscViewerDestroy(&viewer));

     82:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     83:   PetscCall(TSSetDM(ts, user.da));
     84:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
     85:   PetscCall(TSSetType(ts, TSARKIMEX));
     86:   PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
     87:   /*  PetscCall(TSSetIJacobian(ts,NULL,NULL,IJacobian,&user));  */
     88:   PetscCall(TSSetApplicationContext(ts, &user));
     89:   PetscCall(TSSetTimeStep(ts, .005));
     90:   PetscCall(TSSetFromOptions(ts));
     91:   PetscCall(TSSetPostStep(ts, PostStep));
     92:   PetscCall(TSSolve(ts, x));

     94:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer));
     95:   PetscCall(VecView(x, viewer));
     96:   PetscCall(PetscViewerDestroy(&viewer));

     98:   PetscCall(VecDestroy(&x));
     99:   PetscCall(DMDestroy(&user.da));
    100:   PetscCall(TSDestroy(&ts));
    101:   PetscCall(PetscFinalize());
    102:   return 0;
    103: }

    105: PetscErrorCode PostStep(TS ts)
    106: {
    107:   Vec          X, gc;
    108:   AppCtx      *user;
    109:   PetscScalar  sum = 0, asum;
    110:   PetscReal    t, **p;
    111:   DMDACoor2d **coors;
    112:   DM           cda;
    113:   PetscInt     i, j, xs, ys, xm, ym;

    115:   PetscFunctionBegin;
    116:   PetscCall(TSGetApplicationContext(ts, &user));
    117:   PetscCall(TSGetTime(ts, &t));
    118:   PetscCall(TSGetSolution(ts, &X));

    120:   PetscCall(DMGetCoordinateDM(user->da, &cda));
    121:   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
    122:   PetscCall(DMGetCoordinates(user->da, &gc));
    123:   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
    124:   PetscCall(DMDAVecGetArrayRead(user->da, X, &p));
    125:   for (i = xs; i < xs + xm; i++) {
    126:     for (j = ys; j < ys + ym; j++) {
    127:       if (coors[j][i].y < 5) sum += p[j][i];
    128:     }
    129:   }
    130:   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
    131:   PetscCall(DMDAVecRestoreArrayRead(user->da, X, &p));
    132:   PetscCallMPI(MPIU_Allreduce(&sum, &asum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
    133:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %f = %f\n", (double)t, (double)asum));
    134:   if (sum < 1.0e-2) {
    135:     PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER));
    136:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Exiting TS as the integral of PDF is almost zero\n"));
    137:   }
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: PetscErrorCode ini_bou(Vec X, AppCtx *user)
    142: {
    143:   DM            cda;
    144:   DMDACoor2d  **coors;
    145:   PetscScalar **p;
    146:   Vec           gc;
    147:   PetscInt      i, j;
    148:   PetscInt      xs, ys, xm, ym, M, N;
    149:   PetscScalar   xi, yi;
    150:   PetscScalar   sigmax = user->sigmax, sigmay = user->sigmay;
    151:   PetscScalar   rho = user->rho;
    152:   PetscScalar   muy = user->muy, mux;
    153:   PetscMPIInt   rank;
    154:   PetscScalar   sum;

    156:   PetscFunctionBeginUser;
    157:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    158:   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    159:   user->dx = (user->xmax - user->xmin) / (M - 1);
    160:   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));

    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;
    176:         yi = coors[j][i].y;
    177:         if ((xi == mux) && (yi == muy)) p[j][i] = 1.0;
    178:       }
    179:     }
    180:   } else {
    181:     /* Change PM_min accordingly */
    182:     user->PM_min = user->Pmax * PetscSinScalar(mux);
    183:     for (i = xs; i < xs + xm; i++) {
    184:       for (j = ys; j < ys + ym; j++) {
    185:         xi = coors[j][i].x;
    186:         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(PETSC_SUCCESS);
    196: }

    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;

    203:   PetscFunctionBegin;
    204:   f    = (y - user->ws);
    205:   fpos = PetscMax(f, 0);
    206:   fneg = PetscMin(f, 0);
    207:   if (user->st_width == 1) {
    208:     *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx;
    209:   } else if (user->st_width == 2) {
    210:     *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);
    211:   } else if (user->st_width == 3) {
    212:     *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);
    213:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: /* Second advection term */
    218: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
    219: {
    220:   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(PETSC_SUCCESS);
    234: }

    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(PETSC_SUCCESS);
    248: }

    250: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
    251: {
    252:   AppCtx       *user = (AppCtx *)ctx;
    253:   DM            cda;
    254:   DMDACoor2d  **coors;
    255:   PetscScalar **p, **f, **pdot;
    256:   PetscInt      i, j;
    257:   PetscInt      xs, ys, xm, ym, M, N;
    258:   Vec           localX, gc, localXdot;
    259:   PetscScalar   p_adv1 = 0.0, p_adv2 = 0.0, p_diff; /* initialize to prevent incorrect compiler warnings */

    261:   PetscFunctionBeginUser;
    262:   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    263:   PetscCall(DMGetCoordinateDM(user->da, &cda));
    264:   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));

    266:   PetscCall(DMGetLocalVector(user->da, &localX));
    267:   PetscCall(DMGetLocalVector(user->da, &localXdot));

    269:   PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
    270:   PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));
    271:   PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot));
    272:   PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot));

    274:   PetscCall(DMGetCoordinatesLocal(user->da, &gc));

    276:   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
    277:   PetscCall(DMDAVecGetArrayRead(user->da, localX, &p));
    278:   PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot));
    279:   PetscCall(DMDAVecGetArray(user->da, F, &f));

    281:   user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda));
    282:   for (i = xs; i < xs + xm; i++) {
    283:     for (j = ys; j < ys + ym; j++) {
    284:       PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
    285:       PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user));
    286:       PetscCall(diffuse(p, i, j, t, &p_diff, user));
    287:       f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
    288:     }
    289:   }
    290:   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
    291:   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
    292:   PetscCall(DMRestoreLocalVector(user->da, &localX));
    293:   PetscCall(DMRestoreLocalVector(user->da, &localXdot));
    294:   PetscCall(DMDAVecRestoreArray(user->da, F, &f));
    295:   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
    296:   PetscFunctionReturn(PETSC_SUCCESS);
    297: }

    299: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx)
    300: {
    301:   AppCtx      *user = (AppCtx *)ctx;
    302:   DM           cda;
    303:   DMDACoor2d **coors;
    304:   PetscInt     i, j;
    305:   PetscInt     xs, ys, xm, ym, M, N;
    306:   Vec          gc;
    307:   PetscScalar  val[5], xi, yi;
    308:   MatStencil   row, col[5];
    309:   PetscScalar  c1, c3, c5, c1pos, c1neg, c3pos, c3neg;

    311:   PetscFunctionBeginUser;
    312:   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    313:   PetscCall(DMGetCoordinateDM(user->da, &cda));
    314:   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));

    316:   PetscCall(DMGetCoordinatesLocal(user->da, &gc));
    317:   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
    318:   for (i = xs; i < xs + xm; i++) {
    319:     for (j = ys; j < ys + ym; j++) {
    320:       PetscInt nc = 0;
    321:       xi          = coors[j][i].x;
    322:       yi          = coors[j][i].y;
    323:       row.i       = i;
    324:       row.j       = j;
    325:       c1          = (yi - user->ws) / user->dx;
    326:       c1pos       = PetscMax(c1, 0);
    327:       c1neg       = PetscMin(c1, 0);
    328:       c3          = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / user->dy;
    329:       c3pos       = PetscMax(c3, 0);
    330:       c3neg       = PetscMin(c3, 0);
    331:       c5          = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
    332:       col[nc].i   = i - 1;
    333:       col[nc].j   = j;
    334:       val[nc++]   = c1pos;
    335:       col[nc].i   = i + 1;
    336:       col[nc].j   = j;
    337:       val[nc++]   = -c1neg;
    338:       col[nc].i   = i;
    339:       col[nc].j   = j - 1;
    340:       val[nc++]   = c3pos + c5;
    341:       col[nc].i   = i;
    342:       col[nc].j   = j + 1;
    343:       val[nc++]   = -c3neg + c5;
    344:       col[nc].i   = i;
    345:       col[nc].j   = j;
    346:       val[nc++]   = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
    347:       PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
    348:     }
    349:   }
    350:   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));

    352:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    353:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    354:   if (J != Jpre) {
    355:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    356:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    357:   }
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: PetscErrorCode Parameter_settings(AppCtx *user)
    362: {
    363:   PetscBool flg;

    365:   PetscFunctionBeginUser;
    366:   /* Set default parameters */
    367:   user->ws             = 1.0;
    368:   user->H              = 5.0;
    369:   user->Pmax           = 2.1;
    370:   user->PM_min         = 1.0;
    371:   user->lambda         = 0.1;
    372:   user->q              = 1.0;
    373:   user->mux            = PetscAsinScalar(user->PM_min / user->Pmax);
    374:   user->sigmax         = 0.1;
    375:   user->sigmay         = 0.1;
    376:   user->rho            = 0.0;
    377:   user->xmin           = -PETSC_PI;
    378:   user->xmax           = PETSC_PI;
    379:   user->bx             = DM_BOUNDARY_PERIODIC;
    380:   user->by             = DM_BOUNDARY_MIRROR;
    381:   user->nonoiseinitial = PETSC_FALSE;

    383:   /*
    384:      ymin of -3 seems to let the unstable solution move up and leave a zero in its wake
    385:      with an ymin of -1 the wake is never exactly zero
    386:   */
    387:   user->ymin     = -3.0;
    388:   user->ymax     = 10.0;
    389:   user->st_width = 1;

    391:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
    392:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
    393:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
    394:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
    395:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
    396:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
    397:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
    398:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg));
    399:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
    400:   if (flg == 0) user->muy = user->ws;
    401:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg));
    402:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg));
    403:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
    404:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
    405:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
    406:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
    407:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg));
    408:   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg));
    409:   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg));
    410:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonoiseinitial", &user->nonoiseinitial, &flg));
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: /*TEST

    416:    build:
    417:       requires: !complex !single

    419:    test:
    420:       args: -ts_max_steps 2
    421:       localrunfiles: petscopt_ex7

    423:    test:
    424:       suffix: 2
    425:       args: -ts_max_steps 2 -snes_mf_operator
    426:       output_file: output/ex7_1.out
    427:       localrunfiles: petscopt_ex7
    428:       timeoutfactor: 2

    430:    test:
    431:       suffix: 3
    432:       args: -ts_max_steps 2 -snes_mf -pc_type none
    433:       output_file: output/ex7_1.out
    434:       localrunfiles: petscopt_ex7
    435:       timeoutfactor: 2

    437: TEST*/