Actual source code: ex31.c

petsc-master 2016-08-26
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 example is WRONG. It, for example, creates the global vector user->sol_n.rho but then accesses values in it using the element numbering
 32:   which is based on local (ghosted) numbering.

 34:   This example does not produce any meaningful results because it is incomplete; the systems produced may not be reasonable.

 36: */

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

 40:  #include <petscdm.h>
 41:  #include <petscdmda.h>
 42:  #include <petscksp.h>

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

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

 65: extern PetscErrorCode CreateStructures(DM,UserContext*);
 66: extern PetscErrorCode DestroyStructures(DM,UserContext*);
 67: extern PetscErrorCode ComputePredictor(DM,UserContext*);
 68: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
 69: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
 70: extern PetscErrorCode ComputeCorrector(DM,Vec,Vec);

 74: int main(int argc,char **argv)
 75: {
 76:   KSP            ksp;
 77:   DM             da;
 78:   Vec            x, xNew;
 79:   UserContext    user;

 82:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 83:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 84:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 85:   DMDASetElementType(da,DMDA_ELEMENT_P1);
 86:   DMSetApplicationContext(da, &user);
 87:   KSPSetDM(ksp, da);

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

 96:   CreateStructures(da, &user);
 97:   ComputePredictor(da, &user);

 99:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
100:   KSPSetComputeOperators(ksp,ComputeMatrix,&user);
101:   KSPSetFromOptions(ksp);
102:   KSPSolve(ksp, NULL, NULL);

104:   KSPGetSolution(ksp, &x);
105:   VecDuplicate(x, &xNew);
106:   ComputeCorrector(da, x, xNew);
107:   VecDestroy(&xNew);

109:   DestroyStructures(da, &user);
110:   DMDestroy(&da);
111:   KSPDestroy(&ksp);
112:   PetscFinalize();
113:   return ierr;
114: }

118: PetscErrorCode CreateStructures(DM da, UserContext *user)
119: {
120:   const PetscInt *necon;
121:   PetscInt       ne,nc;

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:   DMCreateGlobalVector(da, &user->sol_n.t);
135:   VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho);
136:   VecSetSizes(user->sol_phi.rho, ne, PETSC_DECIDE);
137:   VecSetType(user->sol_phi.rho,VECMPI);
138:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_u);
139:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_v);
140:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.u);
141:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.v);
142:   DMCreateGlobalVector(da, &user->sol_np1.rho);
143:   DMCreateGlobalVector(da, &user->sol_np1.rho_u);
144:   DMCreateGlobalVector(da, &user->sol_np1.rho_v);
145:   DMCreateGlobalVector(da, &user->sol_np1.rho_e);
146:   DMCreateGlobalVector(da, &user->sol_np1.p);
147:   DMCreateGlobalVector(da, &user->sol_np1.u);
148:   DMCreateGlobalVector(da, &user->sol_np1.v);
149:   DMCreateGlobalVector(da, &user->mu);
150:   DMCreateGlobalVector(da, &user->kappa);
151:   return(0);
152: }

156: PetscErrorCode DestroyStructures(DM da, UserContext   *user)
157: {

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

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

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

223: /* This is equation 32,

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

227: which is really simple for linear elements

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

231: where

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

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

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

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

324: /*
325: The element stiffness matrix for the identity in linear elements is

327:   1  /2 1 1\
328:   -  |1 2 1|
329:   12 \1 1 2/

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

355:   PetscObjectGetComm((PetscObject) da, &comm);
356:   DMSetMatType(da,MATAIJ);
357:   DMCreateMatrix(da, &mat);
358:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
359:   MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
360:   DMGetGlobalVector(da, &rhs_u);
361:   DMGetGlobalVector(da, &rhs_v);
362:   KSPCreate(comm, &ksp);
363:   KSPSetFromOptions(ksp);

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

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

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

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

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

457:   1  /2 1 1\
458:   -  |1 2 1|
459:   12 \1 1 2/

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

487:   PetscObjectGetComm((PetscObject) da, &comm);
488:   DMSetMatType(da,MATAIJ);
489:   DMCreateMatrix(da, &mat);
490:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
491:   MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
492:   DMGetGlobalVector(da, &rhs_m);
493:   DMGetGlobalVector(da, &rhs_e);
494:   KSPCreate(comm, &ksp);
495:   KSPSetFromOptions(ksp);

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

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

546:       mu     = 0.0;
547:       kappa  = 0.0;
548:       tau_xx = 0.0;
549:       tau_xy = 0.0;
550:       tau_yy = 0.0;
551:       for (k = 0; k < 3; k++) {
552:         mu     += mu_n[idx[k]];
553:         kappa  += kappa_n[idx[k]];
554:         tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
555:         tau_xy +=     psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
556:         tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
557:         q_x    += psi_x[k]*t_n[idx[k]];
558:         q_y    += psi_y[k]*t_n[idx[k]];
559:       }
560:       mu          /= 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);
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(DM da, UserContext *user)
609: {
610:   Vec            uOldLocal, uLocal,uOld;
611:   PetscScalar    *pOld;
612:   PetscScalar    *p;

616:   DMGetGlobalVector(da, &uOld);
617:   DMGetLocalVector(da, &uOldLocal);
618:   DMGetLocalVector(da, &uLocal);
619:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
620:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
621:   VecGetArray(uOldLocal, &pOld);
622:   VecGetArray(uLocal,    &p);

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


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

650: /*
651:   We integrate over each cell

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

675:   KSPGetDM(ksp,&da);
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(KSP ksp, Mat J, Mat jac, void *ctx)
729: {
730:   UserContext *user = (UserContext*)ctx;
731:   /* not being used!
732:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
733:                                 0.08333333333, 0.16666666667, 0.08333333333,
734:                                 0.08333333333, 0.08333333333, 0.16666666667};
735:   */
736:   PetscScalar    values[3][3];
737:   PetscInt       idx[3];
738:   PetscScalar    hx, hy, hx2, hy2, area,phi_dt2;
739:   PetscInt       i,mx,my,xm,ym,xs,ys;
740:   PetscInt       ne,nc;
741:   const PetscInt *e;
743:   DM             da;

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

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

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

776: PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u)
777: {
778:   Vec            uOldLocal, uLocal;
779:   PetscScalar    *cOld;
780:   PetscScalar    *c;
781:   PetscInt       i,ne,nc;
782:   const PetscInt *e;

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