Actual source code: ex7.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

      9:     Same as ex6.c except the user provided functions take input values as a single vector instead of two vectors
     10: */

     13: /*
     14:    f(U,V) = U + V

     16: */
     17: PetscErrorCode f(PetscReal t,Vec UV,Vec F)
     18: {
     19:   PetscErrorCode    ierr;
     20:   const PetscScalar *u,*v;
     21:   PetscScalar       *f;
     22:   PetscInt          n,i;

     25:   VecGetLocalSize(UV,&n);
     26:   n    = n/2;
     27:   VecGetArrayRead(UV,&u);
     28:   v    = u + n;
     29:   VecGetArray(F,&f);
     30:   for (i=0; i<n; i++) f[i] = u[i] + v[i];
     31:   VecRestoreArrayRead(UV,&u);
     32:   VecRestoreArray(F,&f);
     33:   return(0);
     34: }

     36: /*
     37:    F(U,V) = U - V

     39: */
     40: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
     41: {
     42:   PetscErrorCode    ierr;
     43:   const PetscScalar *u,*v;
     44:   PetscScalar       *f;
     45:   PetscInt          n,i;

     48:   VecGetLocalSize(UV,&n);
     49:   n    = n/2;
     50:   VecGetArrayRead(UV,&u);
     51:   v    = u + n;
     52:   VecGetArray(F,&f);
     53:   for (i=0; i<n; i++) f[i] = u[i] - v[i];
     54:   VecRestoreArrayRead(UV,&u);
     55:   VecRestoreArray(F,&f);
     56:   return(0);
     57: }

     59: typedef struct {
     60:   PetscReal      t;
     61:   SNES           snes;
     62:   Vec            UV,V;
     63:   VecScatter     scatterU,scatterV;
     64:   PetscErrorCode (*f)(PetscReal,Vec,Vec);
     65:   PetscErrorCode (*F)(PetscReal,Vec,Vec);
     66: } AppCtx;

     68: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
     69: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);

     71: int main(int argc,char **argv)
     72: {
     74:   AppCtx         ctx;
     75:   TS             ts;
     76:   Vec            tsrhs,U;
     77:   IS             is;
     78:   PetscInt       I;
     79:   PetscMPIInt    rank;

     81:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     82:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     83:   TSCreate(PETSC_COMM_WORLD,&ts);
     84:   TSSetProblemType(ts,TS_NONLINEAR);
     85:   TSSetType(ts,TSEULER);
     86:   TSSetFromOptions(ts);
     87:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);
     88:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
     89:   TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
     90:   ctx.f = f;

     92:   SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
     93:   SNESSetFromOptions(ctx.snes);
     94:   SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
     95:   SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
     96:   ctx.F = F;
     97:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);

     99:   /* Create scatters to move between separate U and V representation and UV representation of solution */
    100:   VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);
    101:   I    = 2*rank;
    102:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
    103:   VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);
    104:   ISDestroy(&is);
    105:   I    = 2*rank + 1;
    106:   ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
    107:   VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);
    108:   ISDestroy(&is);

    110:   VecSet(U,1.0);
    111:   TSSolve(ts,U);

    113:   VecDestroy(&ctx.V);
    114:   VecDestroy(&ctx.UV);
    115:   VecScatterDestroy(&ctx.scatterU);
    116:   VecScatterDestroy(&ctx.scatterV);
    117:   VecDestroy(&tsrhs);
    118:   VecDestroy(&U);
    119:   SNESDestroy(&ctx.snes);
    120:   TSDestroy(&ts);
    121:   PetscFinalize();
    122:   return ierr;
    123: }

    125: /*
    126:    Defines the RHS function that is passed to the time-integrator.

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

    130: */
    131: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
    132: {
    133:   AppCtx         *ctx = (AppCtx*)actx;

    137:   ctx->t = t;
    138:   VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    139:   VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    140:   SNESSolve(ctx->snes,NULL,ctx->V);
    141:   VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    142:   VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    143:   (*ctx->f)(t,ctx->UV,F);
    144:   return(0);
    145: }

    147: /*
    148:    Defines the nonlinear function that is passed to the nonlinear solver

    150: */
    151: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
    152: {
    153:   AppCtx         *ctx = (AppCtx*)actx;

    157:   VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    158:   VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
    159:   (*ctx->F)(ctx->t,ctx->UV,F);
    160:   return(0);
    161: }