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