自学教程：C++ DMGlobalToLocalEnd函数代码示例

/*   FormFunction - Evaluates nonlinear function, F(x).   Input Parameters:.  snes - the SNES context.  x - input vector.  ctx - optional user-defined context, as set by SNESSetFunction()   Output Parameter:.  f - function vector   Note:   The user-defined context can contain any application-specific   data needed for the function evaluation.*/PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx){  ApplicationCtx *user = (ApplicationCtx*) ctx;  DM             da    = user->da;  PetscScalar    *xx,*ff,*FF,d;  PetscErrorCode ierr;  PetscInt       i,M,xs,xm;  Vec            xlocal;  PetscFunctionBeginUser;  ierr = DMGetLocalVector(da,&xlocal);CHKERRQ(ierr);  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().     By placing code between these two statements, computations can     be done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);  /*     Get pointers to vector data.       - The vector xlocal includes ghost point; the vectors x and f do         NOT include ghost points.       - Using DMDAVecGetArray() allows accessing the values using global ordering  */  ierr = DMDAVecGetArray(da,xlocal,&xx);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,f,&ff);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,user->F,&FF);CHKERRQ(ierr);  /*     Get local grid boundaries (for 1-dimensional DMDA):       xs, xm  - starting grid index, width of local grid (no ghost points)  */  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);  /*     Set function values for boundary points; define local interior grid point range:        xsi - starting interior grid index        xei - ending interior grid index  */  if (xs == 0) { /* left boundary */    ff[0] = xx[0];    xs++;xm--;  }  if (xs+xm == M) {  /* right boundary */    ff[xs+xm-1] = xx[xs+xm-1] - 1.0;    xm--;  }  /*     Compute function over locally owned part of the grid (interior points only)  */  d = 1.0/(user->h*user->h);  for (i=xs; i<xs+xm; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];  /*     Restore vectors  */  ierr = DMDAVecRestoreArray(da,xlocal,&xx);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,f,&ff);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,user->F,&FF);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&xlocal);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*   RHSJacobian - User-provided routine to compute the Jacobian of   the nonlinear right-hand-side function of the ODE.   Input Parameters:   ts - the TS context   t - current time   global_in - global input vector   dummy - optional user-defined context, as set by TSetRHSJacobian()   Output Parameters:   AA - Jacobian matrix   BB - optionally different preconditioning matrix   str - flag indicating matrix structure  Notes:  RHSJacobian computes entries for the locally owned part of the Jacobian.   - Currently, all PETSc parallel matrix formats are partitioned by     contiguous chunks of rows across the processors.   - Each processor needs to insert only elements that it owns     locally (but any non-local elements will be sent to the     appropriate processor during matrix assembly).   - Always specify global row and columns of matrix entries when     using MatSetValues().   - Here, we set all entries for a particular row at once.   - Note that MatSetValues() uses 0-based row and column numbers     in Fortran as well as in C.*/PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx){  Mat            B        = *BB;               /* Jacobian matrix */  AppCtx         *appctx  = (AppCtx*)ctx;    /* user-defined application context */  Vec            local_in = appctx->u_local;   /* local ghosted input vector */  DM             da       = appctx->da;        /* distributed array */  PetscScalar    v[3],*localptr,sc;  PetscErrorCode ierr;  PetscInt       i,mstart,mend,mstarts,mends,idx[3],is;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Get ready for local Jacobian computations     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);CHKERRQ(ierr);  /*     Get pointer to vector data  */  ierr = VecGetArray(local_in,&localptr);CHKERRQ(ierr);  /*     Get starting and ending locally owned rows of the matrix  */  ierr   = MatGetOwnershipRange(B,&mstarts,&mends);CHKERRQ(ierr);  mstart = mstarts; mend = mends;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Compute entries for the locally owned part of the Jacobian.      - Currently, all PETSc parallel matrix formats are partitioned by        contiguous chunks of rows across the processors.      - Each processor needs to insert only elements that it owns        locally (but any non-local elements will be sent to the        appropriate processor during matrix assembly).      - Here, we set all entries for a particular row at once.      - We can set matrix entries either using either        MatSetValuesLocal() or MatSetValues().     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  /*     Set matrix rows corresponding to boundary data  */  if (mstart == 0) {    v[0] = 0.0;    ierr = MatSetValues(B,1,&mstart,1,&mstart,v,INSERT_VALUES);CHKERRQ(ierr);    mstart++;  }  if (mend == appctx->m) {    mend--;    v[0] = 0.0;    ierr = MatSetValues(B,1,&mend,1,&mend,v,INSERT_VALUES);CHKERRQ(ierr);  }  /*     Set matrix rows corresponding to interior data.  We construct the     matrix one row at a time.  */  sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));  for (i=mstart; i<mend; i++) {    idx[0] = i-1; idx[1] = i; idx[2] = i+1;    is     = i - mstart + 1;    v[0]   = sc*localptr[is];    v[1]   = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]);    v[2]   = sc*localptr[is];    ierr   = MatSetValues(B,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);  }//.........这里部分代码省略.........

/*   ComputeJacobian - Evaluates Jacobian matrix.   Input Parameters:.  x - input vector.  user - user-defined application context   Output Parameters:.  jac - Jacobian matrix.  flag - flag indicating matrix structure   Notes:   Due to grid point reordering with DMDAs, we must always work   with the local grid points, and then transform them to the new   global numbering with the "ltog" mapping    We cannot work directly with the global numbers for the original   uniprocessor grid!*/PetscErrorCode ComputeJacobian(AppCtx *user,Vec X,Mat jac){  PetscErrorCode         ierr;  Vec                    localX = user->localX;   /* local vector */  const PetscInt         *ltog;                   /* local-to-global mapping */  PetscInt               i,j,row,mx,my,col[5];  PetscInt               xs,ys,xm,ym,gxs,gys,gxm,gym,grow;  PetscScalar            two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;  ISLocalToGlobalMapping ltogm;  mx = user->mx;            my = user->my;            lambda = user->param;  hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);  sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /*     Get pointer to vector data  */  ierr = VecGetArray(localX,&x);CHKERRQ(ierr);  /*     Get local grid boundaries  */  ierr = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  ierr = DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);  /*     Get the global node numbers for all local nodes, including ghost points  */  ierr = DMGetLocalToGlobalMapping(user->da,&ltogm);CHKERRQ(ierr);  ierr = ISLocalToGlobalMappingGetIndices(ltogm,&ltog);CHKERRQ(ierr);  /*     Compute entries for the locally owned part of the Jacobian.      - Currently, all PETSc parallel matrix formats are partitioned by        contiguous chunks of rows across the processors. The "grow"        parameter computed below specifies the global row number        corresponding to each local grid point.      - Each processor needs to insert only elements that it owns        locally (but any non-local elements will be sent to the        appropriate processor during matrix assembly).      - Always specify global row and columns of matrix entries.      - Here, we set all entries for a particular row at once.  */  for (j=ys; j<ys+ym; j++) {    row = (j - gys)*gxm + xs - gxs - 1;    for (i=xs; i<xs+xm; i++) {      row++;      grow = ltog[row];      /* boundary points */      if (i == 0 || j == 0 || i == mx-1 || j == my-1) {        ierr = MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);CHKERRQ(ierr);        continue;      }      /* interior grid points */      v[0] = -hxdhy; col[0] = ltog[row - gxm];      v[1] = -hydhx; col[1] = ltog[row - 1];      v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;      v[3] = -hydhx; col[3] = ltog[row + 1];      v[4] = -hxdhy; col[4] = ltog[row + gxm];      ierr = MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);CHKERRQ(ierr);    }  }  ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&ltog);CHKERRQ(ierr);  /*     Assemble matrix, using the 2-step process:       MatAssemblyBegin(), MatAssemblyEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);//.........这里部分代码省略.........

PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *ptr){  AppCtx*        user=(AppCtx*)ptr;  PetscErrorCode ierr;  PetscInt       i,j,k,kk;  PetscInt       col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;  PetscReal      one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);  PetscReal      hx,hy,hxhy,hxhx,hyhy;  PetscReal      xi,v[5];  PetscReal      ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;  PetscReal      vmiddle, vup, vdown, vleft, vright;  PetscReal      tt,f1,f2;  PetscReal      *x,*g,zero=0.0;  Vec            localX;  nx=user->nx;  ny=user->ny;  hx=two*pi/(nx+1.0);  hy=two*user->b/(ny+1.0);  hxhy=hx*hy;  hxhx=one/(hx*hx);  hyhy=one/(hy*hy);  ierr = DMGetLocalVector(user->dm,&localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = VecSet(G, zero);CHKERRQ(ierr);  /*    Get local grid boundaries  */  ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);  ierr = VecGetArray(localX,&x);CHKERRQ(ierr);  ierr = VecGetArray(G,&g);CHKERRQ(ierr);  for (i=xs; i< xs+xm; i++){    xi=(i+1)*hx;    trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */    trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */    trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */    trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */    trule5=trule1; /* L(i,j-1) */    trule6=trule2; /* U(i,j+1) */    vdown=-(trule5+trule2)*hyhy;    vleft=-hxhx*(trule2+trule4);    vright= -hxhx*(trule1+trule3);    vup=-hyhy*(trule1+trule6);    vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);    for (j=ys; j<ys+ym; j++){      row=(j-gys)*gxm + (i-gxs);       v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;       k=0;       if (j>gys){         v[k]=vdown; col[k]=row - gxm; k++;       }       if (i>gxs){         v[k]= vleft; col[k]=row - 1; k++;       }       v[k]= vmiddle; col[k]=row; k++;       if (i+1 < gxs+gxm){         v[k]= vright; col[k]=row+1; k++;       }       if (j+1 <gys+gym){         v[k]= vup; col[k] = row+gxm; k++;       }       tt=0;       for (kk=0;kk<k;kk++){         tt+=v[kk]*x[col[kk]];       }       row=(j-ys)*xm + (i-xs);       g[row]=tt;     }  }  ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);  ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->dm,&localX);CHKERRQ(ierr);  ierr = VecDot(X,G,&f1);CHKERRQ(ierr);  ierr = VecDot(user->B,X,&f2);CHKERRQ(ierr);  ierr = VecAXPY(G, one, user->B);CHKERRQ(ierr);  *fcn = f1/2.0 + f2;  ierr = PetscLogFlops((91 + 10*ym) * xm);CHKERRQ(ierr);  return 0;//.........这里部分代码省略.........

PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx){  PetscErrorCode ierr;  AppCtx         *user=(AppCtx*)ctx;  DM             da = (DM)user->da;  PetscInt       i,j,Mx,My,xs,ys,xm,ym;  PetscReal      hx,hy,sx,sy;  PetscScalar    u,uxx,uyy,**uarray,**f,**udot;  Vec            localU;  PetscFunctionBeginUser;  ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,                   PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);  hx     = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);  hy     = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);  if (user->nstencilpts == 9 && hx != hy)SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"hx must equal hy when nstencilpts = 9 for this example");  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,localU,&uarray);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,Udot,&udot);CHKERRQ(ierr);  /* Get local grid boundaries */  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      /* Boundary conditions */      if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {        if (user->boundary == 0){ /* Drichlet BC */          f[j][i] = uarray[j][i]; /* F = U */        } else {                  /* Neumann BC */          if (i == 0 && j == 0){              /* SW corner */            f[j][i] = uarray[j][i] - uarray[j+1][i+1];          } else if (i == Mx-1 && j == 0){    /* SE corner */            f[j][i] = uarray[j][i] - uarray[j+1][i-1];          } else if (i == 0 && j == My-1){    /* NW corner */            f[j][i] = uarray[j][i] - uarray[j-1][i+1];          } else if (i == Mx-1 && j == My-1){ /* NE corner */            f[j][i] = uarray[j][i] - uarray[j-1][i-1];          } else if (i == 0){                  /* Left */            f[j][i] = uarray[j][i] - uarray[j][i+1];          } else if (i == Mx-1){               /* Right */            f[j][i] = uarray[j][i] - uarray[j][i-1];          } else if (j == 0) {                 /* Bottom */            f[j][i] = uarray[j][i] - uarray[j+1][i];          } else if (j == My-1){               /* Top */            f[j][i] = uarray[j][i] - uarray[j-1][i];          }        }      } else { /* Interior */        u       = uarray[j][i];        /* 5-point stencil */        uxx     = (-2.0*u + uarray[j][i-1] + uarray[j][i+1]);        uyy     = (-2.0*u + uarray[j-1][i] + uarray[j+1][i]);        if (user->nstencilpts == 9){        /* 9-point stencil: assume hx=hy */          uxx = 2.0*uxx/3.0 + (0.5*(uarray[j-1][i-1]+uarray[j-1][i+1]+uarray[j+1][i-1]+uarray[j+1][i+1]) - 2.0*u)/6.0;          uyy = 2.0*uyy/3.0 + (0.5*(uarray[j-1][i-1]+uarray[j-1][i+1]+uarray[j+1][i-1]+uarray[j+1][i+1]) - 2.0*u)/6.0;        }        f[j][i] = udot[j][i] - (uxx*sx + uyy*sy);      }    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(da,localU,&uarray);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,Udot,&udot);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);  ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr){  AppCtx         *user = (AppCtx*)ptr;  PetscErrorCode ierr;  PetscInt       i,j,mx,my,xs,ys,xm,ym;  PetscScalar    zero = 0.0,one = 1.0;  PetscScalar    hx,hy,hxdhy,hydhx;  PetscScalar    t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;  PetscScalar    tleft,tright,beta;  PetscScalar    **x,**f;  Vec            localX;  DM             da;  PetscFunctionBeginUser;  ierr  = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr  = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  ierr  = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);  hxdhy = hx/hy;               hydhx = hy/hx;  tleft = user->tleft;         tright = user->tright;  beta  = user->beta;  /* Get ghost points */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Evaluate function */  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      t0 = x[j][i];      if (i > 0 && i < mx-1 && j > 0 && j < my-1) {        /* general interior volume */        tw = x[j][i-1];        aw = 0.5*(t0 + tw);        dw = PetscPowScalar(aw,beta);        fw = dw*(t0 - tw);        te = x[j][i+1];        ae = 0.5*(t0 + te);        de = PetscPowScalar(ae,beta);        fe = de*(te - t0);        ts = x[j-1][i];        as = 0.5*(t0 + ts);        ds = PetscPowScalar(as,beta);        fs = ds*(t0 - ts);        tn = x[j+1][i];        an = 0.5*(t0 + tn);        dn = PetscPowScalar(an,beta);        fn = dn*(tn - t0);      } else if (i == 0) {        /* left-hand boundary */        tw = tleft;        aw = 0.5*(t0 + tw);        dw = PetscPowScalar(aw,beta);        fw = dw*(t0 - tw);        te = x[j][i+1];        ae = 0.5*(t0 + te);        de = PetscPowScalar(ae,beta);        fe = de*(te - t0);        if (j > 0) {          ts = x[j-1][i];          as = 0.5*(t0 + ts);          ds = PetscPowScalar(as,beta);          fs = ds*(t0 - ts);        } else fs = zero;        if (j < my-1) {          tn = x[j+1][i];          an = 0.5*(t0 + tn);          dn = PetscPowScalar(an,beta);          fn = dn*(tn - t0);        } else fn = zero;      } else if (i == mx-1) {        /* right-hand boundary */        tw = x[j][i-1];        aw = 0.5*(t0 + tw);        dw = PetscPowScalar(aw,beta);        fw = dw*(t0 - tw);        te = tright;        ae = 0.5*(t0 + te);        de = PetscPowScalar(ae,beta);        fe = de*(te - t0);        if (j > 0) {          ts = x[j-1][i];//.........这里部分代码省略.........

PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da){  PetscErrorCode ierr;  PetscScalar    **uin,**uout;  Vec            UIN, UOUT;  PetscInt       xs,xm,*outindex;  const PetscInt *index;  PetscInt       k,i,l,n,M,cnt=0;  PetscFunctionBeginUser;  ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = DMGetGlobalVector(da,&UIN);CHKERRQ(ierr);  ierr = VecSet(UIN,0.0);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&UOUT);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(da,UIN,&uin);CHKERRQ(ierr);  ierr = ISGetIndices(act,&index);CHKERRQ(ierr);  ierr = ISGetLocalSize(act,&n);CHKERRQ(ierr);  for (k=0; k<n; k++) {    l = index[k]%5;    i = index[k]/5;    uin[i][l]=1.0;  }  printf("Number of active constraints before applying redundancy %d/n",n);  ierr = ISRestoreIndices(act,&index);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(da,UIN,&uin);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(da,UOUT,&uout);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  for (i=xs; i < xs+xm;i++) {    if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;    if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;  }  for (i=xs; i < xs+xm; i++) {    for (l=0; l < 5; l++) {      if (uout[i][l]) cnt++;    }  }  printf("Number of active constraints after applying redundancy %d/n",cnt);  ierr = PetscMalloc(cnt*sizeof(PetscInt),&outindex);CHKERRQ(ierr);  cnt  = 0;  for (i=xs; i < xs+xm;i++) {    for (l=0;l<5;l++) {      if (uout[i][l]) outindex[cnt++] = 5*(i)+l;    }  }  ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(da,UOUT,&uout);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(da,&UIN);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&UOUT);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ComputeJacobian_LS(DM dm, Vec locX, PetscInt cell, PetscScalar CellValues[], void *ctx){  User              user = (User) ctx;  Physics           phys = user->model->physics;  PetscInt          dof = phys->dof;  const PetscScalar *facegeom, *cellgeom,*x;  PetscErrorCode    ierr;  DM                dmFace, dmCell;  DM                dmGrad = user->dmGrad;  PetscInt          fStart, fEnd, face, cStart;  Vec               locGrad, locGradLimiter, Grad;  /*here the localGradLimiter refers to the gradient that has been multiplied by the limiter function.   The locGradLimiter is used to construct the uL and uR, and the locGrad is used to caculate the diffusion term*/  Vec               TempVec; /*a temperal vec for the vector restore*/  PetscFunctionBeginUser;  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);  ierr = DMGetGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);  ierr = VecDuplicate(Grad, &TempVec);CHKERRQ(ierr);  ierr = VecCopy(Grad, TempVec);CHKERRQ(ierr);  /*Backup the original vector and use it to restore the value of dmGrad,    because I do not want to change the values of the cell gradient*/  ierr = VecGetArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);  {    PetscScalar *grad;    ierr = VecGetArray(Grad,&grad);CHKERRQ(ierr);    /* Limit interior gradients. Using cell-based loop because it generalizes better to vector limiters. */      const PetscInt    *faces;      PetscInt          numFaces,f;      PetscReal         *cellPhi; /* Scalar limiter applied to each component separately */      const PetscScalar *cx;      const CellGeom    *cg;      PetscScalar       *cgrad;      PetscInt          i;      ierr = PetscMalloc(phys->dof*sizeof(PetscScalar),&cellPhi);CHKERRQ(ierr);      ierr = DMPlexGetConeSize(dm,cell,&numFaces);CHKERRQ(ierr);      ierr = DMPlexGetCone(dm,cell,&faces);CHKERRQ(ierr);      ierr = DMPlexPointLocalRead(dm,cell,x,&cx);CHKERRQ(ierr);      ierr = DMPlexPointLocalRead(dmCell,cell,cellgeom,&cg);CHKERRQ(ierr);      ierr = DMPlexPointGlobalRef(dmGrad,cell,grad,&cgrad);CHKERRQ(ierr);      /* Limiter will be minimum value over all neighbors */      for (i=0; i<dof; i++) {        cellPhi[i] = PETSC_MAX_REAL;      }      for (f=0; f<numFaces; f++) {        const PetscScalar *ncx;        const CellGeom    *ncg;        const PetscInt    *fcells;        PetscInt          face = faces[f],ncell;        PetscScalar       v[DIM];        PetscBool         ghost;        ierr = IsExteriorGhostFace(dm,face,&ghost);CHKERRQ(ierr);        if (ghost) continue;        ierr  = DMPlexGetSupport(dm,face,&fcells);CHKERRQ(ierr);        ncell = cell == fcells[0] ? fcells[1] : fcells[0];  /*The expression (x ? y : z) has the value of y if x is nonzero, z otherwise */        ierr  = DMPlexPointLocalRead(dm,ncell,x,&ncx);CHKERRQ(ierr);        ierr  = DMPlexPointLocalRead(dmCell,ncell,cellgeom,&ncg);CHKERRQ(ierr);        Waxpy2(-1, cg->centroid, ncg->centroid, v);        for (i=0; i<dof; i++) {          /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */          PetscScalar phi,flim = 0.5 * (ncx[i] - cx[i]) / Dot2(&cgrad[i*DIM],v);          phi        = (*user->LimitGrad)(flim);          cellPhi[i] = PetscMin(cellPhi[i],phi);        }      }      /* Apply limiter to gradient */      for (i=0; i<dof; i++) Scale2(cellPhi[i],&cgrad[i*DIM],&cgrad[i*DIM]);      ierr = PetscFree(cellPhi);CHKERRQ(ierr);    ierr = VecRestoreArray(Grad,&grad);CHKERRQ(ierr);  }  ierr = DMGetLocalVector(dmGrad,&locGradLimiter);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGradLimiter);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGradLimiter);CHKERRQ(ierr);  ierr = VecCopy(TempVec, Grad);CHKERRQ(ierr);/*Restore the vector*/  ierr = DMGetLocalVector(dmGrad,&locGrad);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);  ierr = VecDestroy(&TempVec);CHKERRQ(ierr);  {//.........这里部分代码省略.........

PetscErrorCode SetupJacobian(DM dm, Vec X, Mat jac, Mat B, void *ctx){  User              user = (User) ctx;  Physics           phys = user->model->physics;  PetscSection      section, globalSection;  PetscInt          cStart, cEnd, c;//  PetscInt          numCells;  PetscInt          dof = phys->dof;  PetscErrorCode    ierr;  Vec               inLocal;  PetscFunctionBegin;  PetscPrintf(PETSC_COMM_WORLD, "dof = %d/n", dof);  ierr = DMGetLocalVector(user->dm, &inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);  ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);  cEnd = user->cEndInterior;  //numCells = cEnd - cStart;  {    PetscInt         NumOfIndices;    PetscInt         indices[dof];    PetscScalar      *values;    for (c = cStart; c < cEnd; ++c) {      ierr = DMPlexGetIndex(dm, section, globalSection, c, &NumOfIndices, indices);CHKERRQ(ierr);      ierr = PetscMalloc1(NumOfIndices*NumOfIndices, &values);CHKERRQ(ierr);      ierr = PetscMemzero(values, NumOfIndices*NumOfIndices* sizeof(PetscScalar));CHKERRQ(ierr);      if (user->second_order){        ierr = ComputeJacobian_LS(dm, inLocal, c, values, user);CHKERRQ(ierr);      }else{        ierr = ComputeJacobian_Upwind(dm, inLocal, c, values, user);CHKERRQ(ierr);      }      ierr = MatSetValues(jac, NumOfIndices, indices, NumOfIndices, indices, values, INSERT_VALUES);      ierr = PetscFree(values);CHKERRQ(ierr);    }  }  ierr = DMLocalToGlobalBegin(user->dm, inLocal, INSERT_VALUES, X);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(user->dm, inLocal, INSERT_VALUES, X);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->dm, &inLocal);CHKERRQ(ierr);  ierr = MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  {    PetscViewer viewer;    char filename[256];    sprintf(filename,"matJac.m");    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,                              &viewer);CHKERRQ(ierr);    ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);    ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);    ierr = PetscViewerASCIIPrintf(viewer, "/n% -----------------------------/n");CHKERRQ(ierr);    ierr = PetscViewerASCIIPrintf(viewer, "%  Matrix Jacobian: /n% -------------------------/n");CHKERRQ(ierr);    ierr = MatView(jac, viewer);CHKERRQ(ierr);    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

int main(int argc,char **argv){  PetscErrorCode  ierr;  DM              dm;  Vec             vec,vecLocal1,vecLocal2;  PetscScalar     *a,***a1,***a2,expected;  PetscInt        startx,starty,nx,ny,i,j,d,is,js,dof0,dof1,dof2,dofTotal,stencilWidth,Nx,Ny;  DMBoundaryType  boundaryTypex,boundaryTypey;  PetscMPIInt     rank;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  dof0 = 1;  dof1 = 1;  dof2 = 1;  stencilWidth = 2;  ierr = DMStagCreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,4,4,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,DMSTAG_STENCIL_BOX,stencilWidth,NULL,NULL,&dm);CHKERRQ(ierr);  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);  ierr = DMSetUp(dm);CHKERRQ(ierr);  ierr = DMStagGetDOF(dm,&dof0,&dof1,&dof2,NULL);CHKERRQ(ierr);  dofTotal = dof0 + 2*dof1 + dof2;  ierr = DMStagGetStencilWidth(dm,&stencilWidth);CHKERRQ(ierr);  ierr = DMCreateLocalVector(dm,&vecLocal1);CHKERRQ(ierr);  ierr = VecDuplicate(vecLocal1,&vecLocal2);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dm,&vec);CHKERRQ(ierr);  ierr = VecSet(vec,1.0);CHKERRQ(ierr);  ierr = VecSet(vecLocal1,0.0);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);  ierr = DMStagGetCorners(dm,&startx,&starty,NULL,&nx,&ny,NULL,NULL,NULL,NULL);CHKERRQ(ierr);  ierr = DMStagVecGetArrayDOFRead(dm,vecLocal1,&a1);CHKERRQ(ierr);  ierr = DMStagVecGetArrayDOF(dm,vecLocal2,&a2);CHKERRQ(ierr);  for (j=starty; j<starty + ny; ++j) {    for (i=startx; i<startx + nx; ++i) {      for (d=0; d<dofTotal; ++d) {        if (a1[j][i][d] != 1.0) {          PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)/n",rank,a1[j][i][d],1.0);CHKERRQ(ierr);        }        a2[j][i][d] = 0.0;        for (js = -stencilWidth; js <= stencilWidth; ++js) {          for (is = -stencilWidth; is <= stencilWidth; ++is) {            a2[j][i][d] += a1[j+js][i+is][d];          }        }      }    }  }  ierr = DMStagVecRestoreArrayDOFRead(dm,vecLocal1,&a1);CHKERRQ(ierr);  ierr = DMStagVecRestoreArrayDOF(dm,vecLocal2,&a2);CHKERRQ(ierr);  DMLocalToGlobalBegin(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);  DMLocalToGlobalEnd(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);  /* For the all-periodic case, all values are the same . Otherwise, just check the local version */  ierr = DMStagGetBoundaryTypes(dm,&boundaryTypex,&boundaryTypey,NULL);CHKERRQ(ierr);  if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {    ierr = VecGetArray(vec,&a);CHKERRQ(ierr);    expected = 1.0; for(d=0;d<2;++d) expected *= (2*stencilWidth+1);    for (i=0; i<ny*nx*dofTotal; ++i) {      if (a[i] != expected) {        ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)/n",rank,a[i],expected);CHKERRQ(ierr);      }    }    ierr = VecRestoreArray(vec,&a);CHKERRQ(ierr);  } else {    ierr = DMStagVecGetArrayDOFRead(dm,vecLocal2,&a2);CHKERRQ(ierr);    ierr = DMStagGetGlobalSizes(dm,&Nx,&Ny,NULL);CHKERRQ(ierr);    if (stencilWidth > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Non-periodic check implemented assuming stencilWidth = 1");      for (j=starty; j<starty + ny; ++j) {        for (i=startx; i<startx + nx; ++i) {          PetscInt  dd,extra[2];          PetscBool bnd[2];          bnd[0] = (PetscBool)((i == 0 || i == Nx-1) && boundaryTypex != DM_BOUNDARY_PERIODIC);          bnd[1] = (PetscBool)((j == 0 || j == Ny-1) && boundaryTypey != DM_BOUNDARY_PERIODIC);          extra[0] = i == Nx-1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;          extra[1] = j == Ny-1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;          { /* vertices */            PetscScalar expected = 1.0;            for (dd=0;dd<2;++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);            for (d=0; d<dof0; ++d) {              if (a2[j][i][d] != expected) {                ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D)[%D] Unexpected value %g (expecting %g)/n",rank,i,j,d,a2[j][i][d],expected);CHKERRQ(ierr);              }            }          }          { /* down edges */            PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2*stencilWidth + 1);            expected *= ((bnd[0] ? 1 : 2) * stencilWidth + 1);            for (d=dof0; d<dof0+dof1; ++d) {              if (a2[j][i][d] != expected) {                ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D)[%D] Unexpected value %g (expecting %g)/n",rank,i,j,d,a2[j][i][d],expected);CHKERRQ(ierr);              }            }          }          { /* left edges */            PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2*stencilWidth + 1);            expected *= ((bnd[1] ? 1 : 2) * stencilWidth + 1);//.........这里部分代码省略.........

/*   FormJacobian - Evaluates Jacobian matrix.   Input Parameters:.  snes - SNES context.  X    - input vector.  ptr  - optional user-defined context, as set by SNESSetJacobian()   Output Parameters:.  tH    - Jacobian matrix*/PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *ptr){  AppCtx         *user;  PetscErrorCode ierr;  PetscInt       i,j,k;  PetscInt       mx, my;  MatStencil     row,col[7];  PetscScalar    hx, hy, hydhx, hxdhy;  PetscScalar    f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;  PetscScalar    hl,hr,ht,hb,hc,htl,hbr;  PetscScalar    **x, v[7];  PetscBool      assembled;  PetscInt       xs,xm,ys,ym;  Vec            localX;  DM             da;  PetscFunctionBeginUser;  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr = SNESGetApplicationContext(snes,(void**)&user);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  hx   = 1.0/(mx+1); hy=1.0/(my+1); hydhx=hy/hx; hxdhy=hx/hy;/* Set various matrix options */  ierr = MatAssembled(H,&assembled);CHKERRQ(ierr);  if (assembled) {ierr = MatZeroEntries(H);CHKERRQ(ierr);}  /* Get local vector */  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  /* Get ghost points */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,localX, &x);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  /* Compute Jacobian over the locally owned part of the mesh */  for (j=ys; j< ys+ym; j++) {    for (i=xs; i< xs+xm; i++) {      xc = x[j][i];      xlt=xrb=xl=xr=xb=xt=xc;      /* Left */      if (i==0) {        xl  = user->left[j+1];        xlt = user->left[j+2];      } else xl = x[j][i-1];      /* Bottom */      if (j==0) {        xb  =user->bottom[i+1];        xrb = user->bottom[i+2];      } else xb = x[j-1][i];      /* Right */      if (i+1 == mx) {        xr  =user->right[j+1];        xrb = user->right[j];      } else xr = x[j][i+1];      /* Top */      if (j+1==my) {        xt  =user->top[i+1];        xlt = user->top[i];      } else xt = x[j+1][i];      /* Top left */      if (i>0 && j+1<my) xlt = x[j+1][i-1];      /* Bottom right */      if (j>0 && i+1<mx) xrb = x[j-1][i+1];      d1 = (xc-xl)/hx;      d2 = (xc-xr)/hx;      d3 = (xc-xt)/hy;      d4 = (xc-xb)/hy;      d5 = (xrb-xr)/hy;      d6 = (xrb-xb)/hx;      d7 = (xlt-xl)/hy;      d8 = (xlt-xt)/hx;      f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);      f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);      f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);      f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);      f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);      f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);//.........这里部分代码省略.........

