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