Actual source code: ex30.c

  1: static char help[] = "Tests Sliced operations\n\n";

 3:  #include petscda.h

  7: int main(int argc,char *argv[])
  8: {
  9:   char mat_type[256]  = "aij";  /* default matrix type */
 11:   MPI_Comm       comm;
 12:   PetscMPIInt    rank,size;
 13:   Sliced         slice;
 14:   PetscInt       i,bs=1,N=5,n,m,rstart,ghosts[2],*d_nnz,*o_nnz,dfill[4]={1,0,0,1},ofill[4]={1,1,1,1};
 15:   PetscReal      alpha=1,K=1,rho0=1,u0=0,sigma=0.2;
 16:   PetscTruth     useblock=PETSC_TRUE;
 17:   PetscScalar    *xx;
 18:   Mat            A;
 19:   Vec            x,b,lf;

 21:   PetscInitialize(&argc,&argv,0,help);
 22:   comm = PETSC_COMM_WORLD;
 23:   MPI_Comm_size(comm,&size);
 24:   MPI_Comm_rank(comm,&rank);

 26:   PetscOptionsBegin(comm,0,"Options for Sliced test",0);
 27:   {
 28:     PetscOptionsInt("-n","Global number of nodes","",N,&N,PETSC_NULL);
 29:     PetscOptionsInt("-bs","Block size (1 or 2)","",bs,&bs,PETSC_NULL);
 30:     if (bs != 1) {
 31:       if (bs != 2) SETERRQ(1,"Block size must be 1 or 2");
 32:       PetscOptionsReal("-alpha","Inverse time step for wave operator","",alpha,&alpha,PETSC_NULL);
 33:       PetscOptionsReal("-K","Bulk modulus of compressibility","",K,&K,PETSC_NULL);
 34:       PetscOptionsReal("-rho0","Reference density","",rho0,&rho0,PETSC_NULL);
 35:       PetscOptionsReal("-u0","Reference velocity","",u0,&u0,PETSC_NULL);
 36:       PetscOptionsReal("-sigma","Width of Gaussian density perturbation","",sigma,&sigma,PETSC_NULL);
 37:       PetscOptionsTruth("-block","Use block matrix assembly","",useblock,&useblock,PETSC_NULL);
 38:     }
 39:     PetscOptionsString("-sliced_mat_type","Matrix type to use (aij or baij)","",mat_type,mat_type,sizeof mat_type,PETSC_NULL);
 40:   }
 41:   PetscOptionsEnd();

 43:   /* Split ownership, set up periodic grid in 1D */
 44:   n = PETSC_DECIDE;
 45:   PetscSplitOwnership(comm,&n,&N);
 46:   rstart = 0;
 47:   MPI_Scan(&n,&rstart,1,MPIU_INT,MPI_SUM,comm);
 48:   rstart -= n;
 49:   ghosts[0] = (N+rstart-1)%N;
 50:   ghosts[1] = (rstart+n)%N;

 52:   SlicedCreate(comm,&slice);
 53:   SlicedSetGhosts(slice,bs,n,2,ghosts);
 54:   PetscMalloc2(n,PetscInt,&d_nnz,n,PetscInt,&o_nnz);
 55:   for (i=0; i<n; i++) {
 56:     if (size > 1 && (i==0 || i==n-1)) {
 57:       d_nnz[i] = 2;
 58:       o_nnz[i] = 1;
 59:     } else {
 60:       d_nnz[i] = 3;
 61:       o_nnz[i] = 0;
 62:     }
 63:   }
 64:   SlicedSetPreallocation(slice,0,d_nnz,0,o_nnz); /* Currently does not copy X_nnz so we can't free them until after SlicedGetMatrix */

 66:   if (!useblock) {SlicedSetBlockFills(slice,dfill,ofill);} /* Irrelevant for baij formats */
 67:   SlicedGetMatrix(slice,mat_type,&A);
 68:   PetscFree2(d_nnz,o_nnz);
 69:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

 71:   SlicedCreateGlobalVector(slice,&x);
 72:   VecDuplicate(x,&b);

 74:   VecGhostGetLocalForm(x,&lf);
 75:   VecGetSize(lf,&m);
 76:   if (m != (n+2)*bs) SETERRQ2(1,"size of local form %D, expected %D",m,(n+2)*bs);
 77:   VecGetArray(lf,&xx);
 78:   for (i=0; i<n; i++) {
 79:     PetscInt row[2],col[9],im,ip;
 80:     PetscScalar v[12];
 81:     const PetscReal xref = 2.0*(rstart+i)/N - 1; /* [-1,1] */
 82:     const PetscReal h = 1.0/N;                   /* grid spacing */
 83:     im = (i==0) ? n : i-1;
 84:     ip = (i==n-1) ? n+1 : i+1;
 85:     switch (bs) {
 86:       case 1:                   /* Laplacian with periodic boundaries */
 87:         col[0] = im;         col[1] = i;        col[2] = ip;
 88:           v[0] = -h;           v[1] = 2*h;        v[2] = -h;
 89:         MatSetValuesLocal(A,1,&i,3,col,v,INSERT_VALUES);
 90:         xx[i] = sin(xref*PETSC_PI);
 91:         break;
 92:       case 2:                   /* Linear acoustic wave operator in variables [rho, u], central differences, periodic, timestep 1/alpha */
 93:         v[0] = -0.5*u0;   v[1] = -0.5*K;      v[2] = alpha; v[3] = 0;       v[4] = 0.5*u0;    v[5] = 0.5*K;
 94:         v[6] = -0.5/rho0; v[7] = -0.5*u0;     v[8] = 0;     v[9] = alpha;   v[10] = 0.5/rho0; v[11] = 0.5*u0;
 95:         if (useblock) {
 96:           row[0] = i; col[0] = im; col[1] = i; col[2] = ip;
 97:           MatSetValuesBlockedLocal(A,1,row,3,col,v,INSERT_VALUES);
 98:         } else {
 99:           row[0] = 2*i; row[1] = 2*i+1;
100:           col[0] = 2*im; col[1] = 2*im+1; col[2] = 2*i; col[3] = 2*ip; col[4] = 2*ip+1;
101:           v[3] = v[4]; v[4] = v[5];                                                     /* pack values in first row */
102:           MatSetValuesLocal(A,1,row,5,col,v,INSERT_VALUES);
103:           col[2] = 2*i+1;
104:           v[8] = v[9]; v[9] = v[10]; v[10] = v[11];                                     /* pack values in second row */
105:           MatSetValuesLocal(A,1,row+1,5,col,v+6,INSERT_VALUES);
106:         }
107:         /* Set current state (gaussian density perturbation) */
108:         xx[2*i] = 0.2*exp(-PetscSqr(xref)/(2*PetscSqr(sigma)));
109:         xx[2*i+1] = 0;
110:         break;
111:       default: SETERRQ1(1,"not implemented for block size %D",bs);
112:     }
113:   }
114:   VecRestoreArray(lf,&xx);
115:   VecGhostRestoreLocalForm(x,&lf);
116:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
117:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

119:   MatMult(A,x,b);
120:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);
121:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
122:   VecView(b,PETSC_VIEWER_STDOUT_WORLD);

124:   /* Update the ghosted values, view the result on rank 0. */
125:   VecGhostUpdateBegin(b,INSERT_VALUES,SCATTER_FORWARD);
126:   VecGhostUpdateEnd(b,INSERT_VALUES,SCATTER_FORWARD);
127:   if (!rank) {
128:     VecGhostGetLocalForm(b,&lf);
129:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"Local form of b on rank 0, last two nodes are ghost nodes\n");
130:     VecView(lf,PETSC_VIEWER_STDOUT_SELF);
131:     VecGhostRestoreLocalForm(b,&lf);
132:   }

134:   SlicedDestroy(slice);
135:   VecDestroy(x);
136:   VecDestroy(b);
137:   MatDestroy(A);
138:   PetscFinalize();
139:   return 0;
140: }