C++ VecGetArray函数代码示例

PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
{
  PetscErrorCode    ierr;
  PetscInt          nele,nen,n,i;
  const PetscInt    *ele;
  Vec               coords, rand1, rand2;
  const PetscScalar *_coords;
  PetscScalar       x[3],y[3];
  PetscInt          idx[3];
  PetscScalar       *xx,*w1,*w2,*u1,*u2,*u3;
  PetscViewer       view_out;

  PetscFunctionBeginUser;
  /* Get ghosted coordinates */
  ierr = DMGetCoordinatesLocal(user->da,&coords);CHKERRQ(ierr);
  ierr = VecDuplicate(user->u1,&rand1);
  ierr = VecDuplicate(user->u1,&rand2);
  ierr = VecSetRandom(rand1,NULL);
  ierr = VecSetRandom(rand2,NULL);

  ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
  ierr = VecGetArrayRead(coords,&_coords);CHKERRQ(ierr);
  ierr = VecGetArray(X,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(user->work1,&w1);
  ierr = VecGetArray(user->work2,&w2);
  ierr = VecGetArray(user->u1,&u1);
  ierr = VecGetArray(user->u2,&u2);
  ierr = VecGetArray(user->u3,&u3);

  /* Get local element info */
  ierr = DMDAGetElements(user->da,&nele,&nen,&ele);CHKERRQ(ierr);
  for (i=0; i < nele; i++) {
    idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
    x[0]   = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
    x[1]   = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
    x[2]   = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

    PetscScalar vals1[3],vals2[3],valsrand[3];
    PetscInt    r;
    for (r=0; r<3; r++) {
      valsrand[r]=5*x[r]*(1-x[r])*y[r]*(1-y[r]);
      if (x[r]>=0.5 && y[r]>=0.5) {
        vals1[r]=0.75;
        vals2[r]=0.0;
      }
      if (x[r]>=0.5 && y[r]<0.5) {
        vals1[r]=0.0;
        vals2[r]=0.0;
      }
      if (x[r]<0.5 && y[r]>=0.5) {
        vals1[r]=0.0;
        vals2[r]=0.75;
      }
      if (x[r]<0.5 && y[r]<0.5) {
        vals1[r]=0.75;
        vals2[r]=0.0;
      }
    }

    ierr = VecSetValues(user->work1,3,idx,vals1,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(user->work2,3,idx,vals2,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(user->work3,3,idx,valsrand,INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = VecAssemblyBegin(user->work1);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work1);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(user->work2);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work2);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(user->work3);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(user->work3);CHKERRQ(ierr);

  ierr = VecAXPY(user->work1,1.0,user->work3);CHKERRQ(ierr);
  ierr = VecAXPY(user->work2,1.0,user->work3);CHKERRQ(ierr);

  for (i=0; i<n/4; i++) {
    xx[4*i] = w1[i];
    if (xx[4*i]>1) xx[4*i]=1;

    xx[4*i+1] = w2[i];
    if (xx[4*i+1]>1) xx[4*i+1]=1;

    if (xx[4*i]+xx[4*i+1]>1) xx[4*i+1] = 1.0 - xx[4*i];

    xx[4*i+2] = 1.0 - xx[4*i] - xx[4*i+1];
    xx[4*i+3] = 0.0;

    u1[i] = xx[4*i];
    u2[i] = xx[4*i+1];
    u3[i] = xx[4*i+2];
  }

  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);CHKERRQ(ierr);
  ierr = VecView(user->u1,view_out);CHKERRQ(ierr);
  ierr = VecView(user->u2,view_out);CHKERRQ(ierr);
  ierr = VecView(user->u3,view_out);CHKERRQ(ierr);
  PetscViewerDestroy(&view_out);

  ierr = DMDARestoreElements(user->da,&nele,&nen,&ele);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(coords,&_coords);CHKERRQ(ierr);
  ierr = VecRestoreArray(X,&xx);CHKERRQ(ierr);
//.........这里部分代码省略.........

PetscInt main(PetscInt argc,char **args)
{
  ptrdiff_t      N0=256,N1=256,N2=256,N3=2,dim[4];
  fftw_plan      bplan,fplan;
  fftw_complex   *out;
  double         *in1,*in2;
  ptrdiff_t      alloc_local,local_n0,local_0_start;
  ptrdiff_t      local_n1,local_1_start;
  PetscInt       i,j,indx,n1;
  PetscInt       size,rank,n,N,*in,N_factor,NM;
  PetscScalar    *data_fin,value1,one=1.57,zero=0.0;
  PetscScalar    a,*x_arr,*y_arr,*z_arr,enorm;
  Vec            fin,fout,fout1,ini,final;
  PetscRandom    rnd;
  PetscErrorCode ierr;
  VecScatter     vecscat,vecscat1;
  IS             indx1,indx2;
  PetscInt       *indx3,k,l,*indx4;
  PetscInt       low,tempindx,tempindx1;


  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
#if defined(PETSC_USE_COMPLEX)
  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
#endif
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);

  PetscRandomCreate(PETSC_COMM_WORLD,&rnd);


  alloc_local = fftw_mpi_local_size_3d_transposed(N0,N1,N2/2+1,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start);

/*    printf("The value alloc_local is %ld from process %d\n",alloc_local,rank);     */
  printf("The value local_n0 is %ld from process %d\n",local_n0,rank);
/*    printf("The value local_0_start is  %ld from process %d\n",local_0_start,rank);*/
/*    printf("The value local_n1 is  %ld from process %d\n",local_n1,rank);          */
/*    printf("The value local_1_start is  %ld from process %d\n",local_1_start,rank);*/

  /* Allocate space for input and output arrays  */

  in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
  in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
  out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);


  N=2*N0*N1*(N2/2+1);N_factor=N0*N1*N2;
  n=2*local_n0*N1*(N2/2+1);n1=local_n1*N0*2*N1;

/*    printf("The value N is  %d from process %d\n",N,rank);   */
/*    printf("The value n is  %d from process %d\n",n,rank);   */
/*    printf("The value n1 is  %d from process %d\n",n1,rank); */
  /* Creating data vector and accompanying array with VeccreateMPIWithArray */
  ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);CHKERRQ(ierr);

/*    VecGetSize(fin,&size); */
/*    printf("The size is %d\n",size); */

  VecSet(fin,one);
  VecSet(fout,zero);
  VecSet(fout1,zero);

  VecAssemblyBegin(fin);
  VecAssemblyEnd(fin);
/*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */


  VecGetArray(fin,&x_arr);
  VecGetArray(fout1,&z_arr);
  VecGetArray(fout,&y_arr);

  fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
  bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);

  fftw_execute(fplan);
  fftw_execute(bplan);

  VecRestoreArray(fin,&x_arr);
  VecRestoreArray(fout1,&z_arr);
  VecRestoreArray(fout,&y_arr);


