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

      2: static char help[] = "Solves 1D heat equation with FEM formulation.\n\
      3: Input arguments are\n\
      4:   -useAlhs: solve Alhs*U' =  (Arhs*U + g) \n\
      5:             otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n";

      7: /*--------------------------------------------------------------------------
      8:   Solves 1D heat equation U_t = U_xx with FEM formulation:
      9:                           Alhs*U' = rhs (= Arhs*U + g)
     10:   We thank Chris Cox <clcox@clemson.edu> for contributing the original code
     11: ----------------------------------------------------------------------------*/

     13:  #include <petscksp.h>
     14:  #include <petscts.h>

     16: /* special variable - max size of all arrays  */
     17: #define num_z 10

     19: /*
     20:    User-defined application context - contains data needed by the
     21:    application-provided call-back routines.
     22: */
     23: typedef struct {
     24:   Mat         Amat;               /* left hand side matrix */
     25:   Vec         ksp_rhs,ksp_sol;    /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */
     26:   int         max_probsz;         /* max size of the problem */
     27:   PetscBool   useAlhs;            /* flag (1 indicates solving Alhs*U' = Arhs*U+g */
     28:   int         nz;                 /* total number of grid points */
     29:   PetscInt    m;                  /* total number of interio grid points */
     30:   Vec         solution;           /* global exact ts solution vector */
     31:   PetscScalar *z;                 /* array of grid points */
     32:   PetscBool   debug;              /* flag (1 indicates activation of debugging printouts) */
     33: } AppCtx;

     35: extern PetscScalar exact(PetscScalar,PetscReal);
     36: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
     37: extern PetscErrorCode Petsc_KSPSolve(AppCtx*);
     38: extern PetscScalar bspl(PetscScalar*,PetscScalar,PetscInt,PetscInt,PetscInt[][2],PetscInt);
     39: extern PetscErrorCode femBg(PetscScalar[][3],PetscScalar*,PetscInt,PetscScalar*,PetscReal);
     40: extern PetscErrorCode femA(AppCtx*,PetscInt,PetscScalar*);
     41: extern PetscErrorCode rhs(AppCtx*,PetscScalar*, PetscInt, PetscScalar*,PetscReal);
     42: extern PetscErrorCode RHSfunction(TS,PetscReal,Vec,Vec,void*);

     44: int main(int argc,char **argv)
     45: {
     46:   PetscInt       i,m,nz,steps,max_steps,k,nphase=1;
     47:   PetscScalar    zInitial,zFinal,val,*z;
     48:   PetscReal      stepsz[4],T,ftime;
     50:   TS             ts;
     51:   SNES           snes;
     52:   Mat            Jmat;
     53:   AppCtx         appctx;   /* user-defined application context */
     54:   Vec            init_sol; /* ts solution vector */
     55:   PetscMPIInt    size;

     57:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     58:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     59:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only");

     61:   PetscOptionsGetInt(NULL,NULL,"-nphase",&nphase,NULL);
     62:   if (nphase > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nphase must be an integer between 1 and 3");

     64:   /* initializations */
     65:   zInitial  = 0.0;
     66:   zFinal    = 1.0;
     67:   T         = 0.014/nphase;
     68:   nz        = num_z;
     69:   m         = nz-2;
     70:   appctx.nz = nz;
     71:   max_steps = (PetscInt)10000;

     73:   appctx.m          = m;
     74:   appctx.max_probsz = nz;
     75:   appctx.debug      = PETSC_FALSE;
     76:   appctx.useAlhs    = PETSC_FALSE;

     78:   PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);
     79:   PetscOptionsHasName(NULL,NULL,"-useAlhs",&appctx.useAlhs);

     81:   /* create vector to hold ts solution */
     82:   /*-----------------------------------*/
     83:   VecCreate(PETSC_COMM_WORLD, &init_sol);
     84:   VecSetSizes(init_sol, PETSC_DECIDE, m);
     85:   VecSetFromOptions(init_sol);

     87:   /* create vector to hold true ts soln for comparison */
     88:   VecDuplicate(init_sol, &appctx.solution);

     90:   /* create LHS matrix Amat */
     91:   /*------------------------*/
     92:   MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat);
     93:   MatSetFromOptions(appctx.Amat);
     94:   MatSetUp(appctx.Amat);
     95:   /* set space grid points - interio points only! */
     96:   PetscMalloc1(nz+1,&z);
     97:   for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1));
     98:   appctx.z = z;
     99:   femA(&appctx,nz,z);

    101:   /* create the jacobian matrix */
    102:   /*----------------------------*/
    103:   MatCreate(PETSC_COMM_WORLD, &Jmat);
    104:   MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m);
    105:   MatSetFromOptions(Jmat);
    106:   MatSetUp(Jmat);

    108:   /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
    109:   VecDuplicate(init_sol,&appctx.ksp_rhs);
    110:   VecDuplicate(init_sol,&appctx.ksp_sol);

    112:   /* set intial guess */
    113:   /*------------------*/
    114:   for (i=0; i<nz-2; i++) {
    115:     val  = exact(z[i+1], 0.0);
    116:     VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES);
    117:   }
    118:   VecAssemblyBegin(init_sol);
    119:   VecAssemblyEnd(init_sol);

    121:   /*create a time-stepping context and set the problem type */
    122:   /*--------------------------------------------------------*/
    123:   TSCreate(PETSC_COMM_WORLD, &ts);
    124:   TSSetProblemType(ts,TS_NONLINEAR);

    126:   /* set time-step method */
    127:   TSSetType(ts,TSCN);

    129:   /* Set optional user-defined monitoring routine */
    130:   TSMonitorSet(ts,Monitor,&appctx,NULL);
    131:   /* set the right hand side of U_t = RHSfunction(U,t) */
    132:   TSSetRHSFunction(ts,NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx);

    134:   if (appctx.useAlhs) {
    135:     /* set the left hand side matrix of Amat*U_t = rhs(U,t) */

    137:     /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
    138:      * Alhs matrix without making a copy.  Either finite difference the entire thing or use analytic Jacobians in both
    139:      * places.
    140:      */
    141:     TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,&appctx);
    142:     TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx);
    143:   }

    145:   /* use petsc to compute the jacobian by finite differences */
    146:   TSGetSNES(ts,&snes);
    147:   SNESSetJacobian(snes,Jmat,Jmat,SNESComputeJacobianDefault,NULL);

    149:   /* get the command line options if there are any and set them */
    150:   TSSetFromOptions(ts);

    152: #if defined(PETSC_HAVE_SUNDIALS)
    153:   {
    154:     TSType    type;
    155:     PetscBool sundialstype=PETSC_FALSE;
    156:     TSGetType(ts,&type);
    157:     PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype);
    158:     if (sundialstype && appctx.useAlhs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type");
    159:   }
    160: #endif
    161:   /* Sets the initial solution */
    162:   TSSetSolution(ts,init_sol);

    164:   stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
    165:   ftime     = 0.0;
    166:   for (k=0; k<nphase; k++) {
    167:     if (nphase > 1) {PetscPrintf(PETSC_COMM_WORLD,"Phase %D initial time %g, stepsz %g, duration: %g\n",k,(double)ftime,(double)stepsz[k],(double)((k+1)*T));}
    168:     TSSetTime(ts,ftime);
    169:     TSSetTimeStep(ts,stepsz[k]);
    170:     TSSetMaxSteps(ts,max_steps);
    171:     TSSetMaxTime(ts,(k+1)*T);
    172:     TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

    174:     /* loop over time steps */
    175:     /*----------------------*/
    176:     TSSolve(ts,init_sol);
    177:     TSGetSolveTime(ts,&ftime);
    178:     TSGetStepNumber(ts,&steps);
    179:     stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */
    180:   }

    182:   /* free space */
    183:   TSDestroy(&ts);
    184:   MatDestroy(&appctx.Amat);
    185:   MatDestroy(&Jmat);
    186:   VecDestroy(&appctx.ksp_rhs);
    187:   VecDestroy(&appctx.ksp_sol);
    188:   VecDestroy(&init_sol);
    189:   VecDestroy(&appctx.solution);
    190:   PetscFree(z);

    192:   PetscFinalize();
    193:   return ierr;
    194: }

    196: /*------------------------------------------------------------------------
    197:   Set exact solution
    198:   u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
    199: --------------------------------------------------------------------------*/
    200: PetscScalar exact(PetscScalar z,PetscReal t)
    201: {
    202:   PetscScalar val, ex1, ex2;

    204:   ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t);
    205:   ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t);
    206:   val = PetscSinScalar(6*PETSC_PI*z)*ex1 + 3.*PetscSinScalar(2*PETSC_PI*z)*ex2;
    207:   return val;
    208: }

    210: /*
    211:    Monitor - User-provided routine to monitor the solution computed at
    212:    each timestep.  This example plots the solution and computes the
    213:    error in two different norms.

    215:    Input Parameters:
    216:    ts     - the timestep context
    217:    step   - the count of the current step (with 0 meaning the
    218:              initial condition)
    219:    time   - the current time
    220:    u      - the solution at this timestep
    221:    ctx    - the user-provided context for this monitoring routine.
    222:             In this case we use the application context which contains
    223:             information about the problem size, workspace and the exact
    224:             solution.
    225: */
    226: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
    227: {
    228:   AppCtx         *appctx = (AppCtx*)ctx;
    230:   PetscInt       i,m=appctx->m;
    231:   PetscReal      norm_2,norm_max,h=1.0/(m+1);
    232:   PetscScalar    *u_exact;

    234:   /* Compute the exact solution */
    235:   VecGetArray(appctx->solution,&u_exact);
    236:   for (i=0; i<m; i++) u_exact[i] = exact(appctx->z[i+1],time);
    237:   VecRestoreArray(appctx->solution,&u_exact);

    239:   /* Print debugging information if desired */
    240:   if (appctx->debug) {
    241:     PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",(double)time);
    242:     VecView(u,PETSC_VIEWER_STDOUT_SELF);
    243:     PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");
    244:     VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
    245:   }

    247:   /* Compute the 2-norm and max-norm of the error */
    248:   VecAXPY(appctx->solution,-1.0,u);
    249:   VecNorm(appctx->solution,NORM_2,&norm_2);

    251:   norm_2 = PetscSqrtReal(h)*norm_2;
    252:   VecNorm(appctx->solution,NORM_MAX,&norm_max);

    254:   PetscPrintf(PETSC_COMM_SELF,"Timestep %D: time = %g, 2-norm error = %6.4f, max norm error = %6.4f\n",step,(double)time,(double)norm_2,(double)norm_max);

    256:   /*
    257:      Print debugging information if desired
    258:   */
    259:   if (appctx->debug) {
    260:     PetscPrintf(PETSC_COMM_SELF,"Error vector\n");
    261:     VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
    262:   }
    263:   return 0;
    264: }

    266: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    267: %%      Function to solve a linear system using KSP                                           %%
    268: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

    270: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
    271: {
    273:   KSP            ksp;
    274:   PC             pc;

    276:   /*create the ksp context and set the operators,that is, associate the system matrix with it*/
    277:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    278:   KSPSetOperators(ksp,obj->Amat,obj->Amat);

    280:   /*get the preconditioner context, set its type and the tolerances*/
    281:   KSPGetPC(ksp,&pc);
    282:   PCSetType(pc,PCLU);
    283:   KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);

    285:   /*get the command line options if there are any and set them*/
    286:   KSPSetFromOptions(ksp);

    288:   /*get the linear system (ksp) solve*/
    289:   KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);

    291:   KSPDestroy(&ksp);
    292:   return 0;
    293: }

    295: /***********************************************************************
    296:  * Function to return value of basis function or derivative of basis   *
    297:  *              function.                                              *
    298:  ***********************************************************************
    299:  *                                                                     *
    300:  *       Arguments:                                                    *
    301:  *         x       = array of xpoints or nodal values                  *
    302:  *         xx      = point at which the basis function is to be        *
    303:  *                     evaluated.                                      *
    304:  *         il      = interval containing xx.                           *
    305:  *         iq      = indicates which of the two basis functions in     *
    306:  *                     interval intrvl should be used                  *
    307:  *         nll     = array containing the endpoints of each interval.  *
    308:  *         id      = If id ~= 2, the value of the basis function       *
    309:  *                     is calculated; if id = 2, the value of the      *
    310:  *                     derivative of the basis function is returned.   *
    311:  ***********************************************************************/

    313: PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)
    314: {
    315:   PetscScalar x1,x2,bfcn;
    316:   PetscInt    i1,i2,iq1,iq2;

    318:   /*** Determine which basis function in interval intrvl is to be used in ***/
    319:   iq1 = iq;
    320:   if (iq1==0) iq2 = 1;
    321:   else iq2 = 0;

    323:   /***  Determine endpoint of the interval intrvl ***/
    324:   i1=nll[il][iq1];
    325:   i2=nll[il][iq2];

    327:   /*** Determine nodal values at the endpoints of the interval intrvl ***/
    328:   x1=x[i1];
    329:   x2=x[i2];
    330:   /* printf("x1=%g\tx2=%g\txx=%g\n",x1,x2,xx); */
    331:   /*** Evaluate basis function ***/
    332:   if (id == 2) bfcn=(1.0)/(x1-x2);
    333:   else bfcn=(xx-x2)/(x1-x2);
    334:   /* printf("bfcn=%g\n",bfcn); */
    335:   return bfcn;
    336: }

    338: /*---------------------------------------------------------
    339:   Function called by rhs function to get B and g
    340: ---------------------------------------------------------*/
    341: PetscErrorCode femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t)
    342: {
    343:   PetscInt    i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq;
    344:   PetscInt    nli[num_z][2],indx[num_z];
    345:   PetscScalar dd,dl,zip,zipq,zz,b_z,bb_z,bij;
    346:   PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3];

    348:   /*  initializing everything - btri and f are initialized in rhs.c  */
    349:   for (i=0; i < nz; i++) {
    350:     nli[i][0]   = 0;
    351:     nli[i][1]   = 0;
    352:     indx[i]     = 0;
    353:     zquad[i][0] = 0.0;
    354:     zquad[i][1] = 0.0;
    355:     zquad[i][2] = 0.0;
    356:     dlen[i]     = 0.0;
    357:   } /*end for (i)*/

    359:   /*  quadrature weights  */
    360:   qdwt[0] = 1.0/6.0;
    361:   qdwt[1] = 4.0/6.0;
    362:   qdwt[2] = 1.0/6.0;

    364:   /* 1st and last nodes have Dirichlet boundary condition -
    365:      set indices there to -1 */

    367:   for (i=0; i < nz-1; i++) indx[i] = i-1;
    368:   indx[nz-1] = -1;

    370:   ipq = 0;
    371:   for (il=0; il < nz-1; il++) {
    372:     ip           = ipq;
    373:     ipq          = ip+1;
    374:     zip          = z[ip];
    375:     zipq         = z[ipq];
    376:     dl           = zipq-zip;
    377:     zquad[il][0] = zip;
    378:     zquad[il][1] = (0.5)*(zip+zipq);
    379:     zquad[il][2] = zipq;
    380:     dlen[il]     = PetscAbsScalar(dl);
    381:     nli[il][0]   = ip;
    382:     nli[il][1]   = ipq;
    383:   }

    385:   for (il=0; il < nz-1; il++) {
    386:     for (iquad=0; iquad < 3; iquad++) {
    387:       dd = (dlen[il])*(qdwt[iquad]);
    388:       zz = zquad[il][iquad];

    390:       for (iq=0; iq < 2; iq++) {
    391:         ip  = nli[il][iq];
    392:         b_z = bspl(z,zz,il,iq,nli,2);
    393:         i   = indx[ip];

    395:         if (i > -1) {
    396:           for (iqq=0; iqq < 2; iqq++) {
    397:             ipp  = nli[il][iqq];
    398:             bb_z = bspl(z,zz,il,iqq,nli,2);
    399:             j    = indx[ipp];
    400:             bij  = -b_z*bb_z;

    402:             if (j > -1) {
    403:               jj = 1+j-i;
    404:               btri[i][jj] += bij*dd;
    405:             } else {
    406:               f[i] += bij*dd*exact(z[ipp], t);
    407:               /* f[i] += 0.0; */
    408:               /* if (il==0 && j==-1) { */
    409:               /* f[i] += bij*dd*exact(zz,t); */
    410:               /* }*/ /*end if*/
    411:             } /*end else*/
    412:           } /*end for (iqq)*/
    413:         } /*end if (i>0)*/
    414:       } /*end for (iq)*/
    415:     } /*end for (iquad)*/
    416:   } /*end for (il)*/
    417:   return 0;
    418: }

    420: PetscErrorCode femA(AppCtx *obj,PetscInt nz,PetscScalar *z)
    421: {
    422:   PetscInt       i,j,il,ip,ipp,ipq,iq,iquad,iqq;
    423:   PetscInt       nli[num_z][2],indx[num_z];
    424:   PetscScalar    dd,dl,zip,zipq,zz,bb,bbb,aij;
    425:   PetscScalar    rquad[num_z][3],dlen[num_z],qdwt[3],add_term;

    428:   /*  initializing everything  */

    430:   for (i=0; i < nz; i++) {
    431:     nli[i][0]   = 0;
    432:     nli[i][1]   = 0;
    433:     indx[i]     = 0;
    434:     rquad[i][0] = 0.0;
    435:     rquad[i][1] = 0.0;
    436:     rquad[i][2] = 0.0;
    437:     dlen[i]     = 0.0;
    438:   } /*end for (i)*/

    440:   /*  quadrature weights  */
    441:   qdwt[0] = 1.0/6.0;
    442:   qdwt[1] = 4.0/6.0;
    443:   qdwt[2] = 1.0/6.0;

    445:   /* 1st and last nodes have Dirichlet boundary condition -
    446:      set indices there to -1 */

    448:   for (i=0; i < nz-1; i++) indx[i]=i-1;
    449:   indx[nz-1]=-1;

    451:   ipq = 0;

    453:   for (il=0; il < nz-1; il++) {
    454:     ip           = ipq;
    455:     ipq          = ip+1;
    456:     zip          = z[ip];
    457:     zipq         = z[ipq];
    458:     dl           = zipq-zip;
    459:     rquad[il][0] = zip;
    460:     rquad[il][1] = (0.5)*(zip+zipq);
    461:     rquad[il][2] = zipq;
    462:     dlen[il]     = PetscAbsScalar(dl);
    463:     nli[il][0]   = ip;
    464:     nli[il][1]   = ipq;
    465:   } /*end for (il)*/

    467:   for (il=0; il < nz-1; il++) {
    468:     for (iquad=0; iquad < 3; iquad++) {
    469:       dd = (dlen[il])*(qdwt[iquad]);
    470:       zz = rquad[il][iquad];

    472:       for (iq=0; iq < 2; iq++) {
    473:         ip = nli[il][iq];
    474:         bb = bspl(z,zz,il,iq,nli,1);
    475:         i = indx[ip];
    476:         if (i > -1) {
    477:           for (iqq=0; iqq < 2; iqq++) {
    478:             ipp = nli[il][iqq];
    479:             bbb = bspl(z,zz,il,iqq,nli,1);
    480:             j = indx[ipp];
    481:             aij = bb*bbb;
    482:             if (j > -1) {
    483:               add_term = aij*dd;
    484:               MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES);
    485:             }/*endif*/
    486:           } /*end for (iqq)*/
    487:         } /*end if (i>0)*/
    488:       } /*end for (iq)*/
    489:     } /*end for (iquad)*/
    490:   } /*end for (il)*/
    491:   MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY);
    492:   MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY);
    493:   return 0;
    494: }

    496: /*---------------------------------------------------------
    497:         Function to fill the rhs vector with
    498:         By + g values ****
    499: ---------------------------------------------------------*/
    500: PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
    501: {
    502:   PetscInt       i,j,js,je,jj;
    503:   PetscScalar    val,g[num_z],btri[num_z][3],add_term;

    506:   for (i=0; i < nz-2; i++) {
    507:     for (j=0; j <= 2; j++) btri[i][j]=0.0;
    508:     g[i] = 0.0;
    509:   }

    511:   /*  call femBg to set the tri-diagonal b matrix and vector g  */
    512:   femBg(btri,g,nz,z,t);

    514:   /*  setting the entries of the right hand side vector  */
    515:   for (i=0; i < nz-2; i++) {
    516:     val = 0.0;
    517:     js  = 0;
    518:     if (i == 0) js = 1;
    519:     je = 2;
    520:     if (i == nz-2) je = 1;

    522:     for (jj=js; jj <= je; jj++) {
    523:       j    = i+jj-1;
    524:       val += (btri[i][jj])*(y[j]);
    525:     }
    526:     add_term = val + g[i];
    527:     VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);
    528:   }
    529:   VecAssemblyBegin(obj->ksp_rhs);
    530:   VecAssemblyEnd(obj->ksp_rhs);

    532:   /*  return to main driver function  */
    533:   return 0;
    534: }

    536: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    537: %%   Function to form the right hand side of the time-stepping problem.                       %%
    538: %% -------------------------------------------------------------------------------------------%%
    539:   if (useAlhs):
    540:     globalout = By+g
    541:   else if (!useAlhs):
    542:     globalout = f(y,t)=Ainv(By+g),
    543:       in which the ksp solver to transform the problem A*ydot=By+g
    544:       to the problem ydot=f(y,t)=inv(A)*(By+g)
    545: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

    547: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
    548: {
    549:   PetscErrorCode    ierr;
    550:   AppCtx            *obj = (AppCtx*)ctx;
    551:   PetscScalar       soln[num_z];
    552:   const PetscScalar *soln_ptr;
    553:   PetscInt          i,nz=obj->nz;
    554:   PetscReal         time;

    556:   /* get the previous solution to compute updated system */
    557:   VecGetArrayRead(globalin,&soln_ptr);
    558:   for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i];
    559:   VecRestoreArrayRead(globalin,&soln_ptr);
    560:   soln[num_z-1] = 0.0;
    561:   soln[num_z-2] = 0.0;

    563:   /* clear out the matrix and rhs for ksp to keep things straight */
    564:   VecSet(obj->ksp_rhs,(PetscScalar)0.0);

    566:   time = t;
    567:   /* get the updated system */
    568:   rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */

    570:   /* do a ksp solve to get the rhs for the ts problem */
    571:   if (obj->useAlhs) {
    572:     /* ksp_sol = ksp_rhs */
    573:     VecCopy(obj->ksp_rhs,globalout);
    574:   } else {
    575:     /* ksp_sol = inv(Amat)*ksp_rhs */
    576:     Petsc_KSPSolve(obj);
    577:     VecCopy(obj->ksp_sol,globalout);
    578:   }
    579:   return 0;
    580: }

    582: /*TEST

    584:     build:
    585:       requires: !complex

    587:     test:
    588:       suffix: euler
    589:       output_file: output/ex3.out

    591:     test:
    592:       suffix: 2
    593:       args:   -useAlhs
    594:       output_file: output/ex3.out
    595:       TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant

    597: TEST*/