/*  FormGradient - Evaluates gradient of f.    Input Parameters:.   snes  - the SNES context.   X     - input vector.   ptr   - optional user-defined context, as set by SNESSetFunction()    Output Parameters:.   G - vector containing the newly evaluated gradient*/PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr){  AppCtx      *user;  int         ierr;  PetscInt    i,j;  PetscInt    mx, my;  PetscScalar hx,hy, hydhx, hxdhy;  PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;  PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;  PetscScalar **g, **x;  PetscInt    xs,xm,ys,ym;  Vec         localX;  DM          da;  PetscFunctionBeginUser;  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr = SNESGetApplicationContext(snes,(void**)&user);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  hx   = 1.0/(mx+1);hy=1.0/(my+1); hydhx=hy/hx; hxdhy=hx/hy;  ierr = VecSet(G,0.0);CHKERRQ(ierr);  /* Get local vector */  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  /* Get ghost points */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /* Get pointer to local vector data */  ierr = DMDAVecGetArray(da,localX, &x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,G, &g);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  /* Compute function over the locally owned part of the mesh */  for (j=ys; j < ys+ym; j++) {    for (i=xs; i< xs+xm; i++) {      xc = x[j][i];      xlt=xrb=xl=xr=xb=xt=xc;      if (i==0) { /* left side */        xl  = user->left[j+1];        xlt = user->left[j+2];      } else xl = x[j][i-1];      if (j==0) { /* bottom side */        xb  = user->bottom[i+1];        xrb = user->bottom[i+2];      } else xb = x[j-1][i];      if (i+1 == mx) { /* right side */        xr  = user->right[j+1];        xrb = user->right[j];      } else xr = x[j][i+1];      if (j+1==0+my) { /* top side */        xt  = user->top[i+1];        xlt = user->top[i];      } else xt = x[j+1][i];      if (i>0 && j+1<my) xlt = x[j+1][i-1]; /* left top side */      if (j>0 && i+1<mx) xrb = x[j-1][i+1]; /* right bottom */      d1 = (xc-xl);      d2 = (xc-xr);      d3 = (xc-xt);      d4 = (xc-xb);      d5 = (xr-xrb);      d6 = (xrb-xb);      d7 = (xlt-xl);      d8 = (xt-xlt);      df1dxc = d1*hydhx;      df2dxc = (d1*hydhx + d4*hxdhy);      df3dxc = d3*hxdhy;      df4dxc = (d2*hydhx + d3*hxdhy);      df5dxc = d2*hydhx;      df6dxc = d4*hxdhy;      d1 /= hx;      d2 /= hx;      d3 /= hy;      d4 /= hy;      d5 /= hy;      d6 /= hx;      d7 /= hy;      d8 /= hx;      f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);      f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);      f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);//.........这里部分代码省略.........

