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 /* 12c4762a1bSJed Brown f(U,V) = U + V 13c4762a1bSJed Brown */ 14c4762a1bSJed Brown PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F) 15c4762a1bSJed Brown { 16c4762a1bSJed Brown PetscErrorCode ierr; 17c4762a1bSJed Brown 18c4762a1bSJed Brown PetscFunctionBeginUser; 19c4762a1bSJed Brown ierr = VecWAXPY(F,1.0,U,V);CHKERRQ(ierr); 20c4762a1bSJed Brown PetscFunctionReturn(0); 21c4762a1bSJed Brown } 22c4762a1bSJed Brown 23c4762a1bSJed Brown /* 24c4762a1bSJed Brown F(U,V) = U - V 25c4762a1bSJed Brown */ 26c4762a1bSJed Brown PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F) 27c4762a1bSJed Brown { 28c4762a1bSJed Brown PetscErrorCode ierr; 29c4762a1bSJed Brown 30c4762a1bSJed Brown PetscFunctionBeginUser; 31c4762a1bSJed Brown ierr = VecWAXPY(F,-1.0,V,U);CHKERRQ(ierr); 32c4762a1bSJed Brown PetscFunctionReturn(0); 33c4762a1bSJed Brown } 34c4762a1bSJed Brown 35c4762a1bSJed Brown typedef struct { 36c4762a1bSJed Brown PetscReal t; 37c4762a1bSJed Brown SNES snes; 38c4762a1bSJed Brown Vec U,V; 39c4762a1bSJed Brown PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec); 40c4762a1bSJed Brown PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec); 41c4762a1bSJed Brown } AppCtx; 42c4762a1bSJed Brown 43c4762a1bSJed Brown extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*); 44c4762a1bSJed Brown extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*); 45c4762a1bSJed Brown 46c4762a1bSJed Brown int main(int argc,char **argv) 47c4762a1bSJed Brown { 48c4762a1bSJed Brown PetscErrorCode ierr; 49c4762a1bSJed Brown AppCtx ctx; 50c4762a1bSJed Brown TS ts; 51c4762a1bSJed Brown Vec tsrhs,U; 52c4762a1bSJed Brown 53c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 54c4762a1bSJed Brown ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 55c4762a1bSJed Brown ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 56c4762a1bSJed Brown ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 57c4762a1bSJed Brown ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 58c4762a1bSJed Brown ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);CHKERRQ(ierr); 59c4762a1bSJed Brown ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);CHKERRQ(ierr); 60c4762a1bSJed Brown ierr = TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);CHKERRQ(ierr); 61*303a5415SBarry Smith ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 62c4762a1bSJed Brown ctx.f = f; 63c4762a1bSJed Brown 64c4762a1bSJed Brown ierr = SNESCreate(PETSC_COMM_WORLD,&ctx.snes);CHKERRQ(ierr); 65c4762a1bSJed Brown ierr = SNESSetFromOptions(ctx.snes);CHKERRQ(ierr); 66c4762a1bSJed Brown ierr = SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);CHKERRQ(ierr); 67c4762a1bSJed Brown ierr = SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);CHKERRQ(ierr); 68c4762a1bSJed Brown ctx.F = F; 69c4762a1bSJed Brown ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);CHKERRQ(ierr); 70c4762a1bSJed Brown 71c4762a1bSJed Brown ierr = VecSet(U,1.0);CHKERRQ(ierr); 72c4762a1bSJed Brown ierr = TSSolve(ts,U);CHKERRQ(ierr); 73c4762a1bSJed Brown 74c4762a1bSJed Brown ierr = VecDestroy(&ctx.V);CHKERRQ(ierr); 75c4762a1bSJed Brown ierr = VecDestroy(&tsrhs);CHKERRQ(ierr); 76c4762a1bSJed Brown ierr = VecDestroy(&U);CHKERRQ(ierr); 77c4762a1bSJed Brown ierr = SNESDestroy(&ctx.snes);CHKERRQ(ierr); 78c4762a1bSJed Brown ierr = TSDestroy(&ts);CHKERRQ(ierr); 79c4762a1bSJed Brown ierr = PetscFinalize(); 80c4762a1bSJed Brown return ierr; 81c4762a1bSJed Brown } 82c4762a1bSJed Brown 83c4762a1bSJed Brown /* 84c4762a1bSJed Brown Defines the RHS function that is passed to the time-integrator. 85c4762a1bSJed Brown 86c4762a1bSJed Brown Solves F(U,V) for V and then computes f(U,V) 87c4762a1bSJed Brown */ 88c4762a1bSJed Brown PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 89c4762a1bSJed Brown { 90c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 91c4762a1bSJed Brown PetscErrorCode ierr; 92c4762a1bSJed Brown 93c4762a1bSJed Brown PetscFunctionBeginUser; 94c4762a1bSJed Brown ctx->t = t; 95c4762a1bSJed Brown ctx->U = U; 96c4762a1bSJed Brown ierr = SNESSolve(ctx->snes,NULL,ctx->V);CHKERRQ(ierr); 97c4762a1bSJed Brown ierr = (*ctx->f)(t,U,ctx->V,F);CHKERRQ(ierr); 98c4762a1bSJed Brown PetscFunctionReturn(0); 99c4762a1bSJed Brown } 100c4762a1bSJed Brown 101c4762a1bSJed Brown /* 102c4762a1bSJed Brown Defines the nonlinear function that is passed to the nonlinear solver 103c4762a1bSJed Brown */ 104c4762a1bSJed Brown PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx) 105c4762a1bSJed Brown { 106c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 107c4762a1bSJed Brown PetscErrorCode ierr; 108c4762a1bSJed Brown 109c4762a1bSJed Brown PetscFunctionBeginUser; 110c4762a1bSJed Brown ierr = (*ctx->F)(ctx->t,ctx->U,V,F);CHKERRQ(ierr); 111c4762a1bSJed Brown PetscFunctionReturn(0); 112c4762a1bSJed Brown } 113c4762a1bSJed Brown 114*303a5415SBarry Smith /*TEST 115c4762a1bSJed Brown 116*303a5415SBarry Smith test: 117*303a5415SBarry Smith args: -ts_monitor -ts_view 118*303a5415SBarry Smith 119*303a5415SBarry Smith TEST*/ 120