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 9*3d5a8a6aSBarry Smith Same as ex6.c and ex7.c except calls the TSROSW integrator on the entire DAE 10c4762a1bSJed Brown */ 11c4762a1bSJed Brown 12c4762a1bSJed Brown 13c4762a1bSJed Brown /* 14c4762a1bSJed Brown f(U,V) = U + V 15c4762a1bSJed Brown 16c4762a1bSJed Brown */ 17c4762a1bSJed Brown PetscErrorCode f(PetscReal t,Vec UV,Vec F) 18c4762a1bSJed Brown { 19c4762a1bSJed Brown PetscErrorCode ierr; 20c4762a1bSJed Brown const PetscScalar *u,*v; 21c4762a1bSJed Brown PetscScalar *f; 22c4762a1bSJed Brown PetscInt n,i; 23c4762a1bSJed Brown 24c4762a1bSJed Brown PetscFunctionBeginUser; 25c4762a1bSJed Brown ierr = VecGetLocalSize(UV,&n);CHKERRQ(ierr); 26c4762a1bSJed Brown n = n/2; 27c4762a1bSJed Brown ierr = VecGetArrayRead(UV,&u);CHKERRQ(ierr); 28c4762a1bSJed Brown v = u + n; 29f7f07198SBarry Smith ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 30c4762a1bSJed Brown for (i=0; i<n; i++) f[i] = u[i] + v[i]; 31c4762a1bSJed Brown ierr = VecRestoreArrayRead(UV,&u);CHKERRQ(ierr); 32f7f07198SBarry Smith ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 33c4762a1bSJed Brown PetscFunctionReturn(0); 34c4762a1bSJed Brown } 35c4762a1bSJed Brown 36c4762a1bSJed Brown /* 37c4762a1bSJed Brown F(U,V) = U - V 38c4762a1bSJed Brown 39c4762a1bSJed Brown */ 40c4762a1bSJed Brown PetscErrorCode F(PetscReal t,Vec UV,Vec F) 41c4762a1bSJed Brown { 42c4762a1bSJed Brown PetscErrorCode ierr; 43c4762a1bSJed Brown const PetscScalar *u,*v; 44c4762a1bSJed Brown PetscScalar *f; 45c4762a1bSJed Brown PetscInt n,i; 46c4762a1bSJed Brown 47c4762a1bSJed Brown PetscFunctionBeginUser; 48c4762a1bSJed Brown ierr = VecGetLocalSize(UV,&n);CHKERRQ(ierr); 49c4762a1bSJed Brown n = n/2; 50c4762a1bSJed Brown ierr = VecGetArrayRead(UV,&u);CHKERRQ(ierr); 51c4762a1bSJed Brown v = u + n; 52f7f07198SBarry Smith ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 53c4762a1bSJed Brown f = f + n; 54c4762a1bSJed Brown for (i=0; i<n; i++) f[i] = u[i] - v[i]; 55c4762a1bSJed Brown f = f - n; 56c4762a1bSJed Brown ierr = VecRestoreArrayRead(UV,&u);CHKERRQ(ierr); 57f7f07198SBarry Smith ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 58c4762a1bSJed Brown PetscFunctionReturn(0); 59c4762a1bSJed Brown } 60c4762a1bSJed Brown 61c4762a1bSJed Brown typedef struct { 62c4762a1bSJed Brown PetscErrorCode (*f)(PetscReal,Vec,Vec); 63c4762a1bSJed Brown PetscErrorCode (*F)(PetscReal,Vec,Vec); 64c4762a1bSJed Brown } AppCtx; 65c4762a1bSJed Brown 66c4762a1bSJed Brown extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*); 67c4762a1bSJed Brown extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*); 68c4762a1bSJed Brown 69c4762a1bSJed Brown int main(int argc,char **argv) 70c4762a1bSJed Brown { 71c4762a1bSJed Brown PetscErrorCode ierr; 72c4762a1bSJed Brown AppCtx ctx; 73c4762a1bSJed Brown TS ts; 74c4762a1bSJed Brown Vec tsrhs,UV; 75c4762a1bSJed Brown 76c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 77c4762a1bSJed Brown ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 78c4762a1bSJed Brown ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 79c4762a1bSJed Brown ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 80c4762a1bSJed Brown ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 81c4762a1bSJed Brown ierr = VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);CHKERRQ(ierr); 82f7f07198SBarry Smith ierr = VecDuplicate(tsrhs,&UV);CHKERRQ(ierr); 83c4762a1bSJed Brown ierr = TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);CHKERRQ(ierr); 84c4762a1bSJed Brown ierr = TSSetIFunction(ts,NULL,TSFunctionI,&ctx);CHKERRQ(ierr); 85f7f07198SBarry Smith ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 86c4762a1bSJed Brown ctx.f = f; 87c4762a1bSJed Brown ctx.F = F; 88c4762a1bSJed Brown 89c4762a1bSJed Brown ierr = VecSet(UV,1.0);CHKERRQ(ierr); 90c4762a1bSJed Brown ierr = TSSolve(ts,UV);CHKERRQ(ierr); 91c4762a1bSJed Brown ierr = VecDestroy(&tsrhs);CHKERRQ(ierr); 92c4762a1bSJed Brown ierr = VecDestroy(&UV);CHKERRQ(ierr); 93c4762a1bSJed Brown ierr = TSDestroy(&ts);CHKERRQ(ierr); 94c4762a1bSJed Brown ierr = PetscFinalize(); 95c4762a1bSJed Brown return ierr; 96c4762a1bSJed Brown } 97c4762a1bSJed Brown 98c4762a1bSJed Brown /* 99c4762a1bSJed Brown Defines the RHS function that is passed to the time-integrator. 100c4762a1bSJed Brown */ 101c4762a1bSJed Brown PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 102c4762a1bSJed Brown { 103c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 104c4762a1bSJed Brown PetscErrorCode ierr; 105c4762a1bSJed Brown 106c4762a1bSJed Brown PetscFunctionBeginUser; 107c4762a1bSJed Brown ierr = VecSet(F,0.0);CHKERRQ(ierr); 108c4762a1bSJed Brown ierr = (*ctx->f)(t,UV,F);CHKERRQ(ierr); 109c4762a1bSJed Brown PetscFunctionReturn(0); 110c4762a1bSJed Brown } 111c4762a1bSJed Brown 112c4762a1bSJed Brown /* 113c4762a1bSJed Brown Defines the nonlinear function that is passed to the time-integrator 114c4762a1bSJed Brown */ 115c4762a1bSJed Brown PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 116c4762a1bSJed Brown { 117c4762a1bSJed Brown AppCtx *ctx = (AppCtx*)actx; 118c4762a1bSJed Brown PetscErrorCode ierr; 119c4762a1bSJed Brown 120c4762a1bSJed Brown PetscFunctionBeginUser; 121c4762a1bSJed Brown ierr = VecCopy(UVdot,F);CHKERRQ(ierr); 122c4762a1bSJed Brown ierr = (*ctx->F)(t,UV,F);CHKERRQ(ierr); 123c4762a1bSJed Brown PetscFunctionReturn(0); 124c4762a1bSJed Brown } 125f7f07198SBarry Smith 126f7f07198SBarry Smith /*TEST 127f7f07198SBarry Smith 128f7f07198SBarry Smith test: 129f7f07198SBarry Smith args: -ts_view 130f7f07198SBarry Smith 131f7f07198SBarry Smith test: 132f7f07198SBarry Smith suffix: 2 133f7f07198SBarry Smith args: -snes_lag_jacobian 2 -ts_view 134f7f07198SBarry Smith 135f7f07198SBarry Smith TEST*/ 136f7f07198SBarry Smith 137