Actual source code: ex31.c
petsc-3.7.5 2017-01-01
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: }