Actual source code: ex30.c
petsc-3.4.2 2013-07-02

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

      3: #include <petscdmsliced.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:   DM             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:   PetscBool      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 DMSliced test",0);
     27:   {
     28:     PetscOptionsInt("-n","Global number of nodes","",N,&N,NULL);
     29:     PetscOptionsInt("-bs","Block size (1 or 2)","",bs,&bs,NULL);
     30:     if (bs != 1) {
     31:       if (bs != 2) SETERRQ(PETSC_COMM_WORLD,1,"Block size must be 1 or 2");
     32:       PetscOptionsReal("-alpha","Inverse time step for wave operator","",alpha,&alpha,NULL);
     33:       PetscOptionsReal("-K","Bulk modulus of compressibility","",K,&K,NULL);
     34:       PetscOptionsReal("-rho0","Reference density","",rho0,&rho0,NULL);
     35:       PetscOptionsReal("-u0","Reference velocity","",u0,&u0,NULL);
     36:       PetscOptionsReal("-sigma","Width of Gaussian density perturbation","",sigma,&sigma,NULL);
     37:       PetscOptionsBool("-block","Use block matrix assembly","",useblock,&useblock,NULL);
     38:     }
     39:     PetscOptionsString("-sliced_mat_type","Matrix type to use (aij or baij)","",mat_type,mat_type,sizeof(mat_type),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:   PetscMalloc2(n,PetscInt,&d_nnz,n,PetscInt,&o_nnz);
     53:   for (i=0; i<n; i++) {
     54:     if (size > 1 && (i==0 || i==n-1)) {
     55:       d_nnz[i] = 2;
     56:       o_nnz[i] = 1;
     57:     } else {
     58:       d_nnz[i] = 3;
     59:       o_nnz[i] = 0;
     60:     }
     61:   }
     62:   DMSlicedCreate(comm,bs,n,2,ghosts,d_nnz,o_nnz,&slice); /* Currently does not copy X_nnz so we can't free them until after DMSlicedGetMatrix */

     64:   if (!useblock) {DMSlicedSetBlockFills(slice,dfill,ofill);} /* Irrelevant for baij formats */
     65:   DMCreateMatrix(slice,mat_type,&A);
     66:   PetscFree2(d_nnz,o_nnz);
     67:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

     69:   DMCreateGlobalVector(slice,&x);
     70:   VecDuplicate(x,&b);

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

    117:   MatMult(A,x,b);
    118:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);
    119:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    120:   VecView(b,PETSC_VIEWER_STDOUT_WORLD);

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

    132:   DMDestroy(&slice);
    133:   VecDestroy(&x);
    134:   VecDestroy(&b);
    135:   MatDestroy(&A);
    136:   PetscFinalize();
    137:   return 0;
    138: }