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

int main(int argc, char **argv){#if !defined(PETSC_USE_COMPLEX)  PetscErrorCode ierr;  Vec            x,yp1,yp2,yp3,yp4,ym1,ym2,ym3,ym4;  PetscReal      *values;  PetscViewer    viewer_in,viewer_outp1,viewer_outp2,viewer_outp3,viewer_outp4;  PetscViewer    viewer_outm1,viewer_outm2,viewer_outm3,viewer_outm4;  DM             daf,dac1,dac2,dac3,dac4,daf1,daf2,daf3,daf4;  Vec            scaling_p1,scaling_p2,scaling_p3,scaling_p4;  Mat            interp_p1,interp_p2,interp_p3,interp_p4,interp_m1,interp_m2,interp_m3,interp_m4;#endif  PetscInitialize(&argc,&argv, (char*)0, help);#if defined(PETSC_USE_COMPLEX)  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for complex numbers");#else  ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,1024,1024,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&daf);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(daf,&x);CHKERRQ(ierr);  ierr = VecGetArray(x,&values);CHKERRQ(ierr);  ierr = DMCoarsen(daf,PETSC_COMM_WORLD,&dac1);CHKERRQ(ierr);  ierr = DMCoarsen(dac1,PETSC_COMM_WORLD,&dac2);CHKERRQ(ierr);  ierr = DMCoarsen(dac2,PETSC_COMM_WORLD,&dac3);CHKERRQ(ierr);  ierr = DMCoarsen(dac3,PETSC_COMM_WORLD,&dac4);CHKERRQ(ierr);  ierr = DMRefine(daf,PETSC_COMM_WORLD,&daf1);CHKERRQ(ierr);  ierr = DMRefine(daf1,PETSC_COMM_WORLD,&daf2);CHKERRQ(ierr);  ierr = DMRefine(daf2,PETSC_COMM_WORLD,&daf3);CHKERRQ(ierr);  ierr = DMRefine(daf3,PETSC_COMM_WORLD,&daf4);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dac1,&yp1);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dac2,&yp2);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dac3,&yp3);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dac4,&yp4);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(daf1,&ym1);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(daf2,&ym2);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(daf3,&ym3);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(daf4,&ym4);CHKERRQ(ierr);  ierr = DMCreateInterpolation(dac1,daf,&interp_p1,&scaling_p1);CHKERRQ(ierr);  ierr = DMCreateInterpolation(dac2,dac1,&interp_p2,&scaling_p2);CHKERRQ(ierr);  ierr = DMCreateInterpolation(dac3,dac2,&interp_p3,&scaling_p3);CHKERRQ(ierr);  ierr = DMCreateInterpolation(dac4,dac3,&interp_p4,&scaling_p4);CHKERRQ(ierr);  ierr = DMCreateInterpolation(daf,daf1,&interp_m1,NULL);CHKERRQ(ierr);  ierr = DMCreateInterpolation(daf1,daf2,&interp_m2,NULL);CHKERRQ(ierr);  ierr = DMCreateInterpolation(daf2,daf3,&interp_m3,NULL);CHKERRQ(ierr);  ierr = DMCreateInterpolation(daf3,daf4,&interp_m4,NULL);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer_in);CHKERRQ(ierr);  ierr = PetscViewerBinaryRead(viewer_in,values,1048576,PETSC_DOUBLE);CHKERRQ(ierr);  ierr = MatRestrict(interp_p1,x,yp1);  ierr = VecPointwiseMult(yp1,yp1,scaling_p1);CHKERRQ(ierr);  ierr = MatRestrict(interp_p2,yp1,yp2);  ierr = VecPointwiseMult(yp2,yp2,scaling_p2);CHKERRQ(ierr);  ierr = MatRestrict(interp_p3,yp2,yp3);  ierr = VecPointwiseMult(yp3,yp3,scaling_p3);CHKERRQ(ierr);  ierr = MatRestrict(interp_p4,yp3,yp4);  ierr = VecPointwiseMult(yp4,yp4,scaling_p4);CHKERRQ(ierr);  ierr = MatRestrict(interp_m1,x,ym1);  ierr = MatRestrict(interp_m2,ym1,ym2);  ierr = MatRestrict(interp_m3,ym2,ym3);  ierr = MatRestrict(interp_m4,ym3,ym4);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_WRITE,&viewer_outp1);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_WRITE,&viewer_outp2);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_WRITE,&viewer_outp3);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_WRITE,&viewer_outp4);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_WRITE,&viewer_outm1);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_WRITE,&viewer_outm2);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_WRITE,&viewer_outm3);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_WRITE,&viewer_outm4);CHKERRQ(ierr);  ierr = VecView(yp1,viewer_outp1);CHKERRQ(ierr);  ierr = VecView(x,viewer_outp1);CHKERRQ(ierr);  ierr = VecView(yp2,viewer_outp2);CHKERRQ(ierr);  ierr = VecView(yp3,viewer_outp3);CHKERRQ(ierr);  ierr = VecView(yp4,viewer_outp4);CHKERRQ(ierr);  ierr = VecView(ym1,viewer_outm1);CHKERRQ(ierr);  ierr = VecView(ym2,viewer_outm2);CHKERRQ(ierr);  ierr = VecView(ym3,viewer_outm3);CHKERRQ(ierr);  ierr = VecView(ym4,viewer_outm4);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_in);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outp1);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outp2);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outp3);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outp4);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outm1);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outm2);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outm3);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer_outm4);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = VecDestroy(&yp1);CHKERRQ(ierr);  ierr = VecDestroy(&yp2);CHKERRQ(ierr);  ierr = VecDestroy(&yp3);CHKERRQ(ierr);  ierr = VecDestroy(&yp4);CHKERRQ(ierr);  ierr = VecDestroy(&ym1);CHKERRQ(ierr);  ierr = VecDestroy(&ym2);CHKERRQ(ierr);//.........这里部分代码省略.........

PetscErrorCode efs_setup(efs *slv, int offset[], int stride[]){	PetscErrorCode ierr;	PetscInt xs, ys, zs, xm, ym, zm;	PCType pc_type;	if (efs_log(slv, EFS_LOG_STATUS)) {		ierr = ef_io_print(slv->comm, "Setting up electric field solver");CHKERRQ(ierr);	}	slv->ts = 0;	if (efs_log(slv, EFS_LOG_RESIDUAL)) {		ierr = PetscOptionsSetValue("-ksp_monitor_short", NULL);CHKERRQ(ierr);	}	ierr = DMDASetFieldName(slv->dm, 0,"potential");CHKERRQ(ierr);	if (slv->grid.nd == 2) {		ierr = DMDAGetCorners(slv->dm, &xs, &ys, 0, &xm, &ym, 0);CHKERRQ(ierr);		slv->dmap = ef_dmap_create_2d(xs - offset[0], ys - offset[1], xm, ym, stride);	} else if (slv->grid.nd == 3) {		ierr = DMDAGetCorners(slv->dm, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr);		slv->dmap = ef_dmap_create_3d(xs - offset[0], ys - offset[1], zs - offset[2], xm, ym, zm, stride);	} else {		SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimmension: %d", slv->grid.nd);	}	ierr = ef_callback_create(&slv->callback);CHKERRQ(ierr);	ierr = KSPCreate(slv->comm, &slv->ksp);CHKERRQ(ierr);	if (efs_log(slv, EFS_LOG_EIGS)) {		ierr = KSPSetComputeEigenvalues(slv->ksp, PETSC_TRUE);CHKERRQ(ierr);	}	ierr = KSPSetDM(slv->ksp, slv->dm);CHKERRQ(ierr);	ierr = KSPGetPC(slv->ksp, &slv->pc);CHKERRQ(ierr);	ierr = PCSetType(slv->pc, PCMG);CHKERRQ(ierr);	if (slv->options.galerkin) {		ierr = PCMGSetGalerkin(slv->pc, PETSC_TRUE);CHKERRQ(ierr);	} else {		ierr = PCMGSetGalerkin(slv->pc, PETSC_FALSE);CHKERRQ(ierr);	}	ierr = KSPSetComputeOperators(slv->ksp, slv->callback->matrix, slv);CHKERRQ(ierr);	ierr = KSPSetComputeRHS(slv->ksp, slv->callback->rhs, slv);CHKERRQ(ierr);	ierr = KSPSetComputeInitialGuess(slv->ksp, slv->callback->guess, slv);CHKERRQ(ierr);	ierr = KSPSetFromOptions(slv->ksp);CHKERRQ(ierr);	ierr = PCGetType(slv->pc, &pc_type);CHKERRQ(ierr);	ierr = PCMGGetLevels(slv->pc, &slv->options.levels);CHKERRQ(ierr);	if (slv->options.levels < 1) slv->options.levels++;	ierr = PCMGGetGalerkin(slv->pc, &slv->options.galerkin);CHKERRQ(ierr);	if (strcmp(pc_type, PCGAMG) == 0	    || strcmp(pc_type, PCHYPRE) == 0) slv->options.galerkin = 1;	if (!slv->options.galerkin) {		slv->levels = (ef_level*) malloc(slv->options.levels*sizeof(ef_level));		// setup callback for transforming rhs on coarse levels	} else {		slv->levels = (ef_level*) malloc(sizeof(ef_level));	}	ierr = ef_fd_create(&slv->fd, EF_FD_STANDARD_O2);CHKERRQ(ierr);	ierr = ef_operator_create(&slv->op, slv->levels, slv->fd, slv->grid.nd);CHKERRQ(ierr);	ierr = ef_boundary_create(&slv->boundary, slv->levels, slv->options.levels,	                          slv->dmap, &slv->state, slv->fd);CHKERRQ(ierr);	slv->op->axisymmetric = slv->options.axisymmetric;	slv->boundary->axisymmetric = slv->options.axisymmetric;	ierr = DMSetMatType(slv->dm, MATAIJ);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].eps);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].g);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].ag);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].gcomp);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].scale);CHKERRQ(ierr);	ierr = DMCreateGlobalVector(slv->dm, &slv->levels[0].nscale);CHKERRQ(ierr);	ierr = VecSet(slv->levels[0].g, 0);CHKERRQ(ierr);	ierr = VecSet(slv->levels[0].gcomp, 1);CHKERRQ(ierr);	ierr = VecSet(slv->levels[0].scale, 1);CHKERRQ(ierr);	ierr = VecSet(slv->levels[0].nscale, 1);CHKERRQ(ierr);	ierr = DMCoarsenHookAdd(slv->dm, slv->callback->coarsen, slv->callback->restrct, slv);CHKERRQ(ierr);	return 0;}

int main(int argc,char **argv){  AppCtx         appctx;                 /* user-defined application context */  PetscErrorCode ierr;  PetscInt       i, xs, xm, ind, j, lenglob;  PetscReal      x, *wrk_ptr1, *wrk_ptr2;  MatNullSpace   nsp;  PetscMPIInt    size;   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Initialize program and set problem parameters     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  PetscFunctionBegin;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  /*initialize parameters */  appctx.param.N    = 10;  /* order of the spectral element */  appctx.param.E    = 10;  /* number of elements */  appctx.param.L    = 4.0;  /* length of the domain */  appctx.param.mu   = 0.01; /* diffusion coefficient */  appctx.initial_dt = 5e-3;  appctx.param.steps = PETSC_MAX_INT;  appctx.param.Tend  = 4;  ierr = PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);CHKERRQ(ierr);  appctx.param.Le = appctx.param.L/appctx.param.E;  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);  if (appctx.param.E % size) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of elements must be divisible by number of processes");  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create GLL data structures     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = PetscGLLCreate(appctx.param.N,PETSCGLL_VIA_LINEARALGEBRA,&appctx.SEMop.gll);CHKERRQ(ierr);  lenglob  = appctx.param.E*(appctx.param.N-1);  /*     Create distributed array (DMDA) to manage parallel grid and vectors     and to set up the ghost point communication pattern.  There are E*(Nl-1)+1     total grid values spread equally among all the processors, except first and last  */  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);CHKERRQ(ierr);  ierr = DMSetFromOptions(appctx.da);CHKERRQ(ierr);  ierr = DMSetUp(appctx.da);CHKERRQ(ierr);   /*     Extract global and local vectors from DMDA; we use these to store the     approximate solution.  Then duplicate these for remaining vectors that     have the same types.  */  ierr = DMCreateGlobalVector(appctx.da,&appctx.dat.curr_sol);CHKERRQ(ierr);  ierr = VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.grid);CHKERRQ(ierr);  ierr = VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.mass);CHKERRQ(ierr);  ierr = DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  ierr = DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);CHKERRQ(ierr);  ierr = DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);CHKERRQ(ierr);    /* Compute function over the locally owned part of the grid */      xs=xs/(appctx.param.N-1);    xm=xm/(appctx.param.N-1);    /*      Build total grid and mass over entire mesh (multi-elemental)   */   for (i=xs; i<xs+xm; i++) {    for (j=0; j<appctx.param.N-1; j++) {      x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;       ind=i*(appctx.param.N-1)+j;      wrk_ptr1[ind]=x;      wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];      if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];    }   }  ierr = DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   Create matrix data structure; set matrix evaluation routine.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);CHKERRQ(ierr);  ierr = DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);CHKERRQ(ierr);  ierr = DMCreateMatrix(appctx.da,&appctx.SEMop.grad);CHKERRQ(ierr);  /*   For linear problems with a time-dependent f(u,t) in the equation   u_t = f(u,t), the user provides the discretized right-hand-side   as a time-dependent matrix.   */  ierr = RHSMatrixLaplaciangllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);CHKERRQ(ierr);  ierr = RHSMatrixAdvectiongllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.grad,appctx.SEMop.grad,&appctx);CHKERRQ(ierr);   /*       For linear problems with a time-dependent f(u,t) in the equation//.........这里部分代码省略.........

/*@    DMDASetUniformCoordinates - Sets a DMDA coordinates to be a uniform grid  Collective on DMDA  Input Parameters:+  da - the distributed array object.  xmin,xmax - extremes in the x direction.  ymin,ymax - extremes in the y direction (use NULL for 1 dimensional problems)-  zmin,zmax - extremes in the z direction (use NULL for 1 or 2 dimensional problems)  Level: beginner.seealso: DMSetCoordinates(), DMGetCoordinates(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()@*/PetscErrorCode  DMDASetUniformCoordinates(DM da,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax){  MPI_Comm         comm;  PetscSection     section;  DM               cda;  DMBoundaryType   bx,by,bz;  Vec              xcoor;  PetscScalar      *coors;  PetscReal        hx,hy,hz_;  PetscInt         i,j,k,M,N,P,istart,isize,jstart,jsize,kstart,ksize,dim,cnt;  PetscErrorCode   ierr;  PetscFunctionBegin;  PetscValidHeaderSpecific(da,DM_CLASSID,1);  ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr);  if (xmax <= xmin) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"xmax must be larger than xmin %g %g",(double)xmin,(double)xmax);  if ((ymax <= ymin) && (dim > 1)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %g %g",(double)ymin,(double)ymax);  if ((zmax <= zmin) && (dim > 2)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %g %g",(double)zmin,(double)zmax);  ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr);  ierr = DMGetDefaultSection(da,&section);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&istart,&jstart,&kstart,&isize,&jsize,&ksize);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da, &cda);CHKERRQ(ierr);  if (section) {    /* This would be better as a vector, but this is compatible */    PetscInt numComp[3]      = {1, 1, 1};    PetscInt numVertexDof[3] = {1, 1, 1};    ierr = DMDASetFieldName(cda, 0, "x");CHKERRQ(ierr);    if (dim > 1) {ierr = DMDASetFieldName(cda, 1, "y");CHKERRQ(ierr);}    if (dim > 2) {ierr = DMDASetFieldName(cda, 2, "z");CHKERRQ(ierr);}    ierr = DMDACreateSection(cda, numComp, numVertexDof, NULL, NULL);CHKERRQ(ierr);  }  ierr = DMCreateGlobalVector(cda, &xcoor);CHKERRQ(ierr);  if (section) {    PetscSection csection;    PetscInt     vStart, vEnd;    ierr = DMGetDefaultGlobalSection(cda,&csection);CHKERRQ(ierr);    ierr = VecGetArray(xcoor,&coors);CHKERRQ(ierr);    ierr = DMDAGetHeightStratum(da, dim, &vStart, &vEnd);CHKERRQ(ierr);    if (bx == DM_BOUNDARY_PERIODIC) hx  = (xmax-xmin)/(M+1);    else                              hx  = (xmax-xmin)/(M ? M : 1);    if (by == DM_BOUNDARY_PERIODIC) hy  = (ymax-ymin)/(N+1);    else                              hy  = (ymax-ymin)/(N ? N : 1);    if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax-zmin)/(P+1);    else                              hz_ = (zmax-zmin)/(P ? P : 1);    switch (dim) {    case 1:      for (i = 0; i < isize+1; ++i) {        PetscInt v = i+vStart, dof, off;        ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr);        ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr);        if (off >= 0) {          coors[off] = xmin + hx*(i+istart);        }      }      break;    case 2:      for (j = 0; j < jsize+1; ++j) {        for (i = 0; i < isize+1; ++i) {          PetscInt v = j*(isize+1)+i+vStart, dof, off;          ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr);          ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr);          if (off >= 0) {            coors[off+0] = xmin + hx*(i+istart);            coors[off+1] = ymin + hy*(j+jstart);          }        }      }      break;    case 3:      for (k = 0; k < ksize+1; ++k) {        for (j = 0; j < jsize+1; ++j) {          for (i = 0; i < isize+1; ++i) {            PetscInt v = (k*(jsize+1)+j)*(isize+1)+i+vStart, dof, off;            ierr = PetscSectionGetDof(csection, v, &dof);CHKERRQ(ierr);            ierr = PetscSectionGetOffset(csection, v, &off);CHKERRQ(ierr);            if (off >= 0) {              coors[off+0] = xmin + hx*(i+istart);              coors[off+1] = ymin + hy*(j+jstart);              coors[off+2] = zmin + hz_*(k+kstart);//.........这里部分代码省略.........

int main(int argc, char **argv){    AppCtx         ctx;    DM             dm;    TS             ts;    Vec            u, r;    PetscReal      t       = 0.0;    PetscReal      L2error = 0.0;    PetscErrorCode ierr;    ierr = PetscInitialize(&argc, &argv, NULL, help);    CHKERRQ(ierr);    ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);    CHKERRQ(ierr);    ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);    CHKERRQ(ierr);    ierr = DMSetApplicationContext(dm, &ctx);    CHKERRQ(ierr);    ierr = PetscMalloc1(1, &ctx.exactFuncs);    CHKERRQ(ierr);    ierr = SetupDiscretization(dm, &ctx);    CHKERRQ(ierr);    ierr = DMCreateGlobalVector(dm, &u);    CHKERRQ(ierr);    ierr = PetscObjectSetName((PetscObject) u, "temperature");    CHKERRQ(ierr);    ierr = VecDuplicate(u, &r);    CHKERRQ(ierr);    ierr = TSCreate(PETSC_COMM_WORLD, &ts);    CHKERRQ(ierr);    ierr = TSSetDM(ts, dm);    CHKERRQ(ierr);    ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);    CHKERRQ(ierr);    ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);    CHKERRQ(ierr);    ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);    CHKERRQ(ierr);    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);    CHKERRQ(ierr);    ierr = TSSetFromOptions(ts);    CHKERRQ(ierr);    ierr = DMProjectFunction(dm, t, ctx.exactFuncs, NULL, INSERT_ALL_VALUES, u);    CHKERRQ(ierr);    ierr = TSSolve(ts, u);    CHKERRQ(ierr);    ierr = TSGetTime(ts, &t);    CHKERRQ(ierr);    ierr = DMComputeL2Diff(dm, t, ctx.exactFuncs, NULL, u, &L2error);    CHKERRQ(ierr);    if (L2error < 1.0e-11) {        ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11/n");        CHKERRQ(ierr);    }    else                   {        ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g/n", L2error);        CHKERRQ(ierr);    }    ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");    CHKERRQ(ierr);    ierr = VecDestroy(&u);    CHKERRQ(ierr);    ierr = VecDestroy(&r);    CHKERRQ(ierr);    ierr = TSDestroy(&ts);    CHKERRQ(ierr);    ierr = DMDestroy(&dm);    CHKERRQ(ierr);    ierr = PetscFree(ctx.exactFuncs);    CHKERRQ(ierr);    ierr = PetscFinalize();    return ierr;}

