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