/*    a = 1.0/(PetscReal)N_factor; */
/*    ierr = VecScale(fout1,a);CHKERRQ(ierr); */
  VecCreate(PETSC_COMM_WORLD,&ini);
  VecCreate(PETSC_COMM_WORLD,&final);
  VecSetSizes(ini,local_n0*N1*N2,N_factor);
  VecSetSizes(final,local_n0*N1*N2,N_factor);
/*    VecSetSizes(ini,PETSC_DECIDE,N_factor); */
/*    VecSetSizes(final,PETSC_DECIDE,N_factor); */
  VecSetFromOptions(ini);
  VecSetFromOptions(final);

  if (N2%2==0) NM=N2+2;
  else NM=N2+1;

  ierr = VecGetOwnershipRange(fin,&low,NULL);
  printf("The local index is %d from %d\n",low,rank);
//.........这里部分代码省略.........

/*@
    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;
  DMDABoundaryType 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",xmin,xmax);
  if ((ymax <= ymin) && (dim > 1)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"ymax must be larger than ymin %G %G",ymin,ymax);
  if ((zmax <= zmin) && (dim > 2)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"zmax must be larger than zmin %G %G",zmin,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 == DMDA_BOUNDARY_PERIODIC) hx  = (xmax-xmin)/(M+1);
    else                              hx  = (xmax-xmin)/(M ? M : 1);
    if (by == DMDA_BOUNDARY_PERIODIC) hy  = (ymax-ymin)/(N+1);
    else                              hy  = (ymax-ymin)/(N ? N : 1);
    if (bz == DMDA_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);
//.........这里部分代码省略.........

/*
   FormFunction - Evaluates the function and corresponding gradient.

   Input Parameters:
   tao - the Tao context
   X   - the input vector
   ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

   Output Parameters:
   f   - the newly evaluated function
*/
PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
{
  AppCtx         *ctx = (AppCtx*)ctx0;
  TS             ts;

  Vec            U;             /* solution will be stored here */
  Mat            A;             /* Jacobian matrix */
  Mat            Jacp;          /* Jacobian matrix */
  PetscErrorCode ierr;
  PetscInt       n = 2;
  PetscReal      ftime;
  PetscInt       steps;
  PetscScalar    *u;
  PetscScalar    *x_ptr,*y_ptr;
  Vec            lambda[1],q,mu[1];

  ierr = VecGetArray(P,&x_ptr);CHKERRQ(ierr);
  ctx->Pm = x_ptr[0];
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Create necessary matrix and vectors
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);

  ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&Jacp);CHKERRQ(ierr);
  ierr = MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr);
  ierr = MatSetFromOptions(Jacp);CHKERRQ(ierr);
  ierr = MatSetUp(Jacp);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 = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);CHKERRQ(ierr);
  ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,ctx);CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set initial conditions
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = VecGetArray(U,&u);CHKERRQ(ierr);
  u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
  u[1] = 1.0;
  ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
  ierr = TSSetSolution(ts,U);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Save trajectory of solution so that TSAdjointSolve() may be used
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set solver options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSetDuration(ts,PETSC_DEFAULT,10.0);CHKERRQ(ierr);
  ierr = TSSetInitialTimeStep(ts,0.0,.01);CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Solve nonlinear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSolve(ts,U);CHKERRQ(ierr);

  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Adjoint model starts here
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = MatCreateVecs(A,&lambda[0],NULL);CHKERRQ(ierr);
  /*   Set initial conditions for the adjoint integration */
  ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr);
  y_ptr[0] = 0.0; y_ptr[1] = 0.0;
  ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr);

  ierr = MatCreateVecs(Jacp,&mu[0],NULL);CHKERRQ(ierr);
  ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr);
  x_ptr[0] = -1.0;
  ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr);
  ierr = TSAdjointSetCostGradients(ts,1,lambda,mu);CHKERRQ(ierr);

  ierr = TSAdjointSetRHSJacobian(ts,Jacp,RHSJacobianP,ctx);CHKERRQ(ierr);
//.........这里部分代码省略.........

void getForces(Vec params, std::vector<Vec> &forces, DA da, timeInfo ti, int numParams) {
#ifdef __DEBUG__
  std::cout << RED"Entering "NRM << __func__ << std::endl;
#endif
  // Generate the force vector based on the current parameters ...
  // F = Sum { B_i a_i}
  
  PetscScalar * pVec;
  VecGetArray(params, &pVec);
  // Clear the Forces
  for (unsigned int i=0; i<forces.size(); i++) {
    if (forces[i] != NULL) {
      VecDestroy(forces[i]);
    }
  }
  forces.clear();

  unsigned int numSteps = (unsigned int)(ceil(( ti.stop - ti.start)/ti.step));
  // create and initialize to 0
  for (unsigned int i=0; i<numSteps+1; i++) {
    Vec tmp;
    DACreateGlobalVector(da, &tmp);
    VecZeroEntries(tmp);
    forces.push_back(tmp);
  }

  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  PetscScalar ***tauArray;

  unsigned int paramTimeSteps = (unsigned int)(ceil(( (double)(numSteps))/ ((double)(2*numParams)) ));
  double acx,acy,acz;

  int x, y, z, m, n, p;
  int mx,my,mz;

  DAGetCorners(da, &x, &y, &z, &m, &n, &p);
  DAGetInfo(da,0, &mx, &my, &mz, 0,0,0,0,0,0,0);

  double hx = 1.0/(mx-1.0);

  for (int b=0; b<numParams; b++) {
    std::vector<Vec> tau;
    unsigned int tBegin = paramTimeSteps*b;
    unsigned int tEnd   = tBegin + numSteps/2; // paramTimeSteps*(b+2);

    // std::cout << "For param " << b << ": Time step range is " << tBegin << " -> " << tEnd << std::endl; 
    for (unsigned int t=0; t<numSteps+1; t++) {
      double newTime = (ti.step*(t-tBegin)*numSteps)/((double)(paramTimeSteps));
      
      if ( (t>=tBegin) && (t<=tEnd)) {
        DAVecGetArray(da, forces[t], &tauArray);
        for (int k = z; k < z + p ; k++) {
          for (int j = y; j < y + n; j++) {
            for (int i = x; i < x + m; i++) {
              acx = (i)*hx; acy = (j)*hx; acz = (k)*hx;
              tauArray[k][j][i] += pVec[b]*sin(M_PI*newTime)*cos(2*M_PI*acx)*cos(2*M_PI*acy)*cos(2*M_PI*acz);
            }
          }
        }
        DAVecRestoreArray ( da, forces[t], &tauArray ) ;
      }
    }
  }
  VecRestoreArray(params, &pVec);

#ifdef __DEBUG__
  // Get the norms of the forces ... just to be safe ..
  double fNorm1, fNorm2;

  for (unsigned int i=0; i<forces.size(); i++) {
    VecNorm(forces[i], NORM_INFINITY, &fNorm1);
    VecNorm(forces[i], NORM_2, &fNorm2);
    PetscPrintf(0, "Force Norms at timestep %d are %g and %g\n", i, fNorm1, fNorm2);
  }
#endif
  
#ifdef __DEBUG__
  std::cout << GRN"Leaving "NRM << __func__ << std::endl;
#endif
}

