Actual source code: ex31.c

petsc-master 2017-11-23
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);

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

 80:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 81:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 82:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 83:   DMSetFromOptions(da);
 84:   DMSetUp(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: }

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

123:   DMDAGetElements(da,&ne,&nc,&necon);
124:   DMDARestoreElements(da,&ne,&nc,&necon);
125:   DMCreateGlobalVector(da, &user->sol_n.rho);
126:   DMCreateGlobalVector(da, &user->sol_n.rho_u);
127:   DMCreateGlobalVector(da, &user->sol_n.rho_v);
128:   DMCreateGlobalVector(da, &user->sol_n.rho_e);
129:   DMCreateGlobalVector(da, &user->sol_n.p);
130:   DMCreateGlobalVector(da, &user->sol_n.u);
131:   DMCreateGlobalVector(da, &user->sol_n.v);
132:   DMCreateGlobalVector(da, &user->sol_n.t);
133:   VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho);
134:   VecSetSizes(user->sol_phi.rho, ne, PETSC_DECIDE);
135:   VecSetType(user->sol_phi.rho,VECMPI);
136:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_u);
137:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_v);
138:   VecDuplicate(user->sol_phi.rho, &user->sol_phi.u);
139:   VecDuplicate(user->sol_phi.rho, &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: }

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

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

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

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

215: /* This is equation 32,

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

219: which is really simple for linear elements

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

223: where

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

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

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

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

314: /*
315: The element stiffness matrix for the identity in linear elements is

317:   1  /2 1 1\
318:   -  |1 2 1|
319:   12 \1 1 2/

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

345:   PetscObjectGetComm((PetscObject) da, &comm);
346:   DMSetMatType(da,MATAIJ);
347:   DMCreateMatrix(da, &mat);
348:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
349:   MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
350:   DMGetGlobalVector(da, &rhs_u);
351:   DMGetGlobalVector(da, &rhs_v);
352:   KSPCreate(comm, &ksp);
353:   KSPSetFromOptions(ksp);

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

422:   VecAssemblyBegin(rhs_u);
423:   VecAssemblyBegin(rhs_v);
424:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
425:   VecAssemblyEnd(rhs_u);
426:   VecAssemblyEnd(rhs_v);
427:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
428:   VecScale(rhs_u,user->dt);
429:   VecScale(rhs_v,user->dt);

431:   KSPSetOperators(ksp, mat, mat);
432:   KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
433:   KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
434:   KSPDestroy(&ksp);
435:   MatDestroy(&mat);
436:   DMRestoreGlobalVector(da, &rhs_u);
437:   DMRestoreGlobalVector(da, &rhs_v);
438:   return(0);
439: }

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

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

445:   1  /2 1 1\
446:   -  |1 2 1|
447:   12 \1 1 2/

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

475:   PetscObjectGetComm((PetscObject) da, &comm);
476:   DMSetMatType(da,MATAIJ);
477:   DMCreateMatrix(da, &mat);
478:   MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
479:   MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
480:   DMGetGlobalVector(da, &rhs_m);
481:   DMGetGlobalVector(da, &rhs_e);
482:   KSPCreate(comm, &ksp);
483:   KSPSetFromOptions(ksp);

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

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

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

575:   VecAssemblyBegin(rhs_m);
576:   VecAssemblyBegin(rhs_e);
577:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
578:   VecAssemblyEnd(rhs_m);
579:   VecAssemblyEnd(rhs_e);
580:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
581:   VecScale(rhs_m, user->dt);
582:   VecScale(rhs_e, user->dt);

584:   KSPSetOperators(ksp, mat, mat);
585:   KSPSolve(ksp, rhs_m, user->sol_np1.rho);
586:   KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
587:   KSPDestroy(&ksp);
588:   MatDestroy(&mat);
589:   DMRestoreGlobalVector(da, &rhs_m);
590:   DMRestoreGlobalVector(da, &rhs_e);
591:   return(0);
592: }

594: PetscErrorCode ComputePredictor(DM da, UserContext *user)
595: {
596:   Vec            uOldLocal, uLocal,uOld;
597:   PetscScalar    *pOld;
598:   PetscScalar    *p;

602:   DMGetGlobalVector(da, &uOld);
603:   DMGetLocalVector(da, &uOldLocal);
604:   DMGetLocalVector(da, &uLocal);
605:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
606:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
607:   VecGetArray(uOldLocal, &pOld);
608:   VecGetArray(uLocal,    &p);

610:   /* Source terms are all zero right now */
611:   CalculateElementVelocity(da, user);
612:   TaylorGalerkinStepI(da, user);
613:   TaylorGalerkinStepIIMomentum(da, user);
614:   TaylorGalerkinStepIIMassEnergy(da, user);
615:   /* Solve equation (9) for \delta(\rho\vu) and (\rho\vu)^* */
616:   /* Solve equation (13) for \delta\rho and \rho^* */
617:   /* Solve equation (15) for \delta(\rho e_t) and (\rho e_t)^* */
618:   /* Apply artifical dissipation */
619:   /* Determine the smoothed explicit pressure, \tilde P and temperature \tilde T using the equation of state */


