1 static char help[] = "Tests DMSLICED operations\n\n"; 2 3 #include <petscdmsliced.h> 4 5 int main(int argc,char *argv[]) 6 { 7 char mat_type[256] = MATAIJ; /* default matrix type */ 8 PetscErrorCode ierr; 9 MPI_Comm comm; 10 PetscMPIInt rank,size; 11 DM slice; 12 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}; 13 PetscReal alpha =1,K=1,rho0=1,u0=0,sigma=0.2; 14 PetscBool useblock=PETSC_TRUE; 15 PetscScalar *xx; 16 Mat A; 17 Vec x,b,lf; 18 19 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 20 comm = PETSC_COMM_WORLD; 21 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 22 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 23 24 ierr = PetscOptionsBegin(comm,0,"Options for DMSliced test",0);CHKERRQ(ierr); 25 { 26 ierr = PetscOptionsInt("-n","Global number of nodes","",N,&N,NULL);CHKERRQ(ierr); 27 ierr = PetscOptionsInt("-bs","Block size (1 or 2)","",bs,&bs,NULL);CHKERRQ(ierr); 28 if (bs != 1) { 29 if (bs != 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Block size must be 1 or 2"); 30 ierr = PetscOptionsReal("-alpha","Inverse time step for wave operator","",alpha,&alpha,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsReal("-K","Bulk modulus of compressibility","",K,&K,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsReal("-rho0","Reference density","",rho0,&rho0,NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsReal("-u0","Reference velocity","",u0,&u0,NULL);CHKERRQ(ierr); 34 ierr = PetscOptionsReal("-sigma","Width of Gaussian density perturbation","",sigma,&sigma,NULL);CHKERRQ(ierr); 35 ierr = PetscOptionsBool("-block","Use block matrix assembly","",useblock,&useblock,NULL);CHKERRQ(ierr); 36 } 37 ierr = PetscOptionsString("-sliced_mat_type","Matrix type to use (aij or baij)","",mat_type,mat_type,sizeof(mat_type),NULL);CHKERRQ(ierr); 38 } 39 ierr = PetscOptionsEnd();CHKERRQ(ierr); 40 41 /* Split ownership, set up periodic grid in 1D */ 42 n = PETSC_DECIDE; 43 ierr = PetscSplitOwnership(comm,&n,&N);CHKERRQ(ierr); 44 rstart = 0; 45 ierr = MPI_Scan(&n,&rstart,1,MPIU_INT,MPI_SUM,comm);CHKERRMPI(ierr); 46 rstart -= n; 47 ghosts[0] = (N+rstart-1)%N; 48 ghosts[1] = (rstart+n)%N; 49 50 ierr = PetscMalloc2(n,&d_nnz,n,&o_nnz);CHKERRQ(ierr); 51 for (i=0; i<n; i++) { 52 if (size > 1 && (i==0 || i==n-1)) { 53 d_nnz[i] = 2; 54 o_nnz[i] = 1; 55 } else { 56 d_nnz[i] = 3; 57 o_nnz[i] = 0; 58 } 59 } 60 ierr = DMSlicedCreate(comm,bs,n,2,ghosts,d_nnz,o_nnz,&slice);CHKERRQ(ierr); /* Currently does not copy X_nnz so we can't free them until after DMSlicedGetMatrix */ 61 62 if (!useblock) {ierr = DMSlicedSetBlockFills(slice,dfill,ofill);CHKERRQ(ierr);} /* Irrelevant for baij formats */ 63 ierr = DMSetMatType(slice,mat_type);CHKERRQ(ierr); 64 ierr = DMCreateMatrix(slice,&A);CHKERRQ(ierr); 65 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 66 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 67 68 ierr = DMCreateGlobalVector(slice,&x);CHKERRQ(ierr); 69 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 70 71 ierr = VecGhostGetLocalForm(x,&lf);CHKERRQ(ierr); 72 ierr = VecGetSize(lf,&m);CHKERRQ(ierr); 73 if (m != (n+2)*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size of local form %D, expected %D",m,(n+2)*bs); 74 ierr = VecGetArray(lf,&xx);CHKERRQ(ierr); 75 for (i=0; i<n; i++) { 76 PetscInt row[2],col[9],im,ip; 77 PetscScalar v[12]; 78 const PetscReal xref = 2.0*(rstart+i)/N - 1; /* [-1,1] */ 79 const PetscReal h = 1.0/N; /* grid spacing */ 80 im = (i==0) ? n : i-1; 81 ip = (i==n-1) ? n+1 : i+1; 82 switch (bs) { 83 case 1: /* Laplacian with periodic boundaries */ 84 col[0] = im; col[1] = i; col[2] = ip; 85 v[0] = -h; v[1] = 2*h; v[2] = -h; 86 ierr = MatSetValuesLocal(A,1,&i,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 87 xx[i] = PetscSinReal(xref*PETSC_PI); 88 break; 89 case 2: /* Linear acoustic wave operator in variables [rho, u], central differences, periodic, timestep 1/alpha */ 90 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; 91 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; 92 if (useblock) { 93 row[0] = i; col[0] = im; col[1] = i; col[2] = ip; 94 ierr = MatSetValuesBlockedLocal(A,1,row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 95 } else { 96 row[0] = 2*i; row[1] = 2*i+1; 97 col[0] = 2*im; col[1] = 2*im+1; col[2] = 2*i; col[3] = 2*ip; col[4] = 2*ip+1; 98 v[3] = v[4]; v[4] = v[5]; /* pack values in first row */ 99 ierr = MatSetValuesLocal(A,1,row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 100 col[2] = 2*i+1; 101 v[8] = v[9]; v[9] = v[10]; v[10] = v[11]; /* pack values in second row */ 102 ierr = MatSetValuesLocal(A,1,row+1,5,col,v+6,INSERT_VALUES);CHKERRQ(ierr); 103 } 104 /* Set current state (gaussian density perturbation) */ 105 xx[2*i] = 0.2*PetscExpReal(-PetscSqr(xref)/(2*PetscSqr(sigma))); 106 xx[2*i+1] = 0; 107 break; 108 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"not implemented for block size %D",bs); 109 } 110 } 111 ierr = VecRestoreArray(lf,&xx);CHKERRQ(ierr); 112 ierr = VecGhostRestoreLocalForm(x,&lf);CHKERRQ(ierr); 113 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 114 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 115 116 ierr = MatMult(A,x,b);CHKERRQ(ierr); 117 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 118 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 119 ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 120 121 /* Update the ghosted values, view the result on rank 0. */ 122 ierr = VecGhostUpdateBegin(b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 123 ierr = VecGhostUpdateEnd(b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 124 if (!rank) { 125 ierr = VecGhostGetLocalForm(b,&lf);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"Local form of b on rank 0, last two nodes are ghost nodes\n");CHKERRQ(ierr); 127 ierr = VecView(lf,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 128 ierr = VecGhostRestoreLocalForm(b,&lf);CHKERRQ(ierr); 129 } 130 131 ierr = DMDestroy(&slice);CHKERRQ(ierr); 132 ierr = VecDestroy(&x);CHKERRQ(ierr); 133 ierr = VecDestroy(&b);CHKERRQ(ierr); 134 ierr = MatDestroy(&A);CHKERRQ(ierr); 135 ierr = PetscFinalize(); 136 return ierr; 137 } 138 139 /*TEST 140 141 test: 142 nsize: 2 143 args: -bs 2 -block 0 -sliced_mat_type baij -alpha 10 -u0 0.1 144 145 test: 146 suffix: 2 147 nsize: 2 148 args: -bs 2 -block 1 -sliced_mat_type aij -alpha 10 -u0 0.1 149 150 test: 151 suffix: 3 152 nsize: 2 153 args: -bs 2 -block 0 -sliced_mat_type aij -alpha 10 -u0 0.1 154 155 TEST*/ 156