1 static char help[] = "Simple wrapper object to solve DAE of the form:\n\ 2 \\dot{U} = f(U,V)\n\ 3 F(U,V) = 0\n\n"; 4 5 #include <petscts.h> 6 7 /* ----------------------------------------------------------------------------*/ 8 9 typedef struct _p_TSDAESimple *TSDAESimple; 10 struct _p_TSDAESimple { 11 MPI_Comm comm; 12 PetscErrorCode (*setfromoptions)(PetscOptionItems*,TSDAESimple); 13 PetscErrorCode (*solve)(TSDAESimple,Vec); 14 PetscErrorCode (*destroy)(TSDAESimple); 15 Vec U,V; 16 PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*); 17 PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*); 18 void *fctx,*Fctx; 19 void *data; 20 }; 21 22 PetscErrorCode TSDAESimpleCreate(MPI_Comm comm,TSDAESimple *tsdae) 23 { 24 PetscErrorCode ierr; 25 26 PetscFunctionBegin; 27 ierr = PetscNew(tsdae);CHKERRQ(ierr); 28 (*tsdae)->comm = comm; 29 PetscFunctionReturn(0); 30 } 31 32 PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx) 33 { 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 tsdae->f = f; 38 tsdae->U = U; 39 ierr = PetscObjectReference((PetscObject)U);CHKERRQ(ierr); 40 tsdae->fctx = ctx; 41 PetscFunctionReturn(0); 42 } 43 44 PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx) 45 { 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 tsdae->F = F; 50 tsdae->V = V; 51 ierr = PetscObjectReference((PetscObject)V);CHKERRQ(ierr); 52 tsdae->Fctx = ctx; 53 PetscFunctionReturn(0); 54 } 55 56 PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae) 57 { 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = (*(*tsdae)->destroy)(*tsdae);CHKERRQ(ierr); 62 ierr = VecDestroy(&(*tsdae)->U);CHKERRQ(ierr); 63 ierr = VecDestroy(&(*tsdae)->V);CHKERRQ(ierr); 64 ierr = PetscFree(*tsdae);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution) 69 { 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = (*tsdae->solve)(tsdae,Usolution);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae) 78 { 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 ierr = (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 /* ----------------------------------------------------------------------------*/ 87 /* 88 Integrates the system by integrating the reduced ODE system and solving the 89 algebraic constraints at each stage with a separate SNES solve. 90 */ 91 92 typedef struct { 93 PetscReal t; 94 TS ts; 95 SNES snes; 96 Vec U; 97 } TSDAESimple_Reduced; 98 99 /* 100 Defines the RHS function that is passed to the time-integrator. 101 102 Solves F(U,V) for V and then computes f(U,V) 103 104 */ 105 PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 106 { 107 TSDAESimple tsdae = (TSDAESimple)actx; 108 TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data; 109 PetscErrorCode ierr; 110 111 PetscFunctionBeginUser; 112 red->t = t; 113 red->U = U; 114 ierr = SNESSolve(red->snes,NULL,tsdae->V);CHKERRQ(ierr); 115 ierr = (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 /* 120 Defines the nonlinear function that is passed to the nonlinear solver 121 122 */ 123 PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx) 124 { 125 TSDAESimple tsdae = (TSDAESimple)actx; 126 TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data; 127 PetscErrorCode ierr; 128 129 PetscFunctionBeginUser; 130 ierr = (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U) 135 { 136 PetscErrorCode ierr; 137 TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data; 138 139 PetscFunctionBegin; 140 ierr = TSSolve(red->ts,U);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae) 145 { 146 PetscErrorCode ierr; 147 TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data; 148 149 PetscFunctionBegin; 150 ierr = TSSetFromOptions(red->ts);CHKERRQ(ierr); 151 ierr = SNESSetFromOptions(red->snes);CHKERRQ(ierr); 152 PetscFunctionReturn(0); 153 } 154 155 PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae) 156 { 157 PetscErrorCode ierr; 158 TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data; 159 160 PetscFunctionBegin; 161 ierr = TSDestroy(&red->ts);CHKERRQ(ierr); 162 ierr = SNESDestroy(&red->snes);CHKERRQ(ierr); 163 ierr = PetscFree(red);CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 167 PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae) 168 { 169 PetscErrorCode ierr; 170 TSDAESimple_Reduced *red; 171 Vec tsrhs; 172 173 PetscFunctionBegin; 174 ierr = PetscNew(&red);CHKERRQ(ierr); 175 tsdae->data = red; 176 177 tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced; 178 tsdae->solve = TSDAESimpleSolve_Reduced; 179 tsdae->destroy = TSDAESimpleDestroy_Reduced; 180 181 ierr = TSCreate(tsdae->comm,&red->ts);CHKERRQ(ierr); 182 ierr = TSSetProblemType(red->ts,TS_NONLINEAR);CHKERRQ(ierr); 183 ierr = TSSetType(red->ts,TSEULER);CHKERRQ(ierr); 184 ierr = TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 185 ierr = VecDuplicate(tsdae->U,&tsrhs);CHKERRQ(ierr); 186 ierr = TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);CHKERRQ(ierr); 187 ierr = VecDestroy(&tsrhs);CHKERRQ(ierr); 188 189 ierr = SNESCreate(tsdae->comm,&red->snes);CHKERRQ(ierr); 190 ierr = SNESSetOptionsPrefix(red->snes,"tsdaesimple_");CHKERRQ(ierr); 191 ierr = SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);CHKERRQ(ierr); 192 ierr = SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 /* ----------------------------------------------------------------------------*/ 197 198 /* 199 Integrates the system by integrating directly the entire DAE system 200 */ 201 202 typedef struct { 203 TS ts; 204 Vec UV,UF,VF; 205 VecScatter scatterU,scatterV; 206 } TSDAESimple_Full; 207 208 /* 209 Defines the RHS function that is passed to the time-integrator. 210 211 f(U,V) 212 0 213 214 */ 215 PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 216 { 217 TSDAESimple tsdae = (TSDAESimple)actx; 218 TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data; 219 PetscErrorCode ierr; 220 221 PetscFunctionBeginUser; 222 ierr = VecSet(F,0.0);CHKERRQ(ierr); 223 ierr = VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 224 ierr = VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 225 ierr = VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 226 ierr = VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 227 ierr = (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);CHKERRQ(ierr); 228 ierr = VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 229 ierr = VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 233 /* 234 Defines the nonlinear function that is passed to the nonlinear solver 235 236 \dot{U} 237 F(U,V) 238 239 */ 240 PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 241 { 242 TSDAESimple tsdae = (TSDAESimple)actx; 243 TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data; 244 PetscErrorCode ierr; 245 246 PetscFunctionBeginUser; 247 ierr = VecCopy(UVdot,F);CHKERRQ(ierr); 248 ierr = VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 249 ierr = VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 250 ierr = VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 251 ierr = VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 252 ierr = (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);CHKERRQ(ierr); 253 ierr = VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 254 ierr = VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U) 259 { 260 PetscErrorCode ierr; 261 TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data; 262 263 PetscFunctionBegin; 264 ierr = VecSet(full->UV,1.0);CHKERRQ(ierr); 265 ierr = VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 266 ierr = VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 267 ierr = TSSolve(full->ts,full->UV);CHKERRQ(ierr); 268 ierr = VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 269 ierr = VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae) 274 { 275 PetscErrorCode ierr; 276 TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data; 277 278 PetscFunctionBegin; 279 ierr = TSSetFromOptions(full->ts);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae) 284 { 285 PetscErrorCode ierr; 286 TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data; 287 288 PetscFunctionBegin; 289 ierr = TSDestroy(&full->ts);CHKERRQ(ierr); 290 ierr = VecDestroy(&full->UV);CHKERRQ(ierr); 291 ierr = VecDestroy(&full->UF);CHKERRQ(ierr); 292 ierr = VecDestroy(&full->VF);CHKERRQ(ierr); 293 ierr = VecScatterDestroy(&full->scatterU);CHKERRQ(ierr); 294 ierr = VecScatterDestroy(&full->scatterV);CHKERRQ(ierr); 295 ierr = PetscFree(full);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae) 300 { 301 PetscErrorCode ierr; 302 TSDAESimple_Full *full; 303 Vec tsrhs; 304 PetscInt nU,nV,UVstart; 305 IS is; 306 307 PetscFunctionBegin; 308 ierr = PetscNew(&full);CHKERRQ(ierr); 309 tsdae->data = full; 310 311 tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full; 312 tsdae->solve = TSDAESimpleSolve_Full; 313 tsdae->destroy = TSDAESimpleDestroy_Full; 314 315 ierr = TSCreate(tsdae->comm,&full->ts);CHKERRQ(ierr); 316 ierr = TSSetProblemType(full->ts,TS_NONLINEAR);CHKERRQ(ierr); 317 ierr = TSSetType(full->ts,TSROSW);CHKERRQ(ierr); 318 ierr = TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 319 ierr = VecDuplicate(tsdae->U,&full->UF);CHKERRQ(ierr); 320 ierr = VecDuplicate(tsdae->V,&full->VF);CHKERRQ(ierr); 321 322 ierr = VecGetLocalSize(tsdae->U,&nU);CHKERRQ(ierr); 323 ierr = VecGetLocalSize(tsdae->V,&nV);CHKERRQ(ierr); 324 ierr = VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);CHKERRQ(ierr); 325 ierr = VecDuplicate(tsrhs,&full->UV);CHKERRQ(ierr); 326 327 ierr = VecGetOwnershipRange(tsrhs,&UVstart,NULL);CHKERRQ(ierr); 328 ierr = ISCreateStride(tsdae->comm,nU,UVstart,1,&is);CHKERRQ(ierr); 329 ierr = VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);CHKERRQ(ierr); 330 ierr = ISDestroy(&is);CHKERRQ(ierr); 331 ierr = ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);CHKERRQ(ierr); 332 ierr = VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);CHKERRQ(ierr); 333 ierr = ISDestroy(&is);CHKERRQ(ierr); 334 335 ierr = TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);CHKERRQ(ierr); 336 ierr = TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);CHKERRQ(ierr); 337 ierr = VecDestroy(&tsrhs);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 /* ----------------------------------------------------------------------------*/ 342 343 /* 344 Simple example: f(U,V) = U + V 345 346 */ 347 PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx) 348 { 349 PetscErrorCode ierr; 350 351 PetscFunctionBeginUser; 352 ierr = VecWAXPY(F,1.0,U,V);CHKERRQ(ierr); 353 PetscFunctionReturn(0); 354 } 355 356 /* 357 Simple example: F(U,V) = U - V 358 359 */ 360 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx) 361 { 362 PetscErrorCode ierr; 363 364 PetscFunctionBeginUser; 365 ierr = VecWAXPY(F,-1.0,V,U);CHKERRQ(ierr); 366 PetscFunctionReturn(0); 367 } 368 369 int main(int argc,char **argv) 370 { 371 PetscErrorCode ierr; 372 TSDAESimple tsdae; 373 Vec U,V,Usolution; 374 375 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 376 ierr = TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);CHKERRQ(ierr); 377 378 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);CHKERRQ(ierr); 379 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);CHKERRQ(ierr); 380 ierr = TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);CHKERRQ(ierr); 381 ierr = TSDAESimpleSetIFunction(tsdae,V,F,NULL);CHKERRQ(ierr); 382 383 ierr = VecDuplicate(U,&Usolution);CHKERRQ(ierr); 384 ierr = VecSet(Usolution,1.0);CHKERRQ(ierr); 385 386 /* ierr = TSDAESimpleSetUp_Full(tsdae);CHKERRQ(ierr); */ 387 ierr = TSDAESimpleSetUp_Reduced(tsdae);CHKERRQ(ierr); 388 389 ierr = TSDAESimpleSetFromOptions(tsdae);CHKERRQ(ierr); 390 ierr = TSDAESimpleSolve(tsdae,Usolution);CHKERRQ(ierr); 391 ierr = TSDAESimpleDestroy(&tsdae);CHKERRQ(ierr); 392 393 ierr = VecDestroy(&U);CHKERRQ(ierr); 394 ierr = VecDestroy(&Usolution);CHKERRQ(ierr); 395 ierr = VecDestroy(&V);CHKERRQ(ierr); 396 ierr = PetscFinalize(); 397 return ierr; 398 } 399 400