/*
   FormFunction - Evaluates nonlinear function, F(x).

   Input Parameters:
.  snes - the SNES context
.  X - input vector
.  ptr - optional user-defined context, as set by SNESSetFunction()

   Output Parameter:
.  F - function vector
 */
int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
{
  AppCtx      *user = (AppCtx*)ptr;
  int         ierr,i,j,row,mx,my;
  PetscReal   two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
  PetscScalar u,uxx,uyy,*x,*f;

  /*
      Process 0 has to wait for all other processes to get here
   before proceeding to write in the shared vector
  */
  ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);

  if (user->rank) {
    /*
       All the non-busy processors have to wait here for process 0 to finish
       evaluating the function; otherwise they will start using the vector values
       before they have been computed
    */
    ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);
    return 0;
  }

  mx = user->mx;            my = user->my;            lambda = user->param;
  hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
  sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

  /*
     Get pointers to vector data
  */
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);

  /*
      The next line tells the SGI compiler that x and f contain no overlapping
    regions and thus it can use addition optimizations.
  */
#pragma arl(4)
#pragma distinct (*x,*f)
#pragma no side effects (exp)

  /*
     Compute function over the entire  grid
  */
  for (j=0; j<my; j++) {
    for (i=0; i<mx; i++) {
      row = i + j*mx;
      if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
        f[row] = x[row];
        continue;
      }
      u      = x[row];
      uxx    = (two*u - x[row-1] - x[row+1])*hydhx;
      uyy    = (two*u - x[row-mx] - x[row+mx])*hxdhy;
      f[row] = uxx + uyy - sc*exp(u);
    }
  }

  /*
     Restore vectors
  */
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);

  ierr = PetscLogFlops(11.0*(mx-2)*(my-2))CHKERRQ(ierr);
  ierr = PetscBarrier((PetscObject)X);CHKERRQ(ierr);
  return 0;
}

/**
Compute the gadient of the cell center gradient obtained by the least-square method
*/
PetscErrorCode GradientGradientJacobian(DM dm, Vec locX, PetscScalar elemMat[], void *ctx)
{
  User              user = (User) ctx;
  Physics           phys = user->model->physics;
  PetscInt          dof = phys->dof;
  const PetscScalar *facegeom, *cellgeom,*x;
  PetscErrorCode    ierr;
  DM                dmFace, dmCell;

  DM                dmGrad = user->dmGrad;
  PetscInt          fStart, fEnd, face, cStart;
  Vec               Grad;
  /*here the localGradLimiter refers to the gradient that has been multiplied by the limiter function.
   The locGradLimiter is used to construct the uL and uR, and the locGrad is used to caculate the diffusion term*/
  Vec               TempVec; /*a temperal vec for the vector restore*/

  PetscFunctionBeginUser;

  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);

  ierr = DMGetGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);
  ierr = VecZeroEntries(Grad);CHKERRQ(ierr);
  ierr = VecDuplicate(Grad, &TempVec);CHKERRQ(ierr);
  ierr = VecCopy(Grad, TempVec);CHKERRQ(ierr);

  ierr = VecGetArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
  {
    PetscScalar *grad;
    ierr = VecGetArray(TempVec,&grad);CHKERRQ(ierr);
    /* Reconstruct gradients */
    for (face=fStart; face<fEnd; face++) {
      const PetscInt    *cells;
      const PetscScalar *cx[2];
      const FaceGeom    *fg;
      PetscScalar       *cgrad[2];
      PetscInt          i,j;
      PetscBool         ghost;

      ierr = IsExteriorGhostFace(dm,face,&ghost);CHKERRQ(ierr);
      if (ghost) continue;
      ierr = DMPlexGetSupport(dm,face,&cells);CHKERRQ(ierr);
      ierr = DMPlexPointLocalRead(dmFace,face,facegeom,&fg);CHKERRQ(ierr);
      for (i=0; i<2; i++) {
        ierr = DMPlexPointLocalRead(dm,cells[i],x,&cx[i]);CHKERRQ(ierr);
        ierr = DMPlexPointGlobalRef(dmGrad,cells[i],grad,&cgrad[i]);CHKERRQ(ierr);
      }
      for (i=0; i<dof; i++) {
        PetscScalar delta = cx[1][i] - cx[0][i];
        for (j=0; j<DIM; j++) {
          if (cgrad[0]) cgrad[0][i*DIM+j] += fg->grad[0][j] * delta;
          if (cgrad[1]) cgrad[1][i*DIM+j] -= fg->grad[1][j] * delta;
        }
      }
      for (i=0; i<phys->dof; i++) {
        for (j=0; j<phys->dof; j++) {
          if(cells[0]<user->cEndInterior) elemMat[cells[0]*dof*dof + i*dof + j] -= cells[0]*1.0;
          if(cells[1]<user->cEndInterior) elemMat[cells[1]*dof*dof + i*dof + j] += cells[1]*1.2;
        }
      }
    }
    ierr = VecRestoreArray(TempVec,&grad);CHKERRQ(ierr);
  }
  ierr = DMRestoreGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);

  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(locX,&x);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