int main(int argc,char **argv){  PetscErrorCode ierr;  KSP            ksp;  PC             pc;  Vec            x,b;  DM             da;  Mat            A;  PetscInt       dof=1;  PetscBool      flg;  PetscScalar    zero=0.0;  PetscInitialize(&argc,&argv,(char *)0,help);  ierr = PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDACreate(PETSC_COMM_WORLD,&da);CHKERRQ(ierr);  ierr = DMDASetDim(da,3);CHKERRQ(ierr);  ierr = DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);CHKERRQ(ierr);  ierr = DMDASetStencilType(da,DMDA_STENCIL_STAR);CHKERRQ(ierr);  ierr = DMDASetSizes(da,3,3,3);CHKERRQ(ierr);  ierr = DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);  ierr = DMDASetDof(da,dof);CHKERRQ(ierr);  ierr = DMDASetStencilWidth(da,1);CHKERRQ(ierr);  ierr = DMDASetOwnershipRanges(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);  ierr = DMSetFromOptions(da);CHKERRQ(ierr);  ierr = DMSetUp(da);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&b);CHKERRQ(ierr);  ierr = DMCreateMatrix(da,MATAIJ,&A);CHKERRQ(ierr);  ierr = VecSet(b,zero);CHKERRQ(ierr);  /* Test sbaij matrix */  flg = PETSC_FALSE;  ierr = PetscOptionsGetBool(PETSC_NULL,"-test_sbaij",&flg,PETSC_NULL);CHKERRQ(ierr);  if (flg) {    Mat sA;    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);    ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);    ierr = MatDestroy(&A);CHKERRQ(ierr);    A = sA;  }  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);  ierr = PCSetDM(pc,(DM)da);CHKERRQ(ierr);  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);  /* check final residual */  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);CHKERRQ(ierr);  if (flg){    Vec            b1;    PetscReal      norm;    ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);    ierr = VecDuplicate(b,&b1);CHKERRQ(ierr);    ierr = MatMult(A,x,b1);CHKERRQ(ierr);    ierr = VecAXPY(b1,-1.0,b);CHKERRQ(ierr);    ierr = VecNorm(b1,NORM_2,&norm);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD,"Final residual %g/n",norm);CHKERRQ(ierr);    ierr = VecDestroy(&b1);CHKERRQ(ierr);  }  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = VecDestroy(&b);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){  PetscMPIInt    rank,size;  PetscInt       M = 14,time_steps = 20,w=1,s=1,localsize,j,i,mybase,myend,globalsize;  PetscErrorCode ierr;  DM             da;  Vec            global,local;  PetscScalar    *globalptr,*localptr;  PetscReal      h,k;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr);  /* Set up the array */  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,M,w,s,NULL,&da);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);  /* Make copy of local array for doing updates */  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  /* determine starting point of each processor */  ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr);  /* Initialize the Array */  ierr = VecGetLocalSize (global,&globalsize);CHKERRQ(ierr);  ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr);  for (i=0; i<globalsize; i++) {    j = i + mybase;    globalptr[i] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 4+4;  }  ierr = VecRestoreArray(global,&localptr);CHKERRQ(ierr);  /* Assign Parameters */  h= 1.0/M;  k= h*h/2.2;  ierr = VecGetLocalSize(local,&localsize);CHKERRQ(ierr);  for (j=0; j<time_steps; j++) {    /* Global to Local */    ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);    /*Extract local array */    ierr = VecGetArray(local,&localptr);CHKERRQ(ierr);    ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr);    /* Update Locally - Make array of new values */    /* Note: I don't do anything for the first and last entry */    for (i=1; i< localsize-1; i++) {      globalptr[i-1] = localptr[i] + (k/(h*h)) * (localptr[i+1]-2.0*localptr[i]+localptr[i-1]);    }    ierr = VecRestoreArray (global,&globalptr);CHKERRQ(ierr);    ierr = VecRestoreArray(local,&localptr);CHKERRQ(ierr);    /* View Wave */    /* Set Up Display to Show Heat Graph */#if defined(PETSC_USE_SOCKET_VIEWER)    ierr = VecView(global,PETSC_VIEWER_SOCKET_WORLD);CHKERRQ(ierr);#endif  }  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){  TS             ts;                   /* nonlinear solver */  Vec            u,r;                  /* solution, residual vectors */  Mat            J,Jmf = PETSC_NULL;   /* Jacobian matrices */  PetscInt       maxsteps = 1000;      /* iterations for convergence */  PetscErrorCode ierr;  DM             da;  PetscReal      dt;  AppCtx         user;              /* user-defined work context */  SNES           snes;  PetscInt       Jtype; /* Jacobian type                            0: user provide Jacobian;                            1: slow finite difference;                            2: fd with coloring; */  PetscInitialize(&argc,&argv,(char *)0,help);  /* Initialize user application context */  user.da            = PETSC_NULL;  user.nstencilpts   = 5;  user.c             = -30.0;  user.boundary      = 0; /* 0: Drichlet BC; 1: Neumann BC */  user.viewJacobian  = PETSC_FALSE;  ierr = PetscOptionsGetInt(PETSC_NULL,"-nstencilpts",&user.nstencilpts,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-boundary",&user.boundary,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsHasName(PETSC_NULL,"-viewJacobian",&user.viewJacobian);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  if (user.nstencilpts == 5){    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);  } else if (user.nstencilpts == 9){    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);  } else {    SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"nstencilpts %d is not supported",user.nstencilpts);  }  user.da = da;  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Extract global vectors from DMDA;   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);  ierr = VecDuplicate(u,&r);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);  ierr = TSSetDM(ts,da);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,r,FormIFunction,&user);CHKERRQ(ierr);  ierr = TSSetDuration(ts,maxsteps,1.0);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = FormInitialSolution(u,&user);CHKERRQ(ierr);  ierr = TSSetSolution(ts,u);CHKERRQ(ierr);  dt   = .01;  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   Set Jacobian evaluation routine  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateMatrix(da,MATAIJ,&J);CHKERRQ(ierr);  Jtype = 0;  ierr = PetscOptionsGetInt(PETSC_NULL, "-Jtype",&Jtype,PETSC_NULL);CHKERRQ(ierr);  if (Jtype == 0){ /* use user provided Jacobian evaluation routine */    if (user.nstencilpts != 5) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"user Jacobian routine FormIJacobian() does not support nstencilpts=%D",user.nstencilpts);    ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr);  } else { /* use finite difference Jacobian J as preconditioner and '-snes_mf_operator' for Mat*vec */    ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);    ierr = MatCreateSNESMF(snes,&Jmf);CHKERRQ(ierr);    if (Jtype == 1){ /* slow finite difference J; */      ierr = SNESSetJacobian(snes,Jmf,J,SNESDefaultComputeJacobian,PETSC_NULL);CHKERRQ(ierr);    } else if (Jtype == 2){ /* Use coloring to compute  finite difference J efficiently */      ierr = SNESSetJacobian(snes,Jmf,J,SNESDefaultComputeJacobianColor,0);CHKERRQ(ierr);    } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Jtype is not supported");  }  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   Sets various TS parameters from user options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSolve(ts,u);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = MatDestroy(&Jmf);CHKERRQ(ierr);  ierr = VecDestroy(&u);CHKERRQ(ierr);  ierr = VecDestroy(&r);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);//.........这里部分代码省略.........

