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 PetscFunctionBeginUser; 19 PetscCall(VecWAXPY(F,1.0,U,V)); 20 PetscFunctionReturn(0); 21 } 22 23 /* 24 F(U,V) = U - V 25 26 */ 27 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F) 28 { 29 PetscFunctionBeginUser; 30 PetscCall(VecWAXPY(F,-1.0,V,U)); 31 PetscFunctionReturn(0); 32 } 33 34 typedef struct { 35 Vec U,V; 36 Vec UF,VF; 37 VecScatter scatterU,scatterV; 38 PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec); 39 PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec); 40 } AppCtx; 41 42 extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*); 43 extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*); 44 45 int main(int argc,char **argv) 46 { 47 AppCtx ctx; 48 TS ts; 49 Vec tsrhs,UV; 50 IS is; 51 PetscInt I; 52 PetscMPIInt rank; 53 54 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 55 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 56 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 57 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 58 PetscCall(TSSetType(ts,TSROSW)); 59 PetscCall(TSSetFromOptions(ts)); 60 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs)); 61 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV)); 62 PetscCall(TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx)); 63 PetscCall(TSSetIFunction(ts,NULL,TSFunctionI,&ctx)); 64 ctx.f = f; 65 ctx.F = F; 66 67 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U)); 68 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V)); 69 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF)); 70 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF)); 71 I = 2*rank; 72 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 73 PetscCall(VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU)); 74 PetscCall(ISDestroy(&is)); 75 I = 2*rank + 1; 76 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 77 PetscCall(VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV)); 78 PetscCall(ISDestroy(&is)); 79 80 PetscCall(VecSet(UV,1.0)); 81 PetscCall(TSSolve(ts,UV)); 82 PetscCall(VecDestroy(&tsrhs)); 83 PetscCall(VecDestroy(&UV)); 84 PetscCall(VecDestroy(&ctx.U)); 85 PetscCall(VecDestroy(&ctx.V)); 86 PetscCall(VecDestroy(&ctx.UF)); 87 PetscCall(VecDestroy(&ctx.VF)); 88 PetscCall(VecScatterDestroy(&ctx.scatterU)); 89 PetscCall(VecScatterDestroy(&ctx.scatterV)); 90 PetscCall(TSDestroy(&ts)); 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /* 96 Defines the RHS function that is passed to the time-integrator. 97 98 */ 99 PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 100 { 101 AppCtx *ctx = (AppCtx*)actx; 102 103 PetscFunctionBeginUser; 104 PetscCall(VecSet(F,0.0)); 105 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 106 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 107 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 108 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 109 PetscCall((*ctx->f)(t,ctx->U,ctx->V,ctx->UF)); 110 PetscCall(VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 111 PetscCall(VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 112 PetscFunctionReturn(0); 113 } 114 115 /* 116 Defines the nonlinear function that is passed to the time-integrator 117 118 */ 119 PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 120 { 121 AppCtx *ctx = (AppCtx*)actx; 122 123 PetscFunctionBeginUser; 124 PetscCall(VecCopy(UVdot,F)); 125 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 126 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 127 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 128 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 129 PetscCall((*ctx->F)(t,ctx->U,ctx->V,ctx->VF)); 130 PetscCall(VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 131 PetscCall(VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 132 PetscFunctionReturn(0); 133 } 134