PetscErrorCode ApplyBC(DM dm, PetscReal time, Vec locX, User user)
{
  const char        *name = "Face Sets"; /*Set up in the function DMPlexCreateExodus. is the side set*/
  DM                dmFace;
  IS                idIS;
  const PetscInt    *ids;
  PetscScalar       *x;
  const PetscScalar *facegeom;
  PetscInt          numFS, fs;
  PetscErrorCode    ierr;
  PetscMPIInt       rank;

  PetscFunctionBeginUser;

  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
  ierr = DMPlexGetLabelIdIS(dm, name, &idIS);CHKERRQ(ierr);
 // ISView(idIS, PETSC_VIEWER_STDOUT_SELF);
  if (!idIS) PetscFunctionReturn(0);
  ierr = ISGetLocalSize(idIS, &numFS);CHKERRQ(ierr);
  ierr = ISGetIndices(idIS, &ids);CHKERRQ(ierr);
  ierr = VecGetArrayRead(user->facegeom, &facegeom);CHKERRQ(ierr);
  ierr = VecGetArray(locX, &x);CHKERRQ(ierr);

  for (fs = 0; fs < numFS; ++fs) {
    IS               faceIS;
    const PetscInt   *faces;
    PetscInt         numFaces, f;

    ierr = DMPlexGetStratumIS(dm, name, ids[fs], &faceIS);CHKERRQ(ierr);
    ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
    ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
    for (f = 0; f < numFaces; ++f) {
//      PetscPrintf(PETSC_COMM_SELF, "rank[%d]: ids[%d] = %d, faceIS[%d] = %d, numFaces = %d\n", rank, fs, ids[fs], f, faces[f], numFaces);
      const PetscInt    face = faces[f], *cells;
      const PetscScalar *xI; /*Inner point*/
      PetscScalar       *xG; /*Ghost point*/
      const FaceGeom    *fg;

      ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
      ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
      ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
      ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr);
      if (ids[fs]==1){
        //PetscPrintf(PETSC_COMM_SELF, "Set Inlfow Boundary Condition! \n");
        ierr = BoundaryInflow(time, fg->centroid, fg->normal, xI, xG, user);CHKERRQ(ierr);
//        DM                dmCell;
//        const PetscScalar *cellgeom;
//        const CellGeom    *cgL, *cgR;
//        ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
//        ierr = VecGetArrayRead(user->cellgeom, &cellgeom);CHKERRQ(ierr);
//        ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr);
//        ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr);
//        ierr = PetscPrintf(PETSC_COMM_WORLD,"cells[0] = (%f, %f, %f), cells[1] = (%f, %f, %f)\n",cgL->centroid[0], cgL->centroid[1], cgL->centroid[2],cgR->centroid[0], cgR->centroid[1], cgR->centroid[2]);CHKERRQ(ierr);
      }else if (ids[fs]==2){
        //PetscPrintf(PETSC_COMM_SELF, "Set Outlfow Boundary Condition! \n");
        ierr = BoundaryOutflow(time, fg->centroid, fg->normal, xI, xG, user);CHKERRQ(ierr);
      }else if (ids[fs]==3){
        //PetscPrintf(PETSC_COMM_SELF, "Set Wall Boundary Condition! \n");
        ierr = BoundaryWallflow(time, fg->centroid, fg->normal, xI, xG, user);CHKERRQ(ierr);
      }else {
        SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Wrong type of boundary condition setup!!! \n The set up of the boundary should be: 1 for the inflow, 2 for the outflow, and 3 for the wallflow");
      }
    }
//    PetscPrintf(PETSC_COMM_SELF, " \n");
    ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
    ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
  }
  ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
  ierr = ISRestoreIndices(idIS, &ids);CHKERRQ(ierr);
  ierr = ISDestroy(&idIS);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

