Actual source code: ex6.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";

      3:  #include <petscts.h>

      5: /*
      6:         \dot{U} = f(U,V)
      7:         F(U,V)  = 0
      8: */

     11: /*
     12:    f(U,V) = U + V

     14: */
     15: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
     16: {

     20:   VecWAXPY(F,1.0,U,V);
     21:   return(0);
     22: }

     24: /*
     25:    F(U,V) = U - V

     27: */
     28: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
     29: {

     33:   VecWAXPY(F,-1.0,V,U);
     34:   return(0);
     35: }

     37: typedef struct {
     38:   PetscReal      t;
     39:   SNES           snes;
     40:   Vec            U,V;
     41:   PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
     42:   PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
     43: } AppCtx;

     45: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
     46: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

     48: int main(int argc,char **argv)
     49: {
     51:   AppCtx         ctx;
     52:   TS             ts;
     53:   Vec            tsrhs,U;

     55:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     56:   TSCreate(PETSC_COMM_WORLD,&ts);
     57:   TSSetProblemType(ts,TS_NONLINEAR);
     58:   TSSetType(ts,TSEULER);
     59:   TSSetFromOptions(ts);
     60:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);
     61:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);
     62:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
     63:   ctx.f = f;

     65:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
     66:   SNESSetFromOptions(ctx.snes);
     67:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
     68:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
     69:   ctx.F = F;
     70:   VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);

     72:   VecSet(U,1.0);
     73:   TSSolve(ts,U);

     75:   VecDestroy(&ctx.V);
     76:   VecDestroy(&tsrhs);
     77:   VecDestroy(&U);
     78:   SNESDestroy(&ctx.snes);
     79:   TSDestroy(&ts);
     80:   PetscFinalize();
     81:   return ierr;
     82: }

     84: /*
     85:    Defines the RHS function that is passed to the time-integrator.

     87:    Solves F(U,V) for V and then computes f(U,V)

     89: */
     90: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
     91: {
     92:   AppCtx         *ctx = (AppCtx*)actx;

     96:   ctx->t = t;
     97:   ctx->U = U;
     98:   SNESSolve(ctx->snes,NULL,ctx->V);
     99:   (*ctx->f)(t,U,ctx->V,F);
    100:   return(0);
    101: }

    103: /*
    104:    Defines the nonlinear function that is passed to the nonlinear solver

    106: */
    107: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
    108: {
    109:   AppCtx         *ctx = (AppCtx*)actx;

    113:   (*ctx->F)(ctx->t,ctx->U,V,F);
    114:   return(0);
    115: }