int main(int argc,char **argv){  TS             ts;                  /* nonlinear solver */  Vec            U;                   /* solution, residual vectors */  PetscErrorCode ierr;  DM             da;  AppCtx         appctx;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Initialize program     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  PetscInitialize(&argc,&argv,(char*)0,help);  appctx.epsilon = 1.0e-3;  appctx.delta   = 1.0;  appctx.alpha   = 10.0;  appctx.beta    = 4.0;  appctx.gamma   = 1.0;  appctx.kappa   = .75;  appctx.lambda  = 1.0;  appctx.mu      = 100.;  appctx.cstar   = .2;  appctx.upwind  = PETSC_TRUE;  ierr = PetscOptionsGetScalar(NULL,"-delta",&appctx.delta,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-upwind",&appctx.upwind,NULL);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,-8,2,1,NULL,&da);CHKERRQ(ierr);  ierr = DMDASetFieldName(da,0,"rho");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,1,"c");CHKERRQ(ierr);  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Extract global vectors from DMDA; then duplicate for remaining     vectors that are the same types   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr);  ierr = TSSetDM(ts,da);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,IFunction,&appctx);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = InitialConditions(da,U);CHKERRQ(ierr);  ierr = TSSetSolution(ts,U);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set solver options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetInitialTimeStep(ts,0.0,.0001);CHKERRQ(ierr);  ierr = TSSetDuration(ts,PETSC_DEFAULT,1.0);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSolve(ts,U);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.  All PETSc objects should be destroyed when they     are no longer needed.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = VecDestroy(&U);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  PetscFunctionReturn(0);}

/*     DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set*/PetscErrorCode  DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2){  PetscErrorCode          ierr;  PetscContainer          isnes;  DM_SNESVI               *dmsnesvi1;  Vec                     finemarked,coarsemarked;  IS                      inactive;  VecScatter              inject;  const PetscInt          *index;  PetscInt                n,k,cnt = 0,rstart,*coarseindex;  PetscScalar             *marked;  PetscFunctionBegin;  ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr);  if (!isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed VI data structure is missing");  ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr);  /* get the original coarsen */  ierr = (*dmsnesvi1->coarsen)(dm1,comm,dm2);CHKERRQ(ierr);  /* not sure why this extra reference is needed, but without the dm2 disappears too early */  ierr = PetscObjectReference((PetscObject)*dm2);CHKERRQ(ierr);  /* need to set back global vectors in order to use the original injection */  ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr);  dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;  ierr = DMCreateGlobalVector(dm1,&finemarked);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(*dm2,&coarsemarked);CHKERRQ(ierr);  /*     fill finemarked with locations of inactive points  */  ierr = ISGetIndices(dmsnesvi1->inactive,&index);CHKERRQ(ierr);  ierr = ISGetLocalSize(dmsnesvi1->inactive,&n);CHKERRQ(ierr);  ierr = VecSet(finemarked,0.0);CHKERRQ(ierr);  for (k=0;k<n;k++){      ierr = VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);CHKERRQ(ierr);  }  ierr = VecAssemblyBegin(finemarked);CHKERRQ(ierr);  ierr = VecAssemblyEnd(finemarked);CHKERRQ(ierr);  ierr = DMCreateInjection(*dm2,dm1,&inject);CHKERRQ(ierr);  ierr = VecScatterBegin(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);  ierr = VecScatterEnd(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);  ierr = VecScatterDestroy(&inject);CHKERRQ(ierr);  /*     create index set list of coarse inactive points from coarsemarked  */  ierr = VecGetLocalSize(coarsemarked,&n);CHKERRQ(ierr);  ierr = VecGetOwnershipRange(coarsemarked,&rstart,PETSC_NULL);CHKERRQ(ierr);  ierr = VecGetArray(coarsemarked,&marked);CHKERRQ(ierr);  for (k=0; k<n; k++) {    if (marked[k] != 0.0) cnt++;  }  ierr = PetscMalloc(cnt*sizeof(PetscInt),&coarseindex);CHKERRQ(ierr);  cnt  = 0;  for (k=0; k<n; k++) {    if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;  }  ierr = VecRestoreArray(coarsemarked,&marked);CHKERRQ(ierr);  ierr = ISCreateGeneral(((PetscObject)coarsemarked)->comm,cnt,coarseindex,PETSC_OWN_POINTER,&inactive);CHKERRQ(ierr);  ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr);  dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;  ierr = DMSetVI(*dm2,inactive);CHKERRQ(ierr);  ierr = VecDestroy(&finemarked);CHKERRQ(ierr);  ierr = VecDestroy(&coarsemarked);CHKERRQ(ierr);  ierr = ISDestroy(&inactive);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc,char **argv){  TS             ts;                  /* nonlinear solver */  Vec            U;                   /* solution, residual vectors */  Mat            J;                   /* Jacobian matrix */  PetscInt       maxsteps = 1000;  PetscErrorCode ierr;  DM             da;  AppCtx         user;  PetscInt       i;  char           Name[16];  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Initialize program     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  PetscInitialize(&argc,&argv,(char*)0,help);  user.N = 1;  ierr   = PetscOptionsGetInt(NULL,"-N",&user.N,NULL);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMDACreate1d(PETSC_COMM_WORLD, DMDA_BOUNDARY_MIRROR,-8,user.N,1,NULL,&da);CHKERRQ(ierr);  for (i=0; i<user.N; i++) {    ierr = PetscSNPrintf(Name,16,"Void size %d",(int)(i+1));    ierr = DMDASetFieldName(da,i,Name);CHKERRQ(ierr);  }  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   Extract global vectors from DMDA; then duplicate for remaining     vectors that are the same types   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr);  ierr = DMCreateMatrix(da,MATAIJ,&J);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);  ierr = TSSetDM(ts,da);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);  ierr = TSSetIJacobian(ts,J,J,IJacobian,&user);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = InitialConditions(da,U);CHKERRQ(ierr);  ierr = TSSetSolution(ts,U);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set solver options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetInitialTimeStep(ts,0.0,.001);CHKERRQ(ierr);  ierr = TSSetDuration(ts,maxsteps,1.0);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSolve(ts,U);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.  All PETSc objects should be destroyed when they     are no longer needed.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = VecDestroy(&U);CHKERRQ(ierr);  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  PetscFunctionReturn(0);}

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);//.........这里部分代码省略.........

int main(int argc, char **argv){  MPI_Comm       comm;  DM             dm;  Vec            v, nv, rv, coord;  PetscBool      test_read = PETSC_FALSE, verbose = PETSC_FALSE, flg;  PetscViewer    hdf5Viewer;  PetscInt       dim         = 2;  PetscInt       numFields   = 1;  PetscInt       numBC       = 0;  PetscInt       numComp[1]  = {2};  PetscInt       numDof[3]   = {2, 0, 0};  PetscInt       bcFields[1] = {0};  IS             bcPoints[1] = {NULL};  PetscSection   section;  PetscReal      norm;  PetscErrorCode ierr;  ierr = PetscInitialize(&argc, &argv, (char *) 0, help);CHKERRQ(ierr);  comm = PETSC_COMM_WORLD;  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");CHKERRQ(ierr);  ierr = PetscOptionsBool("-test_read","Test reading from the HDF5 file","",PETSC_FALSE,&test_read,NULL);CHKERRQ(ierr);  ierr = PetscOptionsBool("-verbose","print the Vecs","",PETSC_FALSE,&verbose,NULL);CHKERRQ(ierr);  ierr = PetscOptionsInt("-dim","the dimension of the problem","",2,&dim,NULL);CHKERRQ(ierr);  ierr = PetscOptionsEnd();  ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, &dm);CHKERRQ(ierr);  ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);  numDof[0] = dim;  ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, &section);CHKERRQ(ierr);  ierr = DMSetDefaultSection(dm, section);CHKERRQ(ierr);  ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);  ierr = DMSetUseNatural(dm, PETSC_TRUE);CHKERRQ(ierr);  {    DM dmDist;    ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr);    if (dmDist) {      ierr = DMDestroy(&dm);CHKERRQ(ierr);      dm   = dmDist;    }  }  ierr = DMCreateGlobalVector(dm, &v);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) v, "V");CHKERRQ(ierr);  ierr = DMGetCoordinates(dm, &coord);CHKERRQ(ierr);  ierr = VecCopy(coord, v);CHKERRQ(ierr);  if (verbose) {    PetscInt size, bs;    ierr = VecGetSize(v, &size);CHKERRQ(ierr);    ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d/tblock size=%d/n", size, bs);CHKERRQ(ierr);    ierr = VecView(v, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  ierr = DMCreateGlobalVector(dm, &nv);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) nv, "NV");CHKERRQ(ierr);  ierr = DMPlexGlobalToNaturalBegin(dm, v, nv);CHKERRQ(ierr);  ierr = DMPlexGlobalToNaturalEnd(dm, v, nv);CHKERRQ(ierr);  if (verbose) {    PetscInt size, bs;    ierr = VecGetSize(nv, &size);CHKERRQ(ierr);    ierr = VecGetBlockSize(nv, &bs);CHKERRQ(ierr);    ierr = PetscPrintf(PETSC_COMM_WORLD, "====  V in natural ordering. size==%d/tblock size=%d/n", size, bs);CHKERRQ(ierr);    ierr = VecView(nv, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  ierr = VecViewFromOptions(v, NULL, "-global_vec_view");CHKERRQ(ierr);  if (test_read) {    ierr = DMCreateGlobalVector(dm, &rv);CHKERRQ(ierr);    ierr = PetscObjectSetName((PetscObject) rv, "V");CHKERRQ(ierr);    /* Test native read */    ierr = PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);CHKERRQ(ierr);    ierr = PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);CHKERRQ(ierr);    ierr = VecLoad(rv, hdf5Viewer);CHKERRQ(ierr);    ierr = PetscViewerDestroy(&hdf5Viewer);CHKERRQ(ierr);    if (verbose) {      PetscInt size, bs;      ierr = VecGetSize(rv, &size);CHKERRQ(ierr);      ierr = VecGetBlockSize(rv, &bs);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d/tblock size=%d/n", size, bs);CHKERRQ(ierr);      ierr = VecView(rv, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    }    ierr = VecEqual(rv, v, &flg);CHKERRQ(ierr);    if (flg) {      ierr = PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal/n");CHKERRQ(ierr);    } else {      ierr = PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal/n/n");CHKERRQ(ierr);      ierr = VecAXPY(rv, -1.0, v);CHKERRQ(ierr);      ierr = VecNorm(rv, NORM_INFINITY, &norm);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g/n", (double) norm);CHKERRQ(ierr);      ierr = VecView(rv, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);    }//.........这里部分代码省略.........

