Actual source code: ex31.c
petsc-3.4.5 2014-06-29
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 uses multigrid to solve the linear system
32: */
34: static char help[] = "Solves 2D compressible Euler using multigrid.\n\n";
36: #include <petscdmda.h>
37: #include <petscksp.h>
39: typedef struct {
40: Vec rho; /* The mass solution \rho */
41: Vec rho_u; /* The x-momentum solution \rho u */
42: Vec rho_v; /* The y-momentum solution \rho v */
43: Vec rho_e; /* The energy solution \rho e_t */
44: Vec p; /* The pressure solution P */
45: Vec t; /* The temperature solution T */
46: Vec u; /* The x-velocity solution u */
47: Vec v; /* The y-velocity solution v */
48: } SolutionContext;
50: typedef struct {
51: SolutionContext sol_n; /* The solution at time t^n */
52: SolutionContext sol_phi; /* The element-averaged solution at time t^{n+\phi} */
53: SolutionContext sol_np1; /* The solution at time t^{n+1} */
54: Vec mu; /* The dynamic viscosity \mu(T) at time n */
55: Vec kappa; /* The thermal conductivity \kappa(T) at time n */
56: PetscScalar phi; /* The time weighting parameter */
57: PetscScalar dt; /* The timestep \Delta t */
58: } UserContext;
60: extern PetscErrorCode CreateStructures(DM,UserContext*);
61: extern PetscErrorCode DestroyStructures(DM,UserContext*);
62: extern PetscErrorCode ComputePredictor(DM,UserContext*);
63: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
64: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
65: extern PetscErrorCode ComputeCorrector(DM,Vec,Vec);
69: int main(int argc,char **argv)
70: {
71: KSP ksp;
72: DM da;
73: Vec x, xNew;
74: UserContext user;
77: PetscInitialize(&argc,&argv,(char*)0,help);
79: KSPCreate(PETSC_COMM_WORLD,&ksp);
80: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
81: DMSetApplicationContext(da, &user);
82: KSPSetDM(ksp, da);
84: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PCICE", "DM");
85: user.phi = 0.5;
86: PetscOptionsScalar("-phi", "The time weighting parameter", "ex31.c", user.phi, &user.phi, NULL);
87: user.dt = 0.1;
88: PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, NULL);
89: PetscOptionsEnd();
91: CreateStructures(da, &user);
92: ComputePredictor(da, &user);
94: KSPSetComputeRHS(ksp,ComputeRHS,&user);
95: KSPSetComputeOperators(ksp,ComputeMatrix,&user);
96: KSPSetFromOptions(ksp);
97: KSPSolve(ksp, NULL, NULL);
99: KSPGetSolution(ksp, &x);
100: VecDuplicate(x, &xNew);
101: ComputeCorrector(da, x, xNew);
102: VecDestroy(&xNew);
104: DestroyStructures(da, &user);
105: DMDestroy(&da);
106: KSPDestroy(&ksp);
107: PetscFinalize();
108: return 0;
109: }
113: PetscErrorCode CreateStructures(DM da, UserContext *user)
114: {
115: const PetscInt *necon;
116: PetscInt ne,nc;
120: DMDAGetElements(da,&ne,&nc,&necon);
121: DMDARestoreElements(da,&ne,&nc,&necon);
122: DMCreateGlobalVector(da, &user->sol_n.rho);
123: DMCreateGlobalVector(da, &user->sol_n.rho_u);
124: DMCreateGlobalVector(da, &user->sol_n.rho_v);
125: DMCreateGlobalVector(da, &user->sol_n.rho_e);
126: DMCreateGlobalVector(da, &user->sol_n.p);
127: DMCreateGlobalVector(da, &user->sol_n.u);
128: DMCreateGlobalVector(da, &user->sol_n.v);
129: DMCreateGlobalVector(da, &user->sol_n.t);
130: VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho);
131: VecSetSizes(user->sol_phi.rho, ne, PETSC_DECIDE);
132: VecSetType(user->sol_phi.rho,VECMPI);
133: VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_u);
134: VecDuplicate(user->sol_phi.rho, &user->sol_phi.rho_v);
135: VecDuplicate(user->sol_phi.rho, &user->sol_phi.u);
136: VecDuplicate(user->sol_phi.rho, &user->sol_phi.v);
137: DMCreateGlobalVector(da, &user->sol_np1.rho);
138: DMCreateGlobalVector(da, &user->sol_np1.rho_u);
139: DMCreateGlobalVector(da, &user->sol_np1.rho_v);
140: DMCreateGlobalVector(da, &user->sol_np1.rho_e);
141: DMCreateGlobalVector(da, &user->sol_np1.p);
142: DMCreateGlobalVector(da, &user->sol_np1.u);
143: DMCreateGlobalVector(da, &user->sol_np1.v);
144: DMCreateGlobalVector(da, &user->mu);
145: DMCreateGlobalVector(da, &user->kappa);
146: return(0);
147: }
151: PetscErrorCode DestroyStructures(DM da, UserContext *user)
152: {
156: VecDestroy(&user->sol_n.rho);
157: VecDestroy(&user->sol_n.rho_u);
158: VecDestroy(&user->sol_n.rho_v);
159: VecDestroy(&user->sol_n.rho_e);
160: VecDestroy(&user->sol_n.p);
161: VecDestroy(&user->sol_n.u);
162: VecDestroy(&user->sol_n.v);
163: VecDestroy(&user->sol_n.t);
164: VecDestroy(&user->sol_phi.rho);
165: VecDestroy(&user->sol_phi.rho_u);
166: VecDestroy(&user->sol_phi.rho_v);
167: VecDestroy(&user->sol_phi.u);
168: VecDestroy(&user->sol_phi.v);
169: VecDestroy(&user->sol_np1.rho);
170: VecDestroy(&user->sol_np1.rho_u);
171: VecDestroy(&user->sol_np1.rho_v);
172: VecDestroy(&user->sol_np1.rho_e);
173: VecDestroy(&user->sol_np1.p);
174: VecDestroy(&user->sol_np1.u);
175: VecDestroy(&user->sol_np1.v);
176: VecDestroy(&user->mu);
177: VecDestroy(&user->kappa);
178: return(0);
179: }
183: /* Average the velocity (u,v) at time t^n over each element for time n+\phi */
184: PetscErrorCode CalculateElementVelocity(DM da, UserContext *user)
185: {
186: PetscScalar *u_n, *v_n;
187: PetscScalar *u_phi, *v_phi;
188: const PetscInt *necon;
189: PetscInt j, e, ne, nc;
193: DMDAGetElements(da, &ne, &nc, &necon);
194: VecGetArray(user->sol_n.u, &u_n);
195: VecGetArray(user->sol_n.v, &v_n);
196: PetscMalloc(ne*sizeof(PetscScalar),&u_phi);
197: PetscMalloc(ne*sizeof(PetscScalar),&v_phi);
198: for (e = 0; e < ne; e++) {
199: u_phi[e] = 0.0;
200: v_phi[e] = 0.0;
201: for (j = 0; j < 3; j++) {
202: u_phi[e] += u_n[necon[3*e+j]];
203: v_phi[e] += v_n[necon[3*e+j]];
204: }
205: u_phi[e] /= 3.0;
206: v_phi[e] /= 3.0;
207: }
208: PetscFree(u_phi);
209: PetscFree(v_phi);
210: DMDARestoreElements(da, &ne,&nc, &necon);
211: VecRestoreArray(user->sol_n.u, &u_n);
212: VecRestoreArray(user->sol_n.v, &v_n);
213: return(0);
214: }
218: /* This is equation 32,
220: 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
222: which is really simple for linear elements
224: 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
226: where
228: U^{n+\phi}_E = {\rho \rho u \rho v}^{n+\phi}_E
230: and the x and y components of the convective fluxes F are
232: f^n = {\rho u \rho u^2 \rho uv}^n g^n = {\rho v \rho uv \rho v^2}^n
233: */
234: PetscErrorCode TaylorGalerkinStepI(DM da, UserContext *user)
235: {
236: PetscScalar phi_dt = user->phi*user->dt;
237: PetscScalar *u_n, *v_n;
238: PetscScalar *rho_n, *rho_u_n, *rho_v_n;
239: PetscScalar *rho_phi, *rho_u_phi, *rho_v_phi;
240: PetscScalar Fx_x, Fy_y;
241: PetscScalar psi_x[3], psi_y[3];
242: PetscInt idx[3];
243: PetscReal hx, hy;
244: const PetscInt *necon;
245: PetscInt j, e, ne, nc, mx, my;
249: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
250: hx = 1.0 / (PetscReal)(mx-1);
251: hy = 1.0 / (PetscReal)(my-1);
252: VecSet(user->sol_phi.rho,0.0);
253: VecSet(user->sol_phi.rho_u,0.0);
254: VecSet(user->sol_phi.rho_v,0.0);
255: VecGetArray(user->sol_n.u, &u_n);
256: VecGetArray(user->sol_n.v, &v_n);
257: VecGetArray(user->sol_n.rho, &rho_n);
258: VecGetArray(user->sol_n.rho_u, &rho_u_n);
259: VecGetArray(user->sol_n.rho_v, &rho_v_n);
260: VecGetArray(user->sol_phi.rho, &rho_phi);
261: VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
262: VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
263: DMDAGetElements(da, &ne, &nc, &necon);
264: for (e = 0; e < ne; e++) {
265: /* Average the existing fields over the element */
266: for (j = 0; j < 3; j++) {
267: idx[j] = necon[3*e+j];
268: rho_phi[e] += rho_n[idx[j]];
269: rho_u_phi[e] += rho_u_n[idx[j]];
270: rho_v_phi[e] += rho_v_n[idx[j]];
271: }
272: rho_phi[e] /= 3.0;
273: rho_u_phi[e] /= 3.0;
274: rho_v_phi[e] /= 3.0;
275: /* Get basis function deriatives (we need the orientation of the element here) */
276: if (idx[1] > idx[0]) {
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: } else {
280: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
281: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
282: }
283: /* Determine the convective fluxes for \rho^{n+\phi} */
284: Fx_x = 0.0; Fy_y = 0.0;
285: for (j = 0; j < 3; j++) {
286: Fx_x += psi_x[j]*rho_u_n[idx[j]];
287: Fy_y += psi_y[j]*rho_v_n[idx[j]];
288: }
289: rho_phi[e] -= phi_dt*(Fx_x + Fy_y);
290: /* Determine the convective fluxes for (\rho u)^{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]]*u_n[idx[j]];
294: Fy_y += psi_y[j]*rho_v_n[idx[j]]*u_n[idx[j]];
295: }
296: rho_u_phi[e] -= phi_dt*(Fx_x + Fy_y);
297: /* Determine the convective fluxes for (\rho v)^{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]]*v_n[idx[j]];
301: Fy_y += psi_y[j]*rho_v_n[idx[j]]*v_n[idx[j]];
302: }
303: rho_v_phi[e] -= phi_dt*(Fx_x + Fy_y);
304: }
305: DMDARestoreElements(da, &ne, &nc, &necon);
306: VecRestoreArray(user->sol_n.u, &u_n);
307: VecRestoreArray(user->sol_n.v, &v_n);
308: VecRestoreArray(user->sol_n.rho, &rho_n);
309: VecRestoreArray(user->sol_n.rho_u, &rho_u_n);
310: VecRestoreArray(user->sol_n.rho_v, &rho_v_n);
311: VecRestoreArray(user->sol_phi.rho, &rho_phi);
312: VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
313: VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
314: return(0);
315: }
319: /*
320: The element stiffness matrix for the identity in linear elements is
322: 1 /2 1 1\
323: - |1 2 1|
324: 12 \1 1 2/
326: no matter what the shape of the triangle. */
327: PetscErrorCode TaylorGalerkinStepIIMomentum(DM da, UserContext *user)
328: {
329: MPI_Comm comm;
330: KSP ksp;
331: Mat mat;
332: Vec rhs_u, rhs_v;
333: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
334: 0.08333333333, 0.16666666667, 0.08333333333,
335: 0.08333333333, 0.08333333333, 0.16666666667};
336: PetscScalar *u_n, *v_n, *mu_n;
337: PetscScalar *u_phi, *v_phi;
338: PetscScalar *rho_u_phi, *rho_v_phi;
339: PetscInt idx[3];
340: PetscScalar values_u[3];
341: PetscScalar values_v[3];
342: PetscScalar psi_x[3], psi_y[3];
343: PetscScalar mu, tau_xx, tau_xy, tau_yy;
344: PetscReal hx, hy, area;
345: const PetscInt *necon;
346: PetscInt j, k, e, ne, nc, mx, my;
350: PetscObjectGetComm((PetscObject) da, &comm);
351: DMCreateMatrix(da, MATAIJ, &mat);
352: MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
353: DMGetGlobalVector(da, &rhs_u);
354: DMGetGlobalVector(da, &rhs_v);
355: KSPCreate(comm, &ksp);
356: KSPSetFromOptions(ksp);
358: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
359: hx = 1.0 / (PetscReal)(mx-1);
360: hy = 1.0 / (PetscReal)(my-1);
361: area = 0.5*hx*hy;
362: VecGetArray(user->sol_n.u, &u_n);
363: VecGetArray(user->sol_n.v, &v_n);
364: VecGetArray(user->mu, &mu_n);
365: VecGetArray(user->sol_phi.u, &u_phi);
366: VecGetArray(user->sol_phi.v, &v_phi);
367: VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
368: VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
369: DMDAGetElements(da, &ne, &nc, &necon);
370: for (e = 0; e < ne; e++) {
371: for (j = 0; j < 3; j++) {
372: idx[j] = necon[3*e+j];
373: values_u[j] = 0.0;
374: values_v[j] = 0.0;
375: }
376: /* Get basis function deriatives (we need the orientation of the element here) */
377: if (idx[1] > idx[0]) {
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: } else {
381: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
382: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
383: }
384: /* <\nabla\psi, F^{n+\phi}_e>: Divergence of the element-averaged convective fluxes */
385: for (j = 0; j < 3; j++) {
386: values_u[j] += psi_x[j]*rho_u_phi[e]*u_phi[e] + psi_y[j]*rho_u_phi[e]*v_phi[e];
387: values_v[j] += psi_x[j]*rho_v_phi[e]*u_phi[e] + psi_y[j]*rho_v_phi[e]*v_phi[e];
388: }
389: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
390: for (j = 0; j < 3; j++) {
391: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
392: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
393: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
394: mu = 0.0;
395: tau_xx = 0.0;
396: tau_xy = 0.0;
397: tau_yy = 0.0;
398: for (k = 0; k < 3; k++) {
399: mu += mu_n[idx[k]];
400: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
401: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
402: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
403: }
404: mu /= 3.0;
405: tau_xx *= (2.0/3.0)*mu;
406: tau_xy *= mu;
407: tau_yy *= (2.0/3.0)*mu;
408: values_u[j] -= area*(psi_x[j]*tau_xx + psi_y[j]*tau_xy);
409: values_v[j] -= area*(psi_x[j]*tau_xy + psi_y[j]*tau_yy);
410: }
411: /* Accumulate to global structures */
412: VecSetValuesLocal(rhs_u, 3, idx, values_u, ADD_VALUES);
413: VecSetValuesLocal(rhs_v, 3, idx, values_v, ADD_VALUES);
414: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
415: }
416: DMDARestoreElements(da, &ne, &nc, &necon);
417: VecRestoreArray(user->sol_n.u, &u_n);
418: VecRestoreArray(user->sol_n.v, &v_n);
419: VecRestoreArray(user->mu, &mu_n);
420: VecRestoreArray(user->sol_phi.u, &u_phi);
421: VecRestoreArray(user->sol_phi.v, &v_phi);
422: VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
423: VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
425: VecAssemblyBegin(rhs_u);
426: VecAssemblyBegin(rhs_v);
427: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
428: VecAssemblyEnd(rhs_u);
429: VecAssemblyEnd(rhs_v);
430: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
431: VecScale(rhs_u,user->dt);
432: VecScale(rhs_v,user->dt);
434: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
435: KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
436: KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
437: KSPDestroy(&ksp);
438: MatDestroy(&mat);
439: DMRestoreGlobalVector(da, &rhs_u);
440: DMRestoreGlobalVector(da, &rhs_v);
441: return(0);
442: }
446: /* Notice that this requires the previous momentum solution.
448: The element stiffness matrix for the identity in linear elements is
450: 1 /2 1 1\
451: - |1 2 1|
452: 12 \1 1 2/
454: no matter what the shape of the triangle. */
455: PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
456: {
457: MPI_Comm comm;
458: Mat mat;
459: Vec rhs_m, rhs_e;
460: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
461: 0.08333333333, 0.16666666667, 0.08333333333,
462: 0.08333333333, 0.08333333333, 0.16666666667};
463: PetscScalar *u_n, *v_n, *p_n, *t_n, *mu_n, *kappa_n;
464: PetscScalar *rho_n, *rho_u_n, *rho_v_n, *rho_e_n;
465: PetscScalar *u_phi, *v_phi;
466: PetscScalar *rho_u_np1, *rho_v_np1;
467: PetscInt idx[3];
468: PetscScalar psi_x[3], psi_y[3];
469: PetscScalar values_m[3];
470: PetscScalar values_e[3];
471: PetscScalar phi = user->phi;
472: PetscScalar mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
473: PetscReal hx, hy, area;
474: KSP ksp;
475: const PetscInt *necon;
476: PetscInt j, k, e, ne, nc, mx, my;
480: PetscObjectGetComm((PetscObject) da, &comm);
481: DMCreateMatrix(da, MATAIJ, &mat);
482: MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
483: DMGetGlobalVector(da, &rhs_m);
484: DMGetGlobalVector(da, &rhs_e);
485: KSPCreate(comm, &ksp);
486: KSPSetFromOptions(ksp);
488: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
489: hx = 1.0 / (PetscReal)(mx-1);
490: hy = 1.0 / (PetscReal)(my-1);
491: area = 0.5*hx*hy;
492: VecGetArray(user->sol_n.u, &u_n);
493: VecGetArray(user->sol_n.v, &v_n);
494: VecGetArray(user->sol_n.p, &p_n);
495: VecGetArray(user->sol_n.t, &t_n);
496: VecGetArray(user->mu, &mu_n);
497: VecGetArray(user->kappa, &kappa_n);
498: VecGetArray(user->sol_n.rho, &rho_n);
499: VecGetArray(user->sol_n.rho_u, &rho_u_n);
500: VecGetArray(user->sol_n.rho_v, &rho_v_n);
501: VecGetArray(user->sol_n.rho_e, &rho_e_n);
502: VecGetArray(user->sol_phi.u, &u_phi);
503: VecGetArray(user->sol_phi.v, &v_phi);
504: VecGetArray(user->sol_np1.rho_u, &rho_u_np1);
505: VecGetArray(user->sol_np1.rho_v, &rho_v_np1);
506: DMDAGetElements(da, &ne, &nc, &necon);
507: for (e = 0; e < ne; e++) {
508: for (j = 0; j < 3; j++) {
509: idx[j] = necon[3*e+j];
510: values_m[j] = 0.0;
511: values_e[j] = 0.0;
512: }
513: /* Get basis function deriatives (we need the orientation of the element here) */
514: if (idx[1] > idx[0]) {
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: } else {
518: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
519: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
520: }
521: /* <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
522: for (j = 0; j < 3; j++) {
523: 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;
524: values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
525: }
526: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
527: for (j = 0; j < 3; j++) {
528: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
529: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
530: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
531: /* q_x = -\kappa(T) {\partial T\over\partial x} */
532: /* q_y = -\kappa(T) {\partial T\over\partial y} */
534: /* above code commeted out - causing ininitialized variables. */
535: q_x =0; q_y =0;
537: mu = 0.0;
538: kappa = 0.0;
539: tau_xx = 0.0;
540: tau_xy = 0.0;
541: tau_yy = 0.0;
542: for (k = 0; k < 3; k++) {
543: mu += mu_n[idx[k]];
544: kappa += kappa_n[idx[k]];
545: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
546: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
547: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
548: q_x += psi_x[k]*t_n[idx[k]];
549: q_y += psi_y[k]*t_n[idx[k]];
550: }
551: mu /= 3.0;
552: kappa /= 3.0;
553: tau_xx *= (2.0/3.0)*mu;
554: tau_xy *= mu;
555: tau_yy *= (2.0/3.0)*mu;
556: 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));
557: }
558: /* Accumulate to global structures */
559: VecSetValuesLocal(rhs_m, 3, idx, values_m, ADD_VALUES);
560: VecSetValuesLocal(rhs_e, 3, idx, values_e, ADD_VALUES);
561: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
562: }
563: DMDARestoreElements(da, &ne, &nc, &necon);
564: VecRestoreArray(user->sol_n.u, &u_n);
565: VecRestoreArray(user->sol_n.v, &v_n);
566: VecRestoreArray(user->sol_n.p, &p_n);
567: VecRestoreArray(user->sol_n.t, &t_n);
568: VecRestoreArray(user->mu, &mu_n);
569: VecRestoreArray(user->kappa, &kappa_n);
570: VecRestoreArray(user->sol_n.rho, &rho_n);
571: VecRestoreArray(user->sol_n.rho_u, &rho_u_n);
572: VecRestoreArray(user->sol_n.rho_v, &rho_v_n);
573: VecRestoreArray(user->sol_n.rho_e, &rho_e_n);
574: VecRestoreArray(user->sol_phi.u, &u_phi);
575: VecRestoreArray(user->sol_phi.v, &v_phi);
576: VecRestoreArray(user->sol_np1.rho_u, &rho_u_np1);
577: VecRestoreArray(user->sol_np1.rho_v, &rho_v_np1);
579: VecAssemblyBegin(rhs_m);
580: VecAssemblyBegin(rhs_e);
581: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
582: VecAssemblyEnd(rhs_m);
583: VecAssemblyEnd(rhs_e);
584: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
585: VecScale(rhs_m, user->dt);
586: VecScale(rhs_e, user->dt);
588: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
589: KSPSolve(ksp, rhs_m, user->sol_np1.rho);
590: KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
591: KSPDestroy(&ksp);
592: MatDestroy(&mat);
593: DMRestoreGlobalVector(da, &rhs_m);
594: DMRestoreGlobalVector(da, &rhs_e);
595: return(0);
596: }
600: PetscErrorCode ComputePredictor(DM da, UserContext *user)
601: {
602: Vec uOldLocal, uLocal,uOld;
603: PetscScalar *pOld;
604: PetscScalar *p;
608: DMGetGlobalVector(da, &uOld);
609: DMGetLocalVector(da, &uOldLocal);
610: DMGetLocalVector(da, &uLocal);
611: DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
612: DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
613: VecGetArray(uOldLocal, &pOld);
614: VecGetArray(uLocal, &p);
616: /* Source terms are all zero right now */
617: CalculateElementVelocity(da, user);
618: TaylorGalerkinStepI(da, user);
619: TaylorGalerkinStepIIMomentum(da, user);
620: TaylorGalerkinStepIIMassEnergy(da, user);
621: /* Solve equation (9) for \delta(\rho\vu) and (\rho\vu)^* */
622: /* Solve equation (13) for \delta\rho and \rho^* */
623: /* Solve equation (15) for \delta(\rho e_t) and (\rho e_t)^* */
624: /* Apply artifical dissipation */
625: /* Determine the smoothed explicit pressure, \tilde P and temperature \tilde T using the equation of state */
628: VecRestoreArray(uOldLocal, &pOld);
629: VecRestoreArray(uLocal, &p);
630: #if 0
631: DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
632: DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
633: DMRestoreLocalVector(da, &uOldLocal);
634: DMRestoreLocalVector(da, &uLocal);
635: #endif
636: return(0);
637: }
641: /*
642: We integrate over each cell
644: (i, j+1)----(i+1, j+1)
645: | \ |
646: | \ |
647: | \ |
648: | \ |
649: | \ |
650: | \ |
651: | \ |
652: (i, j)----(i+1, j)
653: */
654: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
655: {
656: UserContext *user = (UserContext*)ctx;
657: PetscScalar phi = user->phi;
658: PetscScalar *array;
659: PetscInt ne,nc,i;
660: const PetscInt *e;
662: Vec blocal;
663: DM da;
666: KSPGetDM(ksp,&da);
667: /* access a local vector with room for the ghost points */
668: DMGetLocalVector(da,&blocal);
669: VecGetArray(blocal, (PetscScalar**) &array);
671: /* access the list of elements on this processor and loop over them */
672: DMDAGetElements(da,&ne,&nc,&e);
673: for (i=0; i<ne; i++) {
675: /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
676: array[e[3*i]] = phi;
677: array[e[3*i+1]] = phi;
678: array[e[3*i+2]] = phi;
679: }
680: VecRestoreArray(blocal, (PetscScalar**) &array);
681: DMDARestoreElements(da,&ne,&nc,&e);
683: /* add our partial sums over all processors into b */
684: DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
685: DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
686: DMRestoreLocalVector(da,&blocal);
687: return(0);
688: }
692: /*
693: We integrate over each cell
695: (i, j+1)----(i+1, j+1)
696: | \ |
697: | \ |
698: | \ |
699: | \ |
700: | \ |
701: | \ |
702: | \ |
703: (i, j)----(i+1, j)
705: However, the element stiffness matrix for the identity in linear elements is
707: 1 /2 1 1\
708: - |1 2 1|
709: 12 \1 1 2/
711: no matter what the shape of the triangle. The Laplacian stiffness matrix is
713: 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)\
714: - |-(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)|
715: 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 /
717: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
718: */
719: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, MatStructure *flag,void *ctx)
720: {
721: UserContext *user = (UserContext*)ctx;
722: /* not being used!
723: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
724: 0.08333333333, 0.16666666667, 0.08333333333,
725: 0.08333333333, 0.08333333333, 0.16666666667};
726: */
727: PetscScalar values[3][3];
728: PetscInt idx[3];
729: PetscScalar hx, hy, hx2, hy2, area,phi_dt2;
730: PetscInt i,mx,my,xm,ym,xs,ys;
731: PetscInt ne,nc;
732: const PetscInt *e;
734: DM da;
737: KSPGetDM(ksp,&da);
738: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
739: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
740: hx = 1.0 / (mx-1);
741: hy = 1.0 / (my-1);
742: area = 0.5*hx*hy;
743: phi_dt2 = user->phi*user->dt*user->dt;
744: hx2 = hx*hx/area*phi_dt2;
745: hy2 = hy*hy/area*phi_dt2;
747: /* initially all elements have identical geometry so all element stiffness are identical */
748: values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
749: values[1][0] = -hy2; values[1][1] = hy2; values[1][2] = 0.0;
750: values[2][0] = -hx2; values[2][1] = 0.0; values[2][2] = hx2;
752: DMDAGetElements(da,&ne,&nc,&e);
753: for (i=0; i<ne; i++) {
754: idx[0] = e[3*i];
755: idx[1] = e[3*i+1];
756: idx[2] = e[3*i+2];
757: MatSetValuesLocal(jac,3,idx,3,idx,(PetscScalar*)values,ADD_VALUES);
758: }
759: DMDARestoreElements(da,&ne,&nc,&e);
760: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
761: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
762: return(0);
763: }
767: PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u)
768: {
769: Vec uOldLocal, uLocal;
770: PetscScalar *cOld;
771: PetscScalar *c;
772: PetscInt i,ne,nc;
773: const PetscInt *e;
777: VecSet(u,0.0);
778: DMGetLocalVector(da, &uOldLocal);
779: DMGetLocalVector(da, &uLocal);
780: VecSet(uLocal,0.0);
781: DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
782: DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
783: VecGetArray(uOldLocal, &cOld);
784: VecGetArray(uLocal, &c);
786: /* access the list of elements on this processor and loop over them */
787: DMDAGetElements(da,&ne,&nc,&e);
788: for (i=0; i<ne; i++) {
790: /* this is nonsense, but copy each nodal value*/
791: c[e[3*i]] = cOld[e[3*i]];
792: c[e[3*i+1]] = cOld[e[3*i+1]];
793: c[e[3*i+2]] = cOld[e[3*i+2]];
794: }
795: DMDARestoreElements(da,&ne,&nc,&e);
796: VecRestoreArray(uOldLocal, &cOld);
797: VecRestoreArray(uLocal, &c);
798: DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
799: DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
800: DMRestoreLocalVector(da, &uOldLocal);
801: DMRestoreLocalVector(da, &uLocal);
802: return(0);
803: }