Actual source code: ex31.c
petsc-3.8.4 2018-03-24
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: }