int main(int argc, char **argv){  SNES           snes;                 /* nonlinear solver */  DM             dm;                   /* problem definition */  Vec            u,r;                  /* solution, residual vectors */  Mat            A,J;                  /* Jacobian matrix */  MatNullSpace   nullSpace;            /* May be necessary for pressure */  AppCtx         user;                 /* user-defined work context */  JacActionCtx   userJ;                /* context for Jacobian MF action */  PetscInt       its;                  /* iterations for convergence */  PetscReal      error         = 0.0;  /* L_2 error in the solution */  PetscInt       numComponents = 0, f;  PetscErrorCode ierr;  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);  ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);  ierr = SetupElement(dm, &user);CHKERRQ(ierr);  for (f = 0; f < NUM_FIELDS; ++f) {    PetscInt numComp;    ierr = PetscFEGetNumComponents(user.fe[f], &numComp);CHKERRQ(ierr);    numComponents += numComp;  }  ierr = PetscMalloc(NUM_FIELDS * sizeof(void (*)(const PetscReal[], PetscScalar *)), &user.exactFuncs);CHKERRQ(ierr);  user.fem.bcFuncs = (void (**)(const PetscReal[], PetscScalar *)) user.exactFuncs;  ierr = SetupExactSolution(dm, &user);CHKERRQ(ierr);  ierr = SetupSection(dm, &user);CHKERRQ(ierr);  ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);  ierr = VecDuplicate(u, &r);CHKERRQ(ierr);  ierr = DMSetMatType(dm,MATAIJ);CHKERRQ(ierr);  ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);  if (user.jacobianMF) {    PetscInt M, m, N, n;    ierr = MatGetSize(J, &M, &N);CHKERRQ(ierr);    ierr = MatGetLocalSize(J, &m, &n);CHKERRQ(ierr);    ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr);    ierr = MatSetSizes(A, m, n, M, N);CHKERRQ(ierr);    ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr);    ierr = MatSetUp(A);CHKERRQ(ierr);    ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);CHKERRQ(ierr);    userJ.dm   = dm;    userJ.J    = J;    userJ.user = &user;    ierr = DMCreateLocalVector(dm, &userJ.u);CHKERRQ(ierr);    ierr = DMPlexProjectFunctionLocal(dm, user.fe, user.exactFuncs, INSERT_BC_VALUES, userJ.u);CHKERRQ(ierr);    ierr = MatShellSetContext(A, &userJ);CHKERRQ(ierr);  } else {    A = J;  }  ierr = CreatePressureNullSpace(dm, &user, &nullSpace);CHKERRQ(ierr);  ierr = MatSetNullSpace(J, nullSpace);CHKERRQ(ierr);  if (A != J) {    ierr = MatSetNullSpace(A, nullSpace);CHKERRQ(ierr);  }  ierr = DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);CHKERRQ(ierr);  ierr = DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);CHKERRQ(ierr);  ierr = SNESSetJacobian(snes, A, J, NULL, NULL);CHKERRQ(ierr);  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);  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);//.........这里部分代码省略.........

