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 except the user provided functions take input values as a single vector instead of two vectors 10 */ 11 12 13 /* 14 f(U,V) = U + V 15 16 */ 17 PetscErrorCode f(PetscReal t,Vec UV,Vec F) 18 { 19 PetscErrorCode ierr; 20 const PetscScalar *u,*v; 21 PetscScalar *f; 22 PetscInt n,i; 23 24 PetscFunctionBeginUser; 25 ierr = VecGetLocalSize(UV,&n);CHKERRQ(ierr); 26 n = n/2; 27 ierr = VecGetArrayRead(UV,&u);CHKERRQ(ierr); 28 v = u + n; 29 ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 30 for (i=0; i<n; i++) f[i] = u[i] + v[i]; 31 ierr = VecRestoreArrayRead(UV,&u);CHKERRQ(ierr); 32 ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 /* 37 F(U,V) = U - V 38 */ 39 PetscErrorCode F(PetscReal t,Vec UV,Vec F) 40 { 41 PetscErrorCode ierr; 42 const PetscScalar *u,*v; 43 PetscScalar *f; 44 PetscInt n,i; 45 46 PetscFunctionBeginUser; 47 ierr = VecGetLocalSize(UV,&n);CHKERRQ(ierr); 48 n = n/2; 49 ierr = VecGetArrayRead(UV,&u);CHKERRQ(ierr); 50 v = u + n; 51 ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 52 for (i=0; i<n; i++) f[i] = u[i] - v[i]; 53 ierr = VecRestoreArrayRead(UV,&u);CHKERRQ(ierr); 54 ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 typedef struct { 59 PetscReal t; 60 SNES snes; 61 Vec UV,V; 62 VecScatter scatterU,scatterV; 63 PetscErrorCode (*f)(PetscReal,Vec,Vec); 64 PetscErrorCode (*F)(PetscReal,Vec,Vec); 65 } AppCtx; 66 67 extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*); 68 extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*); 69 70 int main(int argc,char **argv) 71 { 72 PetscErrorCode ierr; 73 AppCtx ctx; 74 TS ts; 75 Vec tsrhs,U; 76 IS is; 77 PetscInt i; 78 PetscMPIInt rank; 79 80 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 81 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 82 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 83 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 84 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 85 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 86 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);CHKERRQ(ierr); 87 ierr = VecDuplicate(tsrhs,&U);CHKERRQ(ierr); 88 ierr = TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);CHKERRQ(ierr); 89 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 90 ctx.f = f; 91 92 ierr = SNESCreate(PETSC_COMM_WORLD,&ctx.snes);CHKERRQ(ierr); 93 ierr = SNESSetFromOptions(ctx.snes);CHKERRQ(ierr); 94 ierr = SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);CHKERRQ(ierr); 95 ierr = SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);CHKERRQ(ierr); 96 ctx.F = F; 97 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);CHKERRQ(ierr); 98 99 /* Create scatters to move between separate U and V representation and UV representation of solution */ 100 ierr = VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);CHKERRQ(ierr); 101 i = 2*rank; 102 ierr = ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 103 ierr = VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);CHKERRQ(ierr); 104 ierr = ISDestroy(&is);CHKERRQ(ierr); 105 i = 2*rank + 1; 106 ierr = ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 107 ierr = VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);CHKERRQ(ierr); 108 ierr = ISDestroy(&is);CHKERRQ(ierr); 109 110 ierr = VecSet(U,1.0);CHKERRQ(ierr); 111 ierr = TSSolve(ts,U);CHKERRQ(ierr); 112 113 ierr = VecDestroy(&ctx.V);CHKERRQ(ierr); 114 ierr = VecDestroy(&ctx.UV);CHKERRQ(ierr); 115 ierr = VecScatterDestroy(&ctx.scatterU);CHKERRQ(ierr); 116 ierr = VecScatterDestroy(&ctx.scatterV);CHKERRQ(ierr); 117 ierr = VecDestroy(&tsrhs);CHKERRQ(ierr); 118 ierr = VecDestroy(&U);CHKERRQ(ierr); 119 ierr = SNESDestroy(&ctx.snes);CHKERRQ(ierr); 120 ierr = TSDestroy(&ts);CHKERRQ(ierr); 121 ierr = PetscFinalize(); 122 return ierr; 123 } 124 125 /* 126 Defines the RHS function that is passed to the time-integrator. 127 128 Solves F(U,V) for V and then computes f(U,V) 129 */ 130 PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 131 { 132 AppCtx *ctx = (AppCtx*)actx; 133 PetscErrorCode ierr; 134 135 PetscFunctionBeginUser; 136 ctx->t = t; 137 ierr = VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 138 ierr = VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 139 ierr = SNESSolve(ctx->snes,NULL,ctx->V);CHKERRQ(ierr); 140 ierr = VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 141 ierr = VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 142 ierr = (*ctx->f)(t,ctx->UV,F);CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 146 /* 147 Defines the nonlinear function that is passed to the nonlinear solver 148 149 */ 150 PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx) 151 { 152 AppCtx *ctx = (AppCtx*)actx; 153 PetscErrorCode ierr; 154 155 PetscFunctionBeginUser; 156 ierr = VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 157 ierr = VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 158 ierr = (*ctx->F)(ctx->t,ctx->UV,F);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 /*TEST 163 164 test: 165 args: -ts_monitor -ts_view 166 167 TEST*/ 168