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 = VecGetArray(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 = VecRestoreArray(F,&f);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 /* 37 F(U,V) = U - V 38 39 */ 40 PetscErrorCode F(PetscReal t,Vec UV,Vec F) 41 { 42 PetscErrorCode ierr; 43 const PetscScalar *u,*v; 44 PetscScalar *f; 45 PetscInt n,i; 46 47 PetscFunctionBeginUser; 48 ierr = VecGetLocalSize(UV,&n);CHKERRQ(ierr); 49 n = n/2; 50 ierr = VecGetArrayRead(UV,&u);CHKERRQ(ierr); 51 v = u + n; 52 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 53 for (i=0; i<n; i++) f[i] = u[i] - v[i]; 54 ierr = VecRestoreArrayRead(UV,&u);CHKERRQ(ierr); 55 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 56 PetscFunctionReturn(0); 57 } 58 59 typedef struct { 60 PetscReal t; 61 SNES snes; 62 Vec UV,V; 63 VecScatter scatterU,scatterV; 64 PetscErrorCode (*f)(PetscReal,Vec,Vec); 65 PetscErrorCode (*F)(PetscReal,Vec,Vec); 66 } AppCtx; 67 68 extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*); 69 extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*); 70 71 int main(int argc,char **argv) 72 { 73 PetscErrorCode ierr; 74 AppCtx ctx; 75 TS ts; 76 Vec tsrhs,U; 77 IS is; 78 PetscInt I; 79 PetscMPIInt rank; 80 81 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 82 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 83 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 84 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 85 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 86 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 87 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);CHKERRQ(ierr); 88 ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);CHKERRQ(ierr); 89 ierr = TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);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 */ 131 PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 132 { 133 AppCtx *ctx = (AppCtx*)actx; 134 PetscErrorCode ierr; 135 136 PetscFunctionBeginUser; 137 ctx->t = t; 138 ierr = VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 139 ierr = VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 140 ierr = SNESSolve(ctx->snes,NULL,ctx->V);CHKERRQ(ierr); 141 ierr = VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 142 ierr = VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 143 ierr = (*ctx->f)(t,ctx->UV,F);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 /* 148 Defines the nonlinear function that is passed to the nonlinear solver 149 150 */ 151 PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx) 152 { 153 AppCtx *ctx = (AppCtx*)actx; 154 PetscErrorCode ierr; 155 156 PetscFunctionBeginUser; 157 ierr = VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 158 ierr = VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 159 ierr = (*ctx->F)(ctx->t,ctx->UV,F);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 164