Actual source code: ex31.c

petsc-master 2014-10-20
Report Typos and Errors
  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 <petscdm.h>
 37: #include <petscdmda.h>
 38: #include <petscksp.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

219: /* This is equation 32,

221:    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

223: which is really simple for linear elements

225:    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

227: where

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

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

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

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

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

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

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

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

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

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

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

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

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

452:   1  /2 1 1\
453:   -  |1 2 1|
454:   12 \1 1 2/

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

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

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

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

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

582:   VecAssemblyBegin(rhs_m);
583:   VecAssemblyBegin(rhs_e);
584:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
585:   VecAssemblyEnd(rhs_m);
586:   VecAssemblyEnd(rhs_e);
587:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
588:   VecScale(rhs_m, user->dt);
589:   VecScale(rhs_e, user->dt);

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

603: PetscErrorCode ComputePredictor(DM da, UserContext *user)
604: {
605:   Vec            uOldLocal, uLocal,uOld;
606:   PetscScalar    *pOld;
607:   PetscScalar    *p;

611:   DMGetGlobalVector(da, &uOld);
612:   DMGetLocalVector(da, &uOldLocal);
613:   DMGetLocalVector(da, &uLocal);
614:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
615:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
616:   VecGetArray(uOldLocal, &pOld);
617:   VecGetArray(uLocal,    &p);

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


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

644: /*
645:   We integrate over each cell

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

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

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

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

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

695: /*
696:   We integrate over each cell

698:   (i, j+1)----(i+1, j+1)
699:       | \         |
700:       |  \        |
701:       |   \       |
702:       |    \      |
703:       |     \     |
704:       |      \    |
705:       |       \   |
706:   (i,   j)----(i+1, j)

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

710:   1  /2 1 1\
711:   -  |1 2 1|
712:   12 \1 1 2/

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

716:   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)\
717:   -  |-(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)|
718:   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         /

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

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

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

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

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

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

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

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