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