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