622:   VecRestoreArray(uOldLocal, &pOld);
623:   VecRestoreArray(uLocal,    &p);
624: #if 0
625:   DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
626:   DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
627: #endif
628:   DMRestoreLocalVector(da, &uOldLocal);
629:   DMRestoreLocalVector(da, &uLocal);
630:   DMRestoreGlobalVector(da, &uOld);
631:   return(0);
632: }

634: /*
635:   We integrate over each cell

637:   (i, j+1)----(i+1, j+1)
638:       | \         |
639:       |  \        |
640:       |   \       |
641:       |    \      |
642:       |     \     |
643:       |      \    |
644:       |       \   |
645:   (i,   j)----(i+1, j)
646: */
647: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
648: {
649:   UserContext    *user = (UserContext*)ctx;
650:   PetscScalar    phi   = user->phi;
651:   PetscScalar    *array;
652:   PetscInt       ne,nc,i;
653:   const PetscInt *e;
655:   Vec            blocal;
656:   DM             da;

659:   KSPGetDM(ksp,&da);
660:   /* access a local vector with room for the ghost points */
661:   DMGetLocalVector(da,&blocal);
662:   VecGetArray(blocal, (PetscScalar**) &array);

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

668:     /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
669:     array[e[3*i]]   = phi;
670:     array[e[3*i+1]] = phi;
671:     array[e[3*i+2]] = phi;
672:   }
673:   VecRestoreArray(blocal, (PetscScalar**) &array);
674:   DMDARestoreElements(da,&ne,&nc,&e);

676:   /* add our partial sums over all processors into b */
677:   DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
678:   DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
679:   DMRestoreLocalVector(da,&blocal);
680:   return(0);
681: }

683: /*
684:   We integrate over each cell

686:   (i, j+1)----(i+1, j+1)
687:       | \         |
688:       |  \        |
689:       |   \       |
690:       |    \      |
691:       |     \     |
692:       |      \    |
693:       |       \   |
694:   (i,   j)----(i+1, j)

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

698:   1  /2 1 1\
699:   -  |1 2 1|
700:   12 \1 1 2/

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

704:   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)\
705:   -  |-(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)|
706:   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         /

708: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
709: */
710: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
711: {
712:   UserContext *user = (UserContext*)ctx;
713:   /* not being used!
714:   PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
715:                                 0.08333333333, 0.16666666667, 0.08333333333,
716:                                 0.08333333333, 0.08333333333, 0.16666666667};
717:   */
718:   PetscScalar    values[3][3];
719:   PetscInt       idx[3];
720:   PetscScalar    hx, hy, hx2, hy2, area,phi_dt2;
721:   PetscInt       i,mx,my,xm,ym,xs,ys;
722:   PetscInt       ne,nc;
723:   const PetscInt *e;
725:   DM             da;

728:   KSPGetDM(ksp,&da);
729:   DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
730:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
731:   hx      = 1.0 / (mx-1);
732:   hy      = 1.0 / (my-1);
733:   area    = 0.5*hx*hy;
734:   phi_dt2 = user->phi*user->dt*user->dt;
735:   hx2     = hx*hx/area*phi_dt2;
736:   hy2     = hy*hy/area*phi_dt2;

738:   /* initially all elements have identical geometry so all element stiffness are identical */
739:   values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
740:   values[1][0] = -hy2;      values[1][1] = hy2;  values[1][2] = 0.0;
741:   values[2][0] = -hx2;      values[2][1] = 0.0;  values[2][2] = hx2;

743:   DMDAGetElements(da,&ne,&nc,&e);
744:   for (i=0; i<ne; i++) {
745:     idx[0] = e[3*i];
746:     idx[1] = e[3*i+1];
747:     idx[2] = e[3*i+2];
748:     MatSetValuesLocal(jac,3,idx,3,idx,(PetscScalar*)values,ADD_VALUES);
749:   }
750:   DMDARestoreElements(da,&ne,&nc,&e);
751:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
752:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
753:   return(0);
754: }

756: PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u)
757: {
758:   Vec            uOldLocal, uLocal;
759:   PetscScalar    *cOld;
760:   PetscScalar    *c;
761:   PetscInt       i,ne,nc;
762:   const PetscInt *e;

766:   VecSet(u,0.0);
767:   DMGetLocalVector(da, &uOldLocal);
768:   DMGetLocalVector(da, &uLocal);
769:   VecSet(uLocal,0.0);
770:   DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
771:   DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
772:   VecGetArray(uOldLocal, &cOld);
773:   VecGetArray(uLocal,    &c);

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

779:     /* this is nonsense, but copy each nodal value*/
780:     c[e[3*i]]   = cOld[e[3*i]];
781:     c[e[3*i+1]] = cOld[e[3*i+1]];
782:     c[e[3*i+2]] = cOld[e[3*i+2]];
783:   }
784:   DMDARestoreElements(da,&ne,&nc,&e);
785:   VecRestoreArray(uOldLocal, &cOld);
786:   VecRestoreArray(uLocal,    &c);
787:   DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
788:   DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
789:   DMRestoreLocalVector(da, &uOldLocal);
790:   DMRestoreLocalVector(da, &uLocal);
791:   return(0);
792: }