Actual source code: ex30.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Tests DMSLICED operations\n\n";
3: #include <petscdmsliced.h>
5: int main(int argc,char *argv[])
6: {
7: char mat_type[256] = MATAIJ; /* default matrix type */
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;
19: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
20: comm = PETSC_COMM_WORLD;
21: MPI_Comm_size(comm,&size);
22: MPI_Comm_rank(comm,&rank);
24: PetscOptionsBegin(comm,0,"Options for DMSliced test",0);
25: {
26: PetscOptionsInt("-n","Global number of nodes","",N,&N,NULL);
27: PetscOptionsInt("-bs","Block size (1 or 2)","",bs,&bs,NULL);
28: if (bs != 1) {
29: if (bs != 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Block size must be 1 or 2");
30: PetscOptionsReal("-alpha","Inverse time step for wave operator","",alpha,&alpha,NULL);
31: PetscOptionsReal("-K","Bulk modulus of compressibility","",K,&K,NULL);
32: PetscOptionsReal("-rho0","Reference density","",rho0,&rho0,NULL);
33: PetscOptionsReal("-u0","Reference velocity","",u0,&u0,NULL);
34: PetscOptionsReal("-sigma","Width of Gaussian density perturbation","",sigma,&sigma,NULL);
35: PetscOptionsBool("-block","Use block matrix assembly","",useblock,&useblock,NULL);
36: }
37: PetscOptionsString("-sliced_mat_type","Matrix type to use (aij or baij)","",mat_type,mat_type,sizeof(mat_type),NULL);
38: }
39: PetscOptionsEnd();
41: /* Split ownership, set up periodic grid in 1D */
42: n = PETSC_DECIDE;
43: PetscSplitOwnership(comm,&n,&N);
44: rstart = 0;
45: MPI_Scan(&n,&rstart,1,MPIU_INT,MPI_SUM,comm);
46: rstart -= n;
47: ghosts[0] = (N+rstart-1)%N;
48: ghosts[1] = (rstart+n)%N;
50: PetscMalloc2(n,&d_nnz,n,&o_nnz);
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: 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 */
62: if (!useblock) {DMSlicedSetBlockFills(slice,dfill,ofill);} /* Irrelevant for baij formats */
63: DMSetMatType(slice,mat_type);
64: DMCreateMatrix(slice,&A);
65: PetscFree2(d_nnz,o_nnz);
66: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
68: DMCreateGlobalVector(slice,&x);
69: VecDuplicate(x,&b);
71: VecGhostGetLocalForm(x,&lf);
72: VecGetSize(lf,&m);
73: if (m != (n+2)*bs) SETERRQ2(PETSC_COMM_SELF,1,"size of local form %D, expected %D",m,(n+2)*bs);
74: VecGetArray(lf,&xx);
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: MatSetValuesLocal(A,1,&i,3,col,v,INSERT_VALUES);
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: MatSetValuesBlockedLocal(A,1,row,3,col,v,INSERT_VALUES);
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: MatSetValuesLocal(A,1,row,5,col,v,INSERT_VALUES);
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: MatSetValuesLocal(A,1,row+1,5,col,v+6,INSERT_VALUES);
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,1,"not implemented for block size %D",bs);
109: }
110: }
111: VecRestoreArray(lf,&xx);
112: VecGhostRestoreLocalForm(x,&lf);
113: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
116: MatMult(A,x,b);
117: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
118: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
119: VecView(b,PETSC_VIEWER_STDOUT_WORLD);
121: /* Update the ghosted values, view the result on rank 0. */
122: VecGhostUpdateBegin(b,INSERT_VALUES,SCATTER_FORWARD);
123: VecGhostUpdateEnd(b,INSERT_VALUES,SCATTER_FORWARD);
124: if (!rank) {
125: VecGhostGetLocalForm(b,&lf);
126: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"Local form of b on rank 0, last two nodes are ghost nodes\n");
127: VecView(lf,PETSC_VIEWER_STDOUT_SELF);
128: VecGhostRestoreLocalForm(b,&lf);
129: }
131: DMDestroy(&slice);
132: VecDestroy(&x);
133: VecDestroy(&b);
134: MatDestroy(&A);
135: PetscFinalize();
136: return ierr;
137: }
140: /*TEST
142: test:
143: nsize: 2
144: args: -bs 2 -block 0 -sliced_mat_type baij -alpha 10 -u0 0.1
146: test:
147: suffix: 2
148: nsize: 2
149: args: -bs 2 -block 1 -sliced_mat_type aij -alpha 10 -u0 0.1
151: test:
152: suffix: 3
153: nsize: 2
154: args: -bs 2 -block 0 -sliced_mat_type aij -alpha 10 -u0 0.1
156: TEST*/