1 #include <petscdmda.h> /*I "petscdmda.h" I*/ 2 #include <petsc-private/dmimpl.h> 3 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 4 5 /* This structure holds the user-provided DMDA callbacks */ 6 typedef struct { 7 PetscErrorCode (*ifunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*,void*); 8 PetscErrorCode (*rhsfunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*); 9 PetscErrorCode (*ijacobianlocal)(DMDALocalInfo*,PetscReal,void*,void*,PetscReal,Mat,Mat,MatStructure*,void*); 10 PetscErrorCode (*rhsjacobianlocal)(DMDALocalInfo*,PetscReal,void*,Mat,Mat,MatStructure*,void*); 11 void *ifunctionlocalctx; 12 void *ijacobianlocalctx; 13 void *rhsfunctionlocalctx; 14 void *rhsjacobianlocalctx; 15 InsertMode ifunctionlocalimode; 16 InsertMode rhsfunctionlocalimode; 17 } DMTS_DA; 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "DMTSDestroy_DMDA" 21 static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm) 22 { 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = PetscFree(sdm->data);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "DMTSDuplicate_DMDA" 32 static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm) 33 { 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 ierr = PetscNewLog(sdm,DMTS_DA,&sdm->data);CHKERRQ(ierr); 38 if (oldsdm->data) {ierr = PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));CHKERRQ(ierr);} 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "DMDATSGetContext" 44 static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats) 45 { 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 *dmdats = PETSC_NULL; 50 if (!sdm->data) { 51 ierr = PetscNewLog(dm,DMTS_DA,&sdm->data);CHKERRQ(ierr); 52 sdm->ops->destroy = DMTSDestroy_DMDA; 53 sdm->ops->duplicate = DMTSDuplicate_DMDA; 54 } 55 *dmdats = (DMTS_DA*)sdm->data; 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "TSComputeIFunction_DMDA" 61 /* 62 This function should eventually replace: 63 DMDAComputeFunction() and DMDAComputeFunction1() 64 */ 65 static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx) 66 { 67 PetscErrorCode ierr; 68 DM dm; 69 DMTS_DA *dmdats = (DMTS_DA*)ctx; 70 DMDALocalInfo info; 71 Vec Xloc; 72 void *x,*f,*xdot; 73 74 PetscFunctionBegin; 75 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 76 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 77 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 78 if (!dmdats->ifunctionlocal) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Corrupt context"); 79 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 80 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 81 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 82 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 83 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 84 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 85 ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr); 86 switch (dmdats->ifunctionlocalimode) { 87 case INSERT_VALUES: { 88 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr); 89 CHKMEMQ; 90 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr); 91 CHKMEMQ; 92 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr); 93 } break; 94 case ADD_VALUES: { 95 Vec Floc; 96 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 97 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 98 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr); 99 CHKMEMQ; 100 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr); 101 CHKMEMQ; 102 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr); 103 ierr = VecZeroEntries(F);CHKERRQ(ierr); 104 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 105 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 106 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 107 } break; 108 default: SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode); 109 } 110 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 111 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 112 ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "TSComputeIJacobian_DMDA" 118 static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *mstr,void *ctx) 119 { 120 PetscErrorCode ierr; 121 DM dm; 122 DMTS_DA *dmdats = (DMTS_DA*)ctx; 123 DMDALocalInfo info; 124 Vec Xloc; 125 void *x,*xdot; 126 127 PetscFunctionBegin; 128 if (!dmdats->ifunctionlocal) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Corrupt context"); 129 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 130 131 if (dmdats->ijacobianlocal) { 132 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 133 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 134 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 135 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 136 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 137 ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr); 138 CHKMEMQ; 139 ierr = (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,*A,*B,mstr,dmdats->ijacobianlocalctx);CHKERRQ(ierr); 140 CHKMEMQ; 141 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 142 ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr); 143 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 144 } else { 145 SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"TSComputeIJacobian_DMDA() called without calling DMDATSSetIJacobian()"); 146 } 147 /* This will be redundant if the user called both, but it's too common to forget. */ 148 if (*A != *B) { 149 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 } 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "TSComputeRHSFunction_DMDA" 157 /* 158 This function should eventually replace: 159 DMDAComputeFunction() and DMDAComputeFunction1() 160 */ 161 static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx) 162 { 163 PetscErrorCode ierr; 164 DM dm; 165 DMTS_DA *dmdats = (DMTS_DA*)ctx; 166 DMDALocalInfo info; 167 Vec Xloc; 168 void *x,*f; 169 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 172 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 173 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 174 if (!dmdats->rhsfunctionlocal) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Corrupt context"); 175 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 176 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 177 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 178 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 179 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 180 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 181 switch (dmdats->rhsfunctionlocalimode) { 182 case INSERT_VALUES: { 183 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr); 184 CHKMEMQ; 185 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr); 186 CHKMEMQ; 187 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr); 188 } break; 189 case ADD_VALUES: { 190 Vec Floc; 191 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 192 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 193 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr); 194 CHKMEMQ; 195 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr); 196 CHKMEMQ; 197 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr); 198 ierr = VecZeroEntries(F);CHKERRQ(ierr); 199 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 200 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 201 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 202 } break; 203 default: SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode); 204 } 205 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 206 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "TSComputeRHSJacobian_DMDA" 212 static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx) 213 { 214 PetscErrorCode ierr; 215 DM dm; 216 DMTS_DA *dmdats = (DMTS_DA*)ctx; 217 DMDALocalInfo info; 218 Vec Xloc; 219 void *x; 220 221 PetscFunctionBegin; 222 if (!dmdats->rhsfunctionlocal) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Corrupt context"); 223 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 224 225 if (dmdats->rhsjacobianlocal) { 226 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 227 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 228 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 229 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 230 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 231 CHKMEMQ; 232 ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr); 233 CHKMEMQ; 234 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 235 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 236 } else { 237 SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()"); 238 } 239 /* This will be redundant if the user called both, but it's too common to forget. */ 240 if (*A != *B) { 241 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 242 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 243 } 244 PetscFunctionReturn(0); 245 } 246 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "DMDATSSetRHSFunctionLocal" 250 /*@C 251 DMDATSSetRHSFunctionLocal - set a local residual evaluation function 252 253 Logically Collective 254 255 Input Arguments: 256 + dm - DM to associate callback with 257 . imode - insert mode for the residual 258 . func - local residual evaluation 259 - ctx - optional context for local residual evaluation 260 261 Calling sequence for func: 262 263 $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx) 264 265 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 266 . t - time at which to evaluate residual 267 . x - array of local state information 268 . f - output array of local residual information 269 - ctx - optional user context 270 271 Level: beginner 272 273 .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal() 274 @*/ 275 PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx) 276 { 277 PetscErrorCode ierr; 278 DMTS sdm; 279 DMTS_DA *dmdats; 280 281 PetscFunctionBegin; 282 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 283 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 284 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 285 dmdats->rhsfunctionlocalimode = imode; 286 dmdats->rhsfunctionlocal = func; 287 dmdats->rhsfunctionlocalctx = ctx; 288 ierr = DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);CHKERRQ(ierr); 289 PetscFunctionReturn(0); 290 } 291 292 #undef __FUNCT__ 293 #define __FUNCT__ "DMDATSSetRHSJacobianLocal" 294 /*@C 295 DMDATSSetRHSJacobianLocal - set a local residual evaluation function 296 297 Logically Collective 298 299 Input Arguments: 300 + dm - DM to associate callback with 301 . func - local RHS Jacobian evaluation routine 302 - ctx - optional context for local jacobian evaluation 303 304 Calling sequence for func: 305 306 $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,MatStructure *flg,void *ctx); 307 308 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 309 . t - time at which to evaluate residual 310 . x - array of local state information 311 . J - Jacobian matrix 312 . B - preconditioner matrix; often same as J 313 . flg - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators()) 314 - ctx - optional context passed above 315 316 Level: beginner 317 318 .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal() 319 @*/ 320 PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx) 321 { 322 PetscErrorCode ierr; 323 DMTS sdm; 324 DMTS_DA *dmdats; 325 326 PetscFunctionBegin; 327 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 328 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 329 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 330 dmdats->rhsjacobianlocal = func; 331 dmdats->rhsjacobianlocalctx = ctx; 332 ierr = DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);CHKERRQ(ierr); 333 PetscFunctionReturn(0); 334 } 335 336 337 #undef __FUNCT__ 338 #define __FUNCT__ "DMDATSSetIFunctionLocal" 339 /*@C 340 DMDATSSetIFunctionLocal - set a local residual evaluation function 341 342 Logically Collective 343 344 Input Arguments: 345 + dm - DM to associate callback with 346 . func - local residual evaluation 347 - ctx - optional context for local residual evaluation 348 349 Calling sequence for func: 350 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 351 . t - time at which to evaluate residual 352 . x - array of local state information 353 . xdot - array of local time derivative information 354 . f - output array of local function evaluation information 355 - ctx - optional context passed above 356 357 Level: beginner 358 359 .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal() 360 @*/ 361 PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx) 362 { 363 PetscErrorCode ierr; 364 DMTS sdm; 365 DMTS_DA *dmdats; 366 367 PetscFunctionBegin; 368 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 369 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 370 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 371 dmdats->ifunctionlocalimode = imode; 372 dmdats->ifunctionlocal = func; 373 dmdats->ifunctionlocalctx = ctx; 374 ierr = DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "DMDATSSetIJacobianLocal" 380 /*@C 381 DMDATSSetIJacobianLocal - set a local residual evaluation function 382 383 Logically Collective 384 385 Input Arguments: 386 + dm - DM to associate callback with 387 . func - local residual evaluation 388 - ctx - optional context for local residual evaluation 389 390 Calling sequence for func: 391 392 $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,Mat J,Mat B,MatStructure *flg,void *ctx); 393 394 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 395 . t - time at which to evaluate the jacobian 396 . x - array of local state information 397 . xdot - time derivative at this state 398 . J - Jacobian matrix 399 . B - preconditioner matrix; often same as J 400 . flg - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators()) 401 - ctx - optional context passed above 402 403 Level: beginner 404 405 .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal() 406 @*/ 407 PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx) 408 { 409 PetscErrorCode ierr; 410 DMTS sdm; 411 DMTS_DA *dmdats; 412 413 PetscFunctionBegin; 414 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 415 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 416 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 417 dmdats->ijacobianlocal = func; 418 dmdats->ijacobianlocalctx = ctx; 419 ierr = DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);CHKERRQ(ierr); 420 PetscFunctionReturn(0); 421 } 422