Actual source code: ex6.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;
      5:    x_t = (y - ws)  y_t = (ws/2H)*(Pm - Pmax*sin(x))

      7:    Boundary conditions: -bc_type 0 => Zero dirichlet boundary
      8:                         -bc_type 1 => Steady state boundary condition
      9:    Steady state boundary condition found by setting p_t = 0
     10: */

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

     16: /*
     17:    User-defined data structures and routines
     18: */
     19: typedef struct {
     20:   PetscScalar ws;         /* Synchronous speed */
     21:   PetscScalar H;          /* Inertia constant */
     22:   PetscScalar D;          /* Damping constant */
     23:   PetscScalar Pmax;       /* Maximum power output of generator */
     24:   PetscScalar PM_min;     /* Mean mechanical power input */
     25:   PetscScalar lambda;     /* correlation time */
     26:   PetscScalar q;          /* noise strength */
     27:   PetscScalar mux;        /* Initial average angle */
     28:   PetscScalar sigmax;     /* Standard deviation of initial angle */
     29:   PetscScalar muy;        /* Average speed */
     30:   PetscScalar sigmay;     /* standard deviation of initial speed */
     31:   PetscScalar rho;        /* Cross-correlation coefficient */
     32:   PetscScalar t0;         /* Initial time */
     33:   PetscScalar tmax;       /* Final time */
     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:   PetscInt    bc;         /* Boundary conditions */
     41:   PetscScalar disper_coe; /* Dispersion coefficient */
     42:   DM          da;
     43: } AppCtx;

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

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

     59:   PetscFunctionBeginUser;
     60:   PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex6", help));
     61:   /* Get physics and time parameters */
     62:   PetscCall(Parameter_settings(&user));
     63:   /* Create a 2D DA with dof = 1 */
     64:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.da));
     65:   PetscCall(DMSetFromOptions(user.da));
     66:   PetscCall(DMSetUp(user.da));
     67:   /* Set x and y coordinates */
     68:   PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0.0, 1.0));

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

     73:   PetscCall(ini_bou(x, &user));
     74:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer));
     75:   PetscCall(VecView(x, viewer));
     76:   PetscCall(PetscViewerDestroy(&viewer));

     78:   /* Get Jacobian matrix nonzero structure from the da */
     79:   PetscCall(DMSetMatType(user.da, MATAIJ));
     80:   PetscCall(DMCreateMatrix(user.da, &J));

     82:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     83:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
     84:   PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
     85:   PetscCall(TSSetIJacobian(ts, J, J, IJacobian, &user));
     86:   PetscCall(TSSetApplicationContext(ts, &user));
     87:   PetscCall(TSSetMaxTime(ts, user.tmax));
     88:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
     89:   PetscCall(TSSetTime(ts, user.t0));
     90:   PetscCall(TSSetTimeStep(ts, .005));
     91:   PetscCall(TSSetFromOptions(ts));
     92:   PetscCall(TSSetPostStep(ts, PostStep));
     93:   PetscCall(TSSolve(ts, x));

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

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

    107: PetscErrorCode PostStep(TS ts)
    108: {
    109:   Vec         X;
    110:   AppCtx     *user;
    111:   PetscScalar sum;
    112:   PetscReal   t;

    114:   PetscFunctionBegin;
    115:   PetscCall(TSGetApplicationContext(ts, &user));
    116:   PetscCall(TSGetTime(ts, &t));
    117:   PetscCall(TSGetSolution(ts, &X));
    118:   PetscCall(VecSum(X, &sum));
    119:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %3.2f = %3.6f\n", (double)t, (double)(sum * user->dx * user->dy)));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: PetscErrorCode ini_bou(Vec X, AppCtx *user)
    124: {
    125:   DM            cda;
    126:   DMDACoor2d  **coors;
    127:   PetscScalar **p;
    128:   Vec           gc;
    129:   PetscInt      i, j;
    130:   PetscInt      xs, ys, xm, ym, M, N;
    131:   PetscScalar   xi, yi;
    132:   PetscScalar   sigmax = user->sigmax, sigmay = user->sigmay;
    133:   PetscScalar   rho = user->rho;
    134:   PetscScalar   mux = user->mux, muy = user->muy;
    135:   PetscMPIInt   rank;

    137:   PetscFunctionBeginUser;
    138:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    139:   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    140:   user->dx = (user->xmax - user->xmin) / (M - 1);
    141:   user->dy = (user->ymax - user->ymin) / (N - 1);
    142:   PetscCall(DMGetCoordinateDM(user->da, &cda));
    143:   PetscCall(DMGetCoordinates(user->da, &gc));
    144:   PetscCall(DMDAVecGetArray(cda, gc, &coors));
    145:   PetscCall(DMDAVecGetArray(user->da, X, &p));
    146:   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
    147:   for (i = xs; i < xs + xm; i++) {
    148:     for (j = ys; j < ys + ym; j++) {
    149:       xi = coors[j][i].x;
    150:       yi = coors[j][i].y;
    151:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) p[j][i] = 0.0;
    152:       else
    153:         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)));
    154:     }
    155:   }
    156:   /*  p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */

    158:   PetscCall(DMDAVecRestoreArray(cda, gc, &coors));
    159:   PetscCall(DMDAVecRestoreArray(user->da, X, &p));
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: /* First advection term */
    164: PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
    165: {
    166:   PetscScalar f; /* v1,v2,v3,v4,v5,s1,s2,s3; */

    168:   PetscFunctionBegin;
    169:   /*  if (i > 2 && i < M-2) {
    170:     v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx;
    171:     v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx;
    172:     v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx;
    173:     v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx;
    174:     v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx;

    176:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
    177:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
    178:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

    180:     *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3;
    181:     } else *p1 = 0.0; */
    182:   f   = (y - user->ws);
    183:   *p1 = f * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx);
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: /* Second advection term */
    188: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
    189: {
    190:   PetscScalar f; /* v1,v2,v3,v4,v5,s1,s2,s3; */

    192:   PetscFunctionBegin;
    193:   /*  if (j > 2 && j < N-2) {
    194:     v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy;
    195:     v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy;
    196:     v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy;
    197:     v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy;
    198:     v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy;

    200:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
    201:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
    202:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

    204:     *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3;
    205:     } else *p2 = 0.0; */
    206:   f   = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x));
    207:   *p2 = f * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy);
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /* Diffusion term */
    212: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
    213: {
    214:   PetscFunctionBeginUser;
    215:   *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }

    219: PetscErrorCode BoundaryConditions(PetscScalar **p, DMDACoor2d **coors, PetscInt i, PetscInt j, PetscInt M, PetscInt N, PetscScalar **f, AppCtx *user)
    220: {
    221:   PetscScalar fwc, fthetac;
    222:   PetscScalar w = coors[j][i].y, theta = coors[j][i].x;

    224:   PetscFunctionBeginUser;
    225:   if (user->bc == 0) { /* Natural boundary condition */
    226:     f[j][i] = p[j][i];
    227:   } else { /* Steady state boundary condition */
    228:     fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(theta));
    229:     fwc     = (w * w / 2.0 - user->ws * w);
    230:     if (i == 0 && j == 0) { /* left bottom corner */
    231:       f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / user->dy;
    232:     } else if (i == 0 && j == N - 1) { /* right bottom corner */
    233:       f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / user->dy;
    234:     } else if (i == M - 1 && j == 0) { /* left top corner */
    235:       f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / user->dy;
    236:     } else if (i == M - 1 && j == N - 1) { /* right top corner */
    237:       f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / user->dx + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / user->dy;
    238:     } else if (i == 0) { /* Bottom edge */
    239:       f[j][i] = fwc * (p[j][i + 1] - p[j][i]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy);
    240:     } else if (i == M - 1) { /* Top edge */
    241:       f[j][i] = fwc * (p[j][i] - p[j][i - 1]) / (user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j - 1][i]) / (2 * user->dy);
    242:     } else if (j == 0) { /* Left edge */
    243:       f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j + 1][i] - p[j][i]) / (user->dy);
    244:     } else if (j == N - 1) { /* Right edge */
    245:       f[j][i] = fwc * (p[j][i + 1] - p[j][i - 1]) / (2 * user->dx) + fthetac * p[j][i] - user->disper_coe * (p[j][i] - p[j - 1][i]) / (user->dy);
    246:     }
    247:   }
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }

    251: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
    252: {
    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;

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

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

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

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

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

    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:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
    286:         PetscCall(BoundaryConditions(p, coors, i, j, M, N, f, user));
    287:       } else {
    288:         PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
    289:         PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user));
    290:         PetscCall(diffuse(p, i, j, t, &p_diff, user));
    291:         f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
    292:       }
    293:     }
    294:   }
    295:   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
    296:   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
    297:   PetscCall(DMRestoreLocalVector(user->da, &localX));
    298:   PetscCall(DMRestoreLocalVector(user->da, &localXdot));
    299:   PetscCall(DMDAVecRestoreArray(user->da, F, &f));
    300:   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
    301:   PetscFunctionReturn(PETSC_SUCCESS);
    302: }

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

    316:   PetscFunctionBeginUser;
    317:   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    318:   PetscCall(DMGetCoordinateDM(user->da, &cda));
    319:   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));

    321:   PetscCall(DMGetCoordinatesLocal(user->da, &gc));
    322:   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
    323:   for (i = xs; i < xs + xm; i++) {
    324:     for (j = ys; j < ys + ym; j++) {
    325:       PetscInt nc = 0;
    326:       xi          = coors[j][i].x;
    327:       yi          = coors[j][i].y;
    328:       row.i       = i;
    329:       row.j       = j;
    330:       if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
    331:         if (user->bc == 0) {
    332:           col[nc].i = i;
    333:           col[nc].j = j;
    334:           val[nc++] = 1.0;
    335:         } else {
    336:           PetscScalar fthetac, fwc;
    337:           fthetac = user->ws / (2 * user->H) * (user->PM_min - user->Pmax * PetscSinScalar(xi));
    338:           fwc     = (yi * yi / 2.0 - user->ws * yi);
    339:           if (i == 0 && j == 0) {
    340:             col[nc].i = i + 1;
    341:             col[nc].j = j;
    342:             val[nc++] = fwc / user->dx;
    343:             col[nc].i = i;
    344:             col[nc].j = j + 1;
    345:             val[nc++] = -user->disper_coe / user->dy;
    346:             col[nc].i = i;
    347:             col[nc].j = j;
    348:             val[nc++] = -fwc / user->dx + fthetac + user->disper_coe / user->dy;
    349:           } else if (i == 0 && j == N - 1) {
    350:             col[nc].i = i + 1;
    351:             col[nc].j = j;
    352:             val[nc++] = fwc / user->dx;
    353:             col[nc].i = i;
    354:             col[nc].j = j - 1;
    355:             val[nc++] = user->disper_coe / user->dy;
    356:             col[nc].i = i;
    357:             col[nc].j = j;
    358:             val[nc++] = -fwc / user->dx + fthetac - user->disper_coe / user->dy;
    359:           } else if (i == M - 1 && j == 0) {
    360:             col[nc].i = i - 1;
    361:             col[nc].j = j;
    362:             val[nc++] = -fwc / user->dx;
    363:             col[nc].i = i;
    364:             col[nc].j = j + 1;
    365:             val[nc++] = -user->disper_coe / user->dy;
    366:             col[nc].i = i;
    367:             col[nc].j = j;
    368:             val[nc++] = fwc / user->dx + fthetac + user->disper_coe / user->dy;
    369:           } else if (i == M - 1 && j == N - 1) {
    370:             col[nc].i = i - 1;
    371:             col[nc].j = j;
    372:             val[nc++] = -fwc / user->dx;
    373:             col[nc].i = i;
    374:             col[nc].j = j - 1;
    375:             val[nc++] = user->disper_coe / user->dy;
    376:             col[nc].i = i;
    377:             col[nc].j = j;
    378:             val[nc++] = fwc / user->dx + fthetac - user->disper_coe / user->dy;
    379:           } else if (i == 0) {
    380:             col[nc].i = i + 1;
    381:             col[nc].j = j;
    382:             val[nc++] = fwc / user->dx;
    383:             col[nc].i = i;
    384:             col[nc].j = j + 1;
    385:             val[nc++] = -user->disper_coe / (2 * user->dy);
    386:             col[nc].i = i;
    387:             col[nc].j = j - 1;
    388:             val[nc++] = user->disper_coe / (2 * user->dy);
    389:             col[nc].i = i;
    390:             col[nc].j = j;
    391:             val[nc++] = -fwc / user->dx + fthetac;
    392:           } else if (i == M - 1) {
    393:             col[nc].i = i - 1;
    394:             col[nc].j = j;
    395:             val[nc++] = -fwc / user->dx;
    396:             col[nc].i = i;
    397:             col[nc].j = j + 1;
    398:             val[nc++] = -user->disper_coe / (2 * user->dy);
    399:             col[nc].i = i;
    400:             col[nc].j = j - 1;
    401:             val[nc++] = user->disper_coe / (2 * user->dy);
    402:             col[nc].i = i;
    403:             col[nc].j = j;
    404:             val[nc++] = fwc / user->dx + fthetac;
    405:           } else if (j == 0) {
    406:             col[nc].i = i + 1;
    407:             col[nc].j = j;
    408:             val[nc++] = fwc / (2 * user->dx);
    409:             col[nc].i = i - 1;
    410:             col[nc].j = j;
    411:             val[nc++] = -fwc / (2 * user->dx);
    412:             col[nc].i = i;
    413:             col[nc].j = j + 1;
    414:             val[nc++] = -user->disper_coe / user->dy;
    415:             col[nc].i = i;
    416:             col[nc].j = j;
    417:             val[nc++] = user->disper_coe / user->dy + fthetac;
    418:           } else if (j == N - 1) {
    419:             col[nc].i = i + 1;
    420:             col[nc].j = j;
    421:             val[nc++] = fwc / (2 * user->dx);
    422:             col[nc].i = i - 1;
    423:             col[nc].j = j;
    424:             val[nc++] = -fwc / (2 * user->dx);
    425:             col[nc].i = i;
    426:             col[nc].j = j - 1;
    427:             val[nc++] = user->disper_coe / user->dy;
    428:             col[nc].i = i;
    429:             col[nc].j = j;
    430:             val[nc++] = -user->disper_coe / user->dy + fthetac;
    431:           }
    432:         }
    433:       } else {
    434:         c1        = (yi - user->ws) / (2 * user->dx);
    435:         c3        = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / (2 * user->dy);
    436:         c5        = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
    437:         col[nc].i = i - 1;
    438:         col[nc].j = j;
    439:         val[nc++] = c1;
    440:         col[nc].i = i + 1;
    441:         col[nc].j = j;
    442:         val[nc++] = -c1;
    443:         col[nc].i = i;
    444:         col[nc].j = j - 1;
    445:         val[nc++] = c3 + c5;
    446:         col[nc].i = i;
    447:         col[nc].j = j + 1;
    448:         val[nc++] = -c3 + c5;
    449:         col[nc].i = i;
    450:         col[nc].j = j;
    451:         val[nc++] = -2 * c5 - a;
    452:       }
    453:       PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
    454:     }
    455:   }
    456:   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));

    458:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    459:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    460:   if (J != Jpre) {
    461:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    462:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    463:   }
    464:   PetscFunctionReturn(PETSC_SUCCESS);
    465: }

    467: PetscErrorCode Parameter_settings(AppCtx *user)
    468: {
    469:   PetscBool flg;

    471:   PetscFunctionBeginUser;
    472:   /* Set default parameters */
    473:   user->ws     = 1.0;
    474:   user->H      = 5.0;
    475:   user->Pmax   = 2.1;
    476:   user->PM_min = 1.0;
    477:   user->lambda = 0.1;
    478:   user->q      = 1.0;
    479:   user->mux    = PetscAsinScalar(user->PM_min / user->Pmax);
    480:   user->sigmax = 0.1;
    481:   user->sigmay = 0.1;
    482:   user->rho    = 0.0;
    483:   user->t0     = 0.0;
    484:   user->tmax   = 2.0;
    485:   user->xmin   = -1.0;
    486:   user->xmax   = 10.0;
    487:   user->ymin   = -1.0;
    488:   user->ymax   = 10.0;
    489:   user->bc     = 0;

    491:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
    492:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
    493:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
    494:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
    495:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
    496:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
    497:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
    498:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg));
    499:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
    500:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg));
    501:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg));
    502:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-t0", &user->t0, &flg));
    503:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-tmax", &user->tmax, &flg));
    504:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
    505:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
    506:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
    507:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
    508:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bc_type", &user->bc, &flg));
    509:   user->muy = user->ws;
    510:   PetscFunctionReturn(PETSC_SUCCESS);
    511: }

    513: /*TEST

    515:    build:
    516:       requires: !complex

    518:    test:
    519:       args: -nox -ts_max_steps 2
    520:       localrunfiles: petscopt_ex6
    521:       timeoutfactor: 4

    523: TEST*/