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