int main(int argc,char **argv){  PetscMPIInt      rank;  PetscErrorCode   ierr;  PetscInt         M = 10,N = 8,m = PETSC_DECIDE;  PetscInt         s =2,w=2,n = PETSC_DECIDE,nloc,l,i,j,kk;  PetscInt         Xs,Xm,Ys,Ym,iloc,*iglobal;  const PetscInt   *ltog;  PetscInt         *lx       = NULL,*ly = NULL;  PetscBool        testorder = PETSC_FALSE,flg;  DMBoundaryType   bx        = DM_BOUNDARY_NONE,by= DM_BOUNDARY_NONE;  DM               da;  PetscViewer      viewer;  Vec              local,global;  PetscScalar      value;  DMDAStencilType  st = DMDA_STENCIL_BOX;  AO               ao;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",300,0,400,400,&viewer);CHKERRQ(ierr);  /* Readoptions */  ierr = PetscOptionsGetInt(NULL,NULL,"-NX",&M,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-NY",&N,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-s",&s,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-w",&w,NULL);CHKERRQ(ierr);  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-xperiodic",&flg,NULL);CHKERRQ(ierr); if (flg) bx = DM_BOUNDARY_PERIODIC;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-yperiodic",&flg,NULL);CHKERRQ(ierr); if (flg) by = DM_BOUNDARY_PERIODIC;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-xghosted",&flg,NULL);CHKERRQ(ierr); if (flg) bx = DM_BOUNDARY_GHOSTED;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-yghosted",&flg,NULL);CHKERRQ(ierr); if (flg) by = DM_BOUNDARY_GHOSTED;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-star",&flg,NULL);CHKERRQ(ierr); if (flg) st = DMDA_STENCIL_STAR;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-box",&flg,NULL);CHKERRQ(ierr); if (flg) st = DMDA_STENCIL_BOX;  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-testorder",&testorder,NULL);CHKERRQ(ierr);  /*      Test putting two nodes in x and y on each processor, exact last processor      in x and y gets the rest.  */  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL);CHKERRQ(ierr);  if (flg) {    if (m == PETSC_DECIDE) SETERRQ(PETSC_COMM_WORLD,1,"Must set -m option with -distribute option");    ierr = PetscMalloc1(m,&lx);CHKERRQ(ierr);    for (i=0; i<m-1; i++) { lx[i] = 4;}    lx[m-1] = M - 4*(m-1);    if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_WORLD,1,"Must set -n option with -distribute option");    ierr = PetscMalloc1(n,&ly);CHKERRQ(ierr);    for (i=0; i<n-1; i++) { ly[i] = 2;}    ly[n-1] = N - 2*(n-1);  }  /* Create distributed array and get vectors */  ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,st,M,N,m,n,w,s,lx,ly,&da);CHKERRQ(ierr);  ierr = PetscFree(lx);CHKERRQ(ierr);  ierr = PetscFree(ly);CHKERRQ(ierr);  ierr = DMView(da,viewer);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  /* Set global vector; send ghost points to local vectors */  value = 1;  ierr = VecSet(global,value);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  /* Scale local vectors according to processor rank; pass to global vector */  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  value = rank;  ierr = VecScale(local,value);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(da,local,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da,local,INSERT_VALUES,global);CHKERRQ(ierr);  if (!testorder) { /* turn off printing when testing ordering mappings */    ierr = PetscPrintf(PETSC_COMM_WORLD,"/nGlobal Vectors:/n");CHKERRQ(ierr);    ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD,"/n/n");CHKERRQ(ierr);  }  /* Send ghost points to local vectors */  ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-local_print",&flg,NULL);CHKERRQ(ierr);  if (flg) {    PetscViewer sviewer;    ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"/nLocal Vector: processor %d/n",rank);CHKERRQ(ierr);//.........这里部分代码省略.........

/*   IFunction - Evaluates nonlinear function, F(U).   Input Parameters:.  ts - the TS context.  U - input vector.  ptr - optional user-defined context, as set by SNESSetFunction()   Output Parameter:.  F - function vector */PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr){  DM             da;  PetscErrorCode ierr;  PetscInt       i,c,Mx,xs,xm,N;  PetscReal      hx,sx,x;  PetscScalar    uxx;  PetscScalar    **u,**f,**udot;  Vec            localU;  PetscFunctionBegin;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,&N,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /*    Get pointers to vector data  */  ierr = DMDAVecGetArrayDOF(da,localU,&u);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(da,Udot,&udot);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(da,F,&f);CHKERRQ(ierr);  /*     Get local grid boundaries  */  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  /*     Compute function over the locally owned part of the grid  */  for (i=xs; i<xs+xm; i++) {    x = i*hx;    /*  diffusion term */    for (c=0; c<N; c++) {      uxx     = (-2.0*u[i][c] + u[i-1][c] + u[i+1][c])*sx;      f[i][c] = udot[i][c] - uxx;    }    /* reaction terms */    for (c=0; c<N/3; c++) {      f[i][c]   +=  500*u[i][c]*u[i][c] + 500*u[i][c]*u[i][c+1];      f[i][c+1] += -500*u[i][c]*u[i][c] + 500*u[i][c]*u[i][c+1];      f[i][c+2] -=                        500*u[i][c]*u[i][c+1];    }    /* forcing term */    f[i][0] -= 5*PetscExpScalar((1.0 - x)*(1.0 - x));  }  /*     Restore vectors  */  ierr = DMDAVecRestoreArrayDOF(da,localU,&u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(da,Udot,&udot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc,char **argv){  PetscMPIInt      rank;  PetscInt         M  = 13,s=1,dof=1;  DMDABoundaryType bx = DMDA_BOUNDARY_PERIODIC;  PetscErrorCode   ierr;  DM               da;  PetscViewer      viewer;  Vec              local,global;  PetscScalar      value;  PetscDraw        draw;  PetscBool        flg = PETSC_FALSE;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  /* Create viewers */  ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",280,480,600,200,&viewer);CHKERRQ(ierr);  ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);  ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr);  /* Readoptions */  ierr = PetscOptionsGetInt(NULL,"-M",&M,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetEnum(NULL,"-wrap",DMDABoundaryTypes,(PetscEnum*)&bx,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-dof",&dof,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-s",&s,NULL);CHKERRQ(ierr);  /* Create distributed array and get vectors */  ierr = DMDACreate1d(PETSC_COMM_WORLD,bx,M,dof,s,NULL,&da);CHKERRQ(ierr);  ierr = DMView(da,viewer);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  /* Set global vector; send ghost points to local vectors */  value = 1;  ierr  = VecSet(global,value);CHKERRQ(ierr);  ierr  = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr  = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  /* Scale local vectors according to processor rank; pass to global vector */  ierr  = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  value = rank+1;  ierr  = VecScale(local,value);CHKERRQ(ierr);  ierr  = DMLocalToGlobalBegin(da,local,INSERT_VALUES,global);CHKERRQ(ierr);  ierr  = DMLocalToGlobalEnd(da,local,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = VecView(global,viewer);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"/nGlobal Vector:/n");CHKERRQ(ierr);  ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"/n");CHKERRQ(ierr);  /* Send ghost points to local vectors */  ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-local_print",&flg,NULL);CHKERRQ(ierr);  if (flg) {    PetscViewer            sviewer;    ISLocalToGlobalMapping is;    ierr = PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD,PETSC_TRUE);CHKERRQ(ierr);    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"/nLocal Vector: processor %d/n",rank);CHKERRQ(ierr);    ierr = PetscViewerGetSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);    ierr = VecView(local,sviewer);CHKERRQ(ierr);    ierr = PetscViewerRestoreSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"/nLocal to global mapping: processor %d/n",rank);CHKERRQ(ierr);    ierr = PetscViewerGetSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);    ierr = DMGetLocalToGlobalMapping(da,&is);CHKERRQ(ierr);    ierr = ISLocalToGlobalMappingView(is,sviewer);CHKERRQ(ierr);    ierr = PetscViewerRestoreSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);  }  /* Free memory */  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

