1
2 /*
3 Code for interpolating between grids represented by DMDAs
4 */
5
6 /*
7 For linear elements there are two branches of code to compute the interpolation. They should compute the same results but may not. The "new version" does
8 not work for periodic domains, the old does. Change NEWVERSION to 1 to compile in the new version. Eventually when we are sure the two produce identical results
9 we will remove/merge the new version. Based on current tests, these both produce the same results. We are leaving NEWVERSION for now in the code since some
10 consider it cleaner, but old version is turned on since it handles periodic case.
11 */
12 #define NEWVERSION 0
13
14 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
15
16 /*
17 Since the interpolation uses MATMAIJ for dof > 0 we convert request for non-MATAIJ baseded matrices to MATAIJ.
18 This is a bit of a hack, the reason for it is partially because -dm_mat_type defines the
19 matrix type for both the operator matrices and the interpolation matrices so that users
20 can select matrix types of base MATAIJ for accelerators
21 */
ConvertToAIJ(MatType intype,MatType * outtype)22 static PetscErrorCode ConvertToAIJ(MatType intype,MatType *outtype)
23 {
24 PetscErrorCode ierr;
25 PetscInt i;
26 char const *types[3] = {MATAIJ,MATSEQAIJ,MATMPIAIJ};
27 PetscBool flg;
28
29 PetscFunctionBegin;
30 *outtype = MATAIJ;
31 for (i=0; i<3; i++) {
32 ierr = PetscStrbeginswith(intype,types[i],&flg);CHKERRQ(ierr);
33 if (flg) {
34 *outtype = intype;
35 PetscFunctionReturn(0);
36 }
37 }
38 PetscFunctionReturn(0);
39 }
40
DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat * A)41 PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat *A)
42 {
43 PetscErrorCode ierr;
44 PetscInt i,i_start,m_f,Mx;
45 const PetscInt *idx_f,*idx_c;
46 PetscInt m_ghost,m_ghost_c;
47 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
48 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
49 PetscScalar v[2],x;
50 Mat mat;
51 DMBoundaryType bx;
52 ISLocalToGlobalMapping ltog_f,ltog_c;
53 MatType mattype;
54
55 PetscFunctionBegin;
56 ierr = DMDAGetInfo(dac,NULL,&Mx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&bx,NULL,NULL,NULL);CHKERRQ(ierr);
57 ierr = DMDAGetInfo(daf,NULL,&mx,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
58 if (bx == DM_BOUNDARY_PERIODIC) {
59 ratio = mx/Mx;
60 if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
61 } else {
62 ratio = (mx-1)/(Mx-1);
63 if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
64 }
65
66 ierr = DMDAGetCorners(daf,&i_start,NULL,NULL,&m_f,NULL,NULL);CHKERRQ(ierr);
67 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,NULL,NULL,&m_ghost,NULL,NULL);CHKERRQ(ierr);
68 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
69 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
70
71 ierr = DMDAGetCorners(dac,&i_start_c,NULL,NULL,&m_c,NULL,NULL);CHKERRQ(ierr);
72 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,NULL,NULL,&m_ghost_c,NULL,NULL);CHKERRQ(ierr);
73 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
74 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
75
76 /* create interpolation matrix */
77 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr);
78 #if defined(PETSC_HAVE_CUDA)
79 /*
80 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
81 we don't want the original unconverted matrix copied to the GPU
82 */
83 if (dof > 1) {
84 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
85 }
86 #endif
87 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr);
88 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
89 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
90 ierr = MatSeqAIJSetPreallocation(mat,2,NULL);CHKERRQ(ierr);
91 ierr = MatMPIAIJSetPreallocation(mat,2,NULL,1,NULL);CHKERRQ(ierr);
92
93 /* loop over local fine grid nodes setting interpolation for those*/
94 if (!NEWVERSION) {
95
96 for (i=i_start; i<i_start+m_f; i++) {
97 /* convert to local "natural" numbering and then to PETSc global numbering */
98 row = idx_f[i-i_start_ghost];
99
100 i_c = (i/ratio); /* coarse grid node to left of fine grid node */
101 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
102 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
103
104 /*
105 Only include those interpolation points that are truly
106 nonzero. Note this is very important for final grid lines
107 in x direction; since they have no right neighbor
108 */
109 x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
110 nc = 0;
111 /* one left and below; or we are right on it */
112 col = (i_c-i_start_ghost_c);
113 cols[nc] = idx_c[col];
114 v[nc++] = -x + 1.0;
115 /* one right? */
116 if (i_c*ratio != i) {
117 cols[nc] = idx_c[col+1];
118 v[nc++] = x;
119 }
120 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
121 }
122
123 } else {
124 PetscScalar *xi;
125 PetscInt li,nxi,n;
126 PetscScalar Ni[2];
127
128 /* compute local coordinate arrays */
129 nxi = ratio + 1;
130 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr);
131 for (li=0; li<nxi; li++) {
132 xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
133 }
134
135 for (i=i_start; i<i_start+m_f; i++) {
136 /* convert to local "natural" numbering and then to PETSc global numbering */
137 row = idx_f[(i-i_start_ghost)];
138
139 i_c = (i/ratio); /* coarse grid node to left of fine grid node */
140 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
141 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
142
143 /* remainders */
144 li = i - ratio * (i/ratio);
145 if (i==mx-1) li = nxi-1;
146
147 /* corners */
148 col = (i_c-i_start_ghost_c);
149 cols[0] = idx_c[col];
150 Ni[0] = 1.0;
151 if ((li==0) || (li==nxi-1)) {
152 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr);
153 continue;
154 }
155
156 /* edges + interior */
157 /* remainders */
158 if (i==mx-1) i_c--;
159
160 col = (i_c-i_start_ghost_c);
161 cols[0] = idx_c[col]; /* one left and below; or we are right on it */
162 cols[1] = idx_c[col+1];
163
164 Ni[0] = 0.5*(1.0-xi[li]);
165 Ni[1] = 0.5*(1.0+xi[li]);
166 for (n=0; n<2; n++) {
167 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
168 }
169 ierr = MatSetValues(mat,1,&row,2,cols,Ni,INSERT_VALUES);CHKERRQ(ierr);
170 }
171 ierr = PetscFree(xi);CHKERRQ(ierr);
172 }
173 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
174 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
175 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
176 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
177 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
178 ierr = MatDestroy(&mat);CHKERRQ(ierr);
179 PetscFunctionReturn(0);
180 }
181
DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat * A)182 PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat *A)
183 {
184 PetscErrorCode ierr;
185 PetscInt i,i_start,m_f,Mx;
186 const PetscInt *idx_f,*idx_c;
187 ISLocalToGlobalMapping ltog_f,ltog_c;
188 PetscInt m_ghost,m_ghost_c;
189 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
190 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
191 PetscScalar v[2],x;
192 Mat mat;
193 DMBoundaryType bx;
194 MatType mattype;
195
196 PetscFunctionBegin;
197 ierr = DMDAGetInfo(dac,NULL,&Mx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&bx,NULL,NULL,NULL);CHKERRQ(ierr);
198 ierr = DMDAGetInfo(daf,NULL,&mx,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
199 if (bx == DM_BOUNDARY_PERIODIC) {
200 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
201 ratio = mx/Mx;
202 if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
203 } else {
204 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
205 ratio = (mx-1)/(Mx-1);
206 if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
207 }
208
209 ierr = DMDAGetCorners(daf,&i_start,NULL,NULL,&m_f,NULL,NULL);CHKERRQ(ierr);
210 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,NULL,NULL,&m_ghost,NULL,NULL);CHKERRQ(ierr);
211 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
212 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
213
214 ierr = DMDAGetCorners(dac,&i_start_c,NULL,NULL,&m_c,NULL,NULL);CHKERRQ(ierr);
215 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,NULL,NULL,&m_ghost_c,NULL,NULL);CHKERRQ(ierr);
216 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
217 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
218
219 /* create interpolation matrix */
220 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr);
221 #if defined(PETSC_HAVE_CUDA)
222 /*
223 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
224 we don't want the original unconverted matrix copied to the GPU
225 */
226 if (dof > 1) {
227 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
228 }
229 #endif
230 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr);
231 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
232 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
233 ierr = MatSeqAIJSetPreallocation(mat,2,NULL);CHKERRQ(ierr);
234 ierr = MatMPIAIJSetPreallocation(mat,2,NULL,0,NULL);CHKERRQ(ierr);
235
236 /* loop over local fine grid nodes setting interpolation for those*/
237 for (i=i_start; i<i_start+m_f; i++) {
238 /* convert to local "natural" numbering and then to PETSc global numbering */
239 row = idx_f[(i-i_start_ghost)];
240
241 i_c = (i/ratio); /* coarse grid node to left of fine grid node */
242
243 /*
244 Only include those interpolation points that are truly
245 nonzero. Note this is very important for final grid lines
246 in x direction; since they have no right neighbor
247 */
248 x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
249 nc = 0;
250 /* one left and below; or we are right on it */
251 col = (i_c-i_start_ghost_c);
252 cols[nc] = idx_c[col];
253 v[nc++] = -x + 1.0;
254 /* one right? */
255 if (i_c*ratio != i) {
256 cols[nc] = idx_c[col+1];
257 v[nc++] = x;
258 }
259 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
260 }
261 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
262 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
263 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
264 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
265 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
266 ierr = MatDestroy(&mat);CHKERRQ(ierr);
267 ierr = PetscLogFlops(5.0*m_f);CHKERRQ(ierr);
268 PetscFunctionReturn(0);
269 }
270
DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat * A)271 PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat *A)
272 {
273 PetscErrorCode ierr;
274 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
275 const PetscInt *idx_c,*idx_f;
276 ISLocalToGlobalMapping ltog_f,ltog_c;
277 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
278 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
279 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
280 PetscMPIInt size_c,size_f,rank_f;
281 PetscScalar v[4],x,y;
282 Mat mat;
283 DMBoundaryType bx,by;
284 MatType mattype;
285
286 PetscFunctionBegin;
287 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,NULL,NULL,NULL,NULL,NULL,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr);
288 ierr = DMDAGetInfo(daf,NULL,&mx,&my,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
289 if (bx == DM_BOUNDARY_PERIODIC) {
290 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
291 ratioi = mx/Mx;
292 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
293 } else {
294 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
295 ratioi = (mx-1)/(Mx-1);
296 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
297 }
298 if (by == DM_BOUNDARY_PERIODIC) {
299 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
300 ratioj = my/My;
301 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
302 } else {
303 if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
304 ratioj = (my-1)/(My-1);
305 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
306 }
307
308 ierr = DMDAGetCorners(daf,&i_start,&j_start,NULL,&m_f,&n_f,NULL);CHKERRQ(ierr);
309 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,NULL,&m_ghost,&n_ghost,NULL);CHKERRQ(ierr);
310 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
311 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
312
313 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,NULL,&m_c,&n_c,NULL);CHKERRQ(ierr);
314 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,NULL,&m_ghost_c,&n_ghost_c,NULL);CHKERRQ(ierr);
315 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
316 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
317
318 /*
319 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
320 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
321 processors). It's effective length is hence 4 times its normal length, this is
322 why the col_scale is multiplied by the interpolation matrix column sizes.
323 sol_shift allows each set of 1/4 processors do its own interpolation using ITS
324 copy of the coarse vector. A bit of a hack but you do better.
325
326 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
327 */
328 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr);
329 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr);
330 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr);
331 col_scale = size_f/size_c;
332 col_shift = Mx*My*(rank_f/size_c);
333
334 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr);
335 for (j=j_start; j<j_start+n_f; j++) {
336 for (i=i_start; i<i_start+m_f; i++) {
337 /* convert to local "natural" numbering and then to PETSc global numbering */
338 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
339
340 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
341 j_c = (j/ratioj); /* coarse grid node below fine grid node */
342
343 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
344 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
345 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
346 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
347
348 /*
349 Only include those interpolation points that are truly
350 nonzero. Note this is very important for final grid lines
351 in x and y directions; since they have no right/top neighbors
352 */
353 nc = 0;
354 /* one left and below; or we are right on it */
355 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
356 cols[nc++] = col_shift + idx_c[col];
357 /* one right and below */
358 if (i_c*ratioi != i) cols[nc++] = col_shift + idx_c[col+1];
359 /* one left and above */
360 if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c];
361 /* one right and above */
362 if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)];
363 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr);
364 }
365 }
366 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr);
367 #if defined(PETSC_HAVE_CUDA)
368 /*
369 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
370 we don't want the original unconverted matrix copied to the GPU
371 */
372 if (dof > 1) {
373 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
374 }
375 #endif
376 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr);
377 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
378 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
379 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr);
380 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr);
381 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
382
383 /* loop over local fine grid nodes setting interpolation for those*/
384 if (!NEWVERSION) {
385
386 for (j=j_start; j<j_start+n_f; j++) {
387 for (i=i_start; i<i_start+m_f; i++) {
388 /* convert to local "natural" numbering and then to PETSc global numbering */
389 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
390
391 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
392 j_c = (j/ratioj); /* coarse grid node below fine grid node */
393
394 /*
395 Only include those interpolation points that are truly
396 nonzero. Note this is very important for final grid lines
397 in x and y directions; since they have no right/top neighbors
398 */
399 x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
400 y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);
401
402 nc = 0;
403 /* one left and below; or we are right on it */
404 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
405 cols[nc] = col_shift + idx_c[col];
406 v[nc++] = x*y - x - y + 1.0;
407 /* one right and below */
408 if (i_c*ratioi != i) {
409 cols[nc] = col_shift + idx_c[col+1];
410 v[nc++] = -x*y + x;
411 }
412 /* one left and above */
413 if (j_c*ratioj != j) {
414 cols[nc] = col_shift + idx_c[col+m_ghost_c];
415 v[nc++] = -x*y + y;
416 }
417 /* one right and above */
418 if (j_c*ratioj != j && i_c*ratioi != i) {
419 cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)];
420 v[nc++] = x*y;
421 }
422 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
423 }
424 }
425
426 } else {
427 PetscScalar Ni[4];
428 PetscScalar *xi,*eta;
429 PetscInt li,nxi,lj,neta;
430
431 /* compute local coordinate arrays */
432 nxi = ratioi + 1;
433 neta = ratioj + 1;
434 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr);
435 ierr = PetscMalloc1(neta,&eta);CHKERRQ(ierr);
436 for (li=0; li<nxi; li++) {
437 xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
438 }
439 for (lj=0; lj<neta; lj++) {
440 eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
441 }
442
443 /* loop over local fine grid nodes setting interpolation for those*/
444 for (j=j_start; j<j_start+n_f; j++) {
445 for (i=i_start; i<i_start+m_f; i++) {
446 /* convert to local "natural" numbering and then to PETSc global numbering */
447 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
448
449 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
450 j_c = (j/ratioj); /* coarse grid node below fine grid node */
451
452 /* remainders */
453 li = i - ratioi * (i/ratioi);
454 if (i==mx-1) li = nxi-1;
455 lj = j - ratioj * (j/ratioj);
456 if (j==my-1) lj = neta-1;
457
458 /* corners */
459 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
460 cols[0] = col_shift + idx_c[col]; /* left, below */
461 Ni[0] = 1.0;
462 if ((li==0) || (li==nxi-1)) {
463 if ((lj==0) || (lj==neta-1)) {
464 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr);
465 continue;
466 }
467 }
468
469 /* edges + interior */
470 /* remainders */
471 if (i==mx-1) i_c--;
472 if (j==my-1) j_c--;
473
474 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
475 cols[0] = col_shift + idx_c[col]; /* left, below */
476 cols[1] = col_shift + idx_c[col+1]; /* right, below */
477 cols[2] = col_shift + idx_c[col+m_ghost_c]; /* left, above */
478 cols[3] = col_shift + idx_c[col+(m_ghost_c+1)]; /* right, above */
479
480 Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]);
481 Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]);
482 Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]);
483 Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]);
484
485 nc = 0;
486 if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1;
487 if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1;
488 if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1;
489 if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1;
490
491 ierr = MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);CHKERRQ(ierr);
492 }
493 }
494 ierr = PetscFree(xi);CHKERRQ(ierr);
495 ierr = PetscFree(eta);CHKERRQ(ierr);
496 }
497 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
498 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
499 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
500 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
501 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
502 ierr = MatDestroy(&mat);CHKERRQ(ierr);
503 PetscFunctionReturn(0);
504 }
505
506 /*
507 Contributed by Andrei Draganescu <aidraga@sandia.gov>
508 */
DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat * A)509 PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A)
510 {
511 PetscErrorCode ierr;
512 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
513 const PetscInt *idx_c,*idx_f;
514 ISLocalToGlobalMapping ltog_f,ltog_c;
515 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
516 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
517 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
518 PetscMPIInt size_c,size_f,rank_f;
519 PetscScalar v[4];
520 Mat mat;
521 DMBoundaryType bx,by;
522 MatType mattype;
523
524 PetscFunctionBegin;
525 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,NULL,NULL,NULL,NULL,NULL,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr);
526 ierr = DMDAGetInfo(daf,NULL,&mx,&my,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
527 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
528 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
529 ratioi = mx/Mx;
530 ratioj = my/My;
531 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
532 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
533 if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2");
534 if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2");
535
536 ierr = DMDAGetCorners(daf,&i_start,&j_start,NULL,&m_f,&n_f,NULL);CHKERRQ(ierr);
537 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,NULL,&m_ghost,&n_ghost,NULL);CHKERRQ(ierr);
538 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
539 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
540
541 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,NULL,&m_c,&n_c,NULL);CHKERRQ(ierr);
542 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,NULL,&m_ghost_c,&n_ghost_c,NULL);CHKERRQ(ierr);
543 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
544 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
545
546 /*
547 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
548 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
549 processors). It's effective length is hence 4 times its normal length, this is
550 why the col_scale is multiplied by the interpolation matrix column sizes.
551 sol_shift allows each set of 1/4 processors do its own interpolation using ITS
552 copy of the coarse vector. A bit of a hack but you do better.
553
554 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
555 */
556 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr);
557 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr);
558 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr);
559 col_scale = size_f/size_c;
560 col_shift = Mx*My*(rank_f/size_c);
561
562 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr);
563 for (j=j_start; j<j_start+n_f; j++) {
564 for (i=i_start; i<i_start+m_f; i++) {
565 /* convert to local "natural" numbering and then to PETSc global numbering */
566 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
567
568 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
569 j_c = (j/ratioj); /* coarse grid node below fine grid node */
570
571 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
572 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
573 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
574 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
575
576 /*
577 Only include those interpolation points that are truly
578 nonzero. Note this is very important for final grid lines
579 in x and y directions; since they have no right/top neighbors
580 */
581 nc = 0;
582 /* one left and below; or we are right on it */
583 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
584 cols[nc++] = col_shift + idx_c[col];
585 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr);
586 }
587 }
588 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr);
589 #if defined(PETSC_HAVE_CUDA)
590 /*
591 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
592 we don't want the original unconverted matrix copied to the GPU
593 */
594 if (dof > 1) {
595 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
596 }
597 #endif
598 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr);
599 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
600 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
601 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr);
602 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr);
603 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
604
605 /* loop over local fine grid nodes setting interpolation for those*/
606 for (j=j_start; j<j_start+n_f; j++) {
607 for (i=i_start; i<i_start+m_f; i++) {
608 /* convert to local "natural" numbering and then to PETSc global numbering */
609 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
610
611 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
612 j_c = (j/ratioj); /* coarse grid node below fine grid node */
613 nc = 0;
614 /* one left and below; or we are right on it */
615 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
616 cols[nc] = col_shift + idx_c[col];
617 v[nc++] = 1.0;
618
619 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
620 }
621 }
622 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
623 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
624 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
625 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
626 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
627 ierr = MatDestroy(&mat);CHKERRQ(ierr);
628 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr);
629 PetscFunctionReturn(0);
630 }
631
632 /*
633 Contributed by Jianming Yang <jianming-yang@uiowa.edu>
634 */
DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat * A)635 PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A)
636 {
637 PetscErrorCode ierr;
638 PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,dof;
639 const PetscInt *idx_c,*idx_f;
640 ISLocalToGlobalMapping ltog_f,ltog_c;
641 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz;
642 PetscInt row,col,i_start_ghost,j_start_ghost,l_start_ghost,cols[8],mx,m_c,my,n_c,mz,p_c,ratioi,ratioj,ratiol;
643 PetscInt i_c,j_c,l_c,i_start_c,j_start_c,l_start_c,i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,col_shift,col_scale;
644 PetscMPIInt size_c,size_f,rank_f;
645 PetscScalar v[8];
646 Mat mat;
647 DMBoundaryType bx,by,bz;
648 MatType mattype;
649
650 PetscFunctionBegin;
651 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,&Mz,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);CHKERRQ(ierr);
652 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
653 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
654 if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
655 ierr = DMDAGetInfo(daf,NULL,&mx,&my,&mz,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
656 ratioi = mx/Mx;
657 ratioj = my/My;
658 ratiol = mz/Mz;
659 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
660 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
661 if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z");
662 if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2");
663 if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2");
664 if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2");
665
666 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr);
667 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr);
668 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
669 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
670
671 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr);
672 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr);
673 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
674 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
675
676 /*
677 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
678 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
679 processors). It's effective length is hence 4 times its normal length, this is
680 why the col_scale is multiplied by the interpolation matrix column sizes.
681 sol_shift allows each set of 1/4 processors do its own interpolation using ITS
682 copy of the coarse vector. A bit of a hack but you do better.
683
684 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
685 */
686 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr);
687 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr);
688 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr);
689 col_scale = size_f/size_c;
690 col_shift = Mx*My*Mz*(rank_f/size_c);
691
692 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr);
693 for (l=l_start; l<l_start+p_f; l++) {
694 for (j=j_start; j<j_start+n_f; j++) {
695 for (i=i_start; i<i_start+m_f; i++) {
696 /* convert to local "natural" numbering and then to PETSc global numbering */
697 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
698
699 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
700 j_c = (j/ratioj); /* coarse grid node below fine grid node */
701 l_c = (l/ratiol);
702
703 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
704 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
705 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
706 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
707 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
708 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
709
710 /*
711 Only include those interpolation points that are truly
712 nonzero. Note this is very important for final grid lines
713 in x and y directions; since they have no right/top neighbors
714 */
715 nc = 0;
716 /* one left and below; or we are right on it */
717 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
718 cols[nc++] = col_shift + idx_c[col];
719 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr);
720 }
721 }
722 }
723 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr);
724 #if defined(PETSC_HAVE_CUDA)
725 /*
726 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
727 we don't want the original unconverted matrix copied to the GPU
728 */
729 if (dof > 1) {
730 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
731 }
732 #endif
733 ierr = MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);CHKERRQ(ierr);
734 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
735 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
736 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr);
737 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr);
738 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
739
740 /* loop over local fine grid nodes setting interpolation for those*/
741 for (l=l_start; l<l_start+p_f; l++) {
742 for (j=j_start; j<j_start+n_f; j++) {
743 for (i=i_start; i<i_start+m_f; i++) {
744 /* convert to local "natural" numbering and then to PETSc global numbering */
745 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
746
747 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
748 j_c = (j/ratioj); /* coarse grid node below fine grid node */
749 l_c = (l/ratiol);
750 nc = 0;
751 /* one left and below; or we are right on it */
752 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
753 cols[nc] = col_shift + idx_c[col];
754 v[nc++] = 1.0;
755
756 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
757 }
758 }
759 }
760 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
761 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
762 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
763 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
764 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
765 ierr = MatDestroy(&mat);CHKERRQ(ierr);
766 ierr = PetscLogFlops(13.0*m_f*n_f*p_f);CHKERRQ(ierr);
767 PetscFunctionReturn(0);
768 }
769
DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat * A)770 PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A)
771 {
772 PetscErrorCode ierr;
773 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof,l;
774 const PetscInt *idx_c,*idx_f;
775 ISLocalToGlobalMapping ltog_f,ltog_c;
776 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,Mz,mz;
777 PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
778 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
779 PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
780 PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
781 PetscScalar v[8],x,y,z;
782 Mat mat;
783 DMBoundaryType bx,by,bz;
784 MatType mattype;
785
786 PetscFunctionBegin;
787 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,&Mz,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);CHKERRQ(ierr);
788 ierr = DMDAGetInfo(daf,NULL,&mx,&my,&mz,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
789 if (mx == Mx) {
790 ratioi = 1;
791 } else if (bx == DM_BOUNDARY_PERIODIC) {
792 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
793 ratioi = mx/Mx;
794 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
795 } else {
796 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
797 ratioi = (mx-1)/(Mx-1);
798 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
799 }
800 if (my == My) {
801 ratioj = 1;
802 } else if (by == DM_BOUNDARY_PERIODIC) {
803 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
804 ratioj = my/My;
805 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
806 } else {
807 if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
808 ratioj = (my-1)/(My-1);
809 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
810 }
811 if (mz == Mz) {
812 ratiok = 1;
813 } else if (bz == DM_BOUNDARY_PERIODIC) {
814 if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
815 ratiok = mz/Mz;
816 if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D Mz %D",mz,Mz);
817 } else {
818 if (Mz < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be greater than 1",Mz);
819 ratiok = (mz-1)/(Mz-1);
820 if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
821 }
822
823 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr);
824 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr);
825 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
826 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
827
828 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr);
829 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr);
830 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
831 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
832
833 /* create interpolation matrix, determining exact preallocation */
834 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr);
835 /* loop over local fine grid nodes counting interpolating points */
836 for (l=l_start; l<l_start+p_f; l++) {
837 for (j=j_start; j<j_start+n_f; j++) {
838 for (i=i_start; i<i_start+m_f; i++) {
839 /* convert to local "natural" numbering and then to PETSc global numbering */
840 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
841 i_c = (i/ratioi);
842 j_c = (j/ratioj);
843 l_c = (l/ratiok);
844 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
845 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
846 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
847 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
848 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
849 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
850
851 /*
852 Only include those interpolation points that are truly
853 nonzero. Note this is very important for final grid lines
854 in x and y directions; since they have no right/top neighbors
855 */
856 nc = 0;
857 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
858 cols[nc++] = idx_c[col];
859 if (i_c*ratioi != i) {
860 cols[nc++] = idx_c[col+1];
861 }
862 if (j_c*ratioj != j) {
863 cols[nc++] = idx_c[col+m_ghost_c];
864 }
865 if (l_c*ratiok != l) {
866 cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c];
867 }
868 if (j_c*ratioj != j && i_c*ratioi != i) {
869 cols[nc++] = idx_c[col+(m_ghost_c+1)];
870 }
871 if (j_c*ratioj != j && l_c*ratiok != l) {
872 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
873 }
874 if (i_c*ratioi != i && l_c*ratiok != l) {
875 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
876 }
877 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
878 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
879 }
880 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr);
881 }
882 }
883 }
884 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr);
885 #if defined(PETSC_HAVE_CUDA)
886 /*
887 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
888 we don't want the original unconverted matrix copied to the GPU
889 */
890 if (dof > 1) {
891 ierr = MatBindToCPU(mat,PETSC_TRUE);CHKERRQ(ierr);
892 }
893 #endif
894 ierr = MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);CHKERRQ(ierr);
895 ierr = ConvertToAIJ(dac->mattype,&mattype);CHKERRQ(ierr);
896 ierr = MatSetType(mat,mattype);CHKERRQ(ierr);
897 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr);
898 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr);
899 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
900
901 /* loop over local fine grid nodes setting interpolation for those*/
902 if (!NEWVERSION) {
903
904 for (l=l_start; l<l_start+p_f; l++) {
905 for (j=j_start; j<j_start+n_f; j++) {
906 for (i=i_start; i<i_start+m_f; i++) {
907 /* convert to local "natural" numbering and then to PETSc global numbering */
908 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
909
910 i_c = (i/ratioi);
911 j_c = (j/ratioj);
912 l_c = (l/ratiok);
913
914 /*
915 Only include those interpolation points that are truly
916 nonzero. Note this is very important for final grid lines
917 in x and y directions; since they have no right/top neighbors
918 */
919 x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
920 y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);
921 z = ((PetscReal)(l - l_c*ratiok))/((PetscReal)ratiok);
922
923 nc = 0;
924 /* one left and below; or we are right on it */
925 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
926
927 cols[nc] = idx_c[col];
928 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
929
930 if (i_c*ratioi != i) {
931 cols[nc] = idx_c[col+1];
932 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
933 }
934
935 if (j_c*ratioj != j) {
936 cols[nc] = idx_c[col+m_ghost_c];
937 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
938 }
939
940 if (l_c*ratiok != l) {
941 cols[nc] = idx_c[col+m_ghost_c*n_ghost_c];
942 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
943 }
944
945 if (j_c*ratioj != j && i_c*ratioi != i) {
946 cols[nc] = idx_c[col+(m_ghost_c+1)];
947 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
948 }
949
950 if (j_c*ratioj != j && l_c*ratiok != l) {
951 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
952 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
953 }
954
955 if (i_c*ratioi != i && l_c*ratiok != l) {
956 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
957 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
958 }
959
960 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
961 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
962 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
963 }
964 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
965 }
966 }
967 }
968
969 } else {
970 PetscScalar *xi,*eta,*zeta;
971 PetscInt li,nxi,lj,neta,lk,nzeta,n;
972 PetscScalar Ni[8];
973
974 /* compute local coordinate arrays */
975 nxi = ratioi + 1;
976 neta = ratioj + 1;
977 nzeta = ratiok + 1;
978 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr);
979 ierr = PetscMalloc1(neta,&eta);CHKERRQ(ierr);
980 ierr = PetscMalloc1(nzeta,&zeta);CHKERRQ(ierr);
981 for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
982 for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
983 for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
984
985 for (l=l_start; l<l_start+p_f; l++) {
986 for (j=j_start; j<j_start+n_f; j++) {
987 for (i=i_start; i<i_start+m_f; i++) {
988 /* convert to local "natural" numbering and then to PETSc global numbering */
989 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
990
991 i_c = (i/ratioi);
992 j_c = (j/ratioj);
993 l_c = (l/ratiok);
994
995 /* remainders */
996 li = i - ratioi * (i/ratioi);
997 if (i==mx-1) li = nxi-1;
998 lj = j - ratioj * (j/ratioj);
999 if (j==my-1) lj = neta-1;
1000 lk = l - ratiok * (l/ratiok);
1001 if (l==mz-1) lk = nzeta-1;
1002
1003 /* corners */
1004 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
1005 cols[0] = idx_c[col];
1006 Ni[0] = 1.0;
1007 if ((li==0) || (li==nxi-1)) {
1008 if ((lj==0) || (lj==neta-1)) {
1009 if ((lk==0) || (lk==nzeta-1)) {
1010 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr);
1011 continue;
1012 }
1013 }
1014 }
1015
1016 /* edges + interior */
1017 /* remainders */
1018 if (i==mx-1) i_c--;
1019 if (j==my-1) j_c--;
1020 if (l==mz-1) l_c--;
1021
1022 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
1023 cols[0] = idx_c[col]; /* one left and below; or we are right on it */
1024 cols[1] = idx_c[col+1]; /* one right and below */
1025 cols[2] = idx_c[col+m_ghost_c]; /* one left and above */
1026 cols[3] = idx_c[col+(m_ghost_c+1)]; /* one right and above */
1027
1028 cols[4] = idx_c[col+m_ghost_c*n_ghost_c]; /* one left and below and front; or we are right on it */
1029 cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; /* one right and below, and front */
1030 cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; /* one left and above and front*/
1031 cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; /* one right and above and front */
1032
1033 Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
1034 Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
1035 Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
1036 Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
1037
1038 Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
1039 Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
1040 Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
1041 Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
1042
1043 for (n=0; n<8; n++) {
1044 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
1045 }
1046 ierr = MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr);
1047
1048 }
1049 }
1050 }
1051 ierr = PetscFree(xi);CHKERRQ(ierr);
1052 ierr = PetscFree(eta);CHKERRQ(ierr);
1053 ierr = PetscFree(zeta);CHKERRQ(ierr);
1054 }
1055 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1056 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
1057 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1058 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1059
1060 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr);
1061 ierr = MatDestroy(&mat);CHKERRQ(ierr);
1062 PetscFunctionReturn(0);
1063 }
1064
DMCreateInterpolation_DA(DM dac,DM daf,Mat * A,Vec * scale)1065 PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale)
1066 {
1067 PetscErrorCode ierr;
1068 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1069 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1070 DMDAStencilType stc,stf;
1071 DM_DA *ddc = (DM_DA*)dac->data;
1072
1073 PetscFunctionBegin;
1074 PetscValidHeaderSpecific(dac,DM_CLASSID,1);
1075 PetscValidHeaderSpecific(daf,DM_CLASSID,2);
1076 PetscValidPointer(A,3);
1077 if (scale) PetscValidPointer(scale,4);
1078
1079 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr);
1080 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr);
1081 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1082 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1083 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1084 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1085 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1086 if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1087 if (dimc > 1 && Nc < 2 && Nf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1088 if (dimc > 2 && Pc < 2 && Pf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1089
1090 if (ddc->interptype == DMDA_Q1) {
1091 if (dimc == 1) {
1092 ierr = DMCreateInterpolation_DA_1D_Q1(dac,daf,A);CHKERRQ(ierr);
1093 } else if (dimc == 2) {
1094 ierr = DMCreateInterpolation_DA_2D_Q1(dac,daf,A);CHKERRQ(ierr);
1095 } else if (dimc == 3) {
1096 ierr = DMCreateInterpolation_DA_3D_Q1(dac,daf,A);CHKERRQ(ierr);
1097 } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1098 } else if (ddc->interptype == DMDA_Q0) {
1099 if (dimc == 1) {
1100 ierr = DMCreateInterpolation_DA_1D_Q0(dac,daf,A);CHKERRQ(ierr);
1101 } else if (dimc == 2) {
1102 ierr = DMCreateInterpolation_DA_2D_Q0(dac,daf,A);CHKERRQ(ierr);
1103 } else if (dimc == 3) {
1104 ierr = DMCreateInterpolation_DA_3D_Q0(dac,daf,A);CHKERRQ(ierr);
1105 } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1106 }
1107 if (scale) {
1108 ierr = DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr);
1109 }
1110 PetscFunctionReturn(0);
1111 }
1112
DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter * inject)1113 PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject)
1114 {
1115 PetscErrorCode ierr;
1116 PetscInt i,i_start,m_f,Mx,dof;
1117 const PetscInt *idx_f;
1118 ISLocalToGlobalMapping ltog_f;
1119 PetscInt m_ghost,m_ghost_c;
1120 PetscInt row,i_start_ghost,mx,m_c,nc,ratioi;
1121 PetscInt i_start_c,i_start_ghost_c;
1122 PetscInt *cols;
1123 DMBoundaryType bx;
1124 Vec vecf,vecc;
1125 IS isf;
1126
1127 PetscFunctionBegin;
1128 ierr = DMDAGetInfo(dac,NULL,&Mx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&bx,NULL,NULL,NULL);CHKERRQ(ierr);
1129 ierr = DMDAGetInfo(daf,NULL,&mx,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1130 if (bx == DM_BOUNDARY_PERIODIC) {
1131 ratioi = mx/Mx;
1132 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1133 } else {
1134 ratioi = (mx-1)/(Mx-1);
1135 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1136 }
1137
1138 ierr = DMDAGetCorners(daf,&i_start,NULL,NULL,&m_f,NULL,NULL);CHKERRQ(ierr);
1139 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,NULL,NULL,&m_ghost,NULL,NULL);CHKERRQ(ierr);
1140 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
1141 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1142
1143 ierr = DMDAGetCorners(dac,&i_start_c,NULL,NULL,&m_c,NULL,NULL);CHKERRQ(ierr);
1144 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,NULL,NULL,&m_ghost_c,NULL,NULL);CHKERRQ(ierr);
1145
1146
1147 /* loop over local fine grid nodes setting interpolation for those*/
1148 nc = 0;
1149 ierr = PetscMalloc1(m_f,&cols);CHKERRQ(ierr);
1150
1151
1152 for (i=i_start_c; i<i_start_c+m_c; i++) {
1153 PetscInt i_f = i*ratioi;
1154
1155 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\ni_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1156
1157 row = idx_f[(i_f-i_start_ghost)];
1158 cols[nc++] = row;
1159 }
1160
1161 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1162 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr);
1163 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr);
1164 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr);
1165 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr);
1166 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr);
1167 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr);
1168 ierr = ISDestroy(&isf);CHKERRQ(ierr);
1169 PetscFunctionReturn(0);
1170 }
1171
DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter * inject)1172 PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject)
1173 {
1174 PetscErrorCode ierr;
1175 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
1176 const PetscInt *idx_c,*idx_f;
1177 ISLocalToGlobalMapping ltog_f,ltog_c;
1178 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c;
1179 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
1180 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
1181 PetscInt *cols;
1182 DMBoundaryType bx,by;
1183 Vec vecf,vecc;
1184 IS isf;
1185
1186 PetscFunctionBegin;
1187 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,NULL,NULL,NULL,NULL,NULL,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr);
1188 ierr = DMDAGetInfo(daf,NULL,&mx,&my,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1189 if (bx == DM_BOUNDARY_PERIODIC) {
1190 ratioi = mx/Mx;
1191 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1192 } else {
1193 ratioi = (mx-1)/(Mx-1);
1194 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1195 }
1196 if (by == DM_BOUNDARY_PERIODIC) {
1197 ratioj = my/My;
1198 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
1199 } else {
1200 ratioj = (my-1)/(My-1);
1201 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1202 }
1203
1204 ierr = DMDAGetCorners(daf,&i_start,&j_start,NULL,&m_f,&n_f,NULL);CHKERRQ(ierr);
1205 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,NULL,&m_ghost,&n_ghost,NULL);CHKERRQ(ierr);
1206 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
1207 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1208
1209 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,NULL,&m_c,&n_c,NULL);CHKERRQ(ierr);
1210 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,NULL,&m_ghost_c,&n_ghost_c,NULL);CHKERRQ(ierr);
1211 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
1212 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
1213
1214 /* loop over local fine grid nodes setting interpolation for those*/
1215 nc = 0;
1216 ierr = PetscMalloc1(n_f*m_f,&cols);CHKERRQ(ierr);
1217 for (j=j_start_c; j<j_start_c+n_c; j++) {
1218 for (i=i_start_c; i<i_start_c+m_c; i++) {
1219 PetscInt i_f = i*ratioi,j_f = j*ratioj;
1220 if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1221 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1222 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1223 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1224 row = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1225 cols[nc++] = row;
1226 }
1227 }
1228 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1229 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
1230
1231 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr);
1232 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr);
1233 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr);
1234 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr);
1235 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr);
1236 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr);
1237 ierr = ISDestroy(&isf);CHKERRQ(ierr);
1238 PetscFunctionReturn(0);
1239 }
1240
DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter * inject)1241 PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject)
1242 {
1243 PetscErrorCode ierr;
1244 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz;
1245 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c;
1246 PetscInt i_start_ghost,j_start_ghost,k_start_ghost;
1247 PetscInt mx,my,mz,ratioi,ratioj,ratiok;
1248 PetscInt i_start_c,j_start_c,k_start_c;
1249 PetscInt m_c,n_c,p_c;
1250 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c;
1251 PetscInt row,nc,dof;
1252 const PetscInt *idx_c,*idx_f;
1253 ISLocalToGlobalMapping ltog_f,ltog_c;
1254 PetscInt *cols;
1255 DMBoundaryType bx,by,bz;
1256 Vec vecf,vecc;
1257 IS isf;
1258
1259 PetscFunctionBegin;
1260 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,&Mz,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);CHKERRQ(ierr);
1261 ierr = DMDAGetInfo(daf,NULL,&mx,&my,&mz,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1262
1263 if (bx == DM_BOUNDARY_PERIODIC) {
1264 ratioi = mx/Mx;
1265 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1266 } else {
1267 ratioi = (mx-1)/(Mx-1);
1268 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1269 }
1270 if (by == DM_BOUNDARY_PERIODIC) {
1271 ratioj = my/My;
1272 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
1273 } else {
1274 ratioj = (my-1)/(My-1);
1275 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1276 }
1277 if (bz == DM_BOUNDARY_PERIODIC) {
1278 ratiok = mz/Mz;
1279 if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D My %D",mz,Mz);
1280 } else {
1281 ratiok = (mz-1)/(Mz-1);
1282 if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
1283 }
1284
1285 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr);
1286 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr);
1287 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr);
1288 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1289
1290 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr);
1291 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr);
1292 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr);
1293 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
1294
1295
1296 /* loop over local fine grid nodes setting interpolation for those*/
1297 nc = 0;
1298 ierr = PetscMalloc1(n_f*m_f*p_f,&cols);CHKERRQ(ierr);
1299 for (k=k_start_c; k<k_start_c+p_c; k++) {
1300 for (j=j_start_c; j<j_start_c+n_c; j++) {
1301 for (i=i_start_c; i<i_start_c+m_c; i++) {
1302 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok;
1303 if (k_f < k_start_ghost || k_f >= k_start_ghost+p_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1304 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost);
1305 if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1306 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1307 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1308 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1309 row = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1310 cols[nc++] = row;
1311 }
1312 }
1313 }
1314 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr);
1315 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr);
1316
1317 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr);
1318 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr);
1319 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr);
1320 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr);
1321 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr);
1322 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr);
1323 ierr = ISDestroy(&isf);CHKERRQ(ierr);
1324 PetscFunctionReturn(0);
1325 }
1326
DMCreateInjection_DA(DM dac,DM daf,Mat * mat)1327 PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,Mat *mat)
1328 {
1329 PetscErrorCode ierr;
1330 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1331 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1332 DMDAStencilType stc,stf;
1333 VecScatter inject = NULL;
1334
1335 PetscFunctionBegin;
1336 PetscValidHeaderSpecific(dac,DM_CLASSID,1);
1337 PetscValidHeaderSpecific(daf,DM_CLASSID,2);
1338 PetscValidPointer(mat,3);
1339
1340 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr);
1341 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr);
1342 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1343 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1344 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1345 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1346 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1347 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1348 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1349 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1350
1351 if (dimc == 1) {
1352 ierr = DMCreateInjection_DA_1D(dac,daf,&inject);CHKERRQ(ierr);
1353 } else if (dimc == 2) {
1354 ierr = DMCreateInjection_DA_2D(dac,daf,&inject);CHKERRQ(ierr);
1355 } else if (dimc == 3) {
1356 ierr = DMCreateInjection_DA_3D(dac,daf,&inject);CHKERRQ(ierr);
1357 }
1358 ierr = MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);CHKERRQ(ierr);
1359 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr);
1360 PetscFunctionReturn(0);
1361 }
1362
1363 /*@
1364 DMCreateAggregates - Deprecated, see DMDACreateAggregates.
1365
1366 Level: intermediate
1367 @*/
DMCreateAggregates(DM dac,DM daf,Mat * mat)1368 PetscErrorCode DMCreateAggregates(DM dac,DM daf,Mat *mat)
1369 {
1370 return DMDACreateAggregates(dac,daf,mat);
1371 }
1372
1373 /*@
1374 DMDACreateAggregates - Gets the aggregates that map between
1375 grids associated with two DMDAs.
1376
1377 Collective on dmc
1378
1379 Input Parameters:
1380 + dmc - the coarse grid DMDA
1381 - dmf - the fine grid DMDA
1382
1383 Output Parameters:
1384 . rest - the restriction matrix (transpose of the projection matrix)
1385
1386 Level: intermediate
1387
1388 Note: This routine is not used by PETSc.
1389 It is not clear what its use case is and it may be removed in a future release.
1390 Users should contact petsc-maint@mcs.anl.gov if they plan to use it.
1391
1392 .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
1393 @*/
DMDACreateAggregates(DM dac,DM daf,Mat * rest)1394 PetscErrorCode DMDACreateAggregates(DM dac,DM daf,Mat *rest)
1395 {
1396 PetscErrorCode ierr;
1397 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc;
1398 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1399 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1400 DMDAStencilType stc,stf;
1401 PetscInt i,j,l;
1402 PetscInt i_start,j_start,l_start, m_f,n_f,p_f;
1403 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost;
1404 const PetscInt *idx_f;
1405 PetscInt i_c,j_c,l_c;
1406 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c;
1407 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c;
1408 const PetscInt *idx_c;
1409 PetscInt d;
1410 PetscInt a;
1411 PetscInt max_agg_size;
1412 PetscInt *fine_nodes;
1413 PetscScalar *one_vec;
1414 PetscInt fn_idx;
1415 ISLocalToGlobalMapping ltogmf,ltogmc;
1416
1417 PetscFunctionBegin;
1418 PetscValidHeaderSpecificType(dac,DM_CLASSID,1,DMDA);
1419 PetscValidHeaderSpecificType(daf,DM_CLASSID,2,DMDA);
1420 PetscValidPointer(rest,3);
1421
1422 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr);
1423 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr);
1424 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1425 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1426 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1427 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1428 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1429
1430 if (Mf < Mc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Mc %D, Mf %D", Mc, Mf);
1431 if (Nf < Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Nc %D, Nf %D", Nc, Nf);
1432 if (Pf < Pc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Pc %D, Pf %D", Pc, Pf);
1433
1434 if (Pc < 0) Pc = 1;
1435 if (Pf < 0) Pf = 1;
1436 if (Nc < 0) Nc = 1;
1437 if (Nf < 0) Nf = 1;
1438
1439 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr);
1440 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr);
1441
1442 ierr = DMGetLocalToGlobalMapping(daf,<ogmf);CHKERRQ(ierr);
1443 ierr = ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);CHKERRQ(ierr);
1444
1445 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr);
1446 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr);
1447
1448 ierr = DMGetLocalToGlobalMapping(dac,<ogmc);CHKERRQ(ierr);
1449 ierr = ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);CHKERRQ(ierr);
1450
1451 /*
1452 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1453 for dimension 1 and 2 respectively.
1454 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1455 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate.
1456 Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1457 */
1458
1459 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1);
1460
1461 /* create the matrix that will contain the restriction operator */
1462 ierr = MatCreateAIJ(PetscObjectComm((PetscObject)daf), m_c*n_c*p_c*dofc, m_f*n_f*p_f*doff, Mc*Nc*Pc*dofc, Mf*Nf*Pf*doff,
1463 max_agg_size, NULL, max_agg_size, NULL, rest);CHKERRQ(ierr);
1464
1465 /* store nodes in the fine grid here */
1466 ierr = PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);CHKERRQ(ierr);
1467 for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0;
1468
1469 /* loop over all coarse nodes */
1470 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) {
1471 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) {
1472 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) {
1473 for (d=0; d<dofc; d++) {
1474 /* convert to local "natural" numbering and then to PETSc global numbering */
1475 a = idx_c[dofc*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c))] + d;
1476
1477 fn_idx = 0;
1478 /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1479 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1480 (same for other dimensions)
1481 */
1482 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) {
1483 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) {
1484 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) {
1485 fine_nodes[fn_idx] = idx_f[doff*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))] + d;
1486 fn_idx++;
1487 }
1488 }
1489 }
1490 /* add all these points to one aggregate */
1491 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr);
1492 }
1493 }
1494 }
1495 }
1496 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);CHKERRQ(ierr);
1497 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);CHKERRQ(ierr);
1498 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr);
1499 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1500 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1501 PetscFunctionReturn(0);
1502 }
1503