Actual source code: ex7.c

      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;

      7:    Boundary conditions Neumman using mirror values

      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))

     12: */

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

     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;

     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);

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

     62:   PetscInitialize(&argc,&argv,"petscopt_ex7", help);if (ierr) return ierr;

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

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

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

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

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

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

    106: PetscErrorCode PostStep(TS ts)
    107: {
    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;

    118:   TSGetApplicationContext(ts,&user);
    119:   TSGetTime(ts,&t);
    120:   TSGetSolution(ts,&X);

    122:   DMGetCoordinateDM(user->da,&cda);
    123:   DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
    124:   DMGetCoordinates(user->da,&gc);
    125:   DMDAVecGetArrayRead(cda,gc,&coors);
    126:   DMDAVecGetArrayRead(user->da,X,&p);
    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:   DMDAVecRestoreArrayRead(cda,gc,&coors);
    133:   DMDAVecRestoreArrayRead(user->da,X,&p);
    134:   MPI_Allreduce(&sum,&asum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ts));
    135:   PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %f = %f\n",(double)t,(double)(asum));
    136:   if (sum  < 1.0e-2) {
    137:     TSSetConvergedReason(ts,TS_CONVERGED_USER);
    138:     PetscPrintf(PETSC_COMM_WORLD,"Exiting TS as the integral of PDF is almost zero\n");
    139:   }
    140:   return(0);
    141: }

    143: PetscErrorCode ini_bou(Vec X,AppCtx* user)
    144: {
    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;

    160:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    161:   DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    162:   user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
    163:   DMGetCoordinateDM(user->da,&cda);
    164:   DMGetCoordinates(user->da,&gc);
    165:   DMDAVecGetArray(cda,gc,&coors);
    166:   DMDAVecGetArray(user->da,X,&p);
    167:   DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);

    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:   DMDAVecRestoreArray(cda,gc,&coors);
    194:   DMDAVecRestoreArray(user->da,X,&p);
    195:   VecSum(X,&sum);
    196:   VecScale(X,1.0/sum);
    197:   return(0);
    198: }

    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;
    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:   return(0);
    216: }

    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;
    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:   return(0);
    234: }

    236: /* Diffusion term */
    237: PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
    238: {
    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:   return(0);
    248: }

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

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

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

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

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

    277:   DMDAVecGetArrayRead(cda,gc,&coors);
    278:   DMDAVecGetArrayRead(user->da,localX,&p);
    279:   DMDAVecGetArrayRead(user->da,localXdot,&pdot);
    280:   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:       adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);
    286:       adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);
    287:       diffuse(p,i,j,t,&p_diff,user);
    288:       f[j][i] = -p_adv1 - p_adv2  + p_diff - pdot[j][i];
    289:     }
    290:   }
    291:   DMDAVecRestoreArrayRead(user->da,localX,&p);
    292:   DMDAVecRestoreArrayRead(user->da,localX,&pdot);
    293:   DMRestoreLocalVector(user->da,&localX);
    294:   DMRestoreLocalVector(user->da,&localXdot);
    295:   DMDAVecRestoreArray(user->da,F,&f);
    296:   DMDAVecRestoreArrayRead(cda,gc,&coors);

    298:   return(0);
    299: }

    301: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
    302: {
    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;

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

    319:   DMGetCoordinatesLocal(user->da,&gc);
    320:   DMDAVecGetArrayRead(cda,gc,&coors);
    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:       MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);
    339:     }
    340:   }
    341:   DMDAVecRestoreArrayRead(cda,gc,&coors);

    343:    MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
    344:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
    345:   if (J != Jpre) {
    346:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    347:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    348:   }
    349:   return(0);
    350: }

    354: PetscErrorCode Parameter_settings(AppCtx *user)
    355: {
    357:   PetscBool      flg;

    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;

    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;

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

    409:   return(0);
    410: }

    413: /*TEST

    415:    build:
    416:       requires: !complex !single

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

    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

    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

    436: TEST*/