static PetscErrorCode  KSPSolve_CGS(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscScalar    rho,rhoold,a,s,b;
  Vec            X,B,V,P,R,RP,T,Q,U,AUQ;
  PetscReal      dp = 0.0;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  /* not sure what residual norm it does use, should use for right preconditioning */

  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  X   = ksp->vec_sol;
  B   = ksp->vec_rhs;
  R   = ksp->work[0];
  RP  = ksp->work[1];
  V   = ksp->work[2];
  T   = ksp->work[3];
  Q   = ksp->work[4];
  P   = ksp->work[5];
  U   = ksp->work[6];
  AUQ = V;

  /* Compute initial preconditioned residual */
  ierr = KSPInitialResidual(ksp,X,V,T,R,B);CHKERRQ(ierr);

  /* Test for nothing to do */
  ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
  KSPCheckNorm(ksp,dp);
  if (ksp->normtype == KSP_NORM_NATURAL) dp *= dp;
  ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->its   = 0;
  ksp->rnorm = dp;
  ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  /* Make the initial Rp == R */
  ierr = VecCopy(R,RP);CHKERRQ(ierr);
  /*  added for Fidap */
  /* Penalize Startup - Isaac Hasbani Trick for CGS
     Since most initial conditions result in a mostly 0 residual,
     we change all the 0 values in the vector RP to the maximum.
  */
  if (ksp->normtype == KSP_NORM_NATURAL) {
    PetscReal   vr0max;
    PetscScalar *tmp_RP=0;
    PetscInt    numnp  =0, *max_pos=0;
    ierr = VecMax(RP, max_pos, &vr0max);CHKERRQ(ierr);
    ierr = VecGetArray(RP, &tmp_RP);CHKERRQ(ierr);
    ierr = VecGetLocalSize(RP, &numnp);CHKERRQ(ierr);
    for (i=0; i<numnp; i++) {
      if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max;
    }
    ierr = VecRestoreArray(RP, &tmp_RP);CHKERRQ(ierr);
  }
  /*  end of addition for Fidap */

  /* Set the initial conditions */
  ierr = VecDot(R,RP,&rhoold);CHKERRQ(ierr);        /* rhoold = (r,rp)      */
  ierr = VecCopy(R,U);CHKERRQ(ierr);
  ierr = VecCopy(R,P);CHKERRQ(ierr);
  ierr = KSP_PCApplyBAorAB(ksp,P,V,T);CHKERRQ(ierr);

  i = 0;
  do {

    ierr = VecDot(V,RP,&s);CHKERRQ(ierr);           /* s <- (v,rp)          */
    KSPCheckDot(ksp,s);
    a    = rhoold / s;                               /* a <- rho / s         */
    ierr = VecWAXPY(Q,-a,V,U);CHKERRQ(ierr);      /* q <- u - a v         */
    ierr = VecWAXPY(T,1.0,U,Q);CHKERRQ(ierr);      /* t <- u + q           */
    ierr = VecAXPY(X,a,T);CHKERRQ(ierr);           /* x <- x + a (u + q)   */
    ierr = KSP_PCApplyBAorAB(ksp,T,AUQ,U);CHKERRQ(ierr);
    ierr = VecAXPY(R,-a,AUQ);CHKERRQ(ierr);       /* r <- r - a K (u + q) */
    ierr = VecDot(R,RP,&rho);CHKERRQ(ierr);         /* rho <- (r,rp)        */
    KSPCheckDot(ksp,rho);
    if (ksp->normtype == KSP_NORM_NATURAL) {
      dp = PetscAbsScalar(rho);
    } else {
      ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
      KSPCheckNorm(ksp,dp);
    }

    ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
    ksp->its++;
    ksp->rnorm = dp;
    ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
    ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
    ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr);
    ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
    if (ksp->reason) break;

    b      = rho / rhoold;                           /* b <- rho / rhoold    */
    ierr   = VecWAXPY(U,b,Q,R);CHKERRQ(ierr);       /* u <- r + b q         */
//.........这里部分代码省略.........

PetscErrorCode ComputeJacobian_LS(DM dm, Vec locX, PetscInt cell, PetscScalar CellValues[], void *ctx)
{
  User              user = (User) ctx;
  Physics           phys = user->model->physics;
  PetscInt          dof = phys->dof;
  const PetscScalar *facegeom, *cellgeom,*x;
  PetscErrorCode    ierr;
  DM                dmFace, dmCell;

  DM                dmGrad = user->dmGrad;
  PetscInt          fStart, fEnd, face, cStart;
  Vec               locGrad, locGradLimiter, Grad;
  /*here the localGradLimiter refers to the gradient that has been multiplied by the limiter function.
   The locGradLimiter is used to construct the uL and uR, and the locGrad is used to caculate the diffusion term*/
  Vec               TempVec; /*a temperal vec for the vector restore*/

  PetscFunctionBeginUser;

  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);

  ierr = DMGetGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);
  ierr = VecDuplicate(Grad, &TempVec);CHKERRQ(ierr);
  ierr = VecCopy(Grad, TempVec);CHKERRQ(ierr);
  /*Backup the original vector and use it to restore the value of dmGrad,
    because I do not want to change the values of the cell gradient*/

  ierr = VecGetArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
  {
    PetscScalar *grad;
    ierr = VecGetArray(Grad,&grad);CHKERRQ(ierr);

    /* Limit interior gradients. Using cell-based loop because it generalizes better to vector limiters. */

      const PetscInt    *faces;
      PetscInt          numFaces,f;
      PetscReal         *cellPhi; /* Scalar limiter applied to each component separately */
      const PetscScalar *cx;
      const CellGeom    *cg;
      PetscScalar       *cgrad;
      PetscInt          i;

      ierr = PetscMalloc(phys->dof*sizeof(PetscScalar),&cellPhi);CHKERRQ(ierr);

      ierr = DMPlexGetConeSize(dm,cell,&numFaces);CHKERRQ(ierr);
      ierr = DMPlexGetCone(dm,cell,&faces);CHKERRQ(ierr);
      ierr = DMPlexPointLocalRead(dm,cell,x,&cx);CHKERRQ(ierr);
      ierr = DMPlexPointLocalRead(dmCell,cell,cellgeom,&cg);CHKERRQ(ierr);
      ierr = DMPlexPointGlobalRef(dmGrad,cell,grad,&cgrad);CHKERRQ(ierr);

      /* Limiter will be minimum value over all neighbors */
      for (i=0; i<dof; i++) {
        cellPhi[i] = PETSC_MAX_REAL;
      }
      for (f=0; f<numFaces; f++) {
        const PetscScalar *ncx;
        const CellGeom    *ncg;
        const PetscInt    *fcells;
        PetscInt          face = faces[f],ncell;
        PetscScalar       v[DIM];
        PetscBool         ghost;
        ierr = IsExteriorGhostFace(dm,face,&ghost);CHKERRQ(ierr);
        if (ghost) continue;
        ierr  = DMPlexGetSupport(dm,face,&fcells);CHKERRQ(ierr);
        ncell = cell == fcells[0] ? fcells[1] : fcells[0];  /*The expression (x ? y : z) has the value of y if x is nonzero, z otherwise */
        ierr  = DMPlexPointLocalRead(dm,ncell,x,&ncx);CHKERRQ(ierr);
        ierr  = DMPlexPointLocalRead(dmCell,ncell,cellgeom,&ncg);CHKERRQ(ierr);
        Waxpy2(-1, cg->centroid, ncg->centroid, v);
        for (i=0; i<dof; i++) {
          /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
          PetscScalar phi,flim = 0.5 * (ncx[i] - cx[i]) / Dot2(&cgrad[i*DIM],v);
          phi        = (*user->LimitGrad)(flim);
          cellPhi[i] = PetscMin(cellPhi[i],phi);
        }
      }
      /* Apply limiter to gradient */
      for (i=0; i<dof; i++) Scale2(cellPhi[i],&cgrad[i*DIM],&cgrad[i*DIM]);

      ierr = PetscFree(cellPhi);CHKERRQ(ierr);

    ierr = VecRestoreArray(Grad,&grad);CHKERRQ(ierr);
  }
  ierr = DMGetLocalVector(dmGrad,&locGradLimiter);CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGradLimiter);CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGradLimiter);CHKERRQ(ierr);

  ierr = VecCopy(TempVec, Grad);CHKERRQ(ierr);/*Restore the vector*/

  ierr = DMGetLocalVector(dmGrad,&locGrad);CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);

  ierr = DMRestoreGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);
  ierr = VecDestroy(&TempVec);CHKERRQ(ierr);

  {
//.........这里部分代码省略.........

PetscErrorCode iniCalcBC(PetscScalar tc, PetscInt Iterc, PetscScalar tf, PetscInt Iterf, PetscInt numTracers, Vec *v, Vec *bcc, Vec *bcf)
{

  PetscErrorCode ierr;
  PetscInt itr;
  PetscViewer fd;
  PetscBool flg;
  PetscInt it;
  PetscScalar myTime;
  PetscScalar zero = 0.0;

/*   v[0]=TR1,v[1]=TR2,... */
  ierr = VecGetArray(TR,&localTR);CHKERRQ(ierr);

  for (itr=0; itr<numTracers; itr++) {
    ierr = VecSet(bcc[itr],zero); CHKERRQ(ierr);
    ierr = VecSet(bcf[itr],zero); CHKERRQ(ierr);    
  }
  ierr = VecGetArray(BCc,&localBCc);CHKERRQ(ierr);
  ierr = VecGetArray(BCf,&localBCf);CHKERRQ(ierr);

  ierr = PetscOptionsGetInt(PETSC_NULL,"-inert_gas_id",&gasID,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate inert gas tracer ID with the -inert_gas_id option");  

  ierr = PetscOptionsGetReal(PETSC_NULL,"-biogeochem_deltat",&DeltaT,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate biogeochemical time step in seconds with the -biogeochem_deltat option");  

  ierr = PetscOptionsHasName(PETSC_NULL,"-periodic_biogeochem_forcing",&periodicBiogeochemForcing);CHKERRQ(ierr);

  if (periodicBiogeochemForcing) {    
    ierr=PetscPrintf(PETSC_COMM_WORLD,"Periodic biogeochemical forcing specified\n");CHKERRQ(ierr);
    ierr = iniPeriodicTimer("periodic_biogeochem_", &biogeochemTimer);CHKERRQ(ierr);
  }
  
/* Grid arrays */

/* Forcing fields */  
  ierr = VecDuplicate(BCc,&Tss);CHKERRQ(ierr);
  ierr = VecDuplicate(BCc,&Sss);CHKERRQ(ierr);  
  ierr = VecDuplicate(BCc,&atmosp);CHKERRQ(ierr);

  if (periodicBiogeochemForcing) {    
    Tssp.firstTime = PETSC_TRUE;
    Sssp.firstTime = PETSC_TRUE;
    atmospp.firstTime = PETSC_TRUE;
  } else {
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Tss.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr);
    ierr = VecLoad(Tss,fd);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);    
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Sss.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr);
    ierr = VecLoad(Sss,fd);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);    
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"atmosp.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr);
    ierr = VecLoad(atmosp,fd);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);    
  }  
  
  ierr = VecGetArray(Tss,&localTss);CHKERRQ(ierr);
  ierr = VecGetArray(Sss,&localSss);CHKERRQ(ierr);
  ierr = VecGetArray(atmosp,&localatmosp);CHKERRQ(ierr);
