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 PetscCall(PetscNew(tsdae)); 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 PetscCall(PetscObjectReference((PetscObject)U)); 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 PetscCall(PetscObjectReference((PetscObject)V)); 52 tsdae->Fctx = ctx; 53 PetscFunctionReturn(0); 54 } 55 56 PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae) 57 { 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 PetscCall((*(*tsdae)->destroy)(*tsdae)); 62 PetscCall(VecDestroy(&(*tsdae)->U)); 63 PetscCall(VecDestroy(&(*tsdae)->V)); 64 PetscCall(PetscFree(*tsdae)); 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution) 69 { 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 PetscCall((*tsdae->solve)(tsdae,Usolution)); 74 PetscFunctionReturn(0); 75 } 76 77 PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae) 78 { 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 PetscCall((*tsdae->setfromoptions)(PetscOptionsObject,tsdae)); 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 PetscCall(SNESSolve(red->snes,NULL,tsdae->V)); 115 PetscCall((*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx)); 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 PetscCall((*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx)); 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 PetscCall(TSSolve(red->ts,U)); 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 PetscCall(TSSetFromOptions(red->ts)); 151 PetscCall(SNESSetFromOptions(red->snes)); 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 PetscCall(TSDestroy(&red->ts)); 162 PetscCall(SNESDestroy(&red->snes)); 163 PetscCall(PetscFree(red)); 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 PetscCall(PetscNew(&red)); 175 tsdae->data = red; 176 177 tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced; 178 tsdae->solve = TSDAESimpleSolve_Reduced; 179 tsdae->destroy = TSDAESimpleDestroy_Reduced; 180 181 PetscCall(TSCreate(tsdae->comm,&red->ts)); 182 PetscCall(TSSetProblemType(red->ts,TS_NONLINEAR)); 183 PetscCall(TSSetType(red->ts,TSEULER)); 184 PetscCall(TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER)); 185 PetscCall(VecDuplicate(tsdae->U,&tsrhs)); 186 PetscCall(TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae)); 187 PetscCall(VecDestroy(&tsrhs)); 188 189 PetscCall(SNESCreate(tsdae->comm,&red->snes)); 190 PetscCall(SNESSetOptionsPrefix(red->snes,"tsdaesimple_")); 191 PetscCall(SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae)); 192 PetscCall(SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae)); 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 PetscCall(VecSet(F,0.0)); 223 PetscCall(VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE)); 224 PetscCall(VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE)); 225 PetscCall(VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE)); 226 PetscCall(VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE)); 227 PetscCall((*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx)); 228 PetscCall(VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 229 PetscCall(VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 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 PetscCall(VecCopy(UVdot,F)); 248 PetscCall(VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE)); 249 PetscCall(VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE)); 250 PetscCall(VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE)); 251 PetscCall(VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE)); 252 PetscCall((*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx)); 253 PetscCall(VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 254 PetscCall(VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 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 PetscCall(VecSet(full->UV,1.0)); 265 PetscCall(VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD)); 266 PetscCall(VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD)); 267 PetscCall(TSSolve(full->ts,full->UV)); 268 PetscCall(VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE)); 269 PetscCall(VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE)); 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 PetscCall(TSSetFromOptions(full->ts)); 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 PetscCall(TSDestroy(&full->ts)); 290 PetscCall(VecDestroy(&full->UV)); 291 PetscCall(VecDestroy(&full->UF)); 292 PetscCall(VecDestroy(&full->VF)); 293 PetscCall(VecScatterDestroy(&full->scatterU)); 294 PetscCall(VecScatterDestroy(&full->scatterV)); 295 PetscCall(PetscFree(full)); 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 PetscCall(PetscNew(&full)); 309 tsdae->data = full; 310 311 tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full; 312 tsdae->solve = TSDAESimpleSolve_Full; 313 tsdae->destroy = TSDAESimpleDestroy_Full; 314 315 PetscCall(TSCreate(tsdae->comm,&full->ts)); 316 PetscCall(TSSetProblemType(full->ts,TS_NONLINEAR)); 317 PetscCall(TSSetType(full->ts,TSROSW)); 318 PetscCall(TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER)); 319 PetscCall(VecDuplicate(tsdae->U,&full->UF)); 320 PetscCall(VecDuplicate(tsdae->V,&full->VF)); 321 322 PetscCall(VecGetLocalSize(tsdae->U,&nU)); 323 PetscCall(VecGetLocalSize(tsdae->V,&nV)); 324 PetscCall(VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs)); 325 PetscCall(VecDuplicate(tsrhs,&full->UV)); 326 327 PetscCall(VecGetOwnershipRange(tsrhs,&UVstart,NULL)); 328 PetscCall(ISCreateStride(tsdae->comm,nU,UVstart,1,&is)); 329 PetscCall(VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU)); 330 PetscCall(ISDestroy(&is)); 331 PetscCall(ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is)); 332 PetscCall(VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV)); 333 PetscCall(ISDestroy(&is)); 334 335 PetscCall(TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae)); 336 PetscCall(TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae)); 337 PetscCall(VecDestroy(&tsrhs)); 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 PetscCall(VecWAXPY(F,1.0,U,V)); 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 PetscCall(VecWAXPY(F,-1.0,V,U)); 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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 376 PetscCall(TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae)); 377 378 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U)); 379 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V)); 380 PetscCall(TSDAESimpleSetRHSFunction(tsdae,U,f,NULL)); 381 PetscCall(TSDAESimpleSetIFunction(tsdae,V,F,NULL)); 382 383 PetscCall(VecDuplicate(U,&Usolution)); 384 PetscCall(VecSet(Usolution,1.0)); 385 386 /* PetscCall(TSDAESimpleSetUp_Full(tsdae)); */ 387 PetscCall(TSDAESimpleSetUp_Reduced(tsdae)); 388 389 PetscCall(TSDAESimpleSetFromOptions(tsdae)); 390 PetscCall(TSDAESimpleSolve(tsdae,Usolution)); 391 PetscCall(TSDAESimpleDestroy(&tsdae)); 392 393 PetscCall(VecDestroy(&U)); 394 PetscCall(VecDestroy(&Usolution)); 395 PetscCall(VecDestroy(&V)); 396 PetscCall(PetscFinalize()); 397 return 0; 398 } 399