1 #include <private/dmimpl.h> /*I "petscdm.h" I*/ 2 #include <petscdmredundant.h> /*I "petscdmredundant.h" I*/ 3 #include <petscmat.h> /*I "petscmat.h" I*/ 4 5 typedef struct { 6 PetscInt rank; /* owner */ 7 PetscInt N; /* total number of dofs */ 8 PetscInt n; /* owned number of dofs, n=N on owner, n=0 on non-owners */ 9 } DM_Redundant; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "DMGetMatrix_Redundant" 13 static PetscErrorCode DMGetMatrix_Redundant(DM dm,const MatType mtype,Mat *J) 14 { 15 DM_Redundant *red = (DM_Redundant*)dm->data; 16 PetscErrorCode ierr; 17 ISLocalToGlobalMapping ltog,ltogb; 18 19 PetscFunctionBegin; 20 ierr = MatCreate(((PetscObject)dm)->comm,J);CHKERRQ(ierr); 21 ierr = MatSetSizes(*J,red->n,red->n,red->N,red->N);CHKERRQ(ierr); 22 ierr = MatSetType(*J,mtype);CHKERRQ(ierr); 23 ierr = MatSeqAIJSetPreallocation(*J,red->n,PETSC_NULL);CHKERRQ(ierr); 24 ierr = MatSeqBAIJSetPreallocation(*J,1,red->n,PETSC_NULL);CHKERRQ(ierr); 25 ierr = MatMPIAIJSetPreallocation(*J,red->n,PETSC_NULL,red->N-red->n,PETSC_NULL);CHKERRQ(ierr); 26 ierr = MatMPIBAIJSetPreallocation(*J,1,red->n,PETSC_NULL,red->N-red->n,PETSC_NULL);CHKERRQ(ierr); 27 28 ierr = DMGetLocalToGlobalMapping(dm,<og);CHKERRQ(ierr); 29 ierr = DMGetLocalToGlobalMappingBlock(dm,<ogb);CHKERRQ(ierr); 30 ierr = MatSetLocalToGlobalMapping(*J,ltog,ltog);CHKERRQ(ierr); 31 ierr = MatSetLocalToGlobalMappingBlock(*J,ltogb,ltogb);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "DMDestroy_Redundant" 37 static PetscErrorCode DMDestroy_Redundant(DM dm) 38 { 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PetscObjectComposeFunctionDynamic((PetscObject)dm,"DMRedundantSetSize_C","",PETSC_NULL);CHKERRQ(ierr); 43 ierr = PetscObjectComposeFunctionDynamic((PetscObject)dm,"DMRedundantGetSize_C","",PETSC_NULL);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "DMCreateGlobalVector_Redundant" 49 static PetscErrorCode DMCreateGlobalVector_Redundant(DM dm,Vec *gvec) 50 { 51 PetscErrorCode ierr; 52 DM_Redundant *red = (DM_Redundant*)dm->data; 53 ISLocalToGlobalMapping ltog,ltogb; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 57 PetscValidPointer(gvec,2); 58 *gvec = 0; 59 ierr = VecCreate(((PetscObject)dm)->comm,gvec);CHKERRQ(ierr); 60 ierr = VecSetSizes(*gvec,red->n,red->N);CHKERRQ(ierr); 61 ierr = VecSetType(*gvec,dm->vectype);CHKERRQ(ierr); 62 ierr = DMGetLocalToGlobalMapping(dm,<og);CHKERRQ(ierr); 63 ierr = DMGetLocalToGlobalMappingBlock(dm,<ogb);CHKERRQ(ierr); 64 ierr = VecSetLocalToGlobalMapping(*gvec,ltog);CHKERRQ(ierr); 65 ierr = VecSetLocalToGlobalMappingBlock(*gvec,ltog);CHKERRQ(ierr); 66 ierr = PetscObjectCompose((PetscObject)*gvec,"DM",(PetscObject)dm);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "DMCreateLocalVector_Redundant" 72 static PetscErrorCode DMCreateLocalVector_Redundant(DM dm,Vec *lvec) 73 { 74 PetscErrorCode ierr; 75 DM_Redundant *red = (DM_Redundant*)dm->data; 76 77 PetscFunctionBegin; 78 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 79 PetscValidPointer(lvec,2); 80 *lvec = 0; 81 ierr = VecCreate(PETSC_COMM_SELF,lvec);CHKERRQ(ierr); 82 ierr = VecSetSizes(*lvec,red->N,red->N);CHKERRQ(ierr); 83 ierr = VecSetType(*lvec,dm->vectype);CHKERRQ(ierr); 84 ierr = PetscObjectCompose((PetscObject)*lvec,"DM",(PetscObject)dm);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "DMLocalToGlobalBegin_Redundant" 90 static PetscErrorCode DMLocalToGlobalBegin_Redundant(DM dm,Vec l,InsertMode imode,Vec g) 91 { 92 PetscErrorCode ierr; 93 DM_Redundant *red = (DM_Redundant*)dm->data; 94 const PetscScalar *lv; 95 PetscScalar *gv; 96 PetscMPIInt rank; 97 98 PetscFunctionBegin; 99 ierr = MPI_Comm_rank(((PetscObject)dm)->comm,&rank);CHKERRQ(ierr); 100 ierr = VecGetArrayRead(l,&lv);CHKERRQ(ierr); 101 ierr = VecGetArray(g,&gv);CHKERRQ(ierr); 102 switch (imode) { 103 case ADD_VALUES: 104 case MAX_VALUES: 105 { 106 void *source; 107 PetscScalar *buffer; 108 PetscInt i; 109 if (rank == red->rank) { 110 #if defined (PETSC_HAVE_MPI_IN_PLACE) 111 buffer = gv; 112 source = MPI_IN_PLACE; 113 #else 114 ierr = PetscMalloc(red->N*sizeof(PetscScalar),&buffer);CHKERRQ(ierr); 115 source = buffer; 116 #endif 117 if (imode == ADD_VALUES) for (i=0; i<red->N; i++) buffer[i] = gv[i] + lv[i]; 118 #if !defined(PETSC_USE_COMPLEX) 119 if (imode == MAX_VALUES) for (i=0; i<red->N; i++) buffer[i] = PetscMax(gv[i],lv[i]); 120 #endif 121 } else { 122 source = (void*)lv; 123 } 124 ierr = MPI_Reduce(source,gv,red->N,MPIU_SCALAR,(imode == ADD_VALUES)?MPI_SUM:MPI_MAX,red->rank,((PetscObject)dm)->comm);CHKERRQ(ierr); 125 #if !defined(PETSC_HAVE_MPI_IN_PLACE) 126 if (rank == mine->rank) {ierr = PetscFree(buffer);CHKERRQ(ierr);} 127 #endif 128 } break; 129 case INSERT_VALUES: 130 ierr = PetscMemcpy(gv,lv,red->n*sizeof(PetscScalar));CHKERRQ(ierr); 131 break; 132 default: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"InsertMode not supported"); 133 } 134 ierr = VecRestoreArrayRead(l,&lv);CHKERRQ(ierr); 135 ierr = VecRestoreArray(g,&gv);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "DMLocalToGlobalEnd_Redundant" 141 static PetscErrorCode DMLocalToGlobalEnd_Redundant(DM dm,Vec l,InsertMode imode,Vec g) 142 { 143 PetscFunctionBegin; 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "DMGlobalToLocalBegin_Redundant" 149 static PetscErrorCode DMGlobalToLocalBegin_Redundant(DM dm,Vec g,InsertMode imode,Vec l) 150 { 151 PetscErrorCode ierr; 152 DM_Redundant *red = (DM_Redundant*)dm->data; 153 const PetscScalar *gv; 154 PetscScalar *lv; 155 156 PetscFunctionBegin; 157 ierr = VecGetArrayRead(g,&gv);CHKERRQ(ierr); 158 ierr = VecGetArray(l,&lv);CHKERRQ(ierr); 159 switch (imode) { 160 case INSERT_VALUES: 161 if (red->n) {ierr = PetscMemcpy(lv,gv,red->n*sizeof(PetscScalar));CHKERRQ(ierr);} 162 ierr = MPI_Bcast(lv,red->N,MPIU_SCALAR,red->rank,((PetscObject)dm)->comm);CHKERRQ(ierr); 163 break; 164 default: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"InsertMode not supported"); 165 } 166 ierr = VecRestoreArrayRead(g,&gv);CHKERRQ(ierr); 167 ierr = VecRestoreArray(l,&lv);CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "DMGlobalToLocalEnd_Redundant" 173 static PetscErrorCode DMGlobalToLocalEnd_Redundant(DM dm,Vec g,InsertMode imode,Vec l) 174 { 175 PetscFunctionBegin; 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "DMSetUp_Redundant" 181 static PetscErrorCode DMSetUp_Redundant(DM dm) 182 { 183 PetscErrorCode ierr; 184 DM_Redundant *red = (DM_Redundant*)dm->data; 185 PetscInt i,*globals; 186 187 PetscFunctionBegin; 188 ierr = PetscMalloc(red->N*sizeof(PetscInt),&globals);CHKERRQ(ierr); 189 for (i=0; i<red->N; i++) globals[i] = i; 190 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,red->N,globals,PETSC_OWN_POINTER,&dm->ltogmap);CHKERRQ(ierr); 191 ierr = PetscObjectReference((PetscObject)dm->ltogmap);CHKERRQ(ierr); 192 dm->ltogmapb = dm->ltogmap; 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "DMView_Redundant" 198 static PetscErrorCode DMView_Redundant(DM dm,PetscViewer viewer) 199 { 200 PetscErrorCode ierr; 201 DM_Redundant *red = (DM_Redundant*)dm->data; 202 PetscBool iascii; 203 204 PetscFunctionBegin; 205 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 206 if (iascii) { 207 ierr = PetscViewerASCIIPrintf(viewer,"redundant: rank=%D N=%D\n",red->rank,red->N);CHKERRQ(ierr); 208 } 209 PetscFunctionReturn(0); 210 } 211 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "DMRefine_Redundant" 215 static PetscErrorCode DMRefine_Redundant(DM dmc,MPI_Comm comm,DM *dmf) 216 { 217 PetscErrorCode ierr; 218 PetscMPIInt flag; 219 DM_Redundant *redc = (DM_Redundant*)dmc->data; 220 221 PetscFunctionBegin; 222 ierr = MPI_Comm_compare(((PetscObject)dmc)->comm,comm,&flag); CHKERRQ(ierr); 223 if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(((PetscObject)dmc)->comm,PETSC_ERR_SUP,"cannot change communicators"); 224 ierr = DMRedundantCreate(comm,redc->rank,redc->N,dmf);CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 228 #undef __FUNCT__ 229 #define __FUNCT__ "DMCoarsen_Redundant" 230 static PetscErrorCode DMCoarsen_Redundant(DM dmf,MPI_Comm comm,DM *dmc) 231 { 232 PetscErrorCode ierr; 233 PetscMPIInt flag; 234 DM_Redundant *redf = (DM_Redundant*)dmf->data; 235 236 PetscFunctionBegin; 237 ierr = MPI_Comm_compare(((PetscObject)dmf)->comm,comm,&flag); CHKERRQ(ierr); 238 if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(((PetscObject)dmf)->comm,PETSC_ERR_SUP,"cannot change communicators"); 239 ierr = DMRedundantCreate(comm,redf->rank,redf->N,dmc);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "DMGetInterpolation_Redundant" 245 static PetscErrorCode DMGetInterpolation_Redundant(DM dmc,DM dmf,Mat *P,Vec *scale) 246 { 247 PetscErrorCode ierr; 248 DM_Redundant *redc = (DM_Redundant*)dmc->data; 249 DM_Redundant *redf = (DM_Redundant*)dmf->data; 250 PetscMPIInt flag; 251 PetscInt i,rstart,rend; 252 253 PetscFunctionBegin; 254 ierr = MPI_Comm_compare(((PetscObject)dmc)->comm,((PetscObject)dmf)->comm,&flag); CHKERRQ(ierr); 255 if (flag != MPI_CONGRUENT && flag != MPI_IDENT) SETERRQ(((PetscObject)dmf)->comm,PETSC_ERR_SUP,"cannot change communicators"); 256 if (redc->rank != redf->rank) SETERRQ(((PetscObject)dmf)->comm,PETSC_ERR_ARG_INCOMP,"Owning rank does not match"); 257 if (redc->N != redf->N) SETERRQ(((PetscObject)dmf)->comm,PETSC_ERR_ARG_INCOMP,"Global size does not match"); 258 ierr = MatCreate(((PetscObject)dmc)->comm,P);CHKERRQ(ierr); 259 ierr = MatSetSizes(*P,redc->n,redc->n,redc->N,redc->N);CHKERRQ(ierr); 260 ierr = MatSetType(*P,MATAIJ);CHKERRQ(ierr); 261 ierr = MatSeqAIJSetPreallocation(*P,1,0);CHKERRQ(ierr); 262 ierr = MatMPIAIJSetPreallocation(*P,1,0,0,0);CHKERRQ(ierr); 263 ierr = MatGetOwnershipRange(*P,&rstart,&rend);CHKERRQ(ierr); 264 for (i=rstart; i<rend; i++) {ierr = MatSetValue(*P,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr);} 265 ierr = MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 266 ierr = MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 267 if (scale) {ierr = DMGetInterpolationScale(dmc,dmf,*P,scale);CHKERRQ(ierr);} 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "DMRedundantSetSize" 273 /*@ 274 DMRedundantSetSize - Sets the size of a densely coupled redundant object 275 276 Collective on DM 277 278 Input Parameter: 279 + dm - redundant DM 280 . rank - rank of process to own redundant degrees of freedom 281 - N - total number of redundant degrees of freedom 282 283 Level: advanced 284 285 .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantGetSize() 286 @*/ 287 PetscErrorCode DMRedundantSetSize(DM dm,PetscInt rank,PetscInt N) 288 { 289 PetscErrorCode ierr; 290 291 PetscFunctionBegin; 292 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 293 PetscValidType(dm,1); 294 PetscValidLogicalCollectiveInt(dm,rank,2); 295 PetscValidLogicalCollectiveInt(dm,N,3); 296 ierr = PetscTryMethod(dm,"DMRedundantSetSize_C",(DM,PetscInt,PetscInt),(dm,rank,N));CHKERRQ(ierr); 297 PetscFunctionReturn(0); 298 } 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "DMRedundantGetSize" 302 /*@ 303 DMRedundantGetSize - Gets the size of a densely coupled redundant object 304 305 Not Collective 306 307 Input Parameter: 308 + dm - redundant DM 309 310 Output Parameters: 311 + rank - rank of process to own redundant degrees of freedom (or PETSC_NULL) 312 - N - total number of redundant degrees of freedom (or PETSC_NULL) 313 314 Level: advanced 315 316 .seealso DMDestroy(), DMCreateGlobalVector(), DMRedundantCreate(), DMRedundantSetSize() 317 @*/ 318 PetscErrorCode DMRedundantGetSize(DM dm,PetscInt *rank,PetscInt *N) 319 { 320 PetscErrorCode ierr; 321 322 PetscFunctionBegin; 323 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 324 PetscValidType(dm,1); 325 ierr = PetscUseMethod(dm,"DMRedundantGetSize_C",(DM,PetscInt*,PetscInt*),(dm,rank,N));CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328 329 EXTERN_C_BEGIN 330 #undef __FUNCT__ 331 #define __FUNCT__ "DMRedundantSetSize_Redundant" 332 PetscErrorCode DMRedundantSetSize_Redundant(DM dm,PetscInt rank,PetscInt N) 333 { 334 DM_Redundant *red = (DM_Redundant*)dm->data; 335 PetscErrorCode ierr; 336 PetscMPIInt myrank; 337 338 PetscFunctionBegin; 339 ierr = MPI_Comm_rank(((PetscObject)dm)->comm,&myrank);CHKERRQ(ierr); 340 red->rank = rank; 341 red->N = N; 342 red->n = (myrank == rank) ? N : 0; 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "DMRedundantGetSize_Redundant" 348 PetscErrorCode DMRedundantGetSize_Redundant(DM dm,PetscInt *rank,PetscInt *N) 349 { 350 DM_Redundant *red = (DM_Redundant*)dm->data; 351 352 PetscFunctionBegin; 353 if (rank) *rank = red->rank; 354 if (N) *N = red->N; 355 PetscFunctionReturn(0); 356 } 357 EXTERN_C_END 358 359 EXTERN_C_BEGIN 360 #undef __FUNCT__ 361 #define __FUNCT__ "DMCreate_Redundant" 362 PetscErrorCode DMCreate_Redundant(DM dm) 363 { 364 PetscErrorCode ierr; 365 DM_Redundant *red; 366 367 PetscFunctionBegin; 368 ierr = PetscNewLog(dm,DM_Redundant,&red);CHKERRQ(ierr); 369 dm->data = red; 370 371 ierr = PetscObjectChangeTypeName((PetscObject)dm,DMREDUNDANT);CHKERRQ(ierr); 372 dm->ops->setup = DMSetUp_Redundant; 373 dm->ops->view = DMView_Redundant; 374 dm->ops->createglobalvector = DMCreateGlobalVector_Redundant; 375 dm->ops->createlocalvector = DMCreateLocalVector_Redundant; 376 dm->ops->getmatrix = DMGetMatrix_Redundant; 377 dm->ops->destroy = DMDestroy_Redundant; 378 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant; 379 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Redundant; 380 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant; 381 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Redundant; 382 dm->ops->refine = DMRefine_Redundant; 383 dm->ops->coarsen = DMCoarsen_Redundant; 384 dm->ops->getinterpolation = DMGetInterpolation_Redundant; 385 ierr = PetscStrallocpy(VECSTANDARD,&dm->vectype);CHKERRQ(ierr); 386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)dm,"DMRedundantSetSize_C","DMRedundantSetSize_Redundant",DMRedundantSetSize_Redundant);CHKERRQ(ierr); 387 ierr = PetscObjectComposeFunctionDynamic((PetscObject)dm,"DMRedundantGetSize_C","DMRedundantGetSize_Redundant",DMRedundantGetSize_Redundant);CHKERRQ(ierr); 388 PetscFunctionReturn(0); 389 } 390 EXTERN_C_END 391 392 #undef __FUNCT__ 393 #define __FUNCT__ "DMRedundantCreate" 394 /*@C 395 DMRedundantCreate - Creates a DM object, used to manage data for dense globally coupled variables 396 397 Collective on MPI_Comm 398 399 Input Parameter: 400 + comm - the processors that will share the global vector 401 . rank - rank to own the redundant values 402 - N - total number of degrees of freedom 403 404 Output Parameters: 405 . red - the redundant DM 406 407 Level: advanced 408 409 .seealso DMDestroy(), DMCreateGlobalVector(), DMGetMatrix(), DMCompositeAddDM() 410 411 @*/ 412 PetscErrorCode DMRedundantCreate(MPI_Comm comm,PetscInt rank,PetscInt N,DM *dm) 413 { 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 PetscValidPointer(dm,2); 418 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 419 ierr = DMSetType(*dm,DMREDUNDANT);CHKERRQ(ierr); 420 ierr = DMRedundantSetSize(*dm,rank,N);CHKERRQ(ierr); 421 ierr = DMSetUp(*dm);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424