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