Actual source code: ex31.c

petsc-master 2016-12-06
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:   DMSetFromOptions(da);
 86:   DMSetUp(da);
 87:   DMDASetElementType(da,DMDA_ELEMENT_P1);
 88:   DMSetApplicationContext(da, &user);
 89:   KSPSetDM(ksp, da);

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

 98:   CreateStructures(da, &user);
 99:   ComputePredictor(da, &user);

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

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

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

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

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

158: PetscErrorCode DestroyStructures(DM da, UserContext   *user)
159: {

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

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

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

225: /* This is equation 32,

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

229: which is really simple for linear elements

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

233: where

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

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

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

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

326: /*
327: The element stiffness matrix for the identity in linear elements is

329:   1  /2 1 1\
330:   -  |1 2 1|
331:   12 \1 1 2/

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

357:   PetscObjectGetComm((PetscObject) da, &comm);
358:   DMSetMatType(da,MATAIJ);
359:   DMCreateMatrix(da, &mat);
360:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
361:   MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
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);
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;

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

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

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

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

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

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

610: PetscErrorCode ComputePredictor(DM da, UserContext *user)
611: {
612:   Vec            uOldLocal, uLocal,uOld;
613:   PetscScalar    *pOld;
614:   PetscScalar    *p;

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: #endif
644:   DMRestoreLocalVector(da, &uOldLocal);
645:   DMRestoreLocalVector(da, &uLocal);
646:   DMRestoreGlobalVector(da, &uOld);
647:   return(0);
648: }

652: /*
653:   We integrate over each cell

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

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

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

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

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

703: /*
704:   We integrate over each cell

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

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

718:   1  /2 1 1\
719:   -  |1 2 1|
720:   12 \1 1 2/

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

724:   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)\
725:   -  |-(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)|
726:   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         /

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

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

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

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

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

788:   VecSet(u,0.0);
789:   DMGetLocalVector(da, &uOldLocal);
790:   DMGetLocalVector(da, &uLocal);
791:   VecSet(uLocal,0.0);
792:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
793:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
794:   VecGetArray(uOldLocal, &cOld);
795:   VecGetArray(uLocal,    &c);

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

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