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