Actual source code: ex31.c

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

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

 12:   1) An explicit predictor step

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

 16:   2) An implicit corrector step

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

 20:   3) An explicit update step

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

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

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

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

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

 35: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

 37: #include <petscdmda.h>
 38: #include <petscksp.h>
 39: #include <petscpcmg.h>
 40: #include <petscdmmg.h>


 50: typedef struct {
 51:   Vec rho;     /* The mass solution \rho */
 52:   Vec rho_u;   /* The x-momentum solution \rho u */
 53:   Vec rho_v;   /* The y-momentum solution \rho v */
 54:   Vec rho_e;   /* The energy solution \rho e_t */
 55:   Vec p;       /* The pressure solution P */
 56:   Vec t;       /* The temperature solution T */
 57:   Vec u;       /* The x-velocity solution u */
 58:   Vec v;       /* The y-velocity solution v */
 59: } SolutionContext;

 61: typedef struct {
 62:   SolutionContext sol_n;   /* The solution at time t^n */
 63:   SolutionContext sol_phi; /* The element-averaged solution at time t^{n+\phi} */
 64:   SolutionContext sol_np1; /* The solution at time t^{n+1} */
 65:   Vec             mu;      /* The dynamic viscosity \mu(T) at time n */
 66:   Vec             kappa;   /* The thermal conductivity \kappa(T) at time n */
 67:   PetscScalar     phi;     /* The time weighting parameter */
 68:   PetscScalar     dt;      /* The timestep \Delta t */
 69: } UserContext;

 73: int main(int argc,char **argv)
 74: {
 75:   DMMG           *dmmg;
 76:   DM             da;
 77:   UserContext    user;
 79:   PetscInt       l;

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

 83:   DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
 84:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 85:   DMMGSetDM(dmmg,(DM)da);
 86:   DMDestroy(&da);
 87:   for (l = 0; l < DMMGGetLevels(dmmg); l++) {
 88:     DMMGSetUser(dmmg,l,&user);
 89:   }

 91:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PCICE", "DMMG");
 92:     user.phi = 0.5;
 93:     PetscOptionsScalar("-phi", "The time weighting parameter", "ex31.c", user.phi, &user.phi, PETSC_NULL);
 94:     user.dt  = 0.1;
 95:     PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, PETSC_NULL);
 96:   PetscOptionsEnd();

 98:   CreateStructures(DMMGGetFine(dmmg));
 99:   ComputeInitialGuess(DMMGGetFine(dmmg));
100:   ComputePredictor(DMMGGetFine(dmmg));

102:   DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
103:   DMMGSetInitialGuess(dmmg, DMMGInitialGuessCurrent);
104:   DMMGSolve(dmmg);

106:   ComputeCorrector(DMMGGetFine(dmmg), DMMGGetx(dmmg), DMMGGetr(dmmg));

108:   DestroyStructures(DMMGGetFine(dmmg));
109:   DMMGDestroy(dmmg);
110:   PetscFinalize();
111:   return 0;
112: }

