1c4762a1bSJed Brown static char help[] = "Solves DAE with integrator only on non-algebraic terms \n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscts.h> 4c4762a1bSJed Brown 5c4762a1bSJed Brown /* 6c4762a1bSJed Brown \dot{U} = f(U,V) 7c4762a1bSJed Brown F(U,V) = 0 8c4762a1bSJed Brown */ 9c4762a1bSJed Brown 10c4762a1bSJed Brown /* 11c4762a1bSJed Brown f(U,V) = U + V 12c4762a1bSJed Brown */ 13d71ae5a4SJacob Faibussowitsch PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F) 14d71ae5a4SJacob Faibussowitsch { 15c4762a1bSJed Brown PetscFunctionBeginUser; 169566063dSJacob Faibussowitsch PetscCall(VecWAXPY(F, 1.0, U, V)); 173ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 18c4762a1bSJed Brown } 19c4762a1bSJed Brown 20c4762a1bSJed Brown /* 21c4762a1bSJed Brown F(U,V) = U - V 22c4762a1bSJed Brown */ 23d71ae5a4SJacob Faibussowitsch PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F) 24d71ae5a4SJacob Faibussowitsch { 25c4762a1bSJed Brown PetscFunctionBeginUser; 269566063dSJacob Faibussowitsch PetscCall(VecWAXPY(F, -1.0, V, U)); 273ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 28c4762a1bSJed Brown } 29c4762a1bSJed Brown 30c4762a1bSJed Brown typedef struct { 31c4762a1bSJed Brown PetscReal t; 32c4762a1bSJed Brown SNES snes; 33c4762a1bSJed Brown Vec U, V; 34c4762a1bSJed Brown PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec); 35c4762a1bSJed Brown PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec); 36c4762a1bSJed Brown } AppCtx; 37c4762a1bSJed Brown 38c4762a1bSJed Brown extern PetscErrorCode TSFunction(TS, PetscReal, Vec, Vec, void *); 39c4762a1bSJed Brown extern PetscErrorCode SNESFunction(SNES, Vec, Vec, void *); 40c4762a1bSJed Brown 41d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 42d71ae5a4SJacob Faibussowitsch { 43c4762a1bSJed Brown AppCtx ctx; 44c4762a1bSJed Brown TS ts; 45c4762a1bSJed Brown Vec tsrhs, U; 46c4762a1bSJed Brown 47327415f7SBarry Smith PetscFunctionBeginUser; 48*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 499566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 509566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 519566063dSJacob Faibussowitsch PetscCall(TSSetType(ts, TSEULER)); 529566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts)); 539566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &tsrhs)); 549566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &U)); 559566063dSJacob Faibussowitsch PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx)); 569566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts, 1.0)); 57c4762a1bSJed Brown ctx.f = f; 58c4762a1bSJed Brown 599566063dSJacob Faibussowitsch PetscCall(SNESCreate(PETSC_COMM_WORLD, &ctx.snes)); 609566063dSJacob Faibussowitsch PetscCall(SNESSetFromOptions(ctx.snes)); 619566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx)); 629566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx)); 63c4762a1bSJed Brown ctx.F = F; 649566063dSJacob Faibussowitsch PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &ctx.V)); 65c4762a1bSJed Brown 669566063dSJacob Faibussowitsch PetscCall(VecSet(U, 1.0)); 679566063dSJacob Faibussowitsch PetscCall(TSSolve(ts, U)); 68c4762a1bSJed Brown 699566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ctx.V)); 709566063dSJacob Faibussowitsch PetscCall(VecDestroy(&tsrhs)); 719566063dSJacob Faibussowitsch PetscCall(VecDestroy(&U)); 729566063dSJacob Faibussowitsch PetscCall(SNESDestroy(&ctx.snes)); 739566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts)); 749566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 75b122ec5aSJacob Faibussowitsch return 0; 76c4762a1bSJed Brown } 77c4762a1bSJed Brown 78c4762a1bSJed Brown /* 79c4762a1bSJed Brown Defines the RHS function that is passed to the time-integrator. 80c4762a1bSJed Brown 81c4762a1bSJed Brown Solves F(U,V) for V and then computes f(U,V) 82c4762a1bSJed Brown */ 83d71ae5a4SJacob Faibussowitsch PetscErrorCode TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx) 84d71ae5a4SJacob Faibussowitsch { 85c4762a1bSJed Brown AppCtx *ctx = (AppCtx *)actx; 86c4762a1bSJed Brown 87c4762a1bSJed Brown PetscFunctionBeginUser; 88c4762a1bSJed Brown ctx->t = t; 89c4762a1bSJed Brown ctx->U = U; 909566063dSJacob Faibussowitsch PetscCall(SNESSolve(ctx->snes, NULL, ctx->V)); 919566063dSJacob Faibussowitsch PetscCall((*ctx->f)(t, U, ctx->V, F)); 923ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 93c4762a1bSJed Brown } 94c4762a1bSJed Brown 95c4762a1bSJed Brown /* 96c4762a1bSJed Brown Defines the nonlinear function that is passed to the nonlinear solver 97c4762a1bSJed Brown */ 98d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESFunction(SNES snes, Vec V, Vec F, void *actx) 99d71ae5a4SJacob Faibussowitsch { 100c4762a1bSJed Brown AppCtx *ctx = (AppCtx *)actx; 101c4762a1bSJed Brown 102c4762a1bSJed Brown PetscFunctionBeginUser; 1039566063dSJacob Faibussowitsch PetscCall((*ctx->F)(ctx->t, ctx->U, V, F)); 1043ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 105c4762a1bSJed Brown } 106c4762a1bSJed Brown 107303a5415SBarry Smith /*TEST 108c4762a1bSJed Brown 109303a5415SBarry Smith test: 110303a5415SBarry Smith args: -ts_monitor -ts_view 111303a5415SBarry Smith 112303a5415SBarry Smith TEST*/ 113