Actual source code: ex29.c

      2: /*#define HAVE_DA_HDF*/

      4: /* solve the equations for the perturbed magnetic field only */
      5: #define EQ

      7: /* turning this on causes instability?!? */
      8: /* #define UPWINDING  */

     10: static char help[] = "Hall MHD with in two dimensions with time stepping and multigrid.\n\n\
     11: -options_file ex29.options\n\
     12: other PETSc options\n\
     13: -resistivity <eta>\n\
     14: -viscosity <nu>\n\
     15: -skin_depth <d_e>\n\
     16: -larmor_radius <rho_s>\n\
     17: -contours : draw contour plots of solution\n\n";

     19: /*T
     20:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
     21:    Concepts: DA^using distributed arrays;
     22:    Concepts: multicomponent
     23:    Processors: n
     24: T*/

     26: /* ------------------------------------------------------------------------

     28:     We thank Kai Germaschewski for contributing the FormFunctionLocal()

     30:   ------------------------------------------------------------------------- */

     32: /*
     33:    Include "petscda.h" so that we can use distributed arrays (DAs).
     34:    Include "petscsnes.h" so that we can use SNES solvers.
     35:    Include "petscmg.h" to control the multigrid solvers.
     36:    Note that these automatically include:
     37:      petsc.h       - base PETSc routines   petscvec.h - vectors
     38:      petscsys.h    - system routines       petscmat.h - matrices
     39:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     40:      petscviewer.h - viewers               petscpc.h  - preconditioners
     41:      petscksp.h   - linear solvers
     42: */
     43:  #include petscsnes.h
     44:  #include petscda.h
     45:  #include petscmg.h
     46:  #include petscdmmg.h

     48: #ifdef HAVE_DA_HDF
     49: PetscInt DAVecHDFOutput(DA,Vec,char*);
     50: #endif

     52: #define D_x(x,m,i,j)  (p5 * (x[(j)][(i)+1].m - x[(j)][(i)-1].m) * dhx)
     53: #define D_xm(x,m,i,j) ((x[(j)][(i)].m - x[(j)][(i)-1].m) * dhx)
     54: #define D_xp(x,m,i,j) ((x[(j)][(i+1)].m - x[(j)][(i)].m) * dhx)
     55: #define D_y(x,m,i,j)  (p5 * (x[(j)+1][(i)].m - x[(j)-1][(i)].m) * dhy)
     56: #define D_ym(x,m,i,j) ((x[(j)][(i)].m - x[(j)-1][(i)].m) * dhy)
     57: #define D_yp(x,m,i,j) ((x[(j)+1][(i)].m - x[(j)][(i)].m) * dhy)
     58: #define D_xx(x,m,i,j) ((x[(j)][(i)+1].m - two*x[(j)][(i)].m + x[(j)][(i)-1].m) * hydhx * dhxdhy)
     59: #define D_yy(x,m,i,j) ((x[(j)+1][(i)].m - two*x[(j)][(i)].m + x[(j)-1][(i)].m) * hxdhy * dhxdhy)
     60: #define Lapl(x,m,i,j) (D_xx(x,m,i,j) + D_yy(x,m,i,j))
     61: #define lx            (2.*PETSC_PI)
     62: #define ly            (4.*PETSC_PI)
     63: #define sqr(a)        ((a)*(a))

     65: /*
     66:    User-defined routines and data structures
     67: */

     69: typedef struct {
     70:   PetscReal      fnorm_ini,dt_ini;
     71:   PetscReal      dt,dt_out;
     72:   PetscReal      ptime;
     73:   PetscReal      max_time;
     74:   PetscReal      fnorm_ratio;
     75:   PetscInt       ires,itstep;
     76:   PetscInt       max_steps,print_freq;
     77:   PetscReal      t;
     78:   PetscScalar    fnorm;

     80:   PetscTruth     ts_monitor;           /* print information about each time step */
     81:   PetscReal      dump_time;            /* time to dump solution to a file */
     82:   PetscViewer    socketviewer;         /* socket to dump solution at each timestep for visualization */
     83: } TstepCtx;

     85: typedef struct {
     86:   PetscScalar phi,psi,U,F;
     87: } Field;

     89: typedef struct {
     90:   PassiveScalar phi,psi,U,F;
     91: } PassiveField;

     93: typedef struct {
     94:   PetscInt     mglevels;
     95:   PetscInt     cycles;           /* number of time steps for integration */
     96:   PassiveReal  nu,eta,d_e,rho_s; /* physical parameters */
     97:   PetscTruth   draw_contours;    /* flag - 1 indicates drawing contours */
     98:   PetscTruth   second_order;
     99:   PetscTruth   PreLoading;
    100: } Parameter;

    102: typedef struct {
    103:   Vec          Xoldold,Xold;
    104:   TstepCtx     *tsCtx;
    105:   Parameter    *param;
    106: } AppCtx;

    119: int main(int argc,char **argv)
    120: {
    122:   DMMG           *dmmg;       /* multilevel grid structure */
    123:   AppCtx         *user;       /* user-defined work context (one for each level) */
    124:   TstepCtx       tsCtx;       /* time-step parameters (one total) */
    125:   Parameter      param;       /* physical parameters (one total) */
    126:   PetscInt       i,m,n,mx,my;
    127:   DA             da;
    128:   PetscReal      dt_ratio;
    129:   PetscInt       dfill[16] = {1,0,1,0,
    130:                               0,1,0,1,
    131:                               0,0,1,1,
    132:                               0,1,1,1};
    133:   PetscInt       ofill[16] = {1,0,0,0,
    134:                               0,1,0,0,
    135:                               1,1,1,1,
    136:                               1,1,1,1};

    138:   PetscInitialize(&argc,&argv,(char *)0,help);

    141:   PreLoadBegin(PETSC_TRUE,"SetUp");

    143:   param.PreLoading = PreLoading;
    144:     DMMGCreate(PETSC_COMM_WORLD,1,&user,&dmmg);
    145:     param.mglevels = DMMGGetLevels(dmmg);

    147:     /*
    148:       Create distributed array multigrid object (DMMG) to manage
    149:       parallel grid and vectors for principal unknowns (x) and
    150:       governing residuals (f)
    151:     */
    152:     DACreate2d(PETSC_COMM_WORLD, DA_XYPERIODIC, DA_STENCIL_STAR, -5, -5,
    153:                       PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da);

    155:     /* overwrite the default sparsity pattern toone specific for
    156:        this code's nonzero structure */
    157:     DASetBlockFills(da,dfill,ofill);

    159:     DMMGSetDM(dmmg,(DM)da);
    160:     DADestroy(da);

    162:     /* default physical parameters */
    163:     param.nu    = 0;
    164:     param.eta   = 1e-3;
    165:     param.d_e   = 0.2;
    166:     param.rho_s = 0;

    168:     PetscOptionsGetReal(PETSC_NULL, "-viscosity", &param.nu,PETSC_NULL);

    170:     PetscOptionsGetReal(PETSC_NULL, "-resistivity", &param.eta,PETSC_NULL);

    172:     PetscOptionsGetReal(PETSC_NULL, "-skin_depth", &param.d_e,PETSC_NULL);

    174:     PetscOptionsGetReal(PETSC_NULL, "-larmor_radius", &param.rho_s,PETSC_NULL);

    176:     PetscOptionsHasName(PETSC_NULL, "-contours", &param.draw_contours);

    178:     PetscOptionsHasName(PETSC_NULL, "-second_order", &param.second_order);

    180:     DASetFieldName(DMMGGetDA(dmmg), 0, "phi");

    182:     DASetFieldName(DMMGGetDA(dmmg), 1, "psi");

    184:     DASetFieldName(DMMGGetDA(dmmg), 2, "U");

    186:     DASetFieldName(DMMGGetDA(dmmg), 3, "F");

    188:     /*======================================================================*/
    189:     /* Initialize stuff related to time stepping */
    190:     /*======================================================================*/

    192:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,0,0,0,0,0,0,0,0);

    194:     tsCtx.fnorm_ini   = 0;
    195:     tsCtx.max_steps   = 1000000;
    196:     tsCtx.max_time    = 1.0e+12;
    197:     /* use for dt = min(dx,dy); multiplied by dt_ratio below */
    198:     tsCtx.dt          = PetscMin(lx/mx,ly/my);
    199:     tsCtx.fnorm_ratio = 1.0e+10;
    200:     tsCtx.t           = 0;
    201:     tsCtx.dt_out      = 10;
    202:     tsCtx.print_freq  = tsCtx.max_steps;
    203:     tsCtx.ts_monitor  = PETSC_FALSE;
    204:     tsCtx.dump_time   = -1.0;

    206:     PetscOptionsGetInt(PETSC_NULL, "-max_st", &tsCtx.max_steps,PETSC_NULL);
    207:     PetscOptionsGetReal(PETSC_NULL, "-ts_rtol", &tsCtx.fnorm_ratio,PETSC_NULL);
    208:     PetscOptionsGetInt(PETSC_NULL, "-print_freq", &tsCtx.print_freq,PETSC_NULL);

    210:     dt_ratio = 1.0;
    211:     PetscOptionsGetReal(PETSC_NULL, "-dt_ratio", &dt_ratio,PETSC_NULL);
    212:     tsCtx.dt *= dt_ratio;

    214:     PetscOptionsHasName(PETSC_NULL, "-ts_monitor", &tsCtx.ts_monitor);
    215:     PetscOptionsGetReal(PETSC_NULL, "-dump", &tsCtx.dump_time,PETSC_NULL);

    218:     tsCtx.socketviewer = 0;
    219: #if defined(PETSC_USE_SOCKET_VIEWER)
    220:     {
    221:       PetscTruth flg;
    222:       PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
    223:       if (flg && !PreLoading) {
    224:         tsCtx.ts_monitor = PETSC_TRUE;
    225:         PetscViewerSocketOpen(PETSC_COMM_WORLD,0,PETSC_DECIDE,&tsCtx.socketviewer);
    226:       }
    227:     }
    228: #endif

    230:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    231:        Create user context, set problem data, create vector data structures.
    232:        Also, compute the initial guess.
    233:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    234:     /* create application context for each level */
    235:     PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
    236:     for (i=0; i<param.mglevels; i++) {
    237:       /* create work vectors to hold previous time-step solution and
    238:          function value */
    239:       VecDuplicate(dmmg[i]->x, &user[i].Xoldold);
    240:       VecDuplicate(dmmg[i]->x, &user[i].Xold);
    241:       user[i].tsCtx = &tsCtx;
    242:       user[i].param = &param;
    243:       DMMGSetUser(dmmg,i,&user[i]);
    244:     }
    245:
    246:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    247:        Create nonlinear solver context
    248:
    249:        Process adiC(20):  AddTSTermLocal AddTSTermLocal2 FormFunctionLocal FormFunctionLocali
    250:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    251:     DMMGSetSNESLocal(dmmg, FormFunctionLocal, 0,ad_FormFunctionLocal, admf_FormFunctionLocal);
    252:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);

    254:     /* attach nullspace to each level of the preconditioner */
    255:     DMMGSetNullSpace(dmmg,PETSC_FALSE,1,CreateNullSpace);

    257:     if (PreLoading) {
    258:       PetscPrintf(PETSC_COMM_WORLD, "# viscosity = %g, resistivity = %g, "
    259:                          "skin_depth # = %g, larmor_radius # = %g\n",
    260:                          param.nu, param.eta, param.d_e, param.rho_s);
    261:       DAGetInfo(DMMGGetDA(dmmg),0,&m,&n,0,0,0,0,0,0,0,0);
    262:       PetscPrintf(PETSC_COMM_WORLD,"Problem size %D by %D\n",m,n);
    263:       PetscPrintf(PETSC_COMM_WORLD,"dx %g dy %g dt %g ratio dt/min(dx,dy) %g\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
    264:     }

    266:
    267:
    268:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    269:        Solve the nonlinear system
    270:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    271:
    272:     PreLoadStage("Solve");

    274:     if (param.draw_contours) {
    275:       VecView(((AppCtx*)DMMGGetUser(dmmg,param.mglevels-1))->Xold,PETSC_VIEWER_DRAW_WORLD);
    276:     }

    278:     Update(dmmg);

    280:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    281:        Free work space.  All PETSc objects should be destroyed when they
    282:        are no longer needed.
    283:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    284:
    285:     for (i=0; i<param.mglevels; i++) {
    286:       VecDestroy(user[i].Xoldold);
    287:       VecDestroy(user[i].Xold);
    288:     }
    289:     PetscFree(user);
    290:     DMMGDestroy(dmmg);

    292:     PreLoadEnd();
    293:
    294:   PetscFinalize();
    295:   return 0;
    296: }

    298: /* ------------------------------------------------------------------- */
    301: /* ------------------------------------------------------------------- */
    302: PetscErrorCode Gnuplot(DA da, Vec X, double mtime)
    303: {
    304:   PetscInt       i,j,xs,ys,xm,ym;
    305:   PetscInt       xints,xinte,yints,yinte;
    307:   Field          **x;
    308:   FILE           *f;
    309:   char           fname[PETSC_MAX_PATH_LEN];
    310:   PetscMPIInt    rank;

    312:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    314:   sprintf(fname, "out-%g-%d.dat", mtime, rank);

    316:   f = fopen(fname, "w");

    318:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

    320:   DAVecGetArray(da,X,&x);

    322:   xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;

    324:   for (j=yints; j<yinte; j++) {
    325:     for (i=xints; i<xinte; i++) {
    326:       PetscFPrintf(PETSC_COMM_WORLD, f,
    327:                           "%D %D %g %g %g %g %g %g\n",
    328:                           i, j, 0.0, 0.0,
    329:                           PetscAbsScalar(x[j][i].U), PetscAbsScalar(x[j][i].F),
    330:                           PetscAbsScalar(x[j][i].phi), PetscAbsScalar(x[j][i].psi));
    331:     }
    332:     PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
    333:   }
    334:   DAVecRestoreArray(da,X,&x);
    335:   fclose(f);
    336:   return 0;
    337: }

    339: /* ------------------------------------------------------------------- */
    342: /* ------------------------------------------------------------------- */
    343: PetscErrorCode Initialize(DMMG *dmmg)
    344: {
    345:   AppCtx         *appCtx = (AppCtx*)DMMGGetUser(dmmg,0);
    346:   Parameter      *param = appCtx->param;
    347:   DA             da;
    349:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
    350:   PetscReal      two = 2.0,one = 1.0;
    351:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
    352:   PetscReal      d_e,rho_s,de2,xx,yy;
    353:   Field          **x, **localx;
    354:   Vec            localX;
    355:   PetscTruth     flg;

    358:   PetscOptionsHasName(0,"-restart",&flg);
    359:   if (flg) {
    360:     PetscViewer viewer;
    361:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",PETSC_FILE_RDONLY,&viewer);
    362:     VecLoadIntoVector(viewer,dmmg[param->mglevels-1]->x);
    363:     PetscViewerDestroy(viewer);
    364:     return(0);
    365:   }

    367:   d_e   = param->d_e;
    368:   rho_s = param->rho_s;
    369:   de2   = sqr(param->d_e);

    371:   da   = (DA)(dmmg[param->mglevels-1]->dm);
    372:   DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);

    374:   dhx   = mx/lx;              dhy = my/ly;
    375:   hx    = one/dhx;             hy = one/dhy;
    376:   hxdhy = hx*dhy;           hydhx = hy*dhx;
    377:   hxhy  = hx*hy;           dhxdhy = dhx*dhy;

    379:   /*
    380:      Get local grid boundaries (for 2-dimensional DA):
    381:        xs, ys   - starting grid indices (no ghost points)
    382:        xm, ym   - widths of local grid (no ghost points)
    383:   */
    384:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

    386:   DAGetLocalVector(da,&localX);
    387:   /*
    388:      Get a pointer to vector data.
    389:        - For default PETSc vectors, VecGetArray() returns a pointer to
    390:          the data array.  Otherwise, the routine is implementation dependent.
    391:        - You MUST call VecRestoreArray() when you no longer need access to
    392:          the array.
    393:   */
    394:   DAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
    395:   DAVecGetArray(da,localX,&localx);

    397:   /*
    398:      Compute initial solution over the locally owned part of the grid
    399:   */
    400:   {
    401:     PetscReal eps = lx/ly;
    402: #if defined(PETSC_HAVE_ERF)
    403:     PetscReal pert = 1e-4;
    404: #endif
    405:     PetscReal k = 1.*eps;
    406:     PetscReal gam;

    408:     if (d_e < rho_s) d_e = rho_s;
    409:     gam = k * d_e;

    411:     for (j=ys-1; j<ys+ym+1; j++) {
    412:       yy = j * hy;
    413:       for (i=xs-1; i<xs+xm+1; i++) {
    414:         xx = i * hx;
    415: #if defined(PETSC_HAVE_ERF)
    416:         if (xx < -PETSC_PI/2) {
    417:           localx[j][i].phi = pert * gam / k * erf((xx + PETSC_PI) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
    418:         } else if (xx < PETSC_PI/2) {
    419:           localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2.0) * d_e)) * (-sin(k*yy));
    420:         } else if (xx < 3*PETSC_PI/2){
    421:           localx[j][i].phi = pert * gam / k * erf((xx - PETSC_PI) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
    422:         } else {
    423:           localx[j][i].phi = - pert * gam / k * erf((xx - 2.*PETSC_PI) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
    424:         }
    425: #endif
    426: #ifdef EQ
    427:         localx[j][i].psi = 0.;
    428: #else
    429:         localx[j][i].psi = cos(xx);
    430: #endif
    431:       }
    432:     }
    433:     for (j=ys; j<ys+ym; j++) {
    434:       for (i=xs; i<xs+xm; i++) {
    435:         x[j][i].U   = Lapl(localx,phi,i,j);
    436:         x[j][i].F   = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
    437:         x[j][i].phi = localx[j][i].phi;
    438:         x[j][i].psi = localx[j][i].psi;
    439:       }
    440:     }
    441:   }
    442:   /*
    443:      Restore vector
    444:   */
    445:   DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
    446:
    447:   DAVecRestoreArray(da,localX,&localx);
    448:
    449:   DARestoreLocalVector(da,&localX);
    450:

    452:   return(0);
    453: }

    457: PetscErrorCode ComputeMaxima(DA da, Vec X, PetscReal t)
    458: {
    460:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
    461:   PetscInt       xints,xinte,yints,yinte;
    462:   Field          **x;
    463:   double         norm[4] = {0,0,0,0};
    464:   double         gnorm[4];
    465:   MPI_Comm       comm;

    469:   DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);

    471:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
    472:

    474:   xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;

    476:   DAVecGetArray(da, X, &x);

    478:   for (j=yints; j<yinte; j++) {
    479:     for (i=xints; i<xinte; i++) {
    480:       norm[0] = PetscMax(norm[0],PetscAbsScalar(x[j][i].U));
    481:       norm[1] = PetscMax(norm[1],PetscAbsScalar(x[j][i].F));
    482:       norm[2] = PetscMax(norm[2],PetscAbsScalar(x[j][i].phi));
    483:       norm[3] = PetscMax(norm[3],PetscAbsScalar(x[j][i].psi));
    484:     }
    485:   }

    487:   DAVecRestoreArray(da,X,&x);

    489:   PetscObjectGetComm((PetscObject)da, &comm);
    490:   MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);

    492:   PetscFPrintf(PETSC_COMM_WORLD, stderr,"%g\t%g\t%g\t%g\t%g\n",t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);

    494:   return(0);
    495: }

    499: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
    500: {
    501:   AppCtx         *user = (AppCtx*)ptr;
    502:   TstepCtx       *tsCtx = user->tsCtx;
    503:   Parameter      *param = user->param;
    505:   PetscInt       xints,xinte,yints,yinte,i,j;
    506:   PassiveReal    hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
    507:   PassiveReal    de2,rhos2,nu,eta,dde2;
    508:   PassiveReal    two = 2.0,one = 1.0,p5 = 0.5;
    509:   PassiveReal    F_eq_x,By_eq;
    510:   PetscScalar    xx;
    511:   PetscScalar    vx,vy,avx,avy,vxp,vxm,vyp,vym;
    512:   PetscScalar    Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;

    515:   de2     = sqr(param->d_e);
    516:   rhos2   = sqr(param->rho_s);
    517:   nu      = param->nu;
    518:   eta     = param->eta;
    519:   dde2    = one/de2;

    521:   /*
    522:      Define mesh intervals ratios for uniform grid.
    523:      [Note: FD formulae below are normalized by multiplying through by
    524:      local volume element to obtain coefficients O(1) in two dimensions.]
    525:   */
    526:   dhx   = info->mx/lx;        dhy   = info->my/ly;
    527:   hx    = one/dhx;             hy   = one/dhy;
    528:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
    529:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

    531:   xints = info->xs; xinte = info->xs+info->xm;
    532:   yints = info->ys; yinte = info->ys+info->ym;

    534:   /* Compute over the interior points */
    535:   for (j=yints; j<yinte; j++) {
    536:     for (i=xints; i<xinte; i++) {
    537: #ifdef EQ
    538:       xx = i * hx;
    539:       F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
    540:       By_eq = sin(PetscAbsScalar(xx));
    541: #else
    542:       F_eq_x = 0.;
    543:       By_eq = 0.;
    544: #endif

    546:       /*
    547:        * convective coefficients for upwinding
    548:        */

    550:       vx  = - D_y(x,phi,i,j);
    551:       vy  =   D_x(x,phi,i,j);
    552:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
    553:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
    554: #ifndef UPWINDING
    555:       vxp = vxm = p5*vx;
    556:       vyp = vym = p5*vy;
    557: #endif

    559:       Bx  =   D_y(x,psi,i,j);
    560:       By  = - D_x(x,psi,i,j) + By_eq;
    561:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
    562:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
    563: #ifndef UPWINDING
    564:       Bxp = Bxm = p5*Bx;
    565:       Byp = Bym = p5*By;
    566: #endif

    568:       /* Lap(phi) - U */
    569:       f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;

    571:       /* psi - d_e^2 * Lap(psi) - F */
    572:       f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;

    574:       /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
    575:          - nu Lap(U) */
    576:       f[j][i].U  =
    577:         ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
    578:           vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
    579:          (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
    580:           Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
    581:          nu * Lapl(x,U,i,j)) * hxhy;
    582:
    583:       /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
    584:          - eta * Lap(psi) */
    585:       f[j][i].F  =
    586:         ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
    587:           vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
    588:          (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
    589:           Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
    590:          eta * Lapl(x,psi,i,j)) * hxhy;
    591:     }
    592:   }

    594:   /* Add time step contribution */
    595:   if (tsCtx->ires) {
    596:     if ((param->second_order) && (tsCtx->itstep > 0)){
    597:       AddTSTermLocal2(info,x,f,user);
    598:     } else {
    599:       AddTSTermLocal(info,x,f,user);
    600:     }
    601:   }

    603:   /*
    604:      Flop count (multiply-adds are counted as 2 operations)
    605:   */
    606:   /*  PetscLogFlops(84*info->ym*info->xm); FIXME */
    607:   return(0);
    608: }

    610: /*---------------------------------------------------------------------*/
    613: PetscErrorCode Update(DMMG *dmmg)
    614: /*---------------------------------------------------------------------*/
    615: {
    616:   AppCtx          *user = (AppCtx *) DMMGGetUser(dmmg,0);
    617:   TstepCtx        *tsCtx = user->tsCtx;
    618:   Parameter       *param = user->param;
    619:   SNES            snes;
    620:   PetscErrorCode  ierr;
    621:   PetscInt        its,lits,i;
    622:   PetscInt        max_steps;
    623:   PetscInt        nfailsCum = 0,nfails = 0;
    624:   static PetscInt ic_out;
    625:   PetscTruth      ts_monitor = (tsCtx->ts_monitor && !param->PreLoading) ? PETSC_TRUE : PETSC_FALSE;

    629:   if (param->PreLoading)
    630:     max_steps = 1;
    631:   else
    632:     max_steps = tsCtx->max_steps;
    633:
    634:   Initialize(dmmg);

    636:   snes = DMMGGetSNES(dmmg);

    638:   for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {

    640:     if ((param->second_order) && (tsCtx->itstep > 0))
    641:     {
    642:       for (i=param->mglevels-1; i>=0 ;i--)
    643:       {
    644:         VecCopy(((AppCtx*)DMMGGetUser(dmmg,i))->Xold,((AppCtx*)DMMGGetUser(dmmg,i))->Xoldold);
    645:       }
    646:     }

    648:     for (i=param->mglevels-1; i>0 ;i--) {
    649:       MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
    650:       VecPointwiseMult(dmmg[i-1]->x,dmmg[i-1]->x,dmmg[i]->Rscale);
    651:       VecCopy(dmmg[i]->x, ((AppCtx*)DMMGGetUser(dmmg,i))->Xold);
    652:     }
    653:     VecCopy(dmmg[0]->x, ((AppCtx*)DMMGGetUser(dmmg,0))->Xold);

    655:     DMMGSolve(dmmg);

    659:     if (tsCtx->itstep == 665000)
    660:     {
    661:       KSP          ksp;
    662:       PC           pc;
    663:       Mat          mat, pmat;
    664:       MatStructure flag;
    665:       PetscViewer  viewer;
    666:       char         file[PETSC_MAX_PATH_LEN];

    668:       PetscStrcpy(file, "matrix");

    670:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, PETSC_FILE_CREATE, &viewer);

    672:       SNESGetKSP(snes, &ksp);

    674:       KSPGetPC(ksp, &pc);

    676:       PCGetOperators(pc, &mat, &pmat, &flag);

    678:       MatView(mat, viewer);

    680:       PetscViewerDestroy(viewer);
    681:       SETERRQ(1,"Done saving Jacobian");
    682:     }

    685:     tsCtx->t += tsCtx->dt;

    687:     /* save restart solution if requested at a particular time, then exit */
    688:     if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
    689:       Vec v = DMMGGetx(dmmg);
    690:       VecView(v,PETSC_VIEWER_BINARY_WORLD);
    691:       SETERRQ1(1,"Saved solution at time %g",tsCtx->t);
    692:     }

    694:     if (ts_monitor)
    695:     {
    696:       SNESGetIterationNumber(snes, &its);
    697:       SNESGetNumberLinearIterations(snes, &lits);
    698:       SNESGetNumberUnsuccessfulSteps(snes, &nfails);
    699:       SNESGetFunctionNorm(snes, &tsCtx->fnorm);

    701:       nfailsCum += nfails;
    702:       if (nfailsCum >= 2)
    703:         SETERRQ(1, "unable to find a newton step");

    705:       PetscPrintf(PETSC_COMM_WORLD,
    706:                          "time step = %D, time = %g, number of nonlinear steps = %D, "
    707:                          "number of linear steps = %D, norm of the function = %g\n",
    708:                          tsCtx->itstep + 1, tsCtx->t, its, lits, PetscAbsScalar(tsCtx->fnorm));

    710:       /* send solution over to Matlab, to be visualized (using ex29.m) */
    711:       if (!param->PreLoading && tsCtx->socketviewer)
    712:       {
    713:         Vec v;
    714:         SNESGetSolution(snes, &v);
    715: #if defined(PETSC_USE_SOCKET_VIEWER)
    716:         VecView(v, tsCtx->socketviewer);
    717: #endif
    718:       }
    719:     }

    721:     if (!param->PreLoading) {
    722:       if (param->draw_contours) {
    723:         VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
    724:       }

    726:       if (1 && ts_monitor) {
    727:         /* compute maxima */
    728:         ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,dmmg[param->mglevels-1]->x,tsCtx->t);
    729:         /* output */
    730:         if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
    731: #ifdef HAVE_DA_HDF
    732:           char fname[PETSC_MAX_PATH_LEN];

    734:           sprintf(fname, "out-%g.hdf", tsCtx->t);
    735:           DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
    736: #else
    737: /*
    738:           Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
    739:
    740: */
    741: #endif
    742:           ic_out = 1;
    743:         }
    744:       }
    745:     }
    746:   } /* End of time step loop */
    747:
    748:   if (!param->PreLoading){
    749:     SNESGetFunctionNorm(snes,&tsCtx->fnorm);
    750:     PetscPrintf(PETSC_COMM_WORLD, "timesteps %D fnorm = %g\n",
    751:                        tsCtx->itstep, PetscAbsScalar(tsCtx->fnorm));
    752:   }

    754:   if (param->PreLoading) {
    755:     tsCtx->fnorm_ini = 0.0;
    756:   }
    757:
    758:   return(0);
    759: }

    761: /*---------------------------------------------------------------------*/
    764: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
    765: /*---------------------------------------------------------------------*/
    766: {
    767:   TstepCtx       *tsCtx = user->tsCtx;
    768:   DA             da = info->da;
    770:   PetscInt       i,j;
    771:   PetscInt       xints,xinte,yints,yinte;
    772:   PassiveReal    hx,hy,dhx,dhy,hxhy;
    773:   PassiveReal    one = 1.0;
    774:   PassiveScalar  dtinv;
    775:   PassiveField   **xold;

    779:   xints = info->xs; xinte = info->xs+info->xm;
    780:   yints = info->ys; yinte = info->ys+info->ym;

    782:   dhx  = info->mx/lx;            dhy = info->my/ly;
    783:   hx   = one/dhx;                 hy = one/dhy;
    784:   hxhy = hx*hy;

    786:   dtinv = hxhy/(tsCtx->dt);

    788:   DAVecGetArray(da,user->Xold,&xold);
    789:   for (j=yints; j<yinte; j++) {
    790:     for (i=xints; i<xinte; i++) {
    791:       f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
    792:       f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
    793:     }
    794:   }
    795:   DAVecRestoreArray(da,user->Xold,&xold);

    797:   return(0);
    798: }

    800: /*---------------------------------------------------------------------*/
    803: PetscErrorCode AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
    804: /*---------------------------------------------------------------------*/
    805: {
    806:   TstepCtx       *tsCtx = user->tsCtx;
    807:   DA             da = info->da;
    809:   PetscInt       i,j,xints,xinte,yints,yinte;
    810:   PassiveReal    hx,hy,dhx,dhy,hxhy;
    811:   PassiveReal    one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
    812:   PassiveScalar  dtinv;
    813:   PassiveField   **xoldold,**xold;

    817:   xints = info->xs; xinte = info->xs+info->xm;
    818:   yints = info->ys; yinte = info->ys+info->ym;

    820:   dhx  = info->mx/lx;            dhy = info->my/ly;
    821:   hx   = one/dhx;                 hy = one/dhy;
    822:   hxhy = hx*hy;

    824:   dtinv = hxhy/(tsCtx->dt);

    826:   DAVecGetArray(da,user->Xoldold,&xoldold);
    827:   DAVecGetArray(da,user->Xold,&xold);
    828:   for (j=yints; j<yinte; j++) {
    829:     for (i=xints; i<xinte; i++) {
    830:       f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
    831:                             p5 * xoldold[j][i].U);
    832:       f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
    833:                             p5 * xoldold[j][i].F);
    834:     }
    835:   }
    836:   DAVecRestoreArray(da,user->Xoldold,&xoldold);
    837:   DAVecRestoreArray(da,user->Xold,&xold);

    839:   return(0);
    840: }

    842: /* Creates the null space of the Jacobian for a particular level */
    845: PetscErrorCode CreateNullSpace(DMMG dmmg,Vec *nulls)
    846: {
    848:   PetscInt       i,N,rstart,rend;
    849:   PetscScalar    scale,*vx;

    852:   VecGetSize(nulls[0],&N);
    853:   scale = 2.0/sqrt((PetscReal)N);
    854:   VecGetArray(nulls[0],&vx);
    855:   VecGetOwnershipRange(nulls[0],&rstart,&rend);
    856:   for (i=rstart; i<rend; i++) {
    857:     if (!(i % 4)) vx[i-rstart] = scale;
    858:     else          vx[i-rstart] = 0.0;
    859:   }
    860:   VecRestoreArray(nulls[0],&vx);
    861:   return(0);
    862: }

    864: /*
    865:     This is an experimental function and can be safely ignored.
    866: */
    867: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
    868: {
    869:   AppCtx         *user = (AppCtx*)ptr;
    870:   TstepCtx       *tsCtx = user->tsCtx;
    871:   Parameter      *param = user->param;
    873:   PetscInt       i,j,c;
    874:   PetscInt       xints,xinte,yints,yinte;
    875:   PassiveReal    hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
    876:   PassiveReal    de2,rhos2,nu,eta,dde2;
    877:   PassiveReal    two = 2.0,one = 1.0,p5 = 0.5;
    878:   PassiveReal    F_eq_x,By_eq,dtinv;
    879:   PetscScalar    xx;
    880:   PetscScalar    vx,vy,avx,avy,vxp,vxm,vyp,vym;
    881:   PetscScalar    Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
    882:   PassiveField   **xold;

    885:   de2     = sqr(param->d_e);
    886:   rhos2   = sqr(param->rho_s);
    887:   nu      = param->nu;
    888:   eta     = param->eta;
    889:   dde2    = one/de2;

    891:   /*
    892:      Define mesh intervals ratios for uniform grid.
    893:      [Note: FD formulae below are normalized by multiplying through by
    894:      local volume element to obtain coefficients O(1) in two dimensions.]
    895:   */
    896:   dhx   = info->mx/lx;        dhy   = info->my/ly;
    897:   hx    = one/dhx;             hy   = one/dhy;
    898:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
    899:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

    901:   xints = info->xs; xinte = info->xs+info->xm;
    902:   yints = info->ys; yinte = info->ys+info->ym;

    905:   i = st->i; j = st->j; c = st->c;

    907: #ifdef EQ
    908:       xx = i * hx;
    909:       F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
    910:       By_eq = sin(PetscAbsScalar(xx));
    911: #else
    912:       F_eq_x = 0.;
    913:       By_eq = 0.;
    914: #endif

    916:       /*
    917:        * convective coefficients for upwinding
    918:        */

    920:       vx  = - D_y(x,phi,i,j);
    921:       vy  =   D_x(x,phi,i,j);
    922:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
    923:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
    924: #ifndef UPWINDING
    925:       vxp = vxm = p5*vx;
    926:       vyp = vym = p5*vy;
    927: #endif

    929:       Bx  =   D_y(x,psi,i,j);
    930:       By  = - D_x(x,psi,i,j) + By_eq;
    931:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
    932:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
    933: #ifndef UPWINDING
    934:       Bxp = Bxm = p5*Bx;
    935:       Byp = Bym = p5*By;
    936: #endif

    938:       DAVecGetArray(info->da,user->Xold,&xold);
    939:       dtinv = hxhy/(tsCtx->dt);
    940:       switch(c) {

    942:         case 0:
    943:           /* Lap(phi) - U */
    944:           *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
    945:           break;

    947:         case 1:
    948:           /* psi - d_e^2 * Lap(psi) - F */
    949:           *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
    950:           break;

    952:         case 2:
    953:           /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
    954:             - nu Lap(U) */
    955:           *f  =
    956:             ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
    957:               vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
    958:              (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
    959:               Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
    960:              nu * Lapl(x,U,i,j)) * hxhy;
    961:           *f += dtinv*(x[j][i].U-xold[j][i].U);
    962:           break;

    964:         case 3:
    965:           /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
    966:            - eta * Lap(psi) */
    967:           *f  =
    968:             ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
    969:              vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
    970:             (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
    971:              Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
    972:             eta * Lapl(x,psi,i,j)) * hxhy;
    973:           *f += dtinv*(x[j][i].F-xold[j][i].F);
    974:           break;
    975:       }
    976:       DAVecRestoreArray(info->da,user->Xold,&xold);

    979:   /*
    980:      Flop count (multiply-adds are counted as 2 operations)
    981:   */
    982:   /*  PetscLogFlops(84*info->ym*info->xm); FIXME */
    983:   return(0);
    984: }