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