1 static char help[] = "Solves DAE with integrator only on non-algebraic terms \n"; 2 3 #include <petscts.h> 4 5 /* 6 \dot{U} = f(U,V) 7 F(U,V) = 0 8 9 Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE 10 */ 11 12 /* 13 f(U,V) = U + V 14 15 */ 16 PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F) 17 { 18 PetscErrorCode ierr; 19 20 PetscFunctionBeginUser; 21 CHKERRQ(VecWAXPY(F,1.0,U,V)); 22 PetscFunctionReturn(0); 23 } 24 25 /* 26 F(U,V) = U - V 27 28 */ 29 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F) 30 { 31 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 34 CHKERRQ(VecWAXPY(F,-1.0,V,U)); 35 PetscFunctionReturn(0); 36 } 37 38 typedef struct { 39 Vec U,V; 40 Vec UF,VF; 41 VecScatter scatterU,scatterV; 42 PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec); 43 PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec); 44 } AppCtx; 45 46 extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*); 47 extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*); 48 49 int main(int argc,char **argv) 50 { 51 PetscErrorCode ierr; 52 AppCtx ctx; 53 TS ts; 54 Vec tsrhs,UV; 55 IS is; 56 PetscInt I; 57 PetscMPIInt rank; 58 59 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 60 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 61 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 62 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 63 CHKERRQ(TSSetType(ts,TSROSW)); 64 CHKERRQ(TSSetFromOptions(ts)); 65 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs)); 66 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV)); 67 CHKERRQ(TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx)); 68 CHKERRQ(TSSetIFunction(ts,NULL,TSFunctionI,&ctx)); 69 ctx.f = f; 70 ctx.F = F; 71 72 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U)); 73 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V)); 74 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF)); 75 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF)); 76 I = 2*rank; 77 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 78 CHKERRQ(VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU)); 79 CHKERRQ(ISDestroy(&is)); 80 I = 2*rank + 1; 81 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 82 CHKERRQ(VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV)); 83 CHKERRQ(ISDestroy(&is)); 84 85 CHKERRQ(VecSet(UV,1.0)); 86 CHKERRQ(TSSolve(ts,UV)); 87 CHKERRQ(VecDestroy(&tsrhs)); 88 CHKERRQ(VecDestroy(&UV)); 89 CHKERRQ(VecDestroy(&ctx.U)); 90 CHKERRQ(VecDestroy(&ctx.V)); 91 CHKERRQ(VecDestroy(&ctx.UF)); 92 CHKERRQ(VecDestroy(&ctx.VF)); 93 CHKERRQ(VecScatterDestroy(&ctx.scatterU)); 94 CHKERRQ(VecScatterDestroy(&ctx.scatterV)); 95 CHKERRQ(TSDestroy(&ts)); 96 ierr = PetscFinalize(); 97 return ierr; 98 } 99 100 /* 101 Defines the RHS function that is passed to the time-integrator. 102 103 */ 104 PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 105 { 106 AppCtx *ctx = (AppCtx*)actx; 107 PetscErrorCode ierr; 108 109 PetscFunctionBeginUser; 110 CHKERRQ(VecSet(F,0.0)); 111 CHKERRQ(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 112 CHKERRQ(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 113 CHKERRQ(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 114 CHKERRQ(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 115 CHKERRQ((*ctx->f)(t,ctx->U,ctx->V,ctx->UF)); 116 CHKERRQ(VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 117 CHKERRQ(VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 118 PetscFunctionReturn(0); 119 } 120 121 /* 122 Defines the nonlinear function that is passed to the time-integrator 123 124 */ 125 PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 126 { 127 AppCtx *ctx = (AppCtx*)actx; 128 PetscErrorCode ierr; 129 130 PetscFunctionBeginUser; 131 CHKERRQ(VecCopy(UVdot,F)); 132 CHKERRQ(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 133 CHKERRQ(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 134 CHKERRQ(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 135 CHKERRQ(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 136 CHKERRQ((*ctx->F)(t,ctx->U,ctx->V,ctx->VF)); 137 CHKERRQ(VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 138 CHKERRQ(VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 139 PetscFunctionReturn(0); 140 } 141