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