Actual source code: beuler.c

      1: #define PETSCTS_DLL

      3: /*
      4:        Code for Timestepping with implicit backwards Euler.
      5: */
     6:  #include src/ts/tsimpl.h

      8: typedef struct {
      9:   Vec  update;      /* work vector where new solution is formed */
     10:   Vec  func;        /* work vector where F(t[i],u[i]) is stored */
     11:   Vec  rhs;         /* work vector for RHS; vec_sol/dt */
     12: } TS_BEuler;

     14: /*------------------------------------------------------------------------------*/

     16: /*
     17:     Version for linear PDE where RHS does not depend on time. Has built a
     18:   single matrix that is to be used for all timesteps.
     19: */
     22: static PetscErrorCode TSStep_BEuler_Linear_Constant_Matrix(TS ts,PetscInt *steps,PetscReal *ptime)
     23: {
     24:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
     25:   Vec            sol = ts->vec_sol,update = beuler->update;
     26:   Vec            rhs = beuler->rhs;
     28:   PetscInt       i,max_steps = ts->max_steps,its;
     29:   PetscScalar    mdt = 1.0/ts->time_step;
     30:   KSP            ksp;

     33:   TSGetKSP(ts,&ksp);
     34:   *steps = -ts->steps;
     35:   TSMonitor(ts,ts->steps,ts->ptime,sol);

     37:   /* set initial guess to be previous solution */
     38:   VecCopy(sol,update);

     40:   for (i=0; i<max_steps; i++) {
     41:     VecCopy(sol,rhs);
     42:     VecScale(rhs,mdt);
     43:     /* apply user-provided boundary conditions (only needed if they are time dependent) */
     44:     TSComputeRHSBoundaryConditions(ts,ts->ptime,rhs);

     46:     ts->ptime += ts->time_step;
     47:     if (ts->ptime > ts->max_time) break;
     48:     KSPSolve(ts->ksp,rhs,update);
     49:     KSPGetIterationNumber(ksp,&its);
     50:     ts->linear_its += its;
     51:     VecCopy(update,sol);
     52:     ts->steps++;
     53:     TSMonitor(ts,ts->steps,ts->ptime,sol);
     54:   }

     56:   *steps += ts->steps;
     57:   *ptime  = ts->ptime;
     58:   return(0);
     59: }

     61: /*
     62:       Version where matrix depends on time
     63: */
     66: static PetscErrorCode TSStep_BEuler_Linear_Variable_Matrix(TS ts,PetscInt *steps,PetscReal *ptime)
     67: {
     68:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
     69:   Vec            sol = ts->vec_sol,update = beuler->update,rhs = beuler->rhs;
     71:   PetscInt       i,max_steps = ts->max_steps,its;
     72:   PetscScalar    mdt = 1.0/ts->time_step;
     73:   MatStructure   str;
     74:   KSP            ksp;

     77:   TSGetKSP(ts,&ksp);
     78:   *steps = -ts->steps;
     79:   TSMonitor(ts,ts->steps,ts->ptime,sol);

     81:   /* set initial guess to be previous solution */
     82:   VecCopy(sol,update);

     84:   for (i=0; i<max_steps; i++) {
     85:     VecCopy(sol,rhs);
     86:     VecScale(rhs,mdt);
     87:     /* apply user-provided boundary conditions (only needed if they are time dependent) */
     88:     TSComputeRHSBoundaryConditions(ts,ts->ptime,rhs);

     90:     ts->ptime += ts->time_step;
     91:     if (ts->ptime > ts->max_time) break;
     92:     /*
     93:         evaluate matrix function
     94:     */
     95:     (*ts->ops->rhsmatrix)(ts,ts->ptime,&ts->A,&ts->B,&str,ts->jacP);
     96:     TSScaleShiftMatrices(ts,ts->A,ts->B,str);
     97:     KSPSetOperators(ts->ksp,ts->A,ts->B,str);
     98:     KSPSolve(ts->ksp,rhs,update);
     99:     KSPGetIterationNumber(ksp,&its);
    100:     ts->linear_its += its;
    101:     VecCopy(update,sol);
    102:     ts->steps++;
    103:     TSMonitor(ts,ts->steps,ts->ptime,sol);
    104:   }

    106:   *steps += ts->steps;
    107:   *ptime  = ts->ptime;
    108:   return(0);
    109: }
    110: /*
    111:     Version for nonlinear PDE.
    112: */
    115: static PetscErrorCode TSStep_BEuler_Nonlinear(TS ts,PetscInt *steps,PetscReal *ptime)
    116: {
    117:   Vec            sol = ts->vec_sol;
    119:   PetscInt       i,max_steps = ts->max_steps,its,lits;
    120:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
    121:
    123:   *steps = -ts->steps;
    124:   TSMonitor(ts,ts->steps,ts->ptime,sol);

    126:   for (i=0; i<max_steps; i++) {
    127:     ts->ptime += ts->time_step;
    128:     if (ts->ptime > ts->max_time) break;
    129:     VecCopy(sol,beuler->update);
    130:     SNESSolve(ts->snes,PETSC_NULL,beuler->update);
    131:     SNESGetNumberLinearIterations(ts->snes,&lits);
    132:     SNESGetIterationNumber(ts->snes,&its);
    133:     ts->nonlinear_its += its; ts->linear_its += lits;
    134:     VecCopy(beuler->update,sol);
    135:     ts->steps++;
    136:     TSMonitor(ts,ts->steps,ts->ptime,sol);
    137:   }

    139:   *steps += ts->steps;
    140:   *ptime  = ts->ptime;
    141:   return(0);
    142: }

    144: /*------------------------------------------------------------*/
    147: static PetscErrorCode TSDestroy_BEuler(TS ts)
    148: {
    149:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

    153:   if (beuler->update) {VecDestroy(beuler->update);}
    154:   if (beuler->func) {VecDestroy(beuler->func);}
    155:   if (beuler->rhs) {VecDestroy(beuler->rhs);}
    156:   PetscFree(beuler);
    157:   return(0);
    158: }

    162: /*
    163:     This defines the nonlinear equation that is to be solved with SNES

    165:               U^{n+1} - dt*F(U^{n+1}) - U^{n}
    166: */
    169: PetscErrorCode TSBEulerFunction(SNES snes,Vec x,Vec y,void *ctx)
    170: {
    171:   TS             ts = (TS) ctx;
    172:   PetscScalar    mdt = 1.0/ts->time_step,*unp1,*un,*Funp1;
    174:   PetscInt       i,n;

    177:   /* apply user-provided function */
    178:   TSComputeRHSFunction(ts,ts->ptime,x,y);
    179:   /* (u^{n+1} - U^{n})/dt - F(u^{n+1}) */
    180:   VecGetArray(ts->vec_sol,&un);
    181:   VecGetArray(x,&unp1);
    182:   VecGetArray(y,&Funp1);
    183:   VecGetLocalSize(x,&n);

    185:   for (i=0; i<n; i++) {
    186:     Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i];
    187:   }
    188:   VecRestoreArray(ts->vec_sol,&un);
    189:   VecRestoreArray(x,&unp1);
    190:   VecRestoreArray(y,&Funp1);
    191:   return(0);
    192: }

    194: /*
    195:    This constructs the Jacobian needed for SNES

    197:              J = I/dt - J_{F}   where J_{F} is the given Jacobian of F.
    198: */
    201: PetscErrorCode TSBEulerJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
    202: {
    203:   TS             ts = (TS) ctx;

    207:   /* construct user's Jacobian */
    208:   TSComputeRHSJacobian(ts,ts->ptime,x,AA,BB,str);

    210:   /* shift and scale Jacobian */
    211:   /* this test is a undesirable hack, we assume that if it is MATMFFD then it is
    212:      obtained from -snes_mf_operator and there is computed directly from the
    213:      FormFunction() SNES is given and therefor does not need to be shifted/scaled
    214:      BUT maybe it could be MATMFFD and does require shift in some other case? */
    215:   TSScaleShiftMatrices(ts,*AA,*BB,*str);
    216:   return(0);
    217: }

    219: /* ------------------------------------------------------------*/
    222: static PetscErrorCode TSSetUp_BEuler_Linear_Constant_Matrix(TS ts)
    223: {
    224:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

    228:   KSPSetFromOptions(ts->ksp);
    229:   VecDuplicate(ts->vec_sol,&beuler->update);
    230:   VecDuplicate(ts->vec_sol,&beuler->rhs);
    231:
    232:   /* build linear system to be solved */
    233:   TSScaleShiftMatrices(ts,ts->A,ts->B,SAME_NONZERO_PATTERN);
    234:   KSPSetOperators(ts->ksp,ts->A,ts->B,SAME_NONZERO_PATTERN);
    235:   return(0);
    236: }

    240: static PetscErrorCode TSSetUp_BEuler_Linear_Variable_Matrix(TS ts)
    241: {
    242:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

    246:   KSPSetFromOptions(ts->ksp);
    247:   VecDuplicate(ts->vec_sol,&beuler->update);
    248:   VecDuplicate(ts->vec_sol,&beuler->rhs);
    249:   return(0);
    250: }

    254: static PetscErrorCode TSSetUp_BEuler_Nonlinear(TS ts)
    255: {
    256:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

    260:   VecDuplicate(ts->vec_sol,&beuler->update);
    261:   VecDuplicate(ts->vec_sol,&beuler->func);
    262:   SNESSetFunction(ts->snes,beuler->func,TSBEulerFunction,ts);
    263:   SNESSetJacobian(ts->snes,ts->A,ts->B,TSBEulerJacobian,ts);
    264:   return(0);
    265: }
    266: /*------------------------------------------------------------*/

    270: static PetscErrorCode TSSetFromOptions_BEuler_Linear(TS ts)
    271: {
    273:   return(0);
    274: }

    278: static PetscErrorCode TSSetFromOptions_BEuler_Nonlinear(TS ts)
    279: {
    281:   return(0);
    282: }

    286: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
    287: {
    289:   return(0);
    290: }

    292: /* ------------------------------------------------------------ */
    293: /*MC
    294:       TS_BEULER - ODE solver using the implicit backward Euler method

    296:   Level: beginner

    298: .seealso:  TSCreate(), TS, TSSetType(), TS_EULER

    300: M*/
    304: PetscErrorCode PETSCTS_DLLEXPORT TSCreate_BEuler(TS ts)
    305: {
    306:   TS_BEuler      *beuler;

    310:   ts->ops->destroy = TSDestroy_BEuler;
    311:   ts->ops->view    = TSView_BEuler;

    313:   if (ts->problem_type == TS_LINEAR) {
    314:     if (!ts->A) {
    315:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set rhs matrix for linear problem");
    316:     }
    317:     if (!ts->ops->rhsmatrix) {
    318:       ts->ops->setup  = TSSetUp_BEuler_Linear_Constant_Matrix;
    319:       ts->ops->step   = TSStep_BEuler_Linear_Constant_Matrix;
    320:     } else {
    321:       ts->ops->setup  = TSSetUp_BEuler_Linear_Variable_Matrix;
    322:       ts->ops->step   = TSStep_BEuler_Linear_Variable_Matrix;
    323:     }
    324:     ts->ops->setfromoptions  = TSSetFromOptions_BEuler_Linear;
    325:     KSPCreate(ts->comm,&ts->ksp);
    326:     KSPSetInitialGuessNonzero(ts->ksp,PETSC_TRUE);
    327:   } else if (ts->problem_type == TS_NONLINEAR) {
    328:     ts->ops->setup           = TSSetUp_BEuler_Nonlinear;
    329:     ts->ops->step            = TSStep_BEuler_Nonlinear;
    330:     ts->ops->setfromoptions  = TSSetFromOptions_BEuler_Nonlinear;
    331:     SNESCreate(ts->comm,&ts->snes);
    332:   } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"No such problem");

    334:   PetscNew(TS_BEuler,&beuler);
    335:   PetscLogObjectMemory(ts,sizeof(TS_BEuler));
    336:   ts->data = (void*)beuler;

    338:   return(0);
    339: }