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 PetscFunctionBeginUser; 55 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 56 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 57 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 58 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 59 PetscCall(TSSetType(ts,TSROSW)); 60 PetscCall(TSSetFromOptions(ts)); 61 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs)); 62 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV)); 63 PetscCall(TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx)); 64 PetscCall(TSSetIFunction(ts,NULL,TSFunctionI,&ctx)); 65 ctx.f = f; 66 ctx.F = F; 67 68 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U)); 69 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V)); 70 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF)); 71 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF)); 72 I = 2*rank; 73 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 74 PetscCall(VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU)); 75 PetscCall(ISDestroy(&is)); 76 I = 2*rank + 1; 77 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 78 PetscCall(VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV)); 79 PetscCall(ISDestroy(&is)); 80 81 PetscCall(VecSet(UV,1.0)); 82 PetscCall(TSSolve(ts,UV)); 83 PetscCall(VecDestroy(&tsrhs)); 84 PetscCall(VecDestroy(&UV)); 85 PetscCall(VecDestroy(&ctx.U)); 86 PetscCall(VecDestroy(&ctx.V)); 87 PetscCall(VecDestroy(&ctx.UF)); 88 PetscCall(VecDestroy(&ctx.VF)); 89 PetscCall(VecScatterDestroy(&ctx.scatterU)); 90 PetscCall(VecScatterDestroy(&ctx.scatterV)); 91 PetscCall(TSDestroy(&ts)); 92 PetscCall(PetscFinalize()); 93 return 0; 94 } 95 96 /* 97 Defines the RHS function that is passed to the time-integrator. 98 99 */ 100 PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 101 { 102 AppCtx *ctx = (AppCtx*)actx; 103 104 PetscFunctionBeginUser; 105 PetscCall(VecSet(F,0.0)); 106 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 107 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 108 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 109 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 110 PetscCall((*ctx->f)(t,ctx->U,ctx->V,ctx->UF)); 111 PetscCall(VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 112 PetscCall(VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 113 PetscFunctionReturn(0); 114 } 115 116 /* 117 Defines the nonlinear function that is passed to the time-integrator 118 119 */ 120 PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 121 { 122 AppCtx *ctx = (AppCtx*)actx; 123 124 PetscFunctionBeginUser; 125 PetscCall(VecCopy(UVdot,F)); 126 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 127 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 128 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 129 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 130 PetscCall((*ctx->F)(t,ctx->U,ctx->V,ctx->VF)); 131 PetscCall(VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 132 PetscCall(VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 133 PetscFunctionReturn(0); 134 } 135