Actual source code: ex31.c

petsc-master 2015-07-30
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);

 84:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 85:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 86:   DMDASetElementType(da,DMDA_ELEMENT_P1);
 87:   DMSetApplicationContext(da, &user);
 88:   KSPSetDM(ksp, da);

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

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

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

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

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

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

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

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

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

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

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

224: /* This is equation 32,

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

228: which is really simple for linear elements

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

232: where

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

547:       mu     = 0.0;
548:       kappa  = 0.0;
549:       tau_xx = 0.0;
550:       tau_xy = 0.0;
551:       tau_yy = 0.0;
552:       for (k = 0; k < 3; k++) {
553:         mu     += mu_n[idx[k]];
554:         kappa  += kappa_n[idx[k]];
555:         tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
556:         tau_xy +=     psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
557:         tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
558:         q_x    += psi_x[k]*t_n[idx[k]];
559:         q_y    += psi_y[k]*t_n[idx[k]];
560:       }
561:       mu          /= 3.0;
562:       kappa       /= 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: }