int main(int argc,char **argv) {  PetscErrorCode ierr;  DM             da, da_after;  SNES           snes;  Vec            u_initial, u;  PoissonCtx     user;  SNESConvergedReason reason;  int            snesits;  double         lflops,flops;  DMDALocalInfo  info;  PetscInitialize(&argc,&argv,NULL,help);  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"el_",                           "elasto-plastic torsion solver options",""); CHKERRQ(ierr);  ierr = PetscOptionsReal("-C","f(x,y)=2C is source term",                          "elasto.c",C,&C,NULL); CHKERRQ(ierr);  ierr = PetscOptionsEnd(); CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD,      DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,      3,3,                       // override with -da_grid_x,_y      PETSC_DECIDE,PETSC_DECIDE, // num of procs in each dim      1,1,NULL,NULL,             // dof = 1 and stencil width = 1      &da);CHKERRQ(ierr);  ierr = DMSetFromOptions(da); CHKERRQ(ierr);  ierr = DMSetUp(da); CHKERRQ(ierr);  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,-1.0,-1.0);CHKERRQ(ierr);  user.cx = 1.0;  user.cy = 1.0;  user.cz = 1.0;  user.g_bdry = &zero;  user.f_rhs = &f_fcn;  user.addctx = NULL;  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);  ierr = SNESSetType(snes,SNESVINEWTONRSLS);CHKERRQ(ierr);  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds);CHKERRQ(ierr);  // reuse residual and jacobian from ch6/  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,             (DMDASNESFunction)Poisson2DFunctionLocal,&user); CHKERRQ(ierr);  ierr = DMDASNESSetJacobianLocal(da,             (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);  // initial iterate is zero  ierr = DMCreateGlobalVector(da,&u_initial);CHKERRQ(ierr);  ierr = VecSet(u_initial,0.0); CHKERRQ(ierr);  /* solve; then get solution and DM after solution*/  ierr = SNESSolve(snes,NULL,u_initial);CHKERRQ(ierr);  ierr = VecDestroy(&u_initial); CHKERRQ(ierr);  ierr = DMDestroy(&da); CHKERRQ(ierr);  ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);  ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr); /* do not destroy u */  /* performance measures */  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);  if (reason <= 0) {      ierr = PetscPrintf(PETSC_COMM_WORLD,          "WARNING: SNES not converged ... use -snes_converged_reason to check/n"); CHKERRQ(ierr);  }  ierr = SNESGetIterationNumber(snes,&snesits); CHKERRQ(ierr);  ierr = PetscGetFlops(&lflops); CHKERRQ(ierr);  ierr = MPI_Allreduce(&lflops,&flops,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,      "done on %4d x %4d grid; total flops = %.3e; SNES iterations %d/n",      info.mx,info.my,flops,snesits); CHKERRQ(ierr);  SNESDestroy(&snes);  return PetscFinalize();}

