当前位置: 首页>>代码示例>>C++>>正文

C++ DMCreateGlobalVector函数代码示例

本文整理汇总了C++中DMCreateGlobalVector函数的典型用法代码示例。如果您正苦于以下问题：C++ DMCreateGlobalVector函数的具体用法？C++ DMCreateGlobalVector怎么用？C++ DMCreateGlobalVector使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。


在下文中一共展示了DMCreateGlobalVector函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: main

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);
//.........这里部分代码省略.........
开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:ex65dm.c

示例2: efs_setup

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;
}
开发者ID:tuxfan，项目名称:ska，代码行数:82，代码来源:ef_solver.c

示例3: main

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
//.........这里部分代码省略.........
开发者ID:firedrakeproject，项目名称:petsc，代码行数:101，代码来源:ex50.c

示例4: direction

/*@
    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);
//.........这里部分代码省略.........
开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:gr1.c

示例5: main

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;
}
开发者ID:petsc，项目名称:petsc，代码行数:78，代码来源:ex45.c

示例6: main

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;
}
开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:74，代码来源:ex35.c

示例7: main

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;
}
开发者ID:masa-ito，项目名称:PETScToPoisson，代码行数:78，代码来源:ex12.c

示例8: main

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);
//.........这里部分代码省略.........
开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:101，代码来源:ex15.c

示例9: main

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);
}
开发者ID:fengyuqi，项目名称:petsc，代码行数:79，代码来源:ex4.c

示例10: DMCoarsen_SNESVI

/*
     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);
}
开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:76，代码来源:virs.c

示例11: main

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);
}
开发者ID:feelpp，项目名称:debian-petsc，代码行数:77，代码来源:ex7.c

示例12: main

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);
//.........这里部分代码省略.........
开发者ID:tom-klotz，项目名称:petsc，代码行数:101，代码来源:ex4.c

示例13: main

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);
    }
//.........这里部分代码省略.........
开发者ID:pombredanne，项目名称:petsc，代码行数:101，代码来源:ex15.c

示例14: main

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);
//.........这里部分代码省略.........
开发者ID:hsahasra，项目名称:petsc-magma-dense-mat，代码行数:101，代码来源:ex62.c

示例15: main

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();
}
开发者ID:bueler，项目名称:p4pdes，代码行数:79，代码来源:elasto.c


注：本文中的DMCreateGlobalVector函数示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。