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 PetscCall(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 PetscCall(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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 60 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 61 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 62 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 63 PetscCall(TSSetType(ts,TSROSW)); 64 PetscCall(TSSetFromOptions(ts)); 65 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs)); 66 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV)); 67 PetscCall(TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx)); 68 PetscCall(TSSetIFunction(ts,NULL,TSFunctionI,&ctx)); 69 ctx.f = f; 70 ctx.F = F; 71 72 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U)); 73 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V)); 74 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF)); 75 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF)); 76 I = 2*rank; 77 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 78 PetscCall(VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU)); 79 PetscCall(ISDestroy(&is)); 80 I = 2*rank + 1; 81 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is)); 82 PetscCall(VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV)); 83 PetscCall(ISDestroy(&is)); 84 85 PetscCall(VecSet(UV,1.0)); 86 PetscCall(TSSolve(ts,UV)); 87 PetscCall(VecDestroy(&tsrhs)); 88 PetscCall(VecDestroy(&UV)); 89 PetscCall(VecDestroy(&ctx.U)); 90 PetscCall(VecDestroy(&ctx.V)); 91 PetscCall(VecDestroy(&ctx.UF)); 92 PetscCall(VecDestroy(&ctx.VF)); 93 PetscCall(VecScatterDestroy(&ctx.scatterU)); 94 PetscCall(VecScatterDestroy(&ctx.scatterV)); 95 PetscCall(TSDestroy(&ts)); 96 PetscCall(PetscFinalize()); 97 return 0; 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 PetscCall(VecSet(F,0.0)); 111 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 112 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 113 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 114 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 115 PetscCall((*ctx->f)(t,ctx->U,ctx->V,ctx->UF)); 116 PetscCall(VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD)); 117 PetscCall(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 PetscCall(VecCopy(UVdot,F)); 132 PetscCall(VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 133 PetscCall(VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE)); 134 PetscCall(VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 135 PetscCall(VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE)); 136 PetscCall((*ctx->F)(t,ctx->U,ctx->V,ctx->VF)); 137 PetscCall(VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 138 PetscCall(VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD)); 139 PetscFunctionReturn(0); 140 } 141