/*   ierr = PetscPrintf(PETSC_COMM_WORLD,"Done reading surface T, S, and atmospheric pressure\n");CHKERRQ(ierr); */

  ierr = PetscOptionsHasName(PETSC_NULL,"-calc_diagnostics",&calcDiagnostics);CHKERRQ(ierr);
  if (calcDiagnostics) {    
/*   Read T and S */
    ierr = VecDuplicate(TR,&Ts);CHKERRQ(ierr);
    ierr = VecDuplicate(TR,&Ss);CHKERRQ(ierr);  
    if (periodicBiogeochemForcing) {    
      Tsp.firstTime = PETSC_TRUE;
      Ssp.firstTime = PETSC_TRUE;

      ierr = interpPeriodicVector(tc,&Ts,biogeochemTimer.cyclePeriod,biogeochemTimer.numPerPeriod,biogeochemTimer.tdp,&Tsp,"Ts_");
      ierr = interpPeriodicVector(tc,&Ss,biogeochemTimer.cyclePeriod,biogeochemTimer.numPerPeriod,biogeochemTimer.tdp,&Ssp,"Ss_");	
      
    } else {
	  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ts.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr);
	  ierr = VecLoad(Ts,fd);CHKERRQ(ierr);
	  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);    
	  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ss.petsc",FILE_MODE_READ,&fd);CHKERRQ(ierr);
	  ierr = VecLoad(Ss,fd);CHKERRQ(ierr);
	  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);    
    }  
    
    ierr = VecGetArray(Ts,&localTs);CHKERRQ(ierr);
    ierr = VecGetArray(Ss,&localSs);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Done reading T/S\n");CHKERRQ(ierr);

    ierr = VecDuplicate(TR,&TReqdiag);CHKERRQ(ierr);
    ierr = VecSet(TReqdiag,zero);CHKERRQ(ierr);
    ierr = VecGetArray(TReqdiag,&localTReqdiag);CHKERRQ(ierr);
  
