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