/*   IFunction - Evaluates nonlinear function, F(U).   Input Parameters:.  ts - the TS context.  U - input vector.  ptr - optional user-defined context, as set by SNESSetFunction()   Output Parameter:.  F - function vector */PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr){  AppCtx         *appctx = (AppCtx*)ptr;  DM             da;  PetscErrorCode ierr;  PetscInt       i,Mx,xs,xm;  PetscReal      hx,sx;  PetscScalar    rho,c,rhoxx,cxx,cx,rhox,kcxrhox;  Field          *u,*f,*udot;  Vec            localU;  PetscFunctionBegin;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /*     Get pointers to vector data  */  ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /*     Get local grid boundaries  */  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  if (!xs) {    f[0].rho = udot[0].rho; /* u[0].rho - 0.0; */    f[0].c   = udot[0].c; /* u[0].c   - 1.0; */    xs++;    xm--;  }  if (xs+xm == Mx) {    f[Mx-1].rho = udot[Mx-1].rho; /* u[Mx-1].rho - 1.0; */    f[Mx-1].c   = udot[Mx-1].c;  /* u[Mx-1].c   - 0.0;  */    xm--;  }  /*     Compute function over the locally owned part of the grid  */  for (i=xs; i<xs+xm; i++) {    rho   = u[i].rho;    rhoxx = (-2.0*rho + u[i-1].rho + u[i+1].rho)*sx;    c     = u[i].c;    cxx   = (-2.0*c + u[i-1].c + u[i+1].c)*sx;    if (!appctx->upwind) {      rhox    = .5*(u[i+1].rho - u[i-1].rho)/hx;      cx      = .5*(u[i+1].c - u[i-1].c)/hx;      kcxrhox = appctx->kappa*(cxx*rho + cx*rhox);    } else {      kcxrhox = appctx->kappa*((u[i+1].c - u[i].c)*u[i+1].rho - (u[i].c - u[i-1].c)*u[i].rho)*sx;    }    f[i].rho = udot[i].rho - appctx->epsilon*rhoxx + kcxrhox  - appctx->mu*PetscAbsScalar(rho)*(1.0 - rho)*PetscMax(0,PetscRealPart(c - appctx->cstar)) + appctx->beta*rho;    f[i].c   = udot[i].c - appctx->delta*cxx + appctx->lambda*c + appctx->alpha*rho*c/(appctx->gamma + c);  }  /*     Restore vectors  */  ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx){  PetscErrorCode ierr;  DM            networkdm;  UserCtx       *User=(UserCtx*)appctx;  Vec           localX;  PetscInt      e;  PetscInt      v,vStart,vEnd,vfrom,vto;  const PetscScalar *xarr;  PetscInt      offsetfrom,offsetto,goffsetfrom,goffsetto;  DMNetworkComponentGenericDataType *arr;  PetscInt      row[2],col[8];  PetscScalar   values[8];  PetscFunctionBegin;  ierr = MatZeroEntries(J);CHKERRQ(ierr);  ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);  ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);  ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr);  ierr = DMNetworkGetComponentDataArray(networkdm,&arr);CHKERRQ(ierr);  for (v=vStart; v < vEnd; v++) {    PetscInt    i,j,offsetd,key;    PetscInt    offset,goffset;    PetscScalar Vm;    PetscScalar Sbase=User->Sbase;    VERTEXDATA  bus;    PetscBool   ghostvtex;    PetscInt    numComps;    ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr);    ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr);    for (j = 0; j < numComps; j++) {      ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr);      ierr = DMNetworkGetVariableGlobalOffset(networkdm,v,&goffset);CHKERRQ(ierr);      ierr = DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);CHKERRQ(ierr);      if (key == 1) {        PetscInt       nconnedges;	const PetscInt *connedges;	bus = (VERTEXDATA)(arr+offsetd);	if (!ghostvtex) {	  /* Handle reference bus constrained dofs */	  if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {	    row[0] = goffset; row[1] = goffset+1;	    col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;	    values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;	    ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);	    break;	  }	  	  Vm = xarr[offset+1];	  	  /* Shunt injections */          row[0] = goffset; row[1] = goffset+1;          col[0] = goffset; col[1] = goffset+1;          values[0] = values[1] = values[2] = values[3] = 0.0;          if (bus->ide != PV_BUS) {            values[1] = 2.0*Vm*bus->gl/Sbase;            values[3] = -2.0*Vm*bus->bl/Sbase;          }          ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);	}	ierr = DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);CHKERRQ(ierr);	for (i=0; i < nconnedges; i++) {	  EDGEDATA       branch;	  VERTEXDATA     busf,bust;	  PetscInt       offsetfd,offsettd,keyf,keyt;          PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;          const PetscInt *cone;          PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;	  e = connedges[i];	  ierr = DMNetworkGetComponentTypeOffset(networkdm,e,0,&key,&offsetd);CHKERRQ(ierr);	  branch = (EDGEDATA)(arr+offsetd);	  if (!branch->status) continue;	  	  Gff = branch->yff[0];	  Bff = branch->yff[1];	  Gft = branch->yft[0];	  Bft = branch->yft[1];	  Gtf = branch->ytf[0];	  Btf = branch->ytf[1];	  Gtt = branch->ytt[0];	  Btt = branch->ytt[1];	  ierr = DMNetworkGetConnectedNodes(networkdm,e,&cone);CHKERRQ(ierr);	  vfrom = cone[0];	  vto   = cone[1];	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);//.........这里部分代码省略.........

PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr){  AppCtx         *user = (AppCtx*)ptr;  PetscErrorCode ierr;  PetscInt       i,j,mx,my,xs,ys,xm,ym;  PetscScalar    one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;  PetscScalar    dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;  PetscScalar    tleft,tright,beta,bm1,coef;  PetscScalar    v[5],**x;  Vec            localX;  MatStencil     col[5],row;  DM             da;  PetscFunctionBeginUser;  ierr  = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr  = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  ierr  = DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);  hxdhy = hx/hy;               hydhx  = hy/hx;  tleft = user->tleft;         tright = user->tright;  beta  = user->beta;          bm1    = user->bm1;          coef = user->coef;  /* Get ghost points */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr);  /* Evaluate Jacobian of function */  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      t0 = x[j][i];      if (i > 0 && i < mx-1 && j > 0 && j < my-1) {        /* general interior volume */        tw = x[j][i-1];        aw = 0.5*(t0 + tw);        bw = PetscPowScalar(aw,bm1);        /* dw = bw * aw */        dw = PetscPowScalar(aw,beta);        gw = coef*bw*(t0 - tw);        te = x[j][i+1];        ae = 0.5*(t0 + te);        be = PetscPowScalar(ae,bm1);        /* de = be * ae; */        de = PetscPowScalar(ae,beta);        ge = coef*be*(te - t0);        ts = x[j-1][i];        as = 0.5*(t0 + ts);        bs = PetscPowScalar(as,bm1);        /* ds = bs * as; */        ds = PetscPowScalar(as,beta);        gs = coef*bs*(t0 - ts);        tn = x[j+1][i];        an = 0.5*(t0 + tn);        bn = PetscPowScalar(an,bm1);        /* dn = bn * an; */        dn = PetscPowScalar(an,beta);        gn = coef*bn*(tn - t0);        v[0] = -hxdhy*(ds - gs);                                       col[0].j = j-1;       col[0].i = i;        v[1] = -hydhx*(dw - gw);                                       col[1].j = j;         col[1].i = i-1;        v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;        v[3] = -hydhx*(de + ge);                                       col[3].j = j;         col[3].i = i+1;        v[4] = -hxdhy*(dn + gn);                                       col[4].j = j+1;       col[4].i = i;        ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);      } else if (i == 0) {        /* left-hand boundary */        tw = tleft;        aw = 0.5*(t0 + tw);        bw = PetscPowScalar(aw,bm1);        /* dw = bw * aw */        dw = PetscPowScalar(aw,beta);        gw = coef*bw*(t0 - tw);        te = x[j][i + 1];        ae = 0.5*(t0 + te);        be = PetscPowScalar(ae,bm1);        /* de = be * ae; */        de = PetscPowScalar(ae,beta);        ge = coef*be*(te - t0);        /* left-hand bottom boundary */        if (j == 0) {          tn = x[j+1][i];          an = 0.5*(t0 + tn);          bn = PetscPowScalar(an,bm1);          /* dn = bn * an; */          dn = PetscPowScalar(an,beta);          gn = coef*bn*(tn - t0);          v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;//.........这里部分代码省略.........

PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx){  PetscErrorCode ierr;  DM             networkdm;  UserCtx       *User=(UserCtx*)appctx;  Vec           localX,localF;  PetscInt      e;  PetscInt      v,vStart,vEnd,vfrom,vto;  const PetscScalar *xarr;  PetscScalar   *farr;  PetscInt      offsetfrom,offsetto,offset;  DMNetworkComponentGenericDataType *arr;  PetscFunctionBegin;  ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);  ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr);  ierr = VecSet(F,0.0);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);  ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);  ierr = VecGetArray(localF,&farr);CHKERRQ(ierr);  ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr);  ierr = DMNetworkGetComponentDataArray(networkdm,&arr);CHKERRQ(ierr);  for (v=vStart; v < vEnd; v++) {    PetscInt    i,j,offsetd,key;    PetscScalar Vm;    PetscScalar Sbase=User->Sbase;    VERTEXDATA  bus=NULL;    GEN         gen;    LOAD        load;    PetscBool   ghostvtex;    PetscInt    numComps;    ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr);    ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr);    ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr);    for (j = 0; j < numComps; j++) {      ierr = DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);CHKERRQ(ierr);      if (key == 1) {        PetscInt       nconnedges;	const PetscInt *connedges;	bus = (VERTEXDATA)(arr+offsetd);	/* Handle reference bus constrained dofs */	if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {	  farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;	  farr[offset+1] = xarr[offset+1] - bus->vm;	  break;	}	if (!ghostvtex) {	  Vm = xarr[offset+1];	  /* Shunt injections */	  farr[offset] += Vm*Vm*bus->gl/Sbase;	  if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;	}	ierr = DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);CHKERRQ(ierr);	for (i=0; i < nconnedges; i++) {	  EDGEDATA       branch;	  PetscInt       keye;          PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;          const PetscInt *cone;          PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;	  e = connedges[i];	  ierr = DMNetworkGetComponentTypeOffset(networkdm,e,0,&keye,&offsetd);CHKERRQ(ierr);	  branch = (EDGEDATA)(arr+offsetd);	  if (!branch->status) continue;	  Gff = branch->yff[0];	  Bff = branch->yff[1];	  Gft = branch->yft[0];	  Bft = branch->yft[1];	  Gtf = branch->ytf[0];	  Btf = branch->ytf[1];	  Gtt = branch->ytt[0];	  Btt = branch->ytt[1];	  ierr = DMNetworkGetConnectedNodes(networkdm,e,&cone);CHKERRQ(ierr);	  vfrom = cone[0];	  vto   = cone[1];	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);	  thetaf = xarr[offsetfrom];	  Vmf     = xarr[offsetfrom+1];	  thetat = xarr[offsetto];	  Vmt     = xarr[offsetto+1];	  thetaft = thetaf - thetat;	  thetatf = thetat - thetaf;//.........这里部分代码省略.........

