Actual source code: ex31.c
petsc-3.4.2 2013-07-02

      1: /*T
      2:    Concepts: KSP^solving a system of linear equations
      3:    Concepts: KSP^semi-implicit
      4:    Processors: n
      5: T*/

      7: /*
      8: This is intended to be a prototypical example of the semi-implicit algorithm for
      9: a PDE. We have three phases:

     11:   1) An explicit predictor step

     13:      u^{k+1/3} = P(u^k)

     15:   2) An implicit corrector step

     17:      \Delta u^{k+2/3} = F(u^{k+1/3})

     19:   3) An explicit update step

     21:      u^{k+1} = C(u^{k+2/3})

     23: We will solve on the unit square with Dirichlet boundary conditions

     25:    u = f(x,y) for x = 0, x = 1, y = 0, y = 1

     27: Although we are using a DMDA, and thus have a structured mesh, we will discretize
     28: the problem with finite elements, splitting each cell of the DMDA into two
     29: triangles.

     31: This uses multigrid to solve the linear system
     32: */

     34: static char help[] = "Solves 2D compressible Euler using multigrid.\n\n";

     36: #include <petscdmda.h>
     37: #include <petscksp.h>

     39: typedef struct {
     40:   Vec rho;     /* The mass solution \rho */
     41:   Vec rho_u;   /* The x-momentum solution \rho u */
     42:   Vec rho_v;   /* The y-momentum solution \rho v */
     43:   Vec rho_e;   /* The energy solution \rho e_t */
     44:   Vec p;       /* The pressure solution P */
     45:   Vec t;       /* The temperature solution T */
     46:   Vec u;       /* The x-velocity solution u */
     47:   Vec v;       /* The y-velocity solution v */
     48: } SolutionContext;

     50: typedef struct {
     51:   SolutionContext sol_n;   /* The solution at time t^n */
     52:   SolutionContext sol_phi; /* The element-averaged solution at time t^{n+\phi} */
     53:   SolutionContext sol_np1; /* The solution at time t^{n+1} */
     54:   Vec             mu;      /* The dynamic viscosity \mu(T) at time n */
     55:   Vec             kappa;   /* The thermal conductivity \kappa(T) at time n */
     56:   PetscScalar     phi;     /* The time weighting parameter */
     57:   PetscScalar     dt;      /* The timestep \Delta t */
     58: } UserContext;

     60: extern PetscErrorCode CreateStructures(DM,UserContext*);
     61: extern PetscErrorCode DestroyStructures(DM,UserContext*);
     62: extern PetscErrorCode ComputePredictor(DM,UserContext*);
     63: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
     64: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
     65: extern PetscErrorCode ComputeCorrector(DM,Vec,Vec);

     69: int main(int argc,char **argv)
     70: {
     71:   KSP            ksp;
     72:   DM             da;
     73:   Vec            x, xNew;
     74:   UserContext    user;

     77:   PetscInitialize(&argc,&argv,(char*)0,help);

     79:   KSPCreate(PETSC_COMM_WORLD,&ksp);
     80:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
     81:   DMSetApplicationContext(da, &user);
     82:   KSPSetDM(ksp, da);

     84:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PCICE", "DM");
     85:   user.phi = 0.5;
     86:   PetscOptionsScalar("-phi", "The time weighting parameter", "ex31.c", user.phi, &user.phi, NULL);
     87:   user.dt  = 0.1;
     88:   PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, NULL);
     89:   PetscOptionsEnd();

     91:   CreateStructures(da, &user);
     92:   ComputePredictor(da, &user);

     94:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
     95:   KSPSetComputeOperators(ksp,ComputeMatrix,&user);
     96:   KSPSetFromOptions(ksp);
     97:   KSPSolve(ksp, NULL, NULL);

     99:   KSPGetSolution(ksp, &x);
    100:   VecDuplicate(x, &xNew);
    101:   ComputeCorrector(da, x, xNew);
    102:   VecDestroy(&xNew);

    104:   DestroyStructures(da, &user);
    105:   DMDestroy(&da);
    106:   KSPDestroy(&ksp);
    107:   PetscFinalize();
    108:   return 0;
    109: }

    113: PetscErrorCode CreateStructures(DM da, UserContext *user)
    114: {
    115:   const PetscInt *necon;
    116:   PetscInt       ne,nc;

    120:   DMDAGetElements(da,&ne,&nc,&necon);
    121:   DMDARestoreElements(da,&ne,&nc,&necon);
    122:   DMCreateGlobalVector(da, &user->sol_n.rho);
    123:   DMCreateGlobalVector(da, &user->sol_n.rho_u);
    124:   DMCreateGlobalVector(da, &user->sol_n.rho_v);
    125:   DMCreateGlobalVector(da, &user->sol_n.rho_e);
    126:   DMCreateGlobalVector(da, &user->sol_n.p);
    127:   DMCreateGlobalVector(da, &user->sol_n.u);
    128:   DMCreateGlobalVector(da, &user->sol_n.v);
    129:   DMCreateGlobalVector(da, &user->sol_n.t);
    130:   VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho);
    131:   VecSetSizes(user->sol_phi.rho, ne, PETSC_DECIDE);
    132:   VecSetType(user->sol_phi.rho,VECMPI);
    133:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_u);
    134:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_v);
    135:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.u);
    136:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.v);
    137:   DMCreateGlobalVector(da, &user->sol_np1.rho);
    138:   DMCreateGlobalVector(da, &user->sol_np1.rho_u);
    139:   DMCreateGlobalVector(da, &user->sol_np1.rho_v);
    140:   DMCreateGlobalVector(da, &user->sol_np1.rho_e);
    141:   DMCreateGlobalVector(da, &user->sol_np1.p);
    142:   DMCreateGlobalVector(da, &user->sol_np1.u);
    143:   DMCreateGlobalVector(da, &user->sol_np1.v);
    144:   DMCreateGlobalVector(da, &user->mu);
    145:   DMCreateGlobalVector(da, &user->kappa);
    146:   return(0);
    147: }

    151: PetscErrorCode DestroyStructures(DM da, UserContext   *user)
    152: {

    156:   VecDestroy(&user->sol_n.rho);
    157:   VecDestroy(&user->sol_n.rho_u);
    158:   VecDestroy(&user->sol_n.rho_v);
    159:   VecDestroy(&user->sol_n.rho_e);
    160:   VecDestroy(&user->sol_n.p);
    161:   VecDestroy(&user->sol_n.u);
    162:   VecDestroy(&user->sol_n.v);
    163:   VecDestroy(&user->sol_n.t);
    164:   VecDestroy(&user->sol_phi.rho);
    165:   VecDestroy(&user->sol_phi.rho_u);
    166:   VecDestroy(&user->sol_phi.rho_v);
    167:   VecDestroy(&user->sol_phi.u);
    168:   VecDestroy(&user->sol_phi.v);
    169:   VecDestroy(&user->sol_np1.rho);
    170:   VecDestroy(&user->sol_np1.rho_u);
    171:   VecDestroy(&user->sol_np1.rho_v);
    172:   VecDestroy(&user->sol_np1.rho_e);
    173:   VecDestroy(&user->sol_np1.p);
    174:   VecDestroy(&user->sol_np1.u);
    175:   VecDestroy(&user->sol_np1.v);
    176:   VecDestroy(&user->mu);
    177:   VecDestroy(&user->kappa);
    178:   return(0);
    179: }

    183: /* Average the velocity (u,v) at time t^n over each element for time n+\phi */
    184: PetscErrorCode CalculateElementVelocity(DM da, UserContext *user)
    185: {
    186:   PetscScalar    *u_n,   *v_n;
    187:   PetscScalar    *u_phi, *v_phi;
    188:   const PetscInt *necon;
    189:   PetscInt       j, e, ne, nc;

    193:   DMDAGetElements(da, &ne, &nc, &necon);
    194:   VecGetArray(user->sol_n.u, &u_n);
    195:   VecGetArray(user->sol_n.v, &v_n);
    196:   PetscMalloc(ne*sizeof(PetscScalar),&u_phi);
    197:   PetscMalloc(ne*sizeof(PetscScalar),&v_phi);
    198:   for (e = 0; e < ne; e++) {
    199:     u_phi[e] = 0.0;
    200:     v_phi[e] = 0.0;
    201:     for (j = 0; j < 3; j++) {
    202:       u_phi[e] += u_n[necon[3*e+j]];
    203:       v_phi[e] += v_n[necon[3*e+j]];
    204:     }
    205:     u_phi[e] /= 3.0;
    206:     v_phi[e] /= 3.0;
    207:   }
    208:   PetscFree(u_phi);
    209:   PetscFree(v_phi);
    210:   DMDARestoreElements(da, &ne,&nc, &necon);
    211:   VecRestoreArray(user->sol_n.u, &u_n);
    212:   VecRestoreArray(user->sol_n.v, &v_n);
    213:   return(0);
    214: }

    218: /* This is equation 32,

    220:    U^{n+\phi}_E = {1\over Vol_E} \left(\int_\Omega [N]{U^n} d\Omega - \phi\Delta t \int_\Omega [\nabla N]\cdot{F^n} d\Omega \right) + \phi\Delta t Q^n

    222: which is really simple for linear elements

    224:    U^{n+\phi}_E = {1\over3} \sum^3_{i=1} U^n_i - \phi\Delta t [\nabla N]\cdot{F^n} + \phi\Delta t Q^n

    226: where

    228:    U^{n+\phi}_E = {\rho  \rho u  \rho v}^{n+\phi}_E

    230: and the x and y components of the convective fluxes F are

    232:    f^n = {\rho u  \rho u^2  \rho uv}^n      g^n = {\rho v  \rho uv  \rho v^2}^n
    233: */
    234: PetscErrorCode TaylorGalerkinStepI(DM da, UserContext *user)
    235: {
    236:   PetscScalar    phi_dt = user->phi*user->dt;
    237:   PetscScalar    *u_n,     *v_n;
    238:   PetscScalar    *rho_n,   *rho_u_n,   *rho_v_n;
    239:   PetscScalar    *rho_phi, *rho_u_phi, *rho_v_phi;
    240:   PetscScalar    Fx_x, Fy_y;
    241:   PetscScalar    psi_x[3], psi_y[3];
    242:   PetscInt       idx[3];
    243:   PetscReal      hx, hy;
    244:   const PetscInt *necon;
    245:   PetscInt       j, e, ne, nc, mx, my;

    249:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
    250:   hx   = 1.0 / (PetscReal)(mx-1);
    251:   hy   = 1.0 / (PetscReal)(my-1);
    252:   VecSet(user->sol_phi.rho,0.0);
    253:   VecSet(user->sol_phi.rho_u,0.0);
    254:   VecSet(user->sol_phi.rho_v,0.0);
    255:   VecGetArray(user->sol_n.u,       &u_n);
    256:   VecGetArray(user->sol_n.v,       &v_n);
    257:   VecGetArray(user->sol_n.rho,     &rho_n);
    258:   VecGetArray(user->sol_n.rho_u,   &rho_u_n);
    259:   VecGetArray(user->sol_n.rho_v,   &rho_v_n);
    260:   VecGetArray(user->sol_phi.rho,   &rho_phi);
    261:   VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
    262:   VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
    263:   DMDAGetElements(da, &ne, &nc, &necon);
    264:   for (e = 0; e < ne; e++) {
    265:     /* Average the existing fields over the element */
    266:     for (j = 0; j < 3; j++) {
    267:       idx[j]        = necon[3*e+j];
    268:       rho_phi[e]   += rho_n[idx[j]];
    269:       rho_u_phi[e] += rho_u_n[idx[j]];
    270:       rho_v_phi[e] += rho_v_n[idx[j]];
    271:     }
    272:     rho_phi[e]   /= 3.0;
    273:     rho_u_phi[e] /= 3.0;
    274:     rho_v_phi[e] /= 3.0;
    275:     /* Get basis function deriatives (we need the orientation of the element here) */
    276:     if (idx[1] > idx[0]) {
    277:       psi_x[0] = -hy; psi_x[1] =  hy; psi_x[2] = 0.0;
    278:       psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] =  hx;
    279:     } else {
    280:       psi_x[0] =  hy; psi_x[1] = -hy; psi_x[2] = 0.0;
    281:       psi_y[0] =  hx; psi_y[1] = 0.0; psi_y[2] = -hx;
    282:     }
    283:     /* Determine the convective fluxes for \rho^{n+\phi} */
    284:     Fx_x = 0.0; Fy_y = 0.0;
    285:     for (j = 0; j < 3; j++) {
    286:       Fx_x += psi_x[j]*rho_u_n[idx[j]];
    287:       Fy_y += psi_y[j]*rho_v_n[idx[j]];
    288:     }
    289:     rho_phi[e] -= phi_dt*(Fx_x + Fy_y);
    290:     /* Determine the convective fluxes for (\rho u)^{n+\phi} */
    291:     Fx_x = 0.0; Fy_y = 0.0;
    292:     for (j = 0; j < 3; j++) {
    293:       Fx_x += psi_x[j]*rho_u_n[idx[j]]*u_n[idx[j]];
    294:       Fy_y += psi_y[j]*rho_v_n[idx[j]]*u_n[idx[j]];
    295:     }
    296:     rho_u_phi[e] -= phi_dt*(Fx_x + Fy_y);
    297:     /* Determine the convective fluxes for (\rho v)^{n+\phi} */
    298:     Fx_x = 0.0; Fy_y = 0.0;
    299:     for (j = 0; j < 3; j++) {
    300:       Fx_x += psi_x[j]*rho_u_n[idx[j]]*v_n[idx[j]];
    301:       Fy_y += psi_y[j]*rho_v_n[idx[j]]*v_n[idx[j]];
    302:     }
    303:     rho_v_phi[e] -= phi_dt*(Fx_x + Fy_y);
    304:   }
    305:   DMDARestoreElements(da, &ne, &nc, &necon);
    306:   VecRestoreArray(user->sol_n.u,       &u_n);
    307:   VecRestoreArray(user->sol_n.v,       &v_n);
    308:   VecRestoreArray(user->sol_n.rho,     &rho_n);
    309:   VecRestoreArray(user->sol_n.rho_u,   &rho_u_n);
    310:   VecRestoreArray(user->sol_n.rho_v,   &rho_v_n);
    311:   VecRestoreArray(user->sol_phi.rho,   &rho_phi);
    312:   VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
    313:   VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
    314:   return(0);
    315: }

    319: /*
    320: The element stiffness matrix for the identity in linear elements is

    322:   1  /2 1 1\
    323:   -  |1 2 1|
    324:   12 \1 1 2/

    326:   no matter what the shape of the triangle. */
    327: PetscErrorCode TaylorGalerkinStepIIMomentum(DM da, UserContext *user)
    328: {
    329:   MPI_Comm       comm;
    330:   KSP            ksp;
    331:   Mat            mat;
    332:   Vec            rhs_u, rhs_v;
    333:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
    334:                                 0.08333333333, 0.16666666667, 0.08333333333,
    335:                                 0.08333333333, 0.08333333333, 0.16666666667};
    336:   PetscScalar    *u_n,       *v_n,      *mu_n;
    337:   PetscScalar    *u_phi,     *v_phi;
    338:   PetscScalar    *rho_u_phi, *rho_v_phi;
    339:   PetscInt       idx[3];
    340:   PetscScalar    values_u[3];
    341:   PetscScalar    values_v[3];
    342:   PetscScalar    psi_x[3], psi_y[3];
    343:   PetscScalar    mu, tau_xx, tau_xy, tau_yy;
    344:   PetscReal      hx, hy, area;
    345:   const PetscInt *necon;
    346:   PetscInt       j, k, e, ne, nc, mx, my;

    350:   PetscObjectGetComm((PetscObject) da, &comm);
    351:   DMCreateMatrix(da, MATAIJ, &mat);
    352:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
    353:   DMGetGlobalVector(da, &rhs_u);
    354:   DMGetGlobalVector(da, &rhs_v);
    355:   KSPCreate(comm, &ksp);
    356:   KSPSetFromOptions(ksp);

    358:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
    359:   hx   = 1.0 / (PetscReal)(mx-1);
    360:   hy   = 1.0 / (PetscReal)(my-1);
    361:   area = 0.5*hx*hy;
    362:   VecGetArray(user->sol_n.u,       &u_n);
    363:   VecGetArray(user->sol_n.v,       &v_n);
    364:   VecGetArray(user->mu,            &mu_n);
    365:   VecGetArray(user->sol_phi.u,     &u_phi);
    366:   VecGetArray(user->sol_phi.v,     &v_phi);
    367:   VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
    368:   VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
    369:   DMDAGetElements(da, &ne, &nc, &necon);
    370:   for (e = 0; e < ne; e++) {
    371:     for (j = 0; j < 3; j++) {
    372:       idx[j]      = necon[3*e+j];
    373:       values_u[j] = 0.0;
    374:       values_v[j] = 0.0;
    375:     }
    376:     /* Get basis function deriatives (we need the orientation of the element here) */
    377:     if (idx[1] > idx[0]) {
    378:       psi_x[0] = -hy; psi_x[1] =  hy; psi_x[2] = 0.0;
    379:       psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] =  hx;
    380:     } else {
    381:       psi_x[0] =  hy; psi_x[1] = -hy; psi_x[2] = 0.0;
    382:       psi_y[0] =  hx; psi_y[1] = 0.0; psi_y[2] = -hx;
    383:     }
    384:     /*  <\nabla\psi, F^{n+\phi}_e>: Divergence of the element-averaged convective fluxes */
    385:     for (j = 0; j < 3; j++) {
    386:       values_u[j] += psi_x[j]*rho_u_phi[e]*u_phi[e] + psi_y[j]*rho_u_phi[e]*v_phi[e];
    387:       values_v[j] += psi_x[j]*rho_v_phi[e]*u_phi[e] + psi_y[j]*rho_v_phi[e]*v_phi[e];
    388:     }
    389:     /*  -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
    390:     for (j = 0; j < 3; j++) {
    391:       /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
    392:       /* \tau_{xy} =     \mu(T) (  {\partial u\over\partial y} + {\partial v\over\partial x}) */
    393:       /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
    394:       mu     = 0.0;
    395:       tau_xx = 0.0;
    396:       tau_xy = 0.0;
    397:       tau_yy = 0.0;
    398:       for (k = 0; k < 3; k++) {
    399:         mu     += mu_n[idx[k]];
    400:         tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
    401:         tau_xy +=     psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
    402:         tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
    403:       }
    404:       mu          /= 3.0;
    405:       tau_xx      *= (2.0/3.0)*mu;
    406:       tau_xy      *= mu;
    407:       tau_yy      *= (2.0/3.0)*mu;
    408:       values_u[j] -= area*(psi_x[j]*tau_xx + psi_y[j]*tau_xy);
    409:       values_v[j] -= area*(psi_x[j]*tau_xy + psi_y[j]*tau_yy);
    410:     }
    411:     /* Accumulate to global structures */
    412:     VecSetValuesLocal(rhs_u, 3, idx, values_u, ADD_VALUES);
    413:     VecSetValuesLocal(rhs_v, 3, idx, values_v, ADD_VALUES);
    414:     MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
    415:   }
    416:   DMDARestoreElements(da, &ne, &nc, &necon);
    417:   VecRestoreArray(user->sol_n.u,       &u_n);
    418:   VecRestoreArray(user->sol_n.v,       &v_n);
    419:   VecRestoreArray(user->mu,            &mu_n);
    420:   VecRestoreArray(user->sol_phi.u,     &u_phi);
    421:   VecRestoreArray(user->sol_phi.v,     &v_phi);
    422:   VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
    423:   VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);

    425:   VecAssemblyBegin(rhs_u);
    426:   VecAssemblyBegin(rhs_v);
    427:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
    428:   VecAssemblyEnd(rhs_u);
    429:   VecAssemblyEnd(rhs_v);
    430:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
    431:   VecScale(rhs_u,user->dt);
    432:   VecScale(rhs_v,user->dt);

    434:   KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
    435:   KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
    436:   KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
    437:   KSPDestroy(&ksp);
    438:   MatDestroy(&mat);
    439:   DMRestoreGlobalVector(da, &rhs_u);
    440:   DMRestoreGlobalVector(da, &rhs_v);
    441:   return(0);
    442: }

    446: /* Notice that this requires the previous momentum solution.

    448: The element stiffness matrix for the identity in linear elements is

    450:   1  /2 1 1\
    451:   -  |1 2 1|
    452:   12 \1 1 2/

    454:   no matter what the shape of the triangle. */
    455: PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
    456: {
    457:   MPI_Comm       comm;
    458:   Mat            mat;
    459:   Vec            rhs_m, rhs_e;
    460:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
    461:                                 0.08333333333, 0.16666666667, 0.08333333333,
    462:                                 0.08333333333, 0.08333333333, 0.16666666667};
    463:   PetscScalar    *u_n,       *v_n,     *p_n,     *t_n,     *mu_n,    *kappa_n;
    464:   PetscScalar    *rho_n,     *rho_u_n, *rho_v_n, *rho_e_n;
    465:   PetscScalar    *u_phi,     *v_phi;
    466:   PetscScalar    *rho_u_np1, *rho_v_np1;
    467:   PetscInt       idx[3];
    468:   PetscScalar    psi_x[3], psi_y[3];
    469:   PetscScalar    values_m[3];
    470:   PetscScalar    values_e[3];
    471:   PetscScalar    phi = user->phi;
    472:   PetscScalar    mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
    473:   PetscReal      hx, hy, area;
    474:   KSP            ksp;
    475:   const PetscInt *necon;
    476:   PetscInt       j, k, e, ne, nc, mx, my;

    480:   PetscObjectGetComm((PetscObject) da, &comm);
    481:   DMCreateMatrix(da, MATAIJ, &mat);
    482:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
    483:   DMGetGlobalVector(da, &rhs_m);
    484:   DMGetGlobalVector(da, &rhs_e);
    485:   KSPCreate(comm, &ksp);
    486:   KSPSetFromOptions(ksp);

    488:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
    489:   hx   = 1.0 / (PetscReal)(mx-1);
    490:   hy   = 1.0 / (PetscReal)(my-1);
    491:   area = 0.5*hx*hy;
    492:   VecGetArray(user->sol_n.u,       &u_n);
    493:   VecGetArray(user->sol_n.v,       &v_n);
    494:   VecGetArray(user->sol_n.p,       &p_n);
    495:   VecGetArray(user->sol_n.t,       &t_n);
    496:   VecGetArray(user->mu,            &mu_n);
    497:   VecGetArray(user->kappa,         &kappa_n);
    498:   VecGetArray(user->sol_n.rho,     &rho_n);
    499:   VecGetArray(user->sol_n.rho_u,   &rho_u_n);
    500:   VecGetArray(user->sol_n.rho_v,   &rho_v_n);
    501:   VecGetArray(user->sol_n.rho_e,   &rho_e_n);
    502:   VecGetArray(user->sol_phi.u,     &u_phi);
    503:   VecGetArray(user->sol_phi.v,     &v_phi);
    504:   VecGetArray(user->sol_np1.rho_u, &rho_u_np1);
    505:   VecGetArray(user->sol_np1.rho_v, &rho_v_np1);
    506:   DMDAGetElements(da, &ne, &nc, &necon);
    507:   for (e = 0; e < ne; e++) {
    508:     for (j = 0; j < 3; j++) {
    509:       idx[j]      = necon[3*e+j];
    510:       values_m[j] = 0.0;
    511:       values_e[j] = 0.0;
    512:     }
    513:     /* Get basis function deriatives (we need the orientation of the element here) */
    514:     if (idx[1] > idx[0]) {
    515:       psi_x[0] = -hy; psi_x[1] =  hy; psi_x[2] = 0.0;
    516:       psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] =  hx;
    517:     } else {
    518:       psi_x[0] =  hy; psi_x[1] = -hy; psi_x[2] = 0.0;
    519:       psi_y[0] =  hx; psi_y[1] = 0.0; psi_y[2] = -hx;
    520:     }
    521:     /*  <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
    522:     for (j = 0; j < 3; j++) {
    523:       values_m[j] += (psi_x[j]*(phi*rho_u_np1[idx[j]] + rho_u_n[idx[j]]) + psi_y[j]*(rho_v_np1[idx[j]] + rho_v_n[idx[j]]))/3.0;
    524:       values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
    525:     }
    526:     /*  -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
    527:     for (j = 0; j < 3; j++) {
    528:       /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
    529:       /* \tau_{xy} =     \mu(T) (  {\partial u\over\partial y} + {\partial v\over\partial x}) */
    530:       /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
    531:       /* q_x       = -\kappa(T) {\partial T\over\partial x} */
    532:       /* q_y       = -\kappa(T) {\partial T\over\partial y} */

    534:       /* above code commeted out - causing ininitialized variables. */
    535:       q_x =0; q_y =0;

    537:       mu     = 0.0;
    538:       kappa  = 0.0;
    539:       tau_xx = 0.0;
    540:       tau_xy = 0.0;
    541:       tau_yy = 0.0;
    542:       for (k = 0; k < 3; k++) {
    543:         mu     += mu_n[idx[k]];
    544:         kappa  += kappa_n[idx[k]];
    545:         tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
    546:         tau_xy +=     psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
    547:         tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
    548:         q_x    += psi_x[k]*t_n[idx[k]];
    549:         q_y    += psi_y[k]*t_n[idx[k]];
    550:       }
    551:       mu          /= 3.0;
    552:       kappa       /= 3.0;
    553:       tau_xx      *= (2.0/3.0)*mu;
    554:       tau_xy      *= mu;
    555:       tau_yy      *= (2.0/3.0)*mu;
    556:       values_e[j] -= area*(psi_x[j]*(u_phi[e]*tau_xx + v_phi[e]*tau_xy + q_x) + psi_y[j]*(u_phi[e]*tau_xy + v_phi[e]*tau_yy + q_y));
    557:     }
    558:     /* Accumulate to global structures */
    559:     VecSetValuesLocal(rhs_m, 3, idx, values_m, ADD_VALUES);
    560:     VecSetValuesLocal(rhs_e, 3, idx, values_e, ADD_VALUES);
    561:     MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
    562:   }
    563:   DMDARestoreElements(da, &ne, &nc, &necon);
    564:   VecRestoreArray(user->sol_n.u,       &u_n);
    565:   VecRestoreArray(user->sol_n.v,       &v_n);
    566:   VecRestoreArray(user->sol_n.p,       &p_n);
    567:   VecRestoreArray(user->sol_n.t,       &t_n);
    568:   VecRestoreArray(user->mu,            &mu_n);
    569:   VecRestoreArray(user->kappa,         &kappa_n);
    570:   VecRestoreArray(user->sol_n.rho,     &rho_n);
    571:   VecRestoreArray(user->sol_n.rho_u,   &rho_u_n);
    572:   VecRestoreArray(user->sol_n.rho_v,   &rho_v_n);
    573:   VecRestoreArray(user->sol_n.rho_e,   &rho_e_n);
    574:   VecRestoreArray(user->sol_phi.u,     &u_phi);
    575:   VecRestoreArray(user->sol_phi.v,     &v_phi);
    576:   VecRestoreArray(user->sol_np1.rho_u, &rho_u_np1);
    577:   VecRestoreArray(user->sol_np1.rho_v, &rho_v_np1);

    579:   VecAssemblyBegin(rhs_m);
    580:   VecAssemblyBegin(rhs_e);
    581:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
    582:   VecAssemblyEnd(rhs_m);
    583:   VecAssemblyEnd(rhs_e);
    584:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
    585:   VecScale(rhs_m, user->dt);
    586:   VecScale(rhs_e, user->dt);

    588:   KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
    589:   KSPSolve(ksp, rhs_m, user->sol_np1.rho);
    590:   KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
    591:   KSPDestroy(&ksp);
    592:   MatDestroy(&mat);
    593:   DMRestoreGlobalVector(da, &rhs_m);
    594:   DMRestoreGlobalVector(da, &rhs_e);
    595:   return(0);
    596: }

    600: PetscErrorCode ComputePredictor(DM da, UserContext *user)
    601: {
    602:   Vec            uOldLocal, uLocal,uOld;
    603:   PetscScalar    *pOld;
    604:   PetscScalar    *p;

    608:   DMGetGlobalVector(da, &uOld);
    609:   DMGetLocalVector(da, &uOldLocal);
    610:   DMGetLocalVector(da, &uLocal);
    611:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
    612:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
    613:   VecGetArray(uOldLocal, &pOld);
    614:   VecGetArray(uLocal,    &p);

    616:   /* Source terms are all zero right now */
    617:   CalculateElementVelocity(da, user);
    618:   TaylorGalerkinStepI(da, user);
    619:   TaylorGalerkinStepIIMomentum(da, user);
    620:   TaylorGalerkinStepIIMassEnergy(da, user);
    621:   /* Solve equation (9) for \delta(\rho\vu) and (\rho\vu)^* */
    622:   /* Solve equation (13) for \delta\rho and \rho^* */
    623:   /* Solve equation (15) for \delta(\rho e_t) and (\rho e_t)^* */
    624:   /* Apply artifical dissipation */
    625:   /* Determine the smoothed explicit pressure, \tilde P and temperature \tilde T using the equation of state */

    628:   VecRestoreArray(uOldLocal, &pOld);
    629:   VecRestoreArray(uLocal,    &p);
    630: #if 0
    631:   DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
    632:   DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
    633:   DMRestoreLocalVector(da, &uOldLocal);
    634:   DMRestoreLocalVector(da, &uLocal);
    635: #endif
    636:   return(0);
    637: }

    641: /*
    642:   We integrate over each cell

    644:   (i, j+1)----(i+1, j+1)
    645:       | \         |
    646:       |  \        |
    647:       |   \       |
    648:       |    \      |
    649:       |     \     |
    650:       |      \    |
    651:       |       \   |
    652:   (i,   j)----(i+1, j)
    653: */
    654: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
    655: {
    656:   UserContext    *user = (UserContext*)ctx;
    657:   PetscScalar    phi   = user->phi;
    658:   PetscScalar    *array;
    659:   PetscInt       ne,nc,i;
    660:   const PetscInt *e;
    662:   Vec            blocal;
    663:   DM             da;

    666:   KSPGetDM(ksp,&da);
    667:   /* access a local vector with room for the ghost points */
    668:   DMGetLocalVector(da,&blocal);
    669:   VecGetArray(blocal, (PetscScalar**) &array);

    671:   /* access the list of elements on this processor and loop over them */
    672:   DMDAGetElements(da,&ne,&nc,&e);
    673:   for (i=0; i<ne; i++) {

    675:     /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
    676:     array[e[3*i]]   = phi;
    677:     array[e[3*i+1]] = phi;
    678:     array[e[3*i+2]] = phi;
    679:   }
    680:   VecRestoreArray(blocal, (PetscScalar**) &array);
    681:   DMDARestoreElements(da,&ne,&nc,&e);

    683:   /* add our partial sums over all processors into b */
    684:   DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
    685:   DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
    686:   DMRestoreLocalVector(da,&blocal);
    687:   return(0);
    688: }

    692: /*
    693:   We integrate over each cell

    695:   (i, j+1)----(i+1, j+1)
    696:       | \         |
    697:       |  \        |
    698:       |   \       |
    699:       |    \      |
    700:       |     \     |
    701:       |      \    |
    702:       |       \   |
    703:   (i,   j)----(i+1, j)

    705: However, the element stiffness matrix for the identity in linear elements is

    707:   1  /2 1 1\
    708:   -  |1 2 1|
    709:   12 \1 1 2/

    711: no matter what the shape of the triangle. The Laplacian stiffness matrix is

    713:   1  /         (x_2 - x_1)^2 + (y_2 - y_1)^2           -(x_2 - x_0)(x_2 - x_1) - (y_2 - y_1)(y_2 - y_0)  (x_1 - x_0)(x_2 - x_1) + (y_1 - y_0)(y_2 - y_1)\
    714:   -  |-(x_2 - x_0)(x_2 - x_1) - (y_2 - y_1)(y_2 - y_0)           (x_2 - x_0)^2 + (y_2 - y_0)^2          -(x_1 - x_0)(x_2 - x_0) - (y_1 - y_0)(y_2 - y_0)|
    715:   A  \ (x_1 - x_0)(x_2 - x_1) + (y_1 - y_0)(y_2 - y_1) -(x_1 - x_0)(x_2 - x_0) - (y_1 - y_0)(y_2 - y_0)           (x_1 - x_0)^2 + (y_1 - y_0)^2         /

    717: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
    718: */
    719: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, MatStructure *flag,void *ctx)
    720: {
    721:   UserContext *user = (UserContext*)ctx;
    722:   /* not being used!
    723:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
    724:                                 0.08333333333, 0.16666666667, 0.08333333333,
    725:                                 0.08333333333, 0.08333333333, 0.16666666667};
    726:   */
    727:   PetscScalar    values[3][3];
    728:   PetscInt       idx[3];
    729:   PetscScalar    hx, hy, hx2, hy2, area,phi_dt2;
    730:   PetscInt       i,mx,my,xm,ym,xs,ys;
    731:   PetscInt       ne,nc;
    732:   const PetscInt *e;
    734:   DM             da;

    737:   KSPGetDM(ksp,&da);
    738:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
    739:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
    740:   hx      = 1.0 / (mx-1);
    741:   hy      = 1.0 / (my-1);
    742:   area    = 0.5*hx*hy;
    743:   phi_dt2 = user->phi*user->dt*user->dt;
    744:   hx2     = hx*hx/area*phi_dt2;
    745:   hy2     = hy*hy/area*phi_dt2;

    747:   /* initially all elements have identical geometry so all element stiffness are identical */
    748:   values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
    749:   values[1][0] = -hy2;      values[1][1] = hy2;  values[1][2] = 0.0;
    750:   values[2][0] = -hx2;      values[2][1] = 0.0;  values[2][2] = hx2;

    752:   DMDAGetElements(da,&ne,&nc,&e);
    753:   for (i=0; i<ne; i++) {
    754:     idx[0] = e[3*i];
    755:     idx[1] = e[3*i+1];
    756:     idx[2] = e[3*i+2];
    757:     MatSetValuesLocal(jac,3,idx,3,idx,(PetscScalar*)values,ADD_VALUES);
    758:   }
    759:   DMDARestoreElements(da,&ne,&nc,&e);
    760:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
    761:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
    762:   return(0);
    763: }

    767: PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u)
    768: {
    769:   Vec            uOldLocal, uLocal;
    770:   PetscScalar    *cOld;
    771:   PetscScalar    *c;
    772:   PetscInt       i,ne,nc;
    773:   const PetscInt *e;

    777:   VecSet(u,0.0);
    778:   DMGetLocalVector(da, &uOldLocal);
    779:   DMGetLocalVector(da, &uLocal);
    780:   VecSet(uLocal,0.0);
    781:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
    782:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
    783:   VecGetArray(uOldLocal, &cOld);
    784:   VecGetArray(uLocal,    &c);

    786:   /* access the list of elements on this processor and loop over them */
    787:   DMDAGetElements(da,&ne,&nc,&e);
    788:   for (i=0; i<ne; i++) {

    790:     /* this is nonsense, but copy each nodal value*/
    791:     c[e[3*i]]   = cOld[e[3*i]];
    792:     c[e[3*i+1]] = cOld[e[3*i+1]];
    793:     c[e[3*i+2]] = cOld[e[3*i+2]];
    794:   }
    795:   DMDARestoreElements(da,&ne,&nc,&e);
    796:   VecRestoreArray(uOldLocal, &cOld);
    797:   VecRestoreArray(uLocal,    &c);
    798:   DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
    799:   DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
    800:   DMRestoreLocalVector(da, &uOldLocal);
    801:   DMRestoreLocalVector(da, &uLocal);
    802:   return(0);
    803: }