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