/* RHSFunction:  evaluates nonlinear function in ODE form X' = F(X,t).   But our case is autonomous, so F = F(X) = F(W,P) and F has no t dependence. */PetscErrorCode RHSFunction(TS ts,PetscReal t_unused,Vec X,Vec F,void *ptr){  PetscErrorCode ierr;  PorousCtx      *user = (PorousCtx*)ptr;  DM             da = user->da;  Vec            dPstag,localX;  PetscInt       i,j,Mx,My,xs,ys,xm,ym;  PetscReal      dx = user->dx, dy = user->dy,                 c1, c2, c3, sig, q, Ymin, Wmin,                 Wij, Weast, Wwest, Wnorth, Wsouth,                 dQx, dQy, divQ,                  dPcentx, dPcenty, WdPsqr,                 Pij, pi, Nn, Prelpow, zz;  WPnode         **wp, **f;  Hbnode         **hb;  staggradnode   **stag;  PetscFunctionBegin;  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,                     PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,                     PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  user->fcncount = user->fcncount + 1;  ierr = checkPositivity(user,X);CHKERRQ(ierr);  c1 = user->Kconst / (user->rhow * user->g);  c2 = c1 / (user->rhow * user->Lfusion);  c3 = c1 / (user->rhoi * user->Lfusion);  sig = user->sigma;  q = 1.0 / sig;  Ymin = user->Ymin;  Wmin = user->Wmin;  /* we will be differencing X = (W,P) */  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&dPstag);CHKERRQ(ierr); /* space for staggered grad */  ierr = DMDAVecGetArray(da,localX,&wp);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,dPstag,&stag);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      if (i == Mx-1) {        stag[j][i].ddxE = 0.0;  /* value will not be referenced */      } else {        stag[j][i].ddxE = (wp[j][i+1].P - wp[j][i].P) / dx;      }      if (j == My-1) {        stag[j][i].ddyN = 0.0;  /* value will not be referenced */      } else {        stag[j][i].ddyN = (wp[j+1][i].P - wp[j][i].P) / dy;      }    }  }  ierr = DMDAVecRestoreArray(da,dPstag,&stag);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,localX,&wp);CHKERRQ(ierr);  /* we will be differencing dPstag = (ddxE,ddyN) */  ierr = DMDALocalToLocalBegin(da,dPstag,INSERT_VALUES,dPstag);CHKERRQ(ierr);  ierr = DMDALocalToLocalEnd(da,dPstag,INSERT_VALUES,dPstag);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,localX,&wp);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,user->geom,&hb);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,dPstag,&stag);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {        f[j][i].W = 0.0;  /* no change at Dirichlet boundary conditions */        f[j][i].P = 0.0;      } else {        Wij     = wp[j][i].W;        Weast   = 0.5 * (wp[j][i+1].W + Wij);        Wwest   = 0.5 * (wp[j][i-1].W + Wij);        Wnorth  = 0.5 * (wp[j+1][i].W + Wij);        Wsouth  = 0.5 * (wp[j-1][i].W + Wij);        dPcentx = (wp[j][i+1].P - wp[j][i-1].P) / dx;        dPcenty = (wp[j+1][i].P - wp[j-1][i].P) / dy;        WdPsqr  = Wij * ( dPcentx * dPcentx + dPcenty * dPcenty );        /* W_t = c1 div (W grad P) + c2 W |grad P|^2 */        dQx       = Weast  * stag[j][i].ddxE - Wwest  * stag[j][i-1].ddxE;        dQy       = Wnorth * stag[j][i].ddyN - Wsouth * stag[j-1][i].ddyN;        divQ      = dQx / dx + dQy / dy;        f[j][i].W = c1 * divQ + user->Cmelt * c2 * WdPsqr;        /* P_t = sigma P/(W+Wmin) [ W_t - Cmelt c3 (P/pi)^q) W |grad P|^2                                    + Creep A max{0,N}^n (W - (P/pi)^q Ymin) ] */        pi        = user->rhoi * user->g * hb[j][i].H;        Pij       = wp[j][i].P;        Nn        = pow(PetscMax(0.0,pi - Pij),user->nglen);        Prelpow   = pow(Pij / pi,q);        zz        = f[j][i].W - user->Cmelt * c3 * Prelpow * WdPsqr                      + user->Creep * user->Aglen * Nn * (Wij - Prelpow * Ymin);        f[j][i].P = sig * (Pij / (Wij + Wmin)) * zz;//.........这里部分代码省略.........

PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ctx){  DMDALocalInfo  info;  Field          **x,**b;  PetscErrorCode ierr;  Vec            localX, localB;  DM             da;  PetscInt       xints,xinte,yints,yinte,i,j,k,l;  PetscInt       max_its,tot_its;  PetscInt       sweeps;  PetscReal      rtol,atol,stol;  PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;  PetscReal      grashof,prandtl,lid;  PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;  PetscScalar    fu, fv, fomega, ftemp;  PetscScalar    dfudu;  PetscScalar    dfvdv;  PetscScalar    dfodu, dfodv, dfodo;  PetscScalar    dftdu, dftdv, dftdt;  PetscScalar    yu=0, yv=0, yo=0, yt=0;  PetscScalar    bjiu, bjiv, bjiomega, bjitemp;  PetscBool      ptconverged;  PetscReal      pfnorm,pfnorm0,pynorm,pxnorm;  AppCtx         *user = (AppCtx*)ctx;  PetscFunctionBeginUser;  grashof = user->grashof;  prandtl = user->prandtl;  lid     = user->lidvelocity;  tot_its = 0;  ierr    = SNESNGSGetTolerances(snes,&rtol,&atol,&stol,&max_its);CHKERRQ(ierr);  ierr    = SNESNGSGetSweeps(snes,&sweeps);CHKERRQ(ierr);  ierr    = SNESGetDM(snes,(DM*)&da);CHKERRQ(ierr);  ierr    = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  if (B) {    ierr = DMGetLocalVector(da,&localB);CHKERRQ(ierr);  }  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().  */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  if (B) {    ierr = DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);CHKERRQ(ierr);  }  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr);  if (B) {    ierr = DMDAVecGetArrayRead(da,localB,&b);CHKERRQ(ierr);  }  /* looks like a combination of the formfunction / formjacobian routines */  dhx   = (PetscReal)(info.mx-1);dhy   = (PetscReal)(info.my-1);  hx    = 1.0/dhx;               hy    = 1.0/dhy;  hxdhy = hx*dhy;                hydhx = hy*dhx;  xints = info.xs; xinte = info.xs+info.xm; yints = info.ys; yinte = info.ys+info.ym;  /* Set the boundary conditions on the momentum equations */  /* Test whether we are on the bottom edge of the global array */  if (yints == 0) {    j     = 0;    yints = yints + 1;    /* bottom edge */    for (i=info.xs; i<info.xs+info.xm; i++) {      if (B) {        bjiu = b[j][i].u;        bjiv = b[j][i].v;      } else {        bjiu = 0.0;        bjiv = 0.0;      }      x[j][i].u = 0.0 + bjiu;      x[j][i].v = 0.0 + bjiv;    }  }  /* Test whether we are on the top edge of the global array */  if (yinte == info.my) {    j     = info.my - 1;    yinte = yinte - 1;    /* top edge */    for (i=info.xs; i<info.xs+info.xm; i++) {      if (B) {        bjiu = b[j][i].u;        bjiv = b[j][i].v;      } else {        bjiu = 0.0;        bjiv = 0.0;      }      x[j][i].u = lid + bjiu;      x[j][i].v = bjiv;    }  }  /* Test whether we are on the left edge of the global array */  if (xints == 0) {    i     = 0;//.........这里部分代码省略.........

/*   RHSFunction - User-provided routine that evalues the right-hand-side   function of the ODE.  This routine is set in the main program by   calling TSSetRHSFunction().  We compute:          global_out = F(global_in)   Input Parameters:   ts         - timesteping context   t          - current time   global_in  - vector containing the current iterate   ctx        - (optional) user-provided context for function evaluation.                In this case we use the appctx defined above.   Output Parameter:   global_out - vector containing the newly evaluated function*/PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx){  AppCtx         *appctx   = (AppCtx*) ctx;     /* user-defined application context */  DM             da        = appctx->da;        /* distributed array */  Vec            local_in  = appctx->u_local;   /* local ghosted input vector */  Vec            localwork = appctx->localwork; /* local ghosted work vector */  PetscErrorCode ierr;  PetscInt       i,localsize;  PetscMPIInt    rank,size;  PetscScalar    *copyptr,*localptr,sc;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Get ready for local function computations     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);CHKERRQ(ierr);  /*      Access directly the values in our local INPUT work array  */  ierr = VecGetArray(local_in,&localptr);CHKERRQ(ierr);  /*      Access directly the values in our local OUTPUT work array  */  ierr = VecGetArray(localwork,&copyptr);CHKERRQ(ierr);  sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));  /*      Evaluate our function on the nodes owned by this processor  */  ierr = VecGetLocalSize(local_in,&localsize);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Compute entries for the locally owned part     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  /*     Handle boundary conditions: This is done by using the boundary condition        u(t,boundary) = g(t,boundary)     for some function g. Now take the derivative with respect to t to obtain        u_{t}(t,boundary) = g_{t}(t,boundary)     In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1             and  u(t,1) = 2t+ 2, so that u_{t}(t,1) = 2  */  ierr = MPI_Comm_rank(appctx->comm,&rank);CHKERRQ(ierr);  ierr = MPI_Comm_size(appctx->comm,&size);CHKERRQ(ierr);  if (!rank) copyptr[0] = 1.0;  if (rank == size-1) copyptr[localsize-1] = (t < .5) ? 2.0 : 0.0;  /*     Handle the interior nodes where the PDE is replace by finite     difference operators.  */  for (i=1; i<localsize-1; i++) copyptr[i] =  localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);  /*     Restore vectors  */  ierr = VecRestoreArray(local_in,&localptr);CHKERRQ(ierr);  ierr = VecRestoreArray(localwork,&copyptr);CHKERRQ(ierr);  /*     Insert values from the local OUTPUT vector into the global     output vector  */  ierr = DMLocalToGlobalBegin(da,localwork,INSERT_VALUES,global_out);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da,localwork,INSERT_VALUES,global_out);CHKERRQ(ierr);  /* Print debugging information if desired */  /*  if (appctx->debug) {     ierr = PetscPrintf(appctx->comm,"RHS function vector/n");CHKERRQ(ierr);     ierr = VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);   } */  return 0;//.........这里部分代码省略.........

