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