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

static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx){  User           user = (User)ctx;  DM             dau,dak;  DMDALocalInfo  infou,infok;  PetscScalar    *u,*k;  PetscScalar    *fu,*fk;  PetscErrorCode ierr;  Vec            Uloc,Kloc,Fu,Fk;  PetscFunctionBeginUser;  ierr = DMCompositeGetEntries(user->pack,&dau,&dak);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dau,&infou);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dak,&infok);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  switch (user->ptype) {  case 0:    ierr = DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,user->Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,F,&fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,F,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,user->Kloc,&k);CHKERRQ(ierr);    break;  case 1:    ierr = DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,F,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,F,&fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  case 2:    ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMCompositeGetAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  }  ierr = DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

	PetscErrorCode TcSolver::generateBC2(){	PetscErrorCode	ierr;	PetscInt	mstart, nstart, m, n, i, j;	PetscReal	**qx, **qy;	PetscReal	**bc22;	Vec		bc2Global;	ierr = VecSet(bc2, 0.0); CHKERRQ(ierr);	ierr = DMCompositeGetAccess(lambdaPack, bc2, &bc2Global, NULL); CHKERRQ(ierr);		ierr = DMDAVecGetArray(uda, qxLocal, &qx); CHKERRQ(ierr);	ierr = DMDAVecGetArray(vda, qyLocal, &qy); CHKERRQ(ierr);	ierr = DMDAVecGetArray(pda, bc2Global, &bc22); CHKERRQ(ierr);	ierr = DMDAGetCorners(pda, &mstart, &nstart, NULL, &m, &n, NULL); CHKERRQ(ierr);	//x-faces		for(j=nstart; j<nstart+n; j++)	{			//-X		if(mstart == 0) //if -X is current process		{			bc22[j][0] -= qx[j][-1];		}		//+X		if(mstart+m-1 == fluid.nx-1)		{			bc22[j][fluid.nx-1] += qx[j][fluid.nx-1];		}	}		//y-faces	for(i=mstart; i<mstart+m; i++)	{		//-Y		if(nstart ==0)		{			bc22[0][i] -= qy[-1][i];		}		//+Y		if(nstart+n-1 == fluid.ny-1)		{			bc22[fluid.ny-1][i] += qy[fluid.ny-1][i];		}	}		ierr = DMDAVecRestoreArray(pda, bc2Global, &bc22); CHKERRQ(ierr);	ierr = DMDAVecRestoreArray(uda, qxLocal, &qx); CHKERRQ(ierr);	ierr = DMDAVecRestoreArray(vda, qyLocal, &qy); CHKERRQ(ierr);	ierr = DMCompositeRestoreAccess(lambdaPack, bc2, &bc2Global, NULL); CHKERRQ(ierr);	return 0;}

PetscErrorCode ini_bou(Vec X,AppCtx* user){  PetscErrorCode ierr;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p;  Vec            gc;  PetscInt       i,j;  PetscInt       xs,ys,xm,ym,M,N;  PetscScalar    xi,yi;  PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;  PetscScalar    rho   =user->rho;  PetscScalar    mux   =user->mux,muy=user->muy;  PetscMPIInt    rank;  PetscScalar    sum;  PetscFunctionBeginUser;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable     muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point     in the y-direction. We only modify mux here  */  mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */  if (user->nonoiseinitial) {    for (i=xs; i < xs+xm; i++) {      for (j=ys; j < ys+ym; j++) {        xi = coors[j][i].x; yi = coors[j][i].y;        if ((xi == mux) && (yi == muy)) {          p[j][i] = 1.0;        }      }    }  } else {    /* Change PM_min accordingly */    user->PM_min = user->Pmax*sin(mux);    for (i=xs; i < xs+xm; i++) {      for (j=ys; j < ys+ym; j++) {        xi = coors[j][i].x; yi = coors[j][i].y;        p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay)));      }    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);  ierr = VecSum(X,&sum);CHKERRQ(ierr);  ierr = VecScale(X,1.0/sum);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*      Evaluates FU = Gradiant(L(w,u,lambda))     This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and   DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).*/PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx){  PetscErrorCode ierr;  PetscInt       xs,xm,i,N;  ULambda        *u_lambda,*fu_lambda;  PetscScalar    d,h,*w,*fw;  Vec            vw,vfw,vu_lambda,vfu_lambda;  DM             packer,red,da;  PetscFunctionBeginUser;  ierr = VecGetDM(U, &packer);CHKERRQ(ierr);  ierr = DMCompositeGetEntries(packer,&red,&da);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeScatter(packer,U,vw,vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = VecGetArray(vw,&w);CHKERRQ(ierr);  ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);  d    = N-1.0;  h    = 1.0/d;  /* derivative of L() w.r.t. w */  if (xs == 0) { /* only first processor computes this */    fw[0] = -2.0*d*u_lambda[0].lambda;  }  /* derivative of L() w.r.t. u */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   fu_lambda[0].lambda   =    h*u_lambda[0].u   + 2.*d*u_lambda[0].lambda   - d*u_lambda[1].lambda;    else if (i == 1)   fu_lambda[1].lambda   = 2.*h*u_lambda[1].u   + 2.*d*u_lambda[1].lambda   - d*u_lambda[2].lambda;    else if (i == N-1) fu_lambda[N-1].lambda =    h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda;    else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda;    else               fu_lambda[i].lambda   = 2.*h*u_lambda[i].u   - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda);  }  /* derivative of L() w.r.t. lambda */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   fu_lambda[0].u   = 2.0*d*(u_lambda[0].u - w[0]);    else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;    else               fu_lambda[i].u   = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h);  }  ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr);  ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);  ierr = PetscLogFlops(13.0*N);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x,*f;  PetscFunctionBegin;  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidHeaderSpecific(F,VEC_CLASSID,3);  if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  switch (dmdasnes->residuallocalimode) {  case INSERT_VALUES: {    ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);    ierr = PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);  } break;  case ADD_VALUES: {    Vec Floc;    ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);    ierr = VecZeroEntries(Floc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = VecZeroEntries(F);CHKERRQ(ierr);    ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);  } break;  default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);  }  ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  if (snes->domainerror) {    ierr = VecSetInf(F);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx){  PetscErrorCode ierr;  AppCtx         *user=(AppCtx*)ctx;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p,**f,**pdot;  PetscInt       i,j;  PetscInt       xs,ys,xm,ym,M,N;  Vec            localX,gc,localXdot;  PetscScalar    p_adv1,p_adv2,p_diff;  PetscFunctionBeginUser;  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMGetLocalVector(user->da,&localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(user->da,&localXdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,localX,&p);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,localXdot,&pdot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,F,&f);CHKERRQ(ierr);  user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda));  for (i=xs; i < xs+xm; i++) {    for (j=ys; j < ys+ym; j++) {      if (i == 0 || j == 0 || i == M-1 || j == N-1) {        ierr = BoundaryConditions(p,coors,i,j,M,N,f,user);CHKERRQ(ierr);      } else {        ierr = adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);CHKERRQ(ierr);        ierr = adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);CHKERRQ(ierr);        ierr = diffuse(p,i,j,t,&p_diff,user);CHKERRQ(ierr);        f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];      }    }  }  ierr = DMDAVecRestoreArray(user->da,localX,&p);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,localX,&pdot);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->da,&localX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->da,&localXdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,F,&f);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr){  User           user = (User)ptr;  DM             da;  DMDALocalInfo  info;  PetscInt       i;  Field          *x,*xdot,*f;  PetscReal      hx;  Vec            Xloc;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)(info.mx-1);  /*     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 = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    if (i == 0) {      f[i].u = hx * (x[i].u - user->uleft);      f[i].v = hx * (x[i].v - user->vleft);    } else if (i == info.mx-1) {      f[i].u = hx * (x[i].u - user->uright);      f[i].v = hx * (x[i].v - user->vright);    } else {      f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;      f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode DiffusiveTerm<2>::initializeFluxes(){  PetscErrorCode ierr;  PetscInt i, j,           // loop indices           m, n,           // local number of nodes along each direction           mstart, nstart; // starting indices  PetscReal x, y; // coordinates of a node  Vec qxGlobal, qyGlobal;  ierr = DMCompositeGetAccess(qPack, q, &qxGlobal, &qyGlobal); CHKERRQ(ierr);  // fluxes in x-direction  PetscReal **qx;  ierr = DMDAVecGetArray(uda, qxGlobal, &qx); CHKERRQ(ierr);  ierr = DMDAGetCorners(uda, &mstart, &nstart, NULL, &m, &n, NULL); CHKERRQ(ierr);  PetscReal u;  // value of u-velocity  for (j=nstart; j<nstart+n; j++)  {    for (i=mstart; i<mstart+m; i++)    {      x = mesh->x[i+1];      y = 0.5*(mesh->y[j]+mesh->y[j+1]);      u = sin(PETSC_PI*x)*sin(PETSC_PI*y) + 1.0;      qx[j][i] = u*mesh->dy[j];    }  }  ierr = DMDAVecRestoreArray(uda, qxGlobal, &qx); CHKERRQ(ierr);  // fluxes in y-direction  PetscReal **qy;  ierr = DMDAVecGetArray(vda, qyGlobal, &qy); CHKERRQ(ierr);  ierr = DMDAGetCorners(vda, &mstart, &nstart, NULL, &m, &n, NULL); CHKERRQ(ierr);  PetscReal v;  // value of v-velocity  for (j=nstart; j<nstart+n; j++)  {    for (i=mstart; i<mstart+m; i++)    {      x = 0.5*(mesh->x[i]+mesh->x[i+1]);      y = mesh->y[j+1];      v = sin(PETSC_PI*x)*sin(PETSC_PI*y);      qy[j][i] = v*mesh->dx[i];    }  }  ierr = DMDAVecRestoreArray(vda, qyGlobal, &qy); CHKERRQ(ierr);  ierr = DMCompositeRestoreAccess(qPack, q, &qxGlobal, &qyGlobal); CHKERRQ(ierr);  return ierr;} // initializeFluxes

PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx){  PetscFunctionBegin;  PetscErrorCode ierr;  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "FormTangent not implemented, use -snes_mf");  DMDALocalInfo    info;  DM               da_dof;  Vec              localU,localUdot; // local versions  Field          **h,**hdot;  /* get the da from the snes */  ierr = TSGetDM(ts,(DM*)&da_dof);CHKERRQ(ierr);  /* handle the vec U */  ierr = DMGetLocalVector(da_dof,&localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da_dof,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da_dof,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /* handle the vec Udot */  ierr = DMGetLocalVector(da_dof,&localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da_dof,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da_dof,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  /* Get the arrays from the vectors */  ierr = DMDAVecGetArray(da_dof,localU,&h);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da_dof,localUdot,&hdot);CHKERRQ(ierr);  /* Grab the local info and call the local tangent routine */  ierr = DMDAGetLocalInfo(da_dof,&info);CHKERRQ(ierr);CHKERRQ(ierr);  ierr = MatZeroEntries(*B);CHKERRQ(ierr); // pre-zero the matrix  ierr = FormTangentLocal(&info,t,h,hdot,shift,B,(AppCtx *) ctx);CHKERRQ(ierr);  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  if (*A != *B) { // then we could be matrix free    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  }  *flag = SAME_NONZERO_PATTERN; /* the sparsity pattern does not change */  ierr = DMDAVecRestoreArray(da_dof,localUdot,&hdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da_dof,localU,&h);CHKERRQ(ierr);  /* Restore the arrays and local vectors */  ierr = DMRestoreLocalVector(da_dof,&localU);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da_dof,&localUdot);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*  This function should eventually replace:    DMDAComputeFunction() and DMDAComputeFunction1() */static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMTS_DA        *dmdats = (DMTS_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x,*f,*xdot;  PetscFunctionBegin;  PetscValidHeaderSpecific(ts,TS_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidHeaderSpecific(F,VEC_CLASSID,3);  if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr);  switch (dmdats->ifunctionlocalimode) {  case INSERT_VALUES: {    ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);  } break;  case ADD_VALUES: {    Vec Floc;    ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);    ierr = VecZeroEntries(Floc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = VecZeroEntries(F);CHKERRQ(ierr);    ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);  } break;  default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);  }  ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr){  User           user = (User)ptr;  DM             da;  Vec            Xloc;  DMDALocalInfo  info;  PetscInt       i,j;  PetscReal      hx;  Field          *x,*f;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)info.mx;  /*     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 = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    const PetscReal *a     = user->a;    PetscReal       u0t[2] = {1. - PetscPowScalar(sin(12*t),4.),0};    for (j=0; j<2; j++) {      if (i == 0)              f[i][j] = a[j]/hx*(1./3*u0t[j] + 0.5*x[i][j] - x[i+1][j] + 1./6*x[i+2][j]);      else if (i == 1)         f[i][j] = a[j]/hx*(-1./12*u0t[j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);      else if (i == info.mx-2) f[i][j] = a[j]/hx*(-1./6*x[i-2][j] + x[i-1][j] - 0.5*x[i][j] - 1./3*x[i+1][j]);      else if (i == info.mx-1) f[i][j] = a[j]/hx*(-x[i][j] + x[i-1][j]);      else                     f[i][j] = a[j]/hx*(-1./12*x[i-2][j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode efs_set_state(efs *slv, double rhs[], double phi[], double eps[],                             double debye[]){	PetscErrorCode ierr;	slv->state.phi = phi;	slv->state.eps = eps;	slv->state.debye = debye;	slv->state.rhs = rhs;	slv->state.sol = NULL;	if (slv->grid.nd == 2) {		PetscScalar **epsvec;		double **eps_state;		PetscInt i, j;		PetscInt xs, ys, xm, ym;		ierr = DMDAGetCorners(slv->dm, &xs, &ys, 0, &xm, &ym, 0);CHKERRQ(ierr);		ierr = ef_dmap_get(slv->dmap, slv->state.eps, &eps_state);CHKERRQ(ierr);		ierr = DMDAVecGetArray(slv->dm, slv->levels[0].eps, &epsvec);CHKERRQ(ierr);		for (j = ys; j < ys + ym; j++) {			for (i = xs; i < xs + xm; i++) {				epsvec[j][i] = eps_state[j][i];			}		}		ierr = ef_dmap_restore(slv->dmap, &eps_state);CHKERRQ(ierr);		ierr = DMDAVecRestoreArray(slv->dm, slv->levels[0].eps, &epsvec);CHKERRQ(ierr);	} else if (slv->grid.nd == 3) {		PetscScalar ***epsvec;		double ***eps_state;		PetscInt i, j, k;		PetscInt xs, ys, zs, xm, ym, zm;		ierr = DMDAGetCorners(slv->dm, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr);		ierr = ef_dmap_get(slv->dmap, slv->state.eps, &eps_state);CHKERRQ(ierr);		ierr = DMDAVecGetArray(slv->dm, slv->levels[0].eps, &epsvec);CHKERRQ(ierr);		for (k = zs; k < zs + zm; k++) {			for (j = ys; j < ys + ym; j++) {				for (i = xs; i < xs + xm; i++) {					epsvec[k][j][i] = eps_state[k][j][i];				}			}		}		ierr = ef_dmap_restore(slv->dmap, &eps_state);CHKERRQ(ierr);		ierr = DMDAVecRestoreArray(slv->dm, slv->levels[0].eps, &epsvec);CHKERRQ(ierr);	}	return 0;}

/*   FormFunctionLocal - Evaluates nonlinear residual, F(x) on local process patch*/PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field *u, Field *f, AppCtx *user){  Vec L;  PetscReal      phi = user->phi;  PetscReal      dt  = user->dt;  PetscReal      dx  = 1.0/(PetscReal)(info->mx-1);  PetscReal      alpha = 2.0;  PetscReal      beta  = 2.0;  PetscReal      kappaWet   = user->kappaWet;  PetscReal      kappaNoWet = user->kappaNoWet;  Field         *uold;  PetscScalar   *Kappa;  PetscInt       i;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMGetGlobalVector(user->cda, &L);CHKERRQ(ierr);  ierr = DMDAVecGetArray(info->da, user->uold,  &uold);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->cda, user->Kappa, &Kappa);CHKERRQ(ierr);  /* Compute residual over the locally owned part of the grid */  for(i = info->xs; i < info->xs+info->xm; ++i) {    if (i == 0) {      f[i].s = u[i].s - user->sl;      f[i].v = u[i].v - user->vl;      f[i].p = u[i].p - user->pl;    } else {      PetscScalar K          = 2*dx/(dx/Kappa[i] + dx/Kappa[i-1]);      PetscReal   lambdaWet  = kappaWet*pow(u[i].s, alpha);      PetscReal   lambda     = lambdaWet + kappaNoWet*pow(1-u[i].s, beta);      PetscReal   lambdaWetL = kappaWet*pow(u[i-1].s, alpha);      PetscReal   lambdaL    = lambdaWetL + kappaNoWet*pow(1-u[i-1].s, beta);      f[i].s = phi*(u[i].s - uold[i].s) + (dt/dx)*((lambdaWet/lambda)*u[i].v - (lambdaWetL/lambdaL)*u[i-1].v);      f[i].v = u[i].v + K*lambda*(u[i].p - u[i-1].p)/dx;      //pxx     = (2.0*u[i].p - u[i-1].p - u[i+1].p)/dx;      f[i].p = u[i].v - u[i-1].v;    }  }  ierr = DMDAVecRestoreArray(info->da, user->uold, &uold);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->cda, user->Kappa, &Kappa);CHKERRQ(ierr);  /* ierr = PetscLogFlops(11.0*info->ym*info->xm);CHKERRQ(ierr); */  ierr = DMRestoreGlobalVector(user->cda, &L);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates1d(DM da){  PetscErrorCode ierr;  PetscInt       i,start,m;  Vec            gc,global;  PetscScalar    *coors;  DM             cda;  PetscFunctionBeginUser;  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&start,0,0,&m,0,0);CHKERRQ(ierr);  for (i=start; i<start+m; i++) {    if (i % 2) {      coors[i] = coors[i-1] + .1*(coors[i+1] - coors[i-1]);    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates3d(DM da){  PetscErrorCode ierr;  PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;  Vec            gc,global;  DMDACoor3d     ***coors;  DM             cda;  PetscFunctionBeginUser;  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);CHKERRQ(ierr);  for (i=mstart; i<mstart+m; i++) {    for (j=nstart; j<nstart+n; j++) {      for (k=pstart; k<pstart+p; k++) {        if (i % 2) {          coors[k][j][i].x = coors[k][j][i-1].x + .1*(coors[k][j][i+1].x - coors[k][j][i-1].x);        }        if (j % 2) {          coors[k][j][i].y = coors[k][j-1][i].y + .3*(coors[k][j+1][i].y - coors[k][j-1][i].y);        }        if (k % 2) {          coors[k][j][i].z = coors[k-1][j][i].z + .4*(coors[k+1][j][i].z - coors[k-1][j][i].z);        }      }    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode checkPositivity(PorousCtx *user,Vec X) {  PetscErrorCode ierr;  PetscInt       i,j,xs,ys,xm,ym;  DM             da = user->da;  WPnode         **wp;  PetscFunctionBegin;  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,X,&wp);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      if (wp[j][i].W < 0) {        SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,          "negative water thickness W at (i,j)=(%d,%d) during function eval %d",          i,j,user->fcncount);      }      if (wp[j][i].P < 0) {        SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,          "negative water pressure P at (i,j)=(%d,%d) during function eval %d",          i,j,user->fcncount);      }    }  }  ierr = DMDAVecRestoreArray(da,X,&wp);CHKERRQ(ierr);  PetscFunctionReturn(0);} 

PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx){  AppCtx         *user;  PetscInt       i,j,xs,ys,xm,ym;  PetscErrorCode ierr;  PetscReal      tleft;  PetscScalar    **x;  DM             da;  PetscFunctionBeginUser;  ierr  = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr  = DMGetApplicationContext(da,&user);CHKERRQ(ierr);  tleft = user->tleft;  /* Get ghost points */  ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  /* Compute initial guess */  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      x[j][i] = tleft;    }  }  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*   Solution - Computes the exact solution at a given time.   Input Parameters:   t - current time   solution - vector in which exact solution will be computed   appctx - user-defined application context   Output Parameter:   solution - vector with the newly computed exact solution*/PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx){  PetscScalar    *u,ex1,ex2,sc1,sc2,h;  PetscErrorCode ierr;  PetscInt       i,mstart,mend,xm,M;  DM             da;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  h    = 1.0/M;  mend = mstart + xm;  /*     Get a pointer to vector data.  */  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);  /*     Simply write the solution directly into the array locations.     Alternatively, we culd use VecSetValues() or VecSetValuesLocal().  */  ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);  ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);  sc1 = PETSC_PI*6.*h;                 sc2 = PETSC_PI*2.*h;  for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(sc1*(PetscReal)i + appctx->a*PETSC_PI*6.*t)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i + appctx->a*PETSC_PI*2.*t)*ex2;  /*     Restore vector  */  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);  return 0;}

