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 */ 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 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 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 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);CHKERRMPI(ierr); 329 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRMPI(ierr); 330 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRMPI(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 */ 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);CHKERRMPI(ierr); 557 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRMPI(ierr); 558 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRMPI(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 */ 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);CHKERRMPI(ierr); 687 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRMPI(ierr); 688 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRMPI(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 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 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 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 /* loop over local fine grid nodes setting interpolation for those*/ 1147 nc = 0; 1148 ierr = PetscMalloc1(m_f,&cols);CHKERRQ(ierr); 1149 1150 for (i=i_start_c; i<i_start_c+m_c; i++) { 1151 PetscInt i_f = i*ratioi; 1152 1153 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); 1154 1155 row = idx_f[(i_f-i_start_ghost)]; 1156 cols[nc++] = row; 1157 } 1158 1159 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1160 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1161 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1162 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1163 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1164 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1165 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1166 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1167 PetscFunctionReturn(0); 1168 } 1169 1170 PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject) 1171 { 1172 PetscErrorCode ierr; 1173 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof; 1174 const PetscInt *idx_c,*idx_f; 1175 ISLocalToGlobalMapping ltog_f,ltog_c; 1176 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c; 1177 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 1178 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 1179 PetscInt *cols; 1180 DMBoundaryType bx,by; 1181 Vec vecf,vecc; 1182 IS isf; 1183 1184 PetscFunctionBegin; 1185 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,NULL,NULL,NULL,NULL,NULL,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr); 1186 ierr = DMDAGetInfo(daf,NULL,&mx,&my,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 1187 if (bx == DM_BOUNDARY_PERIODIC) { 1188 ratioi = mx/Mx; 1189 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); 1190 } else { 1191 ratioi = (mx-1)/(Mx-1); 1192 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); 1193 } 1194 if (by == DM_BOUNDARY_PERIODIC) { 1195 ratioj = my/My; 1196 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); 1197 } else { 1198 ratioj = (my-1)/(My-1); 1199 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); 1200 } 1201 1202 ierr = DMDAGetCorners(daf,&i_start,&j_start,NULL,&m_f,&n_f,NULL);CHKERRQ(ierr); 1203 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,NULL,&m_ghost,&n_ghost,NULL);CHKERRQ(ierr); 1204 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1205 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1206 1207 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,NULL,&m_c,&n_c,NULL);CHKERRQ(ierr); 1208 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,NULL,&m_ghost_c,&n_ghost_c,NULL);CHKERRQ(ierr); 1209 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1210 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1211 1212 /* loop over local fine grid nodes setting interpolation for those*/ 1213 nc = 0; 1214 ierr = PetscMalloc1(n_f*m_f,&cols);CHKERRQ(ierr); 1215 for (j=j_start_c; j<j_start_c+n_c; j++) { 1216 for (i=i_start_c; i<i_start_c+m_c; i++) { 1217 PetscInt i_f = i*ratioi,j_f = j*ratioj; 1218 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\ 1219 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1220 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\ 1221 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1222 row = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1223 cols[nc++] = row; 1224 } 1225 } 1226 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1227 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1228 1229 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1230 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1231 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1232 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1233 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1234 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1235 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1236 PetscFunctionReturn(0); 1237 } 1238 1239 PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject) 1240 { 1241 PetscErrorCode ierr; 1242 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz; 1243 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c; 1244 PetscInt i_start_ghost,j_start_ghost,k_start_ghost; 1245 PetscInt mx,my,mz,ratioi,ratioj,ratiok; 1246 PetscInt i_start_c,j_start_c,k_start_c; 1247 PetscInt m_c,n_c,p_c; 1248 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c; 1249 PetscInt row,nc,dof; 1250 const PetscInt *idx_c,*idx_f; 1251 ISLocalToGlobalMapping ltog_f,ltog_c; 1252 PetscInt *cols; 1253 DMBoundaryType bx,by,bz; 1254 Vec vecf,vecc; 1255 IS isf; 1256 1257 PetscFunctionBegin; 1258 ierr = DMDAGetInfo(dac,NULL,&Mx,&My,&Mz,NULL,NULL,NULL,NULL,NULL,&bx,&by,&bz,NULL);CHKERRQ(ierr); 1259 ierr = DMDAGetInfo(daf,NULL,&mx,&my,&mz,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 1260 1261 if (bx == DM_BOUNDARY_PERIODIC) { 1262 ratioi = mx/Mx; 1263 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); 1264 } else { 1265 ratioi = (mx-1)/(Mx-1); 1266 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); 1267 } 1268 if (by == DM_BOUNDARY_PERIODIC) { 1269 ratioj = my/My; 1270 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); 1271 } else { 1272 ratioj = (my-1)/(My-1); 1273 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); 1274 } 1275 if (bz == DM_BOUNDARY_PERIODIC) { 1276 ratiok = mz/Mz; 1277 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); 1278 } else { 1279 ratiok = (mz-1)/(Mz-1); 1280 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); 1281 } 1282 1283 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1284 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1285 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1286 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1287 1288 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1289 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); 1290 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1291 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1292 1293 /* loop over local fine grid nodes setting interpolation for those*/ 1294 nc = 0; 1295 ierr = PetscMalloc1(n_f*m_f*p_f,&cols);CHKERRQ(ierr); 1296 for (k=k_start_c; k<k_start_c+p_c; k++) { 1297 for (j=j_start_c; j<j_start_c+n_c; j++) { 1298 for (i=i_start_c; i<i_start_c+m_c; i++) { 1299 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok; 1300 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 " 1301 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost); 1302 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 " 1303 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1304 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 " 1305 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1306 row = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1307 cols[nc++] = row; 1308 } 1309 } 1310 } 1311 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1312 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1313 1314 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1315 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1316 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1317 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1318 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1319 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1320 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,Mat *mat) 1325 { 1326 PetscErrorCode ierr; 1327 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1328 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1329 DMDAStencilType stc,stf; 1330 VecScatter inject = NULL; 1331 1332 PetscFunctionBegin; 1333 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1334 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1335 PetscValidPointer(mat,3); 1336 1337 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1338 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1339 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf); 1340 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff); 1341 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf); 1342 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs"); 1343 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs"); 1344 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1345 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1346 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1347 1348 if (dimc == 1) { 1349 ierr = DMCreateInjection_DA_1D(dac,daf,&inject);CHKERRQ(ierr); 1350 } else if (dimc == 2) { 1351 ierr = DMCreateInjection_DA_2D(dac,daf,&inject);CHKERRQ(ierr); 1352 } else if (dimc == 3) { 1353 ierr = DMCreateInjection_DA_3D(dac,daf,&inject);CHKERRQ(ierr); 1354 } 1355 ierr = MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);CHKERRQ(ierr); 1356 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 1357 PetscFunctionReturn(0); 1358 } 1359 1360 /*@ 1361 DMCreateAggregates - Deprecated, see DMDACreateAggregates. 1362 1363 Level: intermediate 1364 @*/ 1365 PetscErrorCode DMCreateAggregates(DM dac,DM daf,Mat *mat) 1366 { 1367 return DMDACreateAggregates(dac,daf,mat); 1368 } 1369 1370 /*@ 1371 DMDACreateAggregates - Gets the aggregates that map between 1372 grids associated with two DMDAs. 1373 1374 Collective on dmc 1375 1376 Input Parameters: 1377 + dmc - the coarse grid DMDA 1378 - dmf - the fine grid DMDA 1379 1380 Output Parameters: 1381 . rest - the restriction matrix (transpose of the projection matrix) 1382 1383 Level: intermediate 1384 1385 Note: This routine is not used by PETSc. 1386 It is not clear what its use case is and it may be removed in a future release. 1387 Users should contact petsc-maint@mcs.anl.gov if they plan to use it. 1388 1389 .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation() 1390 @*/ 1391 PetscErrorCode DMDACreateAggregates(DM dac,DM daf,Mat *rest) 1392 { 1393 PetscErrorCode ierr; 1394 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1395 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1396 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1397 DMDAStencilType stc,stf; 1398 PetscInt i,j,l; 1399 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1400 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1401 const PetscInt *idx_f; 1402 PetscInt i_c,j_c,l_c; 1403 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1404 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1405 const PetscInt *idx_c; 1406 PetscInt d; 1407 PetscInt a; 1408 PetscInt max_agg_size; 1409 PetscInt *fine_nodes; 1410 PetscScalar *one_vec; 1411 PetscInt fn_idx; 1412 ISLocalToGlobalMapping ltogmf,ltogmc; 1413 1414 PetscFunctionBegin; 1415 PetscValidHeaderSpecificType(dac,DM_CLASSID,1,DMDA); 1416 PetscValidHeaderSpecificType(daf,DM_CLASSID,2,DMDA); 1417 PetscValidPointer(rest,3); 1418 1419 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1420 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1421 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf); 1422 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff); 1423 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf); 1424 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs"); 1425 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs"); 1426 1427 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); 1428 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); 1429 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); 1430 1431 if (Pc < 0) Pc = 1; 1432 if (Pf < 0) Pf = 1; 1433 if (Nc < 0) Nc = 1; 1434 if (Nf < 0) Nf = 1; 1435 1436 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1437 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1438 1439 ierr = DMGetLocalToGlobalMapping(daf,<ogmf);CHKERRQ(ierr); 1440 ierr = ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1441 1442 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1443 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); 1444 1445 ierr = DMGetLocalToGlobalMapping(dac,<ogmc);CHKERRQ(ierr); 1446 ierr = ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1447 1448 /* 1449 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1450 for dimension 1 and 2 respectively. 1451 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1452 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1453 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1454 */ 1455 1456 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1457 1458 /* create the matrix that will contain the restriction operator */ 1459 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, 1460 max_agg_size, NULL, max_agg_size, NULL, rest);CHKERRQ(ierr); 1461 1462 /* store nodes in the fine grid here */ 1463 ierr = PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);CHKERRQ(ierr); 1464 for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1465 1466 /* loop over all coarse nodes */ 1467 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1468 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1469 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1470 for (d=0; d<dofc; d++) { 1471 /* convert to local "natural" numbering and then to PETSc global numbering */ 1472 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; 1473 1474 fn_idx = 0; 1475 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1476 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1477 (same for other dimensions) 1478 */ 1479 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1480 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1481 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1482 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; 1483 fn_idx++; 1484 } 1485 } 1486 } 1487 /* add all these points to one aggregate */ 1488 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1489 } 1490 } 1491 } 1492 } 1493 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1494 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1495 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1496 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1497 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1498 PetscFunctionReturn(0); 1499 } 1500