/*Data for diagnostics */
    ierr = iniStepTimer("diag_", Iter0, &diagTimer);CHKERRQ(ierr);
	ierr = PetscPrintf(PETSC_COMM_WORLD,"Diagnostics will be computed starting at (and including) time step: %d\n", diagTimer.startTimeStep);CHKERRQ(ierr);	
	ierr = PetscPrintf(PETSC_COMM_WORLD,"Diagnostics will be computed over %d time steps\n", diagTimer.numTimeSteps);CHKERRQ(ierr);	
    
	ierr = VecDuplicate(TR,&TRsatanomdiag);CHKERRQ(ierr);
	ierr = VecSet(TRsatanomdiag,zero);CHKERRQ(ierr);
	ierr = VecGetArray(TRsatanomdiag,&localTRsatanomdiag);CHKERRQ(ierr);
	ierr = VecDuplicate(TR,&TRsatanomdiagavg);CHKERRQ(ierr);
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  PetscErrorCode  ierr;
  DM              dm;
  Vec             vec,vecLocal1,vecLocal2;
  PetscScalar     *a,***a1,***a2,expected;
  PetscInt        startx,starty,nx,ny,i,j,d,is,js,dof0,dof1,dof2,dofTotal,stencilWidth,Nx,Ny;
  DMBoundaryType  boundaryTypex,boundaryTypey;
  PetscMPIInt     rank;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  dof0 = 1;
  dof1 = 1;
  dof2 = 1;
  stencilWidth = 2;
  ierr = DMStagCreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,4,4,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,DMSTAG_STENCIL_BOX,stencilWidth,NULL,NULL,&dm);CHKERRQ(ierr);
  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
  ierr = DMSetUp(dm);CHKERRQ(ierr);
  ierr = DMStagGetDOF(dm,&dof0,&dof1,&dof2,NULL);CHKERRQ(ierr);
  dofTotal = dof0 + 2*dof1 + dof2;
  ierr = DMStagGetStencilWidth(dm,&stencilWidth);CHKERRQ(ierr);

  ierr = DMCreateLocalVector(dm,&vecLocal1);CHKERRQ(ierr);
  ierr = VecDuplicate(vecLocal1,&vecLocal2);CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(dm,&vec);CHKERRQ(ierr);
  ierr = VecSet(vec,1.0);CHKERRQ(ierr);
  ierr = VecSet(vecLocal1,0.0);CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal1);CHKERRQ(ierr);

  ierr = DMStagGetCorners(dm,&startx,&starty,NULL,&nx,&ny,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = DMStagVecGetArrayDOFRead(dm,vecLocal1,&a1);CHKERRQ(ierr);
  ierr = DMStagVecGetArrayDOF(dm,vecLocal2,&a2);CHKERRQ(ierr);
  for (j=starty; j<starty + ny; ++j) {
    for (i=startx; i<startx + nx; ++i) {
      for (d=0; d<dofTotal; ++d) {
        if (a1[j][i][d] != 1.0) {
          PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a1[j][i][d],1.0);CHKERRQ(ierr);
        }
        a2[j][i][d] = 0.0;
        for (js = -stencilWidth; js <= stencilWidth; ++js) {
          for (is = -stencilWidth; is <= stencilWidth; ++is) {
            a2[j][i][d] += a1[j+js][i+is][d];
          }
        }
      }
    }
  }
  ierr = DMStagVecRestoreArrayDOFRead(dm,vecLocal1,&a1);CHKERRQ(ierr);
  ierr = DMStagVecRestoreArrayDOF(dm,vecLocal2,&a2);CHKERRQ(ierr);

  DMLocalToGlobalBegin(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);
  DMLocalToGlobalEnd(dm,vecLocal2,INSERT_VALUES,vec);CHKERRQ(ierr);

  /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
  ierr = DMStagGetBoundaryTypes(dm,&boundaryTypex,&boundaryTypey,NULL);CHKERRQ(ierr);
  if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
    ierr = VecGetArray(vec,&a);CHKERRQ(ierr);
    expected = 1.0; for(d=0;d<2;++d) expected *= (2*stencilWidth+1);
    for (i=0; i<ny*nx*dofTotal; ++i) {
      if (a[i] != expected) {
        ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a[i],expected);CHKERRQ(ierr);
      }
    }
    ierr = VecRestoreArray(vec,&a);CHKERRQ(ierr);
  } else {
    ierr = DMStagVecGetArrayDOFRead(dm,vecLocal2,&a2);CHKERRQ(ierr);
    ierr = DMStagGetGlobalSizes(dm,&Nx,&Ny,NULL);CHKERRQ(ierr);
    if (stencilWidth > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Non-periodic check implemented assuming stencilWidth = 1");
      for (j=starty; j<starty + ny; ++j) {
        for (i=startx; i<startx + nx; ++i) {
          PetscInt  dd,extra[2];
          PetscBool bnd[2];
          bnd[0] = (PetscBool)((i == 0 || i == Nx-1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
          bnd[1] = (PetscBool)((j == 0 || j == Ny-1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
          extra[0] = i == Nx-1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
          extra[1] = j == Ny-1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
          { /* vertices */
            PetscScalar expected = 1.0;
            for (dd=0;dd<2;++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2*stencilWidth + 1);
            for (d=0; d<dof0; ++d) {
              if (a2[j][i][d] != expected) {
                ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,d,a2[j][i][d],expected);CHKERRQ(ierr);
              }
            }
          }
          { /* down edges */
            PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2*stencilWidth + 1);
            expected *= ((bnd[0] ? 1 : 2) * stencilWidth + 1);
            for (d=dof0; d<dof0+dof1; ++d) {
              if (a2[j][i][d] != expected) {
                ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Element (%D,%D)[%D] Unexpected value %g (expecting %g)\n",rank,i,j,d,a2[j][i][d],expected);CHKERRQ(ierr);
              }
            }
          }
          { /* left edges */
            PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2*stencilWidth + 1);
            expected *= ((bnd[1] ? 1 : 2) * stencilWidth + 1);
//.........这里部分代码省略.........

/*@
   MatNullSpaceCreateRigidBody - create rigid body modes from coordinates

   Collective on Vec

   Input Argument:
.  coords - block of coordinates of each node, must have block size set

   Output Argument:
.  sp - the null space

   Level: advanced

   Notes:
    If you are solving an elasticity problems you should likely use this, in conjunction with ee MatSetNearNullspace(), to provide information that 
           the PCGAMG preconditioner can use to construct a much more efficient preconditioner.

           If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with MatSetNullspace() to
           provide this information to the linear solver so it can handle the null space appropriately in the linear solution.


.seealso: MatNullSpaceCreate(), MatSetNearNullspace(), MatSetNullspace()
@*/
PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
{
  PetscErrorCode    ierr;
  const PetscScalar *x;
  PetscScalar       *v[6],dots[5];
  Vec               vec[6];
  PetscInt          n,N,dim,nmodes,i,j;
  PetscReal         sN;

  PetscFunctionBegin;
  ierr = VecGetBlockSize(coords,&dim);CHKERRQ(ierr);
  ierr = VecGetLocalSize(coords,&n);CHKERRQ(ierr);
  ierr = VecGetSize(coords,&N);CHKERRQ(ierr);
  n   /= dim;
  N   /= dim;
  sN = 1./PetscSqrtReal((PetscReal)N);
  switch (dim) {
  case 1:
    ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);CHKERRQ(ierr);
    break;
  case 2:
  case 3:
    nmodes = (dim == 2) ? 3 : 6;
    ierr   = VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);CHKERRQ(ierr);
    ierr   = VecSetSizes(vec[0],dim*n,dim*N);CHKERRQ(ierr);
    ierr   = VecSetBlockSize(vec[0],dim);CHKERRQ(ierr);
    ierr   = VecSetUp(vec[0]);CHKERRQ(ierr);
    for (i=1; i<nmodes; i++) {ierr = VecDuplicate(vec[0],&vec[i]);CHKERRQ(ierr);}
    for (i=0; i<nmodes; i++) {ierr = VecGetArray(vec[i],&v[i]);CHKERRQ(ierr);}
    ierr = VecGetArrayRead(coords,&x);CHKERRQ(ierr);
    for (i=0; i<n; i++) {
      if (dim == 2) {
        v[0][i*2+0] = sN;
        v[0][i*2+1] = 0.;
        v[1][i*2+0] = 0.;
        v[1][i*2+1] = sN;
        /* Rotations */
        v[2][i*2+0] = -x[i*2+1];
        v[2][i*2+1] = x[i*2+0];
      } else {
        v[0][i*3+0] = sN;
        v[0][i*3+1] = 0.;
        v[0][i*3+2] = 0.;
        v[1][i*3+0] = 0.;
        v[1][i*3+1] = sN;
        v[1][i*3+2] = 0.;
        v[2][i*3+0] = 0.;
        v[2][i*3+1] = 0.;
        v[2][i*3+2] = sN;

        v[3][i*3+0] = x[i*3+1];
        v[3][i*3+1] = -x[i*3+0];
        v[3][i*3+2] = 0.;
        v[4][i*3+0] = 0.;
        v[4][i*3+1] = -x[i*3+2];
        v[4][i*3+2] = x[i*3+1];
        v[5][i*3+0] = x[i*3+2];
        v[5][i*3+1] = 0.;
        v[5][i*3+2] = -x[i*3+0];
      }
    }
    for (i=0; i<nmodes; i++) {ierr = VecRestoreArray(vec[i],&v[i]);CHKERRQ(ierr);}
    ierr = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr);
    for (i=dim; i<nmodes; i++) {
      /* Orthonormalize vec[i] against vec[0:i-1] */
      ierr = VecMDot(vec[i],i,vec,dots);CHKERRQ(ierr);
      for (j=0; j<i; j++) dots[j] *= -1.;
      ierr = VecMAXPY(vec[i],i,dots,vec);CHKERRQ(ierr);
      ierr = VecNormalize(vec[i],NULL);CHKERRQ(ierr);
    }
    ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);CHKERRQ(ierr);
    for (i=0; i<nmodes; i++) {ierr = VecDestroy(&vec[i]);CHKERRQ(ierr);}
  }
  PetscFunctionReturn(0);
}

int main(int argc,char **argv)
{
  TS             ts;            /* ODE integrator */
  Vec            U;             /* solution will be stored here */
  Mat            A;             /* Jacobian matrix */
  PetscErrorCode ierr;
  PetscMPIInt    size;
  PetscInt       n = 2;
  AppCtx         ctx;
  PetscScalar    *u;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Initialize program
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Create necessary matrix and vectors
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);

  ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Set runtime options
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Reaction options","");CHKERRQ(ierr);
  {
    ctx.omega_s = 1.0;
    ierr        = PetscOptionsScalar("-omega_s","","",ctx.omega_s,&ctx.omega_s,NULL);CHKERRQ(ierr);
    ctx.H       = 1.0;
    ierr        = PetscOptionsScalar("-H","","",ctx.H,&ctx.H,NULL);CHKERRQ(ierr);
    ctx.E       = 1.0;
    ierr        = PetscOptionsScalar("-E","","",ctx.E,&ctx.E,NULL);CHKERRQ(ierr);
    ctx.V       = 1.0;
    ierr        = PetscOptionsScalar("-V","","",ctx.V,&ctx.V,NULL);CHKERRQ(ierr);
    ctx.X       = 1.0;
    ierr        = PetscOptionsScalar("-X","","",ctx.X,&ctx.X,NULL);CHKERRQ(ierr);

    ierr = VecGetArray(U,&u);CHKERRQ(ierr);
    u[0] = 1;
    u[1] = .7;
    ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
    ierr = PetscOptionsGetVec(NULL,NULL,"-initial",U,NULL);CHKERRQ(ierr);
  }
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&ctx.rand);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(ctx.rand);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create timestepping solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
  ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr);
  ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);CHKERRQ(ierr);
  ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set initial conditions
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSetSolution(ts,U);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set solver options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSetDuration(ts,100000,2000.0);CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
  ierr = TSSetInitialTimeStep(ts,0.0,.001);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 = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&U);CHKERRQ(ierr);
  ierr = TSDestroy(&ts);CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&ctx.rand);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