116: PetscErrorCode CreateStructures(DMMG dmmg)
117: {
118:   DM              da   = dmmg->dm;
119:   UserContext    *user = (UserContext *) dmmg->user;
120:   const PetscInt *necon;
121:   PetscInt        ne,nc;
122:   PetscErrorCode  ierr;

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

154: PetscErrorCode DestroyStructures(DMMG dmmg)
155: {
156:   UserContext   *user = (UserContext *) dmmg->user;

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

185: PetscErrorCode ComputeInitialGuess(DMMG dmmg)
186: {
188:   return(0);
189: }

193: /* Average the velocity (u,v) at time t^n over each element for time n+\phi */
194: PetscErrorCode CalculateElementVelocity(DM da, UserContext *user)
195: {
196:   PetscScalar    *u_n,   *v_n;
197:   PetscScalar    *u_phi, *v_phi;
198:   const PetscInt *necon;
199:   PetscInt        j, e, ne, nc;
200:   PetscErrorCode  ierr;

203:   DMDAGetElements(da, &ne, &nc, &necon);
204:   VecGetArray(user->sol_n.u, &u_n);
205:   VecGetArray(user->sol_n.v, &v_n);
206:   PetscMalloc(ne*sizeof(PetscScalar),&u_phi);
207:   PetscMalloc(ne*sizeof(PetscScalar),&v_phi);
208:   for(e = 0; e < ne; e++) {
209:     u_phi[e] = 0.0;
210:     v_phi[e] = 0.0;
211:     for(j = 0; j < 3; j++) {
212:       u_phi[e] += u_n[necon[3*e+j]];
213:       v_phi[e] += v_n[necon[3*e+j]];
214:     }
215:     u_phi[e] /= 3.0;
216:     v_phi[e] /= 3.0;
217:   }
218:   PetscFree(u_phi);
219:   PetscFree(v_phi);
220:   DMDARestoreElements(da, &ne,&nc, &necon);
221:   VecRestoreArray(user->sol_n.u, &u_n);
222:   VecRestoreArray(user->sol_n.v, &v_n);
223:   return(0);
224: }

228: /* This is equation 32,

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

232: which is really simple for linear elements

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

236: where

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

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

242:    f^n = {\rho u  \rho u^2  \rho uv}^n      g^n = {\rho v  \rho uv  \rho v^2}^n
243: */
244: PetscErrorCode TaylorGalerkinStepI(DM da, UserContext *user)
245: {
246:   PetscScalar     phi_dt = user->phi*user->dt;
247:   PetscScalar    *u_n,     *v_n;
248:   PetscScalar    *rho_n,   *rho_u_n,   *rho_v_n;
249:   PetscScalar    *rho_phi, *rho_u_phi, *rho_v_phi;
250:   PetscScalar     Fx_x, Fy_y;
251:   PetscScalar     psi_x[3], psi_y[3];
252:   PetscInt        idx[3];
253:   PetscReal       hx, hy;
254:   const PetscInt *necon;
255:   PetscInt        j, e, ne, nc, mx, my;
256:   PetscErrorCode  ierr;

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

329: /*
330: The element stiffness matrix for the identity in linear elements is

332:   1  /2 1 1\
333:   -  |1 2 1|
334:   12 \1 1 2/

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

360:   PetscObjectGetComm((PetscObject) da, &comm);
361:   DMGetMatrix(da, MATAIJ, &mat);
362:   DMGetGlobalVector(da, &rhs_u);
363:   DMGetGlobalVector(da, &rhs_v);
364:   KSPCreate(comm, &ksp);
365:   KSPSetFromOptions(ksp);

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

434:   VecAssemblyBegin(rhs_u);
435:   VecAssemblyBegin(rhs_v);
436:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
437:   VecAssemblyEnd(rhs_u);
438:   VecAssemblyEnd(rhs_v);
439:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
440:   VecScale(rhs_u,user->dt);
441:   VecScale(rhs_v,user->dt);

443:   KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
444:   KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
445:   KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
446:   KSPDestroy(&ksp);
447:   MatDestroy(&mat);
448:   DMRestoreGlobalVector(da, &rhs_u);
449:   DMRestoreGlobalVector(da, &rhs_v);
450:   return(0);
451: }

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

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

459:   1  /2 1 1\
460:   -  |1 2 1|
461:   12 \1 1 2/

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

489:   PetscObjectGetComm((PetscObject) da, &comm);
490:   DMGetMatrix(da, MATAIJ, &mat);
491:   DMGetGlobalVector(da, &rhs_m);
492:   DMGetGlobalVector(da, &rhs_e);
493:   KSPCreate(comm, &ksp);
494:   KSPSetFromOptions(ksp);

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

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

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

587:   VecAssemblyBegin(rhs_m);
588:   VecAssemblyBegin(rhs_e);
589:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
590:   VecAssemblyEnd(rhs_m);
591:   VecAssemblyEnd(rhs_e);
592:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
593:   VecScale(rhs_m, user->dt);
594:   VecScale(rhs_e, user->dt);

596:   KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
597:   KSPSolve(ksp, rhs_m, user->sol_np1.rho);
598:   KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
599:   KSPDestroy(&ksp);
600:   MatDestroy(&mat);
601:   DMRestoreGlobalVector(da, &rhs_m);
602:   DMRestoreGlobalVector(da, &rhs_e);
603:   return(0);
604: }

608: PetscErrorCode ComputePredictor(DMMG dmmg)
609: {
610:   DM             da   = dmmg->dm;
611:   UserContext   *user = (UserContext *) dmmg->user;
612:   Vec            uOldLocal, uLocal,uOld;
613:   PetscScalar   *pOld;
614:   PetscScalar   *p;
616: 
618:   DMGetGlobalVector(da, &uOld);
619:   DMGetLocalVector(da, &uOldLocal);
620:   DMGetLocalVector(da, &uLocal);
621:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
622:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
623:   VecGetArray(uOldLocal, &pOld);
624:   VecGetArray(uLocal,    &p);

626:   /* Source terms are all zero right now */
627:   CalculateElementVelocity(da, user);
628:   TaylorGalerkinStepI(da, user);
629:   TaylorGalerkinStepIIMomentum(da, user);
630:   TaylorGalerkinStepIIMassEnergy(da, user);
631:   /* Solve equation (9) for \delta(\rho\vu) and (\rho\vu)^* */
632:   /* Solve equation (13) for \delta\rho and \rho^* */
633:   /* Solve equation (15) for \delta(\rho e_t) and (\rho e_t)^* */
634:   /* Apply artifical dissipation */
635:   /* Determine the smoothed explicit pressure, \tilde P and temperature \tilde T using the equation of state */


638:   VecRestoreArray(uOldLocal, &pOld);
639:   VecRestoreArray(uLocal,    &p);
640: #if 0
641:   DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
642:   DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
643:   DMRestoreLocalVector(da, &uOldLocal);
644:   DMRestoreLocalVector(da, &uLocal);
645: #endif
646:   return(0);
647: }

651: /*
652:   We integrate over each cell

654:   (i, j+1)----(i+1, j+1)
655:       | \         |
656:       |  \        |
657:       |   \       |
658:       |    \      |
659:       |     \     |
660:       |      \    |
661:       |       \   |
662:   (i,   j)----(i+1, j)
663: */
664: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
665: {
666:   DM             da   = dmmg->dm;
667:   UserContext   *user = (UserContext *) dmmg->user;
668:   PetscScalar    phi  = user->phi;
669:   PetscScalar   *array;
670:   PetscInt       ne,nc,i;
671:   const PetscInt *e;
673:   Vec            blocal;

676:   /* access a local vector with room for the ghost points */
677:   DMGetLocalVector(da,&blocal);
678:   VecGetArray(blocal, (PetscScalar **) &array);

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

684:     /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
685:     array[e[3*i]]   = phi;
686:     array[e[3*i+1]] = phi;
687:     array[e[3*i+2]] = phi;
688:   }
689:   VecRestoreArray(blocal, (PetscScalar **) &array);
690:   DMDARestoreElements(da,&ne,&nc,&e);

692:   /* add our partial sums over all processors into b */
693:   DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
694:   DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
695:   DMRestoreLocalVector(da,&blocal);
696:   return(0);
697: }

701: /*
702:   We integrate over each cell

704:   (i, j+1)----(i+1, j+1)
705:       | \         |
706:       |  \        |
707:       |   \       |
708:       |    \      |
709:       |     \     |
710:       |      \    |
711:       |       \   |
712:   (i,   j)----(i+1, j)

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

716:   1  /2 1 1\
717:   -  |1 2 1|
718:   12 \1 1 2/

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

722:   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)\
723:   -  |-(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)|
724:   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         /

726: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
727: */
728: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J,Mat jac)
729: {
730:   DM             da   =  dmmg->dm;
731:   UserContext   *user = (UserContext *) dmmg->user;
732:   /* not being used!
733:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
734:                                 0.08333333333, 0.16666666667, 0.08333333333,
735:                                 0.08333333333, 0.08333333333, 0.16666666667};
736:   */
737:   PetscScalar    values[3][3];
738:   PetscInt       idx[3];
739:   PetscScalar    hx, hy, hx2, hy2, area,phi_dt2;
740:   PetscInt       i,mx,my,xm,ym,xs,ys;
741:   PetscInt       ne,nc;
742:   const PetscInt *e;

746:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
747:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
748:   hx   = 1.0 / (mx-1);
749:   hy   = 1.0 / (my-1);
750:   area = 0.5*hx*hy;
751:   phi_dt2 = user->phi*user->dt*user->dt;
752:   hx2     = hx*hx/area*phi_dt2;
753:   hy2     = hy*hy/area*phi_dt2;

755:   /* initially all elements have identical geometry so all element stiffness are identical */
756:   values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
757:   values[1][0] = -hy2;      values[1][1] = hy2;  values[1][2] = 0.0;
758:   values[2][0] = -hx2;      values[2][1] = 0.0;  values[2][2] = hx2;

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

775: PetscErrorCode ComputeCorrector(DMMG dmmg, Vec uOld, Vec u)
776: {
777:   DM             da   = dmmg->dm;
778:   Vec            uOldLocal, uLocal;
779:   PetscScalar    *cOld;
780:   PetscScalar    *c;
781:   PetscInt       i,ne,nc;
782:   const PetscInt *e;
784: 
786:   VecSet(u,0.0);
787:   DMGetLocalVector(da, &uOldLocal);
788:   DMGetLocalVector(da, &uLocal);
789:   VecSet(uLocal,0.0);
790:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
791:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
792:   VecGetArray(uOldLocal, &cOld);
793:   VecGetArray(uLocal,    &c);

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

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