int main(int argc,char **argv){  TS                ts;         /* time integrator */  TSAdapt           adapt;  Vec               X;          /* solution vector */  Mat               J;          /* Jacobian matrix */  PetscInt          steps,maxsteps,ncells,xs,xm,i;  PetscErrorCode    ierr;  PetscReal         ftime,dt;  char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp",thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";  struct _User      user;  TSConvergedReason reason;  PetscBool         showsolutions = PETSC_FALSE;  char              **snames,*names;  Vec               lambda;     /* used with TSAdjoint for sensitivities */  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");CHKERRQ(ierr);  ierr = PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);CHKERRQ(ierr);  ierr = PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);CHKERRQ(ierr);  user.pressure = 1.01325e5;    /* Pascal */  ierr = PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);CHKERRQ(ierr);  user.Tini   = 1550;  ierr = PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);CHKERRQ(ierr);  user.diffus = 100;  ierr = PetscOptionsReal("-diffus","Diffusion constant","",user.diffus,&user.diffus,NULL);CHKERRQ(ierr);  ierr = PetscOptionsBool("-draw_solution","Plot the solution for each cell","",showsolutions,&showsolutions,NULL);CHKERRQ(ierr);  user.diffusion = PETSC_TRUE;  ierr = PetscOptionsBool("-diffusion","Have diffusion","",user.diffusion,&user.diffusion,NULL);CHKERRQ(ierr);  user.reactions = PETSC_TRUE;  ierr = PetscOptionsBool("-reactions","Have reactions","",user.reactions,&user.reactions,NULL);CHKERRQ(ierr);  ierr = PetscOptionsEnd();CHKERRQ(ierr);  ierr = TC_initChem(chemfile, thermofile, 0, 1.0);TCCHKERRQ(ierr);  user.Nspec = TC_getNspec();  user.Nreac = TC_getNreac();  ierr    = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,-1,user.Nspec+1,1,NULL,&user.dm);CHKERRQ(ierr);  ierr    = DMDAGetInfo(user.dm,NULL,&ncells,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);  user.dx = 1.0/ncells;  /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */  ierr    = DMDASetUniformCoordinates(user.dm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  /* set the names of each field in the DMDA based on the species name */  ierr = PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);CHKERRQ(ierr);  ierr = PetscStrcpy(names,"Temp");CHKERRQ(ierr);  TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);CHKERRQ(ierr);  ierr = PetscMalloc1((user.Nspec+2),&snames);CHKERRQ(ierr);  for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;  snames[user.Nspec+1] = NULL;  ierr = DMDASetFieldNames(user.dm,(const char * const *)snames);CHKERRQ(ierr);  ierr = PetscFree(snames);CHKERRQ(ierr);  ierr = PetscFree(names);CHKERRQ(ierr);  ierr = DMCreateMatrix(user.dm,&J);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(user.dm,&X);CHKERRQ(ierr);  ierr = PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetDM(ts,user.dm);CHKERRQ(ierr);  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);  ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr);  ierr = TSARKIMEXSetType(ts,TSARKIMEX4);CHKERRQ(ierr);  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);  ierr = TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);CHKERRQ(ierr);  ftime    = 1.0;  maxsteps = 10000;  ierr     = TSSetDuration(ts,maxsteps,ftime);CHKERRQ(ierr);  ierr     = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr);  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);  dt   = 1e-10;                 /* Initial time step */  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);  ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);  ierr = TSAdaptSetStepLimits(adapt,1e-12,1e-4);CHKERRQ(ierr); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */  ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr);            /* Retry step an unlimited number of times */  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Pass information to graphical monitoring routine   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  if (showsolutions) {    ierr = DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);    for (i=xs;i<xs+xm;i++) {      ierr = MonitorCell(ts,&user,i);CHKERRQ(ierr);    }  }  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set runtime options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *///.........这里部分代码省略.........

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){  PetscMPIInt      rank;  PetscInt         M = -10,N = -8;  PetscErrorCode   ierr;  PetscBool        flg = PETSC_FALSE;  DM               da;  PetscViewer      viewer;  Vec              local,global;  PetscScalar      value;  DMBoundaryType   bx    = DM_BOUNDARY_NONE,by = DM_BOUNDARY_NONE;  DMDAStencilType  stype = DMDA_STENCIL_BOX;#if defined(PETSC_HAVE_MATLAB_ENGINE)  PetscViewer      mviewer;  PetscMPIInt      size;#endif  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",300,0,300,300,&viewer);CHKERRQ(ierr);#if defined(PETSC_HAVE_MATLAB_ENGINE)  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);  if (size == 1) {    ierr = PetscViewerMatlabOpen(PETSC_COMM_WORLD,"tmp.mat",FILE_MODE_WRITE,&mviewer);CHKERRQ(ierr);  }#endif  ierr = PetscOptionsGetBool(NULL,NULL,"-star_stencil",&flg,NULL);CHKERRQ(ierr);  if (flg) stype = DMDA_STENCIL_STAR;  /* Create distributed array and get vectors */  ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  value = -3.0;  ierr  = VecSet(global,value);CHKERRQ(ierr);  ierr  = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr  = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr  = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  value = rank+1;  ierr  = VecScale(local,value);CHKERRQ(ierr);  ierr  = DMLocalToGlobalBegin(da,local,ADD_VALUES,global);CHKERRQ(ierr);  ierr  = DMLocalToGlobalEnd(da,local,ADD_VALUES,global);CHKERRQ(ierr);  flg  = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL, "-view_global", &flg,NULL);CHKERRQ(ierr);  if (flg) { /* view global vector in natural ordering */    ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  ierr = DMView(da,viewer);CHKERRQ(ierr);  ierr = VecView(global,viewer);CHKERRQ(ierr);#if defined(PETSC_HAVE_MATLAB_ENGINE)  if (size == 1) {    ierr = DMView(da,mviewer);CHKERRQ(ierr);    ierr = VecView(global,mviewer);CHKERRQ(ierr);  }#endif  /* Free memory */#if defined(PETSC_HAVE_MATLAB_ENGINE)  if (size == 1) {    ierr = PetscViewerDestroy(&mviewer);CHKERRQ(ierr);  }#endif  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return ierr;}

PetscErrorCode CreateStructures(DM da, UserContext *user){  const PetscInt *necon;  PetscInt       ne,nc;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = DMDAGetElements(da,&ne,&nc,&necon);CHKERRQ(ierr);  ierr = DMDARestoreElements(da,&ne,&nc,&necon);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.rho);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.rho_u);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.rho_v);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.rho_e);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.p);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.u);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.v);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_n.t);CHKERRQ(ierr);  ierr = VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho);CHKERRQ(ierr);  ierr = VecSetSizes(user->sol_phi.rho, ne, PETSC_DECIDE);CHKERRQ(ierr);  ierr = VecSetType(user->sol_phi.rho,VECMPI);CHKERRQ(ierr);  ierr = VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_u);CHKERRQ(ierr);  ierr = VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_v);CHKERRQ(ierr);  ierr = VecDuplicate(user->sol_phi.rho, &user->sol_phi.u);CHKERRQ(ierr);  ierr = VecDuplicate(user->sol_phi.rho, &user->sol_phi.v);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.rho);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.rho_u);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.rho_v);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.rho_e);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.p);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.u);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->sol_np1.v);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->mu);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da, &user->kappa);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc, char **argv){  MPI_Comm          comm;  PetscMPIInt       rank;  PetscErrorCode    ierr;  User              user;  PetscLogDouble       v1, v2;  PetscInt          nplot = 0;  char              filename1[2048], fileName[2048];  PetscBool         set = PETSC_FALSE;  PetscInt          steps_output;  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);  comm = PETSC_COMM_WORLD;  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);  ierr = PetscNew(&user);CHKERRQ(ierr);  ierr = PetscNew(&user->algebra);CHKERRQ(ierr);  ierr = PetscNew(&user->model);CHKERRQ(ierr);  ierr = PetscNew(&user->model->physics);CHKERRQ(ierr);  Algebra   algebra = user->algebra;  ierr = LoadOptions(comm, user);CHKERRQ(ierr);  ierr = PetscTime(&v1);CHKERRQ(ierr);  ierr = CreateMesh(comm, user);CHKERRQ(ierr);  ierr = PetscTime(&v2);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,		       "Read and Distribute mesh takes %f sec /n", v2 - v1);CHKERRQ(ierr);  ierr = SetUpLocalSpace(user);CHKERRQ(ierr); //Set up the dofs of each element  ierr = ConstructGeometryFVM(&user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);  ierr = LimiterSetup(user);CHKERRQ(ierr);  if(user->output_solution){  // the output file options    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,0,"Options for output solution",0);CHKERRQ(ierr);    ierr = PetscOptionsString("-solutionfile", "solution file", "AeroSim.c", filename1,filename1, 2048, &set);CHKERRQ(ierr);    if(!set){SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,"please use option -solutionfile to specify solution file name /n");}    ierr = PetscOptionsInt("-steps_output", "the number of time steps between two outputs", "", steps_output, &steps_output, &set);CHKERRQ(ierr);    if(!set){ steps_output = 1;}    ierr = PetscOptionsEnd();CHKERRQ(ierr);  }  if (user->TimeIntegralMethod == EXPLICITMETHOD) {    if(user->myownexplicitmethod){      ierr = PetscPrintf(PETSC_COMM_WORLD,"Using the fully explicit method based on my own routing/n");CHKERRQ(ierr);      user->current_time = user->initial_time;      user->current_step = 1;      ierr = DMCreateGlobalVector(user->dm, &algebra->solution);CHKERRQ(ierr);      ierr = PetscObjectSetName((PetscObject) algebra->solution, "solution");CHKERRQ(ierr);      ierr = SetInitialCondition(user->dm, algebra->solution, user);CHKERRQ(ierr);      ierr = VecDuplicate(algebra->solution, &algebra->fn);CHKERRQ(ierr);      ierr = VecDuplicate(algebra->solution, &algebra->oldsolution);CHKERRQ(ierr);      if(user->Explicit_RK2){        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the second order Runge Kutta method /n");CHKERRQ(ierr);      }else{        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the first order forward Euler method /n");CHKERRQ(ierr);      }      nplot = 0; //the plot step      while(user->current_time < (user->final_time - 0.05 * user->dt)){        user->current_time = user->current_time + user->dt;        ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);        PetscReal fnnorm;        ierr = VecNorm(algebra->fn,NORM_INFINITY,&fnnorm);CHKERRQ(ierr);        if(0){          PetscViewer    viewer;          ierr = OutputVTK(user->dm, "function.vtk", &viewer);CHKERRQ(ierr);          ierr = VecView(algebra->fn, viewer);CHKERRQ(ierr);          ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);          ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with founction norm = %g /n",                                user->current_step, user->current_time, fnnorm);CHKERRQ(ierr);          //break;        }        if(user->Explicit_RK2){          ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);//U^n          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);//U^{(1)}          ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);//f(U^{(1)})          ierr = VecAXPY(algebra->solution, 1.0, algebra->oldsolution);CHKERRQ(ierr);//U^n + U^{(1)}          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);// + dt*f(U^{(1)})          ierr = VecScale(algebra->solution, 0.5);CHKERRQ(ierr);        }else{          ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);        }        {// Monitor the solution and function norms          PetscReal         norm;          PetscLogDouble    space =0;          PetscInt          size;          ierr = VecNorm(algebra->solution,NORM_INFINITY,&norm);CHKERRQ(ierr);          ierr = VecGetSize(algebra->solution, &size);CHKERRQ(ierr);          norm = norm/size;          if (norm>1.e5) {            SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_LIB,            "The norm of the solution is: %f (current time: %f). The explicit method is going to DIVERGE!!!", norm, user->current_time);          }          if (user->current_step%10==0) {            ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with solution norm = %g and founction norm = %g /n",//.........这里部分代码省略.........