//.........这里部分代码省略.........
  CHKERRQ( DAGetInfo(da,0, &mx, &my, &mz, 0,0,0,0,0,0,0) ); 

  double acx,acy,acz;
  double hx = 1.0/((double)Ns);

  // allocate for temporary buffers ...
  // unsigned int elemSize = Ns*Ns*Ns;
  // std::cout << "Elem size is " << elemSize << std::endl;
  // unsigned int nodeSize = (Ns+1)*(Ns+1)*(Ns+1);

  // Now set the activation ...
  unsigned int numSteps = (unsigned int)(ceil(( ti.stop - ti.start)/ti.step));

 // Vec tauVec;

  // PetscScalar ***tauArray;

  unsigned int paramTimeSteps = (unsigned int)(ceil(( (double)(numSteps))/ ((double)(2*numParams)) ));

  /*
  for (int b=0; b<numParams; b++) {
    std::vector<Vec> tau;
    unsigned int tBegin = paramTimeSteps*b;
    unsigned int tEnd   = tBegin + numSteps/2; // paramTimeSteps*(b+2);

    // std::cout << "For param " << b << ": Time step range is " << tBegin << " -> " << tEnd << std::endl; 
    for (unsigned int t=0; t<numSteps+1; t++) {
      double newTime = (dt*(t-tBegin)*numSteps)/((double)(paramTimeSteps));
     // double fff = 0.0;
      CHKERRQ( DACreateGlobalVector(da, &tauVec) );
      CHKERRQ( VecSet( tauVec, 0.0));

      if ( (t>=tBegin) && (t<=tEnd)) {
        CHKERRQ(DAVecGetArray(da, tauVec, &tauArray));
        for (int k = z; k < z + p ; k++) {
          for (int j = y; j < y + n; j++) {
            for (int i = x; i < x + m; i++) {
              acx = (i)*hx; acy = (j)*hx; acz = (k)*hx;
              tauArray[k][j][i] = sin(M_PI*newTime)*cos(2*M_PI*acx)*cos(2*M_PI*acy)*cos(2*M_PI*acz);
            }
          }
        }
        CHKERRQ( DAVecRestoreArray ( da, tauVec, &tauArray ) );
      }
      tau.push_back(tauVec);
    }
    fBasis.push_back(tau);
  }
  */
 // std::cout << "Finished setting basis" << std::endl;

  /*
  // Set initial velocity ...
  CHKERRQ(DAVecGetArray(da, initialVelocity, &solArray));

  for (int k = z; k < z + p ; k++) {
  for (int j = y; j < y + n; j++) {
  for (int i = x; i < x + m; i++) {
  acx = (i)*hx; acy = (j)*hx; acz = (k)*hx;
  solArray[k][j][i] = M_PI*cos(2*M_PI*acx)*cos(2*M_PI*acy)*cos(2*M_PI*acz);
  }
  }
  }
  CHKERRQ( DAVecRestoreArray ( da, initialVelocity, &solArray ) );
  */