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