1c4762a1bSJed Brown static char help[] = "Solves DAE with integrator only on non-algebraic terms \n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscts.h> 4c4762a1bSJed Brown 5c4762a1bSJed Brown /* 6c4762a1bSJed Brown \dot{U} = f(U,V) 7c4762a1bSJed Brown F(U,V) = 0 8c4762a1bSJed Brown */ 9c4762a1bSJed Brown 10c4762a1bSJed Brown /* 11c4762a1bSJed Brown f(U,V) = U + V 12c4762a1bSJed Brown */ 13c4762a1bSJed Brown PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F) 14c4762a1bSJed Brown { 15c4762a1bSJed Brown PetscFunctionBeginUser; 165f80ce2aSJacob Faibussowitsch CHKERRQ(VecWAXPY(F,1.0,U,V)); 17c4762a1bSJed Brown PetscFunctionReturn(0); 18c4762a1bSJed Brown } 19c4762a1bSJed Brown 20c4762a1bSJed Brown /* 21c4762a1bSJed Brown F(U,V) = U - V 22c4762a1bSJed Brown */ 23c4762a1bSJed Brown PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F) 24c4762a1bSJed Brown { 25c4762a1bSJed Brown PetscFunctionBeginUser; 265f80ce2aSJacob Faibussowitsch CHKERRQ(VecWAXPY(F,-1.0,V,U)); 27c4762a1bSJed Brown PetscFunctionReturn(0); 28c4762a1bSJed Brown } 29c4762a1bSJed Brown 30c4762a1bSJed Brown typedef struct { 31c4762a1bSJed Brown PetscReal t; 32c4762a1bSJed Brown SNES snes; 33c4762a1bSJed Brown Vec U,V; 34c4762a1bSJed Brown PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec); 35c4762a1bSJed Brown PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec); 36c4762a1bSJed Brown } AppCtx; 37c4762a1bSJed Brown 38c4762a1bSJed Brown extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*); 39c4762a1bSJed Brown extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*); 40c4762a1bSJed Brown 41c4762a1bSJed Brown int main(int argc,char **argv) 42c4762a1bSJed Brown { 43c4762a1bSJed Brown AppCtx ctx; 44c4762a1bSJed Brown TS ts; 45c4762a1bSJed Brown Vec tsrhs,U; 46c4762a1bSJed Brown 47*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 485f80ce2aSJacob Faibussowitsch CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 495f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 505f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetType(ts,TSEULER)); 515f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetFromOptions(ts)); 525f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs)); 535f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U)); 545f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx)); 555f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxTime(ts,1.0)); 56c4762a1bSJed Brown ctx.f = f; 57c4762a1bSJed Brown 585f80ce2aSJacob Faibussowitsch CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&ctx.snes)); 595f80ce2aSJacob Faibussowitsch CHKERRQ(SNESSetFromOptions(ctx.snes)); 605f80ce2aSJacob Faibussowitsch CHKERRQ(SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx)); 615f80ce2aSJacob Faibussowitsch CHKERRQ(SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx)); 62c4762a1bSJed Brown ctx.F = F; 635f80ce2aSJacob Faibussowitsch CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V)); 64c4762a1bSJed Brown 655f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(U,1.0)); 665f80ce2aSJacob Faibussowitsch CHKERRQ(TSSolve(ts,U)); 67c4762a1bSJed Brown 685f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&ctx.V)); 695f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&tsrhs)); 705f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&U)); 715f80ce2aSJacob Faibussowitsch CHKERRQ(SNESDestroy(&ctx.snes)); 725f80ce2aSJacob Faibussowitsch CHKERRQ(TSDestroy(&ts)); 73*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 74*b122ec5aSJacob Faibussowitsch return 0; 75c4762a1bSJed Brown } 76c4762a1bSJed Brown 77c4762a1bSJed Brown /* 78c4762a1bSJed Brown Defines the RHS function that is passed to the time-integrator. 79c4762a1bSJed Brown 80c4762a1bSJed Brown Solves F(U,V) for V and then computes f(U,V) 81c4762a1bSJed Brown */ 82c4762a1bSJed Brown PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 83c4762a1bSJed Brown { 84c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 85c4762a1bSJed Brown 86c4762a1bSJed Brown PetscFunctionBeginUser; 87c4762a1bSJed Brown ctx->t = t; 88c4762a1bSJed Brown ctx->U = U; 895f80ce2aSJacob Faibussowitsch CHKERRQ(SNESSolve(ctx->snes,NULL,ctx->V)); 905f80ce2aSJacob Faibussowitsch CHKERRQ((*ctx->f)(t,U,ctx->V,F)); 91c4762a1bSJed Brown PetscFunctionReturn(0); 92c4762a1bSJed Brown } 93c4762a1bSJed Brown 94c4762a1bSJed Brown /* 95c4762a1bSJed Brown Defines the nonlinear function that is passed to the nonlinear solver 96c4762a1bSJed Brown */ 97c4762a1bSJed Brown PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx) 98c4762a1bSJed Brown { 99c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 100c4762a1bSJed Brown 101c4762a1bSJed Brown PetscFunctionBeginUser; 1025f80ce2aSJacob Faibussowitsch CHKERRQ((*ctx->F)(ctx->t,ctx->U,V,F)); 103c4762a1bSJed Brown PetscFunctionReturn(0); 104c4762a1bSJed Brown } 105c4762a1bSJed Brown 106303a5415SBarry Smith /*TEST 107c4762a1bSJed Brown 108303a5415SBarry Smith test: 109303a5415SBarry Smith args: -ts_monitor -ts_view 110303a5415SBarry Smith 111303a5415SBarry Smith TEST*/ 112