PetscErrorCode FormInitialGuess(DM da,void *ctx,Vec X){  PetscInt       i,j,Mx,My,xs,ys,xm,ym;  PetscErrorCode ierr;  PetscScalar    **x;  PetscReal      x0,x1;  PetscFunctionBeginUser;  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 = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      x0 = 10.0*(i - 0.5*(Mx-1)) / (Mx-1);      x1 = 10.0*(j - 0.5*(Mx-1)) / (My-1);      if (PetscAbsScalar(x0) < 0.5 && PetscAbsScalar(x1) < 0.5) {        x[j][i] = 1.;      } else {        x[j][i] = 0.0;      }    }  }  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/* FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch */PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,ObsCtx *user) {  PetscErrorCode ierr;  PetscInt       i,j;  PetscReal      dx,dy,uxx,uyy;  PetscReal      **uexact;  /* used for boundary values only */  PetscFunctionBeginUser;  dx = 4.0 / (PetscReal)(info->mx-1);  dy = 4.0 / (PetscReal)(info->my-1);  ierr = DMDAVecGetArray(info->da, user->uexact, &uexact);CHKERRQ(ierr);  for (j=info->ys; j<info->ys+info->ym; j++) {    for (i=info->xs; i<info->xs+info->xm; i++) {      if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {        f[j][i] = 4.0*(x[j][i] - uexact[j][i]);      } else {        uxx     = dy*(x[j][i-1] - 2.0 * x[j][i] + x[j][i+1]) / dx;        uyy     = dx*(x[j-1][i] - 2.0 * x[j][i] + x[j+1][i]) / dy;        f[j][i] = -uxx - uyy;      }    }  }  ierr = DMDAVecRestoreArray(info->da, user->uexact, &uexact);CHKERRQ(ierr);  ierr = PetscLogFlops(10.0*info->ym*info->xm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMTS_DA        *dmdats = (DMTS_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x;  PetscFunctionBegin;  if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  if (dmdats->rhsjacobianlocal) {    ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");  /* This will be redundant if the user called both, but it's too common to forget. */  if (*A != *B) {    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/*   ComputeInitialGuess - Calculates the initial guess   Input Parameters:.  user - user-defined application context.  X - vector for initial guess   Output Parameters:.  X - newly computed initial guess*/PetscErrorCode ComputeInitialGuess(SNES snes, Vec X,void *dummy){  PetscErrorCode ierr;  PetscInt       i,j,mx,my;  DM             da;  AppCtx         *user;  PetscScalar    **x;  PetscInt       xs,xm,ys,ym;  PetscFunctionBeginUser;  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr = SNESGetApplicationContext(snes,(void**)&user);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);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);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  /* Perform local computations */  for (j=ys; j<ys+ym; j++) {    for (i=xs; i< xs+xm; i++) {      x[j][i] = (((j+1.0)*user->bottom[i+1]+(my-j+1.0)*user->top[i+1])/(my+2.0)+((i+1.0)*user->left[j+1]+(mx-i+1.0)*user->right[j+1])/(mx+2.0))/2.0;    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormDiffusionCoefficient(DM da,void *ctx,Vec X){  PetscInt       i,j,Mx,My,xs,ys,xm,ym;  PetscErrorCode ierr;  PetscScalar    **x;  PetscReal      x1,min=0.1;  PetscFunctionBeginUser;  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 = VecSetRandom(X,PETSC_NULL);  ierr = VecMin(X,PETSC_NULL,&min);CHKERRQ(ierr);   */  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      x1 = 10.0*(j - 0.5*(Mx-1)) / (My-1);      x[j][i] += min + (1. - min)*PetscSqr(sin(2.*PETSC_PI*x1));    }  }  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx){  User           user = (User)ctx;  DM             da;  PetscInt       i;  DMDALocalInfo  info;  Field          *x;  PetscReal      hx;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)info.mx;  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    PetscReal r = (i+1)*hx,ik = user->k[1] != 0.0 ? 1.0/user->k[1] : 1.0;    x[i][0] = 1 + user->s[1]*r;    x[i][1] = user->k[0]*ik*x[i][0] + user->s[1]*ik;  }  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ini_bou(Vec X,AppCtx* user){  PetscErrorCode ierr;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p;  Vec            gc;  PetscInt       M,N,I,J;  PetscMPIInt    rank;  PetscFunctionBeginUser;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);  /* Point mass at (mux,muy) */  ierr = PetscPrintf(PETSC_COMM_WORLD,"Original user->mux = %f, user->muy = %f/n",user->mux,user->muy);CHKERRQ(ierr);  ierr = DMDAGetLogicalCoordinate(user->da,user->mux,user->muy,0.0,&I,&J,NULL,&user->mux,&user->muy,NULL);CHKERRQ(ierr);  user->PM_min = user->Pmax*PetscSinScalar(user->mux);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Corrected user->mux = %f, user->muy = %f user->PM_min = %f,user->dx = %f/n",user->mux,user->muy,user->PM_min,user->dx);CHKERRQ(ierr);  if (I > -1 && J > -1) {    p[J][I] = 1.0;  }  ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx){  PetscErrorCode ierr;  PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;  DM             dm;  PetscScalar    Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;  PetscScalar    ***barray;  PetscFunctionBeginUser;  ierr    = KSPGetDM(ksp,&dm);CHKERRQ(ierr);  ierr    = DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  Hx      = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);  HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;  ierr    = DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);  ierr    = DMDAVecGetArray(dm,b,&barray);CHKERRQ(ierr);  for (k=zs; k<zs+zm; k++) {    for (j=ys; j<ys+ym; j++) {      for (i=xs; i<xs+xm; i++) {        if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {          barray[k][j][i] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);        } else {          barray[k][j][i] = Hx*Hy*Hz;        }      }    }  }  ierr = DMDAVecRestoreArray(dm,b,&barray);CHKERRQ(ierr);  PetscFunctionReturn(0);}