Actual source code: ex31.c

petsc-3.4.5 2014-06-29
  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
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 */
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];
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);