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 TSROSW integrator on the entire DAE 10 */ 11 12 /* 13 f(U,V) = U + V 14 15 */ 16 PetscErrorCode f(PetscReal t, Vec UV, Vec F) { 17 const PetscScalar *u, *v; 18 PetscScalar *f; 19 PetscInt n, i; 20 21 PetscFunctionBeginUser; 22 PetscCall(VecGetLocalSize(UV, &n)); 23 n = n / 2; 24 PetscCall(VecGetArrayRead(UV, &u)); 25 v = u + n; 26 PetscCall(VecGetArrayWrite(F, &f)); 27 for (i = 0; i < n; i++) f[i] = u[i] + v[i]; 28 PetscCall(VecRestoreArrayRead(UV, &u)); 29 PetscCall(VecRestoreArrayWrite(F, &f)); 30 PetscFunctionReturn(0); 31 } 32 33 /* 34 F(U,V) = U - V 35 36 */ 37 PetscErrorCode F(PetscReal t, Vec UV, Vec F) { 38 const PetscScalar *u, *v; 39 PetscScalar *f; 40 PetscInt n, i; 41 42 PetscFunctionBeginUser; 43 PetscCall(VecGetLocalSize(UV, &n)); 44 n = n / 2; 45 PetscCall(VecGetArrayRead(UV, &u)); 46 v = u + n; 47 PetscCall(VecGetArrayWrite(F, &f)); 48 f = f + n; 49 for (i = 0; i < n; i++) f[i] = u[i] - v[i]; 50 f = f - n; 51 PetscCall(VecRestoreArrayRead(UV, &u)); 52 PetscCall(VecRestoreArrayWrite(F, &f)); 53 PetscFunctionReturn(0); 54 } 55 56 typedef struct { 57 PetscErrorCode (*f)(PetscReal, Vec, Vec); 58 PetscErrorCode (*F)(PetscReal, Vec, Vec); 59 } AppCtx; 60 61 extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *); 62 extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *); 63 64 int main(int argc, char **argv) { 65 AppCtx ctx; 66 TS ts; 67 Vec tsrhs, UV; 68 69 PetscFunctionBeginUser; 70 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 71 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 72 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 73 PetscCall(TSSetType(ts, TSROSW)); 74 PetscCall(TSSetFromOptions(ts)); 75 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs)); 76 PetscCall(VecDuplicate(tsrhs, &UV)); 77 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx)); 78 PetscCall(TSSetIFunction(ts, NULL, TSFunctionI, &ctx)); 79 PetscCall(TSSetMaxTime(ts, 1.0)); 80 ctx.f = f; 81 ctx.F = F; 82 83 PetscCall(VecSet(UV, 1.0)); 84 PetscCall(TSSolve(ts, UV)); 85 PetscCall(VecDestroy(&tsrhs)); 86 PetscCall(VecDestroy(&UV)); 87 PetscCall(TSDestroy(&ts)); 88 PetscCall(PetscFinalize()); 89 return 0; 90 } 91 92 /* 93 Defines the RHS function that is passed to the time-integrator. 94 */ 95 PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx) { 96 AppCtx *ctx = (AppCtx *)actx; 97 98 PetscFunctionBeginUser; 99 PetscCall(VecSet(F, 0.0)); 100 PetscCall((*ctx->f)(t, UV, F)); 101 PetscFunctionReturn(0); 102 } 103 104 /* 105 Defines the nonlinear function that is passed to the time-integrator 106 */ 107 PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx) { 108 AppCtx *ctx = (AppCtx *)actx; 109 110 PetscFunctionBeginUser; 111 PetscCall(VecCopy(UVdot, F)); 112 PetscCall((*ctx->F)(t, UV, F)); 113 PetscFunctionReturn(0); 114 } 115 116 /*TEST 117 118 test: 119 args: -ts_view 120 121 test: 122 suffix: 2 123 args: -snes_lag_jacobian 2 -ts_view 124 125 TEST*/ 126