/*  FormJacobianAction - Form the global Jacobian action Y = JX from the global input X  Input Parameters:+ mat - The Jacobian shell matrix- X  - Global input vector  Output Parameter:. Y  - Local output vector  Note:  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,  like a GPU, or vectorize on a multicore machine..seealso: FormJacobianActionLocal()*/PetscErrorCode FormJacobianAction(Mat J, Vec X,  Vec Y){  JacActionCtx   *ctx;  DM             dm;  Vec            localX, localY;  PetscInt       N, n;  PetscErrorCode ierr;  PetscFunctionBeginUser;#if 0  /* Needs petscimpl.h */  PetscValidHeaderSpecific(J, MAT_CLASSID, 1);  PetscValidHeaderSpecific(X, VEC_CLASSID, 2);  PetscValidHeaderSpecific(Y, VEC_CLASSID, 3);#endif  ierr = MatShellGetContext(J, &ctx);CHKERRQ(ierr);  dm   = ctx->dm;  /* determine whether X = localX */  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &localY);CHKERRQ(ierr);  ierr = VecGetSize(X, &N);CHKERRQ(ierr);  ierr = VecGetSize(localX, &n);CHKERRQ(ierr);  if (n != N) { /* X != localX */    ierr = VecSet(localX, 0.0);CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);  } else {    ierr   = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);    localX = X;  }  ierr = DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);CHKERRQ(ierr);  if (n != N) {    ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);  }  ierr = VecSet(Y, 0.0);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &localY);CHKERRQ(ierr);  if (0) {    Vec       r;    PetscReal norm;    ierr = VecDuplicate(X, &r);CHKERRQ(ierr);    ierr = MatMult(ctx->J, X, r);CHKERRQ(ierr);    ierr = VecAXPY(r, -1.0, Y);CHKERRQ(ierr);    ierr = VecNorm(r, NORM_2, &norm);CHKERRQ(ierr);    if (norm > 1.0e-8) {      ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:/n");CHKERRQ(ierr);      ierr = VecView(X, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:/n");CHKERRQ(ierr);      ierr = VecView(Y, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "Difference:/n");CHKERRQ(ierr);      ierr = VecView(r, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);      SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);    }    ierr = VecDestroy(&r);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

    PetscErrorCode save_VTK(const char* path, const char* filename, Vec sol, DM dm, std::array<double, 2> h)    {      PetscErrorCode ierr;      DM             dau, dap;      DMDALocalInfo  infou, infop;      Vec            solu, solp;      Vec            locsolu, locsolp;      PetscScalar    ***psolu, **psolp;      ierr = DMCompositeGetEntries(dm, &dau, &dap);CHKERRQ(ierr);      ierr = DMDAGetLocalInfo(dau, &infou);CHKERRQ(ierr);      ierr = DMDAGetLocalInfo(dap, &infop);CHKERRQ(ierr);      ierr = DMCompositeGetAccess(dm,sol,&solu,&solp);CHKERRQ(ierr);      ierr = DMGetLocalVector(dau, &locsolu);CHKERRQ(ierr);      ierr = DMGetLocalVector(dap, &locsolp);CHKERRQ(ierr);      ierr = DMGlobalToLocalBegin(dau, solu, INSERT_VALUES, locsolu);CHKERRQ(ierr);      ierr = DMGlobalToLocalEnd(dau, solu, INSERT_VALUES, locsolu);CHKERRQ(ierr);      ierr = DMGlobalToLocalBegin(dap, solp, INSERT_VALUES, locsolp);CHKERRQ(ierr);      ierr = DMGlobalToLocalEnd(dap, solp, INSERT_VALUES, locsolp);CHKERRQ(ierr);      ierr = DMDAVecGetArrayDOF(dau, locsolu, &psolu);CHKERRQ(ierr);      ierr = DMDAVecGetArray(dap, locsolp, &psolp);CHKERRQ(ierr);            vtkStructuredGrid* pressureDataSet = vtkStructuredGrid::New();      pressureDataSet->SetExtent(infop.gxs, infop.gxs+infop.gxm-1, infop.gys, infop.gys+infop.gym-1, 0, 0);      vtkStructuredGrid* velocityDataSet = vtkStructuredGrid::New();      velocityDataSet->SetExtent(infou.gxs, infou.gxs+infou.gxm-1, infou.gys, infou.gys+infou.gym-1, 0, 0);            vtkPoints* pressurePoints = vtkPoints::New();      vtkPoints* velocityPoints = vtkPoints::New();            vtkDoubleArray* velocity = vtkDoubleArray::New();      velocity->SetNumberOfComponents(3);      velocity->SetName("velocity");            vtkDoubleArray* pressure = vtkDoubleArray::New();      pressure->SetName("pressure");            for (int j=infop.gys; j<infop.gys+infop.gym; j++)        for (int i=infop.gxs; i<infop.gxs+infop.gxm; i++)          {            pressurePoints->InsertNextPoint(2*i*h[0], 2*j*h[1], 0.);            pressure->InsertNextValue(psolp[j][i]);          }      int rank;      MPI_Comm_rank(PETSC_COMM_WORLD, &rank);      for (int j=infou.gys; j<infou.gys+infou.gym; j++)        for (int i=infou.gxs; i<infou.gxs+infou.gxm; i++)          {            velocityPoints->InsertNextPoint(i*h[0], j*h[1], 0.);            velocity->InsertNextTuple3(psolu[j][i][0], psolu[j][i][1], 0.);          }      pressureDataSet->SetPoints(pressurePoints);      pressureDataSet->GetPointData()->SetScalars(pressure);      velocityDataSet->SetPoints(velocityPoints);      velocityDataSet->GetPointData()->SetScalars(velocity);            vtkXMLStructuredGridWriter* pressureDataWriter = vtkXMLStructuredGridWriter::New();      vtkXMLStructuredGridWriter* velocityDataWriter = vtkXMLStructuredGridWriter::New();      std::stringstream op, ov;      int size;//, rank;      MPI_Comm_size(PETSC_COMM_WORLD, &size);      MPI_Comm_rank(PETSC_COMM_WORLD, &rank);      int x[size], y[size], z[size];      int sx[size], sy[size], sz[size];      MPI_Gather(&infop.gxs, 1, MPI_INT, x, 1, MPI_INT, 0, PETSC_COMM_WORLD);      MPI_Gather(&infop.gxm, 1, MPI_INT, sx, 1, MPI_INT, 0, PETSC_COMM_WORLD);      MPI_Gather(&infop.gys, 1, MPI_INT, y, 1, MPI_INT, 0, PETSC_COMM_WORLD);      MPI_Gather(&infop.gym, 1, MPI_INT, sy, 1, MPI_INT, 0, PETSC_COMM_WORLD);      MPI_Gather(&infop.gzs, 1, MPI_INT, z, 1, MPI_INT, 0, PETSC_COMM_WORLD);      MPI_Gather(&infop.gzm, 1, MPI_INT, sz, 1, MPI_INT, 0, PETSC_COMM_WORLD);      op << path << "/" << filename << "_pressure_" << rank << ".vts";      pressureDataWriter->SetFileName(op.str().data());      //dataWriter->SetDataModeToAscii();    #if VTK_MAJOR_VERSION <= 5      pressureDataWriter->SetInput(pressureDataSet);    #else      pressureDataWriter->SetInputData(pressureDataSet);    #endif      pressureDataWriter->Write();      ov << path << "/" << filename << "_velocity_" << rank << ".vts";      velocityDataWriter->SetFileName(ov.str().data());      //dataWriter->SetDataModeToAscii();    #if VTK_MAJOR_VERSION <= 5      velocityDataWriter->SetInput(velocityDataSet);//.........这里部分代码省略.........

//.........这里部分代码省略.........  ierr = DMPlexProjectFunction(dm, user.fe, user.exactFuncs, INSERT_ALL_VALUES, u);CHKERRQ(ierr);  if (user.showInitial) {ierr = DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}  if (user.runType == RUN_FULL) {    ierr = DMPlexProjectFunction(dm, user.fe, user.initialGuess, INSERT_VALUES, u);CHKERRQ(ierr);    if (user.showInitial) {ierr = DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}    if (user.debug) {      ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial guess/n");CHKERRQ(ierr);      ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    }    ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);    ierr = SNESGetIterationNumber(snes, &its);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D/n", its);CHKERRQ(ierr);    ierr = DMPlexComputeL2Diff(dm, user.fe, user.exactFuncs, u, &error);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g/n", error);CHKERRQ(ierr);    if (user.showSolution) {      ierr = PetscPrintf(PETSC_COMM_WORLD, "Solution/n");CHKERRQ(ierr);      ierr = VecChop(u, 3.0e-9);CHKERRQ(ierr);      ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    }  } else {    PetscReal res = 0.0;    /* Check discretization error */    ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial guess/n");CHKERRQ(ierr);    ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    ierr = DMPlexComputeL2Diff(dm, user.fe, user.exactFuncs, u, &error);CHKERRQ(ierr);    if (error >= 1.0e-11) {      ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g/n", error);CHKERRQ(ierr);    } else {      ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11/n", error);CHKERRQ(ierr);    }    /* Check residual */    ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "Initial Residual/n");CHKERRQ(ierr);    ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);    ierr = VecView(r, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g/n", res);CHKERRQ(ierr);    /* Check Jacobian */    {      Vec          b;      MatStructure flag;      PetscBool    isNull;      ierr = SNESComputeJacobian(snes, u, &A, &A, &flag);CHKERRQ(ierr);      ierr = MatNullSpaceTest(nullSpace, J, &isNull);CHKERRQ(ierr);      if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");      ierr = VecDuplicate(u, &b);CHKERRQ(ierr);      ierr = VecSet(r, 0.0);CHKERRQ(ierr);      ierr = SNESComputeFunction(snes, r, b);CHKERRQ(ierr);      ierr = MatMult(A, u, r);CHKERRQ(ierr);      ierr = VecAXPY(r, 1.0, b);CHKERRQ(ierr);      ierr = VecDestroy(&b);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)/n");CHKERRQ(ierr);      ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);      ierr = VecView(r, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);      ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g/n", res);CHKERRQ(ierr);    }  }  if (user.runType == RUN_FULL) {    PetscViewer viewer;    Vec         uLocal;    const char *name;    ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);    ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);    ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);    ierr = PetscViewerFileSetName(viewer, "ex62_sol.vtk");CHKERRQ(ierr);    ierr = DMGetLocalVector(dm, &uLocal);CHKERRQ(ierr);    ierr = PetscObjectGetName((PetscObject) u, &name);CHKERRQ(ierr);    ierr = PetscObjectSetName((PetscObject) uLocal, name);CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uLocal);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uLocal);CHKERRQ(ierr);    ierr = VecView(uLocal, viewer);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm, &uLocal);CHKERRQ(ierr);    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  }  ierr = PetscFree(user.exactFuncs);CHKERRQ(ierr);  ierr = DestroyElement(&user);CHKERRQ(ierr);  ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);  if (user.jacobianMF) {    ierr = VecDestroy(&userJ.u);CHKERRQ(ierr);  }  if (A != J) {    ierr = MatDestroy(&A);CHKERRQ(ierr);  }  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = VecDestroy(&u);CHKERRQ(ierr);  ierr = VecDestroy(&r);CHKERRQ(ierr);  ierr = SNESDestroy(&snes);CHKERRQ(ierr);  ierr = DMDestroy(&dm);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){    PetscErrorCode ierr;    PetscInt       M = -2, N = -3, P = 4,stencil_width = 1, dof = 1,m,n,p,xstart,ystart,zstart,i,j,k,c;    DM             da;    Vec            global,local;    PetscScalar    ****vglobal;    ierr = PetscInitialize(&argc,&argv,(char*)0,help);    CHKERRQ(ierr);    PetscFunctionBeginUser;    ierr = PetscOptionsGetInt(0,"-stencil_width",&stencil_width,0);    CHKERRQ(ierr);    ierr = PetscOptionsGetInt(0,"-dof",&dof,0);    CHKERRQ(ierr);    ierr = DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_MIRROR,DMDA_BOUNDARY_MIRROR,DMDA_BOUNDARY_MIRROR,DMDA_STENCIL_STAR,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da);    CHKERRQ(ierr);    ierr = DMDAGetCorners(da,&xstart,&ystart,&zstart,&m,&n,&p);    CHKERRQ(ierr);    ierr = DMCreateGlobalVector(da,&global);    CHKERRQ(ierr);    ierr = DMDAVecGetArrayDOF(da,global,&vglobal);    CHKERRQ(ierr);    for (k=zstart; k<zstart+p; k++) {        for (j=ystart; j<ystart+n; j++) {            for (i=xstart; i<xstart+m; i++) {                for (c=0; c<dof; c++) {                    vglobal[k][j][i][c] = 1000*k + 100*j + 10*(i+1) + c;                }            }        }    }    ierr = DMDAVecRestoreArrayDOF(da,global,&vglobal);    CHKERRQ(ierr);    ierr = DMCreateLocalVector(da,&local);    CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);    CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);    CHKERRQ(ierr);    ierr = PetscSequentialPhaseBegin(PETSC_COMM_WORLD,1);    CHKERRQ(ierr);    ierr = VecView(local,PETSC_VIEWER_STDOUT_SELF);    CHKERRQ(ierr);    ierr = PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1);    CHKERRQ(ierr);    ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);    CHKERRQ(ierr);    ierr = DMDestroy(&da);    CHKERRQ(ierr);    ierr = VecDestroy(&local);    CHKERRQ(ierr);    ierr = VecDestroy(&global);    CHKERRQ(ierr);    ierr = PetscFinalize();    return 0;}

//------------------------------------------------------------------------------int main(int argc, char *argv[]){   PetscErrorCode ierr;   DM       da;   Vec      ug, ul;   PetscInt i, ibeg, nloc, nx=200;   const PetscInt sw = 3, ndof = 1; // stencil width   PetscMPIInt rank, size;   double cfl = 0.4;   ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);   MPI_Comm_size(PETSC_COMM_WORLD, &size);   ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, nx, ndof, sw, NULL, &da); CHKERRQ(ierr);   ierr = DMSetFromOptions(da); CHKERRQ(ierr);   ierr = DMSetUp(da); CHKERRQ(ierr);   ierr = DMCreateGlobalVector(da, &ug); CHKERRQ(ierr);   ierr = DMDAGetCorners(da, &ibeg, 0, 0, &nloc, 0, 0); CHKERRQ(ierr);   ierr = DMDAGetInfo(da,0,&nx,0,0,0,0,0,0,0,0,0,0,0); CHKERRQ(ierr);   PetscReal dx = (xmax - xmin) / (PetscReal)(nx);   PetscPrintf(PETSC_COMM_WORLD,"nx = %d, dx = %e/n", nx, dx);   for(i=ibeg; i<ibeg+nloc; ++i)   {      PetscReal x = xmin + i*dx;      PetscReal v = initcond(x);      ierr = VecSetValues(ug,1,&i,&v,INSERT_VALUES); CHKERRQ(ierr);   }   ierr = VecAssemblyBegin(ug);  CHKERRQ(ierr);   ierr = VecAssemblyEnd(ug);    CHKERRQ(ierr);   savesol(nx, dx, ug);   // Get local view   ierr = DMGetLocalVector(da, &ul); CHKERRQ(ierr);   PetscInt il, nl;   ierr = DMDAGetGhostCorners(da,&il,0,0,&nl,0,0); CHKERRQ(ierr);   double res[nloc], uold[nloc];   double dt = cfl * dx;   double lam= dt/dx;   double tfinal = 2.0, t = 0.0;   while(t < tfinal)   {      if(t+dt > tfinal)      {         dt = tfinal - t;         lam = dt/dx;      }      for(int rk=0; rk<3; ++rk)      {         ierr = DMGlobalToLocalBegin(da, ug, INSERT_VALUES, ul); CHKERRQ(ierr);         ierr = DMGlobalToLocalEnd(da, ug, INSERT_VALUES, ul); CHKERRQ(ierr);         PetscScalar *u;         ierr = DMDAVecGetArrayRead(da, ul, &u); CHKERRQ(ierr);         PetscScalar *unew;         ierr = DMDAVecGetArray(da, ug, &unew); CHKERRQ(ierr);         if(rk==0)            for(i=ibeg; i<ibeg+nloc; ++i) uold[i-ibeg] = u[i];         for(i=0; i<nloc; ++i)             res[i] = 0.0;         // Loop over faces and compute flux         for(i=0; i<nloc+1; ++i)         {            // face between j-1, j            int j   = il+sw+i;            int jm1 = j-1;            int jm2 = j-2;            int jm3 = j-3;            int jp1 = j+1;            double uleft = weno5(u[jm3],u[jm2],u[jm1],u[j],u[jp1]);            double flux = uleft;            if(i==0)            {               res[i] -= flux;            }            else if(i==nloc)            {               res[i-1] += flux;            }            else            {               res[i]   -= flux;               res[i-1] += flux;            }         }         // Update solution         for(i=ibeg; i<ibeg+nloc; ++i)//.........这里部分代码省略.........