static PetscErrorCode SampleOnGrid(MPI_Comm comm,Op op,const PetscInt M[3],const PetscInt smooth[2],PetscInt nrepeat,PetscLogDouble mintime,PetscLogDouble *memused,PetscLogDouble *memavail,PetscBool monitor) {  PetscErrorCode ierr;  PetscInt pgrid[3],cmax,fedegree,dof,addquadpts,nlevels,M_max,solve_type=0;  PetscMPIInt nranks;  Grid grid;  DM dm;  Vec U,V=NULL,F;  Mat A=NULL;  KSP ksp=NULL;  MG mg=NULL;  const char *solve_types[2] = {"fmg","ksp"};  PetscReal L[3];  PetscBool affine,ksp_only = PETSC_FALSE;#ifdef USE_HPM  char eventname[256];#endif  PetscFunctionBegin;  ierr = PetscOptionsBegin(comm,NULL,"KSP or FMG solver option",NULL);CHKERRQ(ierr);  ierr = PetscOptionsEList("-solve_type","Solve with KSP or FMG","",solve_types,2,solve_types[0],&solve_type,NULL);CHKERRQ(ierr);  if (solve_type) {ksp_only = PETSC_TRUE;}  ierr = PetscOptionsEnd();CHKERRQ(ierr);  ierr = OpGetFEDegree(op,&fedegree);CHKERRQ(ierr);  ierr = OpGetDof(op,&dof);CHKERRQ(ierr);  ierr = OpGetAddQuadPts(op,&addquadpts);CHKERRQ(ierr);  ierr = MPI_Comm_size(comm,&nranks);CHKERRQ(ierr);  ierr = ProcessGridFindSquarest(nranks,pgrid);CHKERRQ(ierr);  // It would make sense to either use a different coarsening criteria (perhaps even specified by the sampler).  On  // large numbers of processes, the coarse grids should be square enough that 192 is a good threshold size.  cmax = 192;  ierr = GridCreate(comm,M,pgrid,cmax,&grid);CHKERRQ(ierr);  ierr = GridGetNumLevels(grid,&nlevels);CHKERRQ(ierr);  ierr = DMCreateFE(grid,fedegree,dof,addquadpts,&dm);CHKERRQ(ierr);  M_max = PetscMax(M[0],PetscMax(M[1],M[2]));  L[0] = M[0]*1./M_max;  L[1] = M[1]*1./M_max;  L[2] = M[2]*1./M_max;  ierr = DMFESetUniformCoordinates(dm,L);CHKERRQ(ierr);  ierr = OpGetAffineOnly(op,&affine);CHKERRQ(ierr);  if (!affine) {ierr = DMCoordDistort(dm,L);CHKERRQ(ierr);}  ierr = DMCreateGlobalVector(dm,&U);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dm,&F);CHKERRQ(ierr);  ierr = OpForcing(op,dm,F);CHKERRQ(ierr);  if (!ksp_only) {    ierr = MGCreate(op,dm,nlevels,&mg);CHKERRQ(ierr);    ierr = MGMonitorSet(mg,monitor);CHKERRQ(ierr);    ierr = MGSetUpPC(mg);CHKERRQ(ierr);  }  else {    ierr = DMCreateGlobalVector(dm,&V);CHKERRQ(ierr);    ierr = OpGetMat(op,dm,&A);CHKERRQ(ierr);    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);  }#ifdef USE_HPM  ierr = PetscSNPrintf(eventname,sizeof eventname,"Solve G[%D %D %D]",M[0],M[1],M[2]);CHKERRQ(ierr);  HPM_Start(eventname);#endif  PetscInt i = 0;  PetscLogDouble sampletime = 0;  while ( (i<nrepeat) || (sampletime < mintime) ) {    PetscLogDouble t0,t1,elapsed,flops,eqs;    ierr = VecZeroEntries(U);CHKERRQ(ierr);    ierr = MPI_Barrier(comm);CHKERRQ(ierr);    ierr = PetscTime(&t0);CHKERRQ(ierr);    flops = petsc_TotalFlops;    if (!ksp_only) {      ierr = MGFCycle(op,mg,smooth[0],smooth[1],F,U);CHKERRQ(ierr);    }    else {      ierr = KSPSolve(ksp,F,V);CHKERRQ(ierr);      ierr = VecAXPY(V,-1.,U);CHKERRQ(ierr);    }    ierr = PetscTime(&t1);CHKERRQ(ierr);    flops = petsc_TotalFlops - flops;    elapsed = t1 - t0;    ierr = MPI_Allreduce(MPI_IN_PLACE,&elapsed,1,MPI_DOUBLE,MPI_MAX,comm);CHKERRQ(ierr);    ierr = MPI_Allreduce(MPI_IN_PLACE,&flops,1,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr);    eqs = (double)(M[0]*fedegree+1)*(M[1]*fedegree+1)*(M[2]*fedegree+1)*dof;    ierr = PetscPrintf(comm,"Q%D G[%5D%5D%5D] P[%3D%3D%3D] %10.3e s  %10f GF  %10f MEq/s/n",fedegree,M[0],M[1],M[2],pgrid[0],pgrid[1],pgrid[2],t1-t0,flops/elapsed*1e-9,eqs/elapsed*1e-6);CHKERRQ(ierr);    i++;    sampletime += elapsed;  }#ifdef USE_HPM  HPM_Stop(eventname);#endif  if (memused) {ierr = MemoryGetUsage(memused,memavail);CHKERRQ(ierr);  }  ierr = MGDestroy(&mg);CHKERRQ(ierr);  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);//.........这里部分代码省略.........

int main(int argc,char **argv){  PetscMPIInt      rank;  PetscInt         M=8,dof=1,stencil_width=1,i,start,end,P=5,N = 6,m=PETSC_DECIDE,n=PETSC_DECIDE,p=PETSC_DECIDE,pt = 0,st = 0;  PetscErrorCode   ierr;  PetscBool        flg = PETSC_FALSE,flg2,flg3;  DMDABoundaryType periodic;  DMDAStencilType  stencil_type;  DM               da;  Vec              local,global,local_copy;  PetscScalar      value;  PetscReal        norm,work;  PetscViewer      viewer;  char             filename[64];  FILE             *file;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-M",&M,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-N",&N,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-dof",&dof,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-stencil_width",&stencil_width,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(NULL,"-periodic",&pt,NULL);CHKERRQ(ierr);  periodic = (DMDABoundaryType) pt;  ierr = PetscOptionsGetInt(NULL,"-stencil_type",&st,NULL);CHKERRQ(ierr);  stencil_type = (DMDAStencilType) st;  ierr = PetscOptionsHasName(NULL,"-2d",&flg2);CHKERRQ(ierr);  ierr = PetscOptionsHasName(NULL,"-3d",&flg3);CHKERRQ(ierr);  if (flg2) {    ierr = DMDACreate2d(PETSC_COMM_WORLD,periodic,periodic,stencil_type,M,N,m,n,dof,stencil_width,                        NULL,NULL,&da);CHKERRQ(ierr);  } else if (flg3) {    ierr = DMDACreate3d(PETSC_COMM_WORLD,periodic,periodic,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,                        NULL,NULL,NULL,&da);CHKERRQ(ierr);  } else {    ierr = DMDACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,NULL,&da);CHKERRQ(ierr);  }  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  ierr = VecDuplicate(local,&local_copy);CHKERRQ(ierr);  /* zero out vectors so that ghostpoints are zero */  value = 0;  ierr  = VecSet(local,value);CHKERRQ(ierr);  ierr  = VecSet(local_copy,value);CHKERRQ(ierr);  ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr);  for (i=start; i<end; i++) {    value = i + 1;    ierr  = VecSetValues(global,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);  }  ierr = VecAssemblyBegin(global);CHKERRQ(ierr);  ierr = VecAssemblyEnd(global);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMDALocalToLocalBegin(da,local,INSERT_VALUES,local_copy);CHKERRQ(ierr);  ierr = DMDALocalToLocalEnd(da,local,INSERT_VALUES,local_copy);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-save",&flg,NULL);CHKERRQ(ierr);  if (flg) {    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);    sprintf(filename,"local.%d",rank);    ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&viewer);CHKERRQ(ierr);    ierr = PetscViewerASCIIGetPointer(viewer,&file);CHKERRQ(ierr);    ierr = VecView(local,viewer);CHKERRQ(ierr);    fprintf(file,"Vector with correct ghost points/n");    ierr = VecView(local_copy,viewer);CHKERRQ(ierr);    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  }  ierr = VecAXPY(local_copy,-1.0,local);CHKERRQ(ierr);  ierr = VecNorm(local_copy,NORM_MAX,&work);CHKERRQ(ierr);  ierr = MPI_Allreduce(&work,&norm,1,MPIU_REAL,MPIU_MAX,PETSC_COMM_WORLD);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of difference %G should be zero/n",norm);CHKERRQ(ierr);  ierr = VecDestroy(&local_copy);CHKERRQ(ierr);  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){  AppCtx         user;                /* user-defined work context */  PetscInt       mx,my;  PetscErrorCode ierr;  MPI_Comm       comm;  DM             da;  Vec            x;  Mat            J = NULL,Jmf = NULL;  MatShellCtx    matshellctx;  PetscInt       mlocal,nlocal;  PC             pc;  KSP            ksp;  PetscBool      errorinmatmult = PETSC_FALSE,errorinpcapply = PETSC_FALSE,errorinpcsetup = PETSC_FALSE;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return(1);  PetscFunctionBeginUser;  ierr = PetscOptionsGetBool(NULL,"-error_in_matmult",&errorinmatmult,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-error_in_pcapply",&errorinpcapply,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-error_in_pcsetup",&errorinpcsetup,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-error_in_domain",&user.errorindomain,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-error_in_domainmf",&user.errorindomainmf,NULL);CHKERRQ(ierr);    comm = PETSC_COMM_WORLD;  ierr = SNESCreate(comm,&user.snes);CHKERRQ(ierr);  /*      Create distributed array object to manage parallel grid and vectors      for principal unknowns (x) and governing residuals (f)  */  ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);CHKERRQ(ierr);  ierr = SNESSetDM(user.snes,da);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,                     PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  /*     Problem parameters (velocity of lid, prandtl, and grashof numbers)  */  user.lidvelocity = 1.0/(mx*my);  user.prandtl     = 1.0;  user.grashof     = 1.0;  ierr = PetscOptionsGetReal(NULL,"-lidvelocity",&user.lidvelocity,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetReal(NULL,"-prandtl",&user.prandtl,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetReal(NULL,"-grashof",&user.grashof,NULL);CHKERRQ(ierr);  ierr = PetscOptionsHasName(NULL,"-contours",&user.draw_contours);CHKERRQ(ierr);  ierr = DMDASetFieldName(da,0,"x_velocity");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,1,"y_velocity");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,2,"Omega");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,3,"temperature");CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create user context, set problem data, create vector data structures.     Also, compute the initial guess.     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create nonlinear solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);CHKERRQ(ierr);  if (errorinmatmult) {    ierr = MatCreateSNESMF(user.snes,&Jmf);CHKERRQ(ierr);    ierr = MatSetFromOptions(Jmf);CHKERRQ(ierr);    ierr = MatGetLocalSize(Jmf,&mlocal,&nlocal);CHKERRQ(ierr);    matshellctx.Jmf = Jmf;    ierr = MatCreateShell(PetscObjectComm((PetscObject)Jmf),mlocal,nlocal,PETSC_DECIDE,PETSC_DECIDE,&matshellctx,&J);CHKERRQ(ierr);    ierr = MatShellSetOperation(J,MATOP_MULT,(void (*)(void))MatMult_MyShell);CHKERRQ(ierr);    ierr = MatShellSetOperation(J,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MyShell);CHKERRQ(ierr);    ierr = SNESSetJacobian(user.snes,J,J,MatMFFDComputeJacobian,NULL);CHKERRQ(ierr);  }  ierr = SNESSetFromOptions(user.snes);CHKERRQ(ierr);  ierr = PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g/n",(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);CHKERRQ(ierr);  if (errorinpcapply) {    ierr = SNESGetKSP(user.snes,&ksp);CHKERRQ(ierr);    ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);    ierr = PCShellSetApply(pc,PCApply_MyShell);CHKERRQ(ierr);  }  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve the nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);  ierr = FormInitialGuess(&user,da,x);CHKERRQ(ierr);  if (errorinpcsetup) {    ierr = SNESSetUp(user.snes);CHKERRQ(ierr);    ierr = SNESSetJacobian(user.snes,NULL,NULL,SNESComputeJacobian_MyShell,NULL);CHKERRQ(ierr);  }  ierr = SNESSolve(user.snes,NULL,x);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -//.........这里部分代码省略.........

int main(int argc,char **argv){  TS                ts;         /* time integrator */  SNES              snes;       /* nonlinear solver */  SNESLineSearch    linesearch; /* line search */  Vec               X;          /* solution, residual vectors */  Mat               J;          /* Jacobian matrix */  PetscInt          steps,maxsteps,mx;  PetscErrorCode    ierr;  DM                da;  PetscReal         ftime,dt;  struct _User      user;       /* user-defined work context */  TSConvergedReason reason;  PetscInitialize(&argc,&argv,(char*)0,help);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-11,2,2,NULL,&da);CHKERRQ(ierr);  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Extract global vectors from DMDA;   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);  /* Initialize user application context */  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","");  {    user.a[0] = 1;           ierr = PetscOptionsReal("-a0","Advection rate 0","",user.a[0],&user.a[0],NULL);CHKERRQ(ierr);    user.a[1] = 0;           ierr = PetscOptionsReal("-a1","Advection rate 1","",user.a[1],&user.a[1],NULL);CHKERRQ(ierr);    user.k[0] = 1e6;         ierr = PetscOptionsReal("-k0","Reaction rate 0","",user.k[0],&user.k[0],NULL);CHKERRQ(ierr);    user.k[1] = 2*user.k[0]; ierr = PetscOptionsReal("-k1","Reaction rate 1","",user.k[1],&user.k[1],NULL);CHKERRQ(ierr);    user.s[0] = 0;           ierr = PetscOptionsReal("-s0","Source 0","",user.s[0],&user.s[0],NULL);CHKERRQ(ierr);    user.s[1] = 1;           ierr = PetscOptionsReal("-s1","Source 1","",user.s[1],&user.s[1],NULL);CHKERRQ(ierr);  }  ierr = PetscOptionsEnd();CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetDM(ts,da);CHKERRQ(ierr);  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr);  ierr = DMCreateMatrix(da,MATAIJ,&J);CHKERRQ(ierr);  ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr);  /* A line search in the nonlinear solve can fail due to ill-conditioning unless an absolute tolerance is set. Since   * this problem is linear, we deactivate the line search. For a linear problem, it is usually recommended to also use   * SNESSetType(snes,SNESKSPONLY). */  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);  ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);  ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);CHKERRQ(ierr);  ftime    = 1.0;  maxsteps = 10000;  ierr     = TSSetDuration(ts,maxsteps,ftime);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr);  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);  ierr = VecGetSize(X,&mx);CHKERRQ(ierr);  dt   = .1 * PetscMax(user.a[0],user.a[1]) / mx; /* Advective CFL, I don't know why it needs so much safety factor. */  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set runtime options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSolve(ts,X);CHKERRQ(ierr);  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps/n",TSConvergedReasons[reason],ftime,steps);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = VecDestroy(&X);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){  PetscErrorCode ierr;  DM             da;                   /* structured grid topology object */  TS             ts;                   /* time-stepping object (contains snes) */  SNES           snes;                 /* Newton solver object */  Vec            X,residual;           /* solution, residual */  Mat            J;                    /* Jacobian matrix */  PetscInt       Mx,My,fsteps,steps;  ISColoring     iscoloring;  PetscReal      tstart,tend,ftime,secperday=3600.0*24.0,Y0;  PetscBool      fdflg = PETSC_FALSE, mfileflg = PETSC_FALSE, optflg = PETSC_FALSE;  char           mfile[PETSC_MAX_PATH_LEN] = "out.m";  MatFDColoring  matfdcoloring;  PorousCtx      user;                 /* user-defined work context */  PetscInitialize(&argc,&argv,(char *)0,help);  ierr = DMDACreate2d(PETSC_COMM_WORLD,             DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, // correct for zero Dirichlet             DMDA_STENCIL_STAR, // nonlinear diffusion but diffusivity                                //   depends on soln W not grad W             -21,-21,           // default to 20x20 grid but override with                                //   -da_grid_x, -da_grid_y (or -da_refine)             PETSC_DECIDE,PETSC_DECIDE, // num of procs in each dim             2,                 // dof = 2:  node = (W,Y)                                //        or node = (P,dPsqr)                                //        or node = (ddxE,ddyN)             1,                 // s = 1 (stencil extends out one cell)             PETSC_NULL,PETSC_NULL, // no specify proc decomposition             &da);CHKERRQ(ierr);  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);  /* get Vecs and Mats for this grid */  ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);  ierr = VecDuplicate(X,&residual);CHKERRQ(ierr);  ierr = VecDuplicate(X,&user.geom);CHKERRQ(ierr);  ierr = DMGetMatrix(da,MATAIJ,&J);CHKERRQ(ierr);  /* set up contexts */  tstart   = 10.0 * secperday; /* 10 days in seconds */  tend     = 30.0 * secperday;  steps    = 20;  Y0       = 1.0;              /* initial value of Y, for computing initial                                  value of P; note Ymin = 0.1 is different */  user.da = da;  ierr = DefaultContext(&user);CHKERRQ(ierr);  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,           "","options to (W,P)-space better hydrology model alt","");CHKERRQ(ierr);  {    ierr = PetscOptionsReal("-alt_sigma","nonlinear power","",                            user.sigma,&user.sigma,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_Ymin",                            "min capacity thickness (esp. in pressure computation)","",                            user.Ymin,&user.Ymin,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_Wmin",                            "min water amount (esp. in pressure computation)","",                            user.Wmin,&user.Wmin,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_Y0",                            "constant initial capacity thickness","",                            Y0,&Y0,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_Cmelt",                            "additional coefficient for amount of melt","",                            user.Cmelt,&user.Cmelt,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_Creep",                            "creep closure coefficient","",                            user.Creep,&user.Creep,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_L","half-width of square region in meters","",                            user.L,&user.L,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsReal("-alt_tstart_days","start time in days","",                            tstart/secperday,&tstart,&optflg);CHKERRQ(ierr);    if (optflg) { tstart *= secperday; }    ierr = PetscOptionsReal("-alt_tend_days","end time in days","",                            tend/secperday,&tend,&optflg);CHKERRQ(ierr);    if (optflg) { tend *= secperday; }    ierr = PetscOptionsInt("-alt_steps","number of timesteps to take","",                           steps,&steps,PETSC_NULL);CHKERRQ(ierr);    ierr = PetscOptionsBool("-alt_converge_check",                            "run silent and check for convergence",                            "",user.run_silent,&user.run_silent,PETSC_NULL);                            CHKERRQ(ierr);    ierr = PetscOptionsString("-mfile",                            "name of Matlab file to write results","",                            mfile,mfile,PETSC_MAX_PATH_LEN,&mfileflg);                            CHKERRQ(ierr);  }  ierr = PetscOptionsEnd();CHKERRQ(ierr);  /* fix remaining parameters */  ierr = DerivedConstants(&user);CHKERRQ(ierr);  ierr = VecStrideSet(user.geom,0,user.H0);CHKERRQ(ierr);  /* H(x,y) = H0 */  ierr = VecStrideSet(user.geom,1,0.0);CHKERRQ(ierr);      /* b(x,y) = 0  */  ierr = DMDASetUniformCoordinates(da,  // square domain              -user.L, user.L, -user.L, user.L, 0.0, 1.0);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);  user.dx = 2.0 * user.L / (Mx-1);//.........这里部分代码省略.........

//.........这里部分代码省略.........        }        if (pfdata->bus[i-vStart].nload) {          for (j=0; j < pfdata->bus[i-vStart].nload; j++) {            ierr = DMNetworkAddComponent(networkdm,i,componentkey[3],&pfdata->load[loadj++]);CHKERRQ(ierr);          }        }        /* Add number of variables */        ierr = DMNetworkAddNumVariables(networkdm,i,2);CHKERRQ(ierr);      }    }    /* Set up DM for use */    ierr = DMSetUp(networkdm);CHKERRQ(ierr);    if (!crank) {      ierr = PetscFree(pfdata->bus);CHKERRQ(ierr);      ierr = PetscFree(pfdata->gen);CHKERRQ(ierr);      ierr = PetscFree(pfdata->branch);CHKERRQ(ierr);      ierr = PetscFree(pfdata->load);CHKERRQ(ierr);      ierr = PetscFree(pfdata);CHKERRQ(ierr);    }        ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);    if (size > 1) {      DM distnetworkdm;      /* Network partitioning and distribution of data */      ierr = DMNetworkDistribute(networkdm,0,&distnetworkdm);CHKERRQ(ierr);      ierr = DMDestroy(&networkdm);CHKERRQ(ierr);      networkdm = distnetworkdm;    }        PetscLogStagePop();    ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr);    ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr);    #if 0    PetscInt numComponents;    EDGEDATA edge;    PetscInt offset,key,kk;    DMNetworkComponentGenericDataType *arr;    VERTEXDATA     bus;    GEN            gen;    LOAD           load;        for (i = eStart; i < eEnd; i++) {      ierr = DMNetworkGetComponentDataArray(networkdm,&arr);CHKERRQ(ierr);      ierr = DMNetworkGetComponentTypeOffset(networkdm,i,0,&key,&offset);CHKERRQ(ierr);      edge = (EDGEDATA)(arr+offset);      ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr);      ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d/n",crank,numComponents,edge->internal_i,edge->internal_j);CHKERRQ(ierr);    }            for (i = vStart; i < vEnd; i++) {      ierr = DMNetworkGetComponentDataArray(networkdm,&arr);CHKERRQ(ierr);      ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr);      for (kk=0; kk < numComponents; kk++) {        ierr = DMNetworkGetComponentTypeOffset(networkdm,i,kk,&key,&offset);CHKERRQ(ierr);        if (key == 1) {          bus = (VERTEXDATA)(arr+offset);          ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d/n",crank,numComponents,bus->internal_i);CHKERRQ(ierr);        } else if (key == 2) {          gen = (GEN)(arr+offset);          ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f/n",crank,gen->pg,gen->qg);CHKERRQ(ierr);        } else if (key == 3) {          load = (LOAD)(arr+offset);          ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f/n",crank,load->pl,load->ql);CHKERRQ(ierr);        }      }    }  #endif      /* Broadcast Sbase to all processors */    ierr = MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);        ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr);    ierr = VecDuplicate(X,&F);CHKERRQ(ierr);        ierr = DMCreateMatrix(networkdm,&J);CHKERRQ(ierr);    ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);        ierr = SetInitialValues(networkdm,X,&User);CHKERRQ(ierr);        /* HOOK UP SOLVER */    ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);    ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr);    ierr = SNESSetFunction(snes,F,FormFunction,&User);CHKERRQ(ierr);    ierr = SNESSetJacobian(snes,J,J,FormJacobian,&User);CHKERRQ(ierr);    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);        ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr);        ierr = VecDestroy(&X);CHKERRQ(ierr);    ierr = VecDestroy(&F);CHKERRQ(ierr);    ierr = MatDestroy(&J);CHKERRQ(ierr);        ierr = SNESDestroy(&snes);CHKERRQ(ierr);    ierr = DMDestroy(&networkdm);CHKERRQ(ierr);  }  ierr = PetscFinalize();  return ierr;}