Actual source code: ex31.c
petsc-3.3-p7 2013-05-11
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>
38: #include <petscpcmg.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,MatStructure*,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, DMDA_BOUNDARY_NONE, DMDA_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, PETSC_NULL);
88: user.dt = 0.1;
89: PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, PETSC_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, PETSC_NULL, PETSC_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;
118: PetscErrorCode ierr;
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;
191: PetscErrorCode ierr;
194: DMDAGetElements(da, &ne, &nc, &necon);
195: VecGetArray(user->sol_n.u, &u_n);
196: VecGetArray(user->sol_n.v, &v_n);
197: PetscMalloc(ne*sizeof(PetscScalar),&u_phi);
198: PetscMalloc(ne*sizeof(PetscScalar),&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;
247: PetscErrorCode ierr;
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;
348: PetscErrorCode ierr;
351: PetscObjectGetComm((PetscObject) da, &comm);
352: DMCreateMatrix(da, MATAIJ, &mat);
353: MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
354: DMGetGlobalVector(da, &rhs_u);
355: DMGetGlobalVector(da, &rhs_v);
356: KSPCreate(comm, &ksp);
357: KSPSetFromOptions(ksp);
359: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
360: hx = 1.0 / (PetscReal)(mx-1);
361: hy = 1.0 / (PetscReal)(my-1);
362: area = 0.5*hx*hy;
363: VecGetArray(user->sol_n.u, &u_n);
364: VecGetArray(user->sol_n.v, &v_n);
365: VecGetArray(user->mu, &mu_n);
366: VecGetArray(user->sol_phi.u, &u_phi);
367: VecGetArray(user->sol_phi.v, &v_phi);
368: VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
369: VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
370: DMDAGetElements(da, &ne, &nc, &necon);
371: for(e = 0; e < ne; e++) {
372: for(j = 0; j < 3; j++) {
373: idx[j] = necon[3*e+j];
374: values_u[j] = 0.0;
375: values_v[j] = 0.0;
376: }
377: /* Get basis function deriatives (we need the orientation of the element here) */
378: if (idx[1] > idx[0]) {
379: psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
380: psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
381: } else {
382: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
383: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
384: }
385: /* <\nabla\psi, F^{n+\phi}_e>: Divergence of the element-averaged convective fluxes */
386: for(j = 0; j < 3; j++) {
387: values_u[j] += psi_x[j]*rho_u_phi[e]*u_phi[e] + psi_y[j]*rho_u_phi[e]*v_phi[e];
388: values_v[j] += psi_x[j]*rho_v_phi[e]*u_phi[e] + psi_y[j]*rho_v_phi[e]*v_phi[e];
389: }
390: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
391: for(j = 0; j < 3; j++) {
392: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
393: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
394: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
395: mu = 0.0;
396: tau_xx = 0.0;
397: tau_xy = 0.0;
398: tau_yy = 0.0;
399: for(k = 0; k < 3; k++) {
400: mu += mu_n[idx[k]];
401: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
402: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
403: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
404: }
405: mu /= 3.0;
406: tau_xx *= (2.0/3.0)*mu;
407: tau_xy *= mu;
408: tau_yy *= (2.0/3.0)*mu;
409: values_u[j] -= area*(psi_x[j]*tau_xx + psi_y[j]*tau_xy);
410: values_v[j] -= area*(psi_x[j]*tau_xy + psi_y[j]*tau_yy);
411: }
412: /* Accumulate to global structures */
413: VecSetValuesLocal(rhs_u, 3, idx, values_u, ADD_VALUES);
414: VecSetValuesLocal(rhs_v, 3, idx, values_v, ADD_VALUES);
415: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
416: }
417: DMDARestoreElements(da, &ne, &nc, &necon);
418: VecRestoreArray(user->sol_n.u, &u_n);
419: VecRestoreArray(user->sol_n.v, &v_n);
420: VecRestoreArray(user->mu, &mu_n);
421: VecRestoreArray(user->sol_phi.u, &u_phi);
422: VecRestoreArray(user->sol_phi.v, &v_phi);
423: VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
424: VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
426: VecAssemblyBegin(rhs_u);
427: VecAssemblyBegin(rhs_v);
428: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
429: VecAssemblyEnd(rhs_u);
430: VecAssemblyEnd(rhs_v);
431: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
432: VecScale(rhs_u,user->dt);
433: VecScale(rhs_v,user->dt);
435: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
436: KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
437: KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
438: KSPDestroy(&ksp);
439: MatDestroy(&mat);
440: DMRestoreGlobalVector(da, &rhs_u);
441: DMRestoreGlobalVector(da, &rhs_v);
442: return(0);
443: }
447: /* Notice that this requires the previous momentum solution.
449: The element stiffness matrix for the identity in linear elements is
451: 1 /2 1 1\
452: - |1 2 1|
453: 12 \1 1 2/
455: no matter what the shape of the triangle. */
456: PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
457: {
458: MPI_Comm comm;
459: Mat mat;
460: Vec rhs_m, rhs_e;
461: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
462: 0.08333333333, 0.16666666667, 0.08333333333,
463: 0.08333333333, 0.08333333333, 0.16666666667};
464: PetscScalar *u_n, *v_n, *p_n, *t_n, *mu_n, *kappa_n;
465: PetscScalar *rho_n, *rho_u_n, *rho_v_n, *rho_e_n;
466: PetscScalar *u_phi, *v_phi;
467: PetscScalar *rho_u_np1, *rho_v_np1;
468: PetscInt idx[3];
469: PetscScalar psi_x[3], psi_y[3];
470: PetscScalar values_m[3];
471: PetscScalar values_e[3];
472: PetscScalar phi = user->phi;
473: PetscScalar mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
474: PetscReal hx, hy, area;
475: KSP ksp;
476: const PetscInt *necon;
477: PetscInt j, k, e, ne, nc, mx, my;
478: PetscErrorCode ierr;
481: PetscObjectGetComm((PetscObject) da, &comm);
482: DMCreateMatrix(da, MATAIJ, &mat);
483: MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
484: DMGetGlobalVector(da, &rhs_m);
485: DMGetGlobalVector(da, &rhs_e);
486: KSPCreate(comm, &ksp);
487: KSPSetFromOptions(ksp);
489: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
490: hx = 1.0 / (PetscReal)(mx-1);
491: hy = 1.0 / (PetscReal)(my-1);
492: area = 0.5*hx*hy;
493: VecGetArray(user->sol_n.u, &u_n);
494: VecGetArray(user->sol_n.v, &v_n);
495: VecGetArray(user->sol_n.p, &p_n);
496: VecGetArray(user->sol_n.t, &t_n);
497: VecGetArray(user->mu, &mu_n);
498: VecGetArray(user->kappa, &kappa_n);
499: VecGetArray(user->sol_n.rho, &rho_n);
500: VecGetArray(user->sol_n.rho_u, &rho_u_n);
501: VecGetArray(user->sol_n.rho_v, &rho_v_n);
502: VecGetArray(user->sol_n.rho_e, &rho_e_n);
503: VecGetArray(user->sol_phi.u, &u_phi);
504: VecGetArray(user->sol_phi.v, &v_phi);
505: VecGetArray(user->sol_np1.rho_u, &rho_u_np1);
506: VecGetArray(user->sol_np1.rho_v, &rho_v_np1);
507: DMDAGetElements(da, &ne, &nc, &necon);
508: for(e = 0; e < ne; e++) {
509: for(j = 0; j < 3; j++) {
510: idx[j] = necon[3*e+j];
511: values_m[j] = 0.0;
512: values_e[j] = 0.0;
513: }
514: /* Get basis function deriatives (we need the orientation of the element here) */
515: if (idx[1] > idx[0]) {
516: psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
517: psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
518: } else {
519: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
520: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
521: }
522: /* <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
523: for(j = 0; j < 3; j++) {
524: 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;
525: values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
526: }
527: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
528: for(j = 0; j < 3; j++) {
529: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
530: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
531: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
532: /* q_x = -\kappa(T) {\partial T\over\partial x} */
533: /* q_y = -\kappa(T) {\partial T\over\partial y} */
535: /* above code commeted out - causing ininitialized variables. */
536: q_x =0; q_y =0;
538: mu = 0.0;
539: kappa = 0.0;
540: tau_xx = 0.0;
541: tau_xy = 0.0;
542: tau_yy = 0.0;
543: for(k = 0; k < 3; k++) {
544: mu += mu_n[idx[k]];
545: kappa += kappa_n[idx[k]];
546: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
547: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
548: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
549: q_x += psi_x[k]*t_n[idx[k]];
550: q_y += psi_y[k]*t_n[idx[k]];
551: }
552: mu /= 3.0;
553: kappa /= 3.0;
554: tau_xx *= (2.0/3.0)*mu;
555: tau_xy *= mu;
556: tau_yy *= (2.0/3.0)*mu;
557: 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));
558: }
559: /* Accumulate to global structures */
560: VecSetValuesLocal(rhs_m, 3, idx, values_m, ADD_VALUES);
561: VecSetValuesLocal(rhs_e, 3, idx, values_e, ADD_VALUES);
562: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
563: }
564: DMDARestoreElements(da, &ne, &nc, &necon);
565: VecRestoreArray(user->sol_n.u, &u_n);
566: VecRestoreArray(user->sol_n.v, &v_n);
567: VecRestoreArray(user->sol_n.p, &p_n);
568: VecRestoreArray(user->sol_n.t, &t_n);
569: VecRestoreArray(user->mu, &mu_n);
570: VecRestoreArray(user->kappa, &kappa_n);
571: VecRestoreArray(user->sol_n.rho, &rho_n);
572: VecRestoreArray(user->sol_n.rho_u, &rho_u_n);
573: VecRestoreArray(user->sol_n.rho_v, &rho_v_n);
574: VecRestoreArray(user->sol_n.rho_e, &rho_e_n);
575: VecRestoreArray(user->sol_phi.u, &u_phi);
576: VecRestoreArray(user->sol_phi.v, &v_phi);
577: VecRestoreArray(user->sol_np1.rho_u, &rho_u_np1);
578: VecRestoreArray(user->sol_np1.rho_v, &rho_v_np1);
580: VecAssemblyBegin(rhs_m);
581: VecAssemblyBegin(rhs_e);
582: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
583: VecAssemblyEnd(rhs_m);
584: VecAssemblyEnd(rhs_e);
585: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
586: VecScale(rhs_m, user->dt);
587: VecScale(rhs_e, user->dt);
589: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
590: KSPSolve(ksp, rhs_m, user->sol_np1.rho);
591: KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
592: KSPDestroy(&ksp);
593: MatDestroy(&mat);
594: DMRestoreGlobalVector(da, &rhs_m);
595: DMRestoreGlobalVector(da, &rhs_e);
596: return(0);
597: }
601: PetscErrorCode ComputePredictor(DM da, UserContext *user)
602: {
603: Vec uOldLocal, uLocal,uOld;
604: PetscScalar *pOld;
605: PetscScalar *p;
609: DMGetGlobalVector(da, &uOld);
610: DMGetLocalVector(da, &uOldLocal);
611: DMGetLocalVector(da, &uLocal);
612: DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
613: DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
614: VecGetArray(uOldLocal, &pOld);
615: VecGetArray(uLocal, &p);
617: /* Source terms are all zero right now */
618: CalculateElementVelocity(da, user);
619: TaylorGalerkinStepI(da, user);
620: TaylorGalerkinStepIIMomentum(da, user);
621: TaylorGalerkinStepIIMassEnergy(da, user);
622: /* Solve equation (9) for \delta(\rho\vu) and (\rho\vu)^* */
623: /* Solve equation (13) for \delta\rho and \rho^* */
624: /* Solve equation (15) for \delta(\rho e_t) and (\rho e_t)^* */
625: /* Apply artifical dissipation */
626: /* Determine the smoothed explicit pressure, \tilde P and temperature \tilde T using the equation of state */
629: VecRestoreArray(uOldLocal, &pOld);
630: VecRestoreArray(uLocal, &p);
631: #if 0
632: DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
633: DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
634: DMRestoreLocalVector(da, &uOldLocal);
635: DMRestoreLocalVector(da, &uLocal);
636: #endif
637: return(0);
638: }
642: /*
643: We integrate over each cell
645: (i, j+1)----(i+1, j+1)
646: | \ |
647: | \ |
648: | \ |
649: | \ |
650: | \ |
651: | \ |
652: | \ |
653: (i, j)----(i+1, j)
654: */
655: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
656: {
657: UserContext *user = (UserContext*)ctx;
658: PetscScalar phi = user->phi;
659: PetscScalar *array;
660: PetscInt ne,nc,i;
661: const PetscInt *e;
663: Vec blocal;
664: DM da;
667: KSPGetDM(ksp,&da);
668: /* access a local vector with room for the ghost points */
669: DMGetLocalVector(da,&blocal);
670: VecGetArray(blocal, (PetscScalar **) &array);
672: /* access the list of elements on this processor and loop over them */
673: DMDAGetElements(da,&ne,&nc,&e);
674: for (i=0; i<ne; i++) {
676: /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
677: array[e[3*i]] = phi;
678: array[e[3*i+1]] = phi;
679: array[e[3*i+2]] = phi;
680: }
681: VecRestoreArray(blocal, (PetscScalar **) &array);
682: DMDARestoreElements(da,&ne,&nc,&e);
684: /* add our partial sums over all processors into b */
685: DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
686: DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
687: DMRestoreLocalVector(da,&blocal);
688: return(0);
689: }
693: /*
694: We integrate over each cell
696: (i, j+1)----(i+1, j+1)
697: | \ |
698: | \ |
699: | \ |
700: | \ |
701: | \ |
702: | \ |
703: | \ |
704: (i, j)----(i+1, j)
706: However, the element stiffness matrix for the identity in linear elements is
708: 1 /2 1 1\
709: - |1 2 1|
710: 12 \1 1 2/
712: no matter what the shape of the triangle. The Laplacian stiffness matrix is
714: 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)\
715: - |-(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)|
716: 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 /
718: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
719: */
720: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, MatStructure *flag,void *ctx)
721: {
722: UserContext *user = (UserContext*)ctx;
723: /* not being used!
724: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
725: 0.08333333333, 0.16666666667, 0.08333333333,
726: 0.08333333333, 0.08333333333, 0.16666666667};
727: */
728: PetscScalar values[3][3];
729: PetscInt idx[3];
730: PetscScalar hx, hy, hx2, hy2, area,phi_dt2;
731: PetscInt i,mx,my,xm,ym,xs,ys;
732: PetscInt ne,nc;
733: const PetscInt *e;
735: DM da;
738: KSPGetDM(ksp,&da);
739: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
740: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
741: hx = 1.0 / (mx-1);
742: hy = 1.0 / (my-1);
743: area = 0.5*hx*hy;
744: phi_dt2 = user->phi*user->dt*user->dt;
745: hx2 = hx*hx/area*phi_dt2;
746: hy2 = hy*hy/area*phi_dt2;
748: /* initially all elements have identical geometry so all element stiffness are identical */
749: values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
750: values[1][0] = -hy2; values[1][1] = hy2; values[1][2] = 0.0;
751: values[2][0] = -hx2; values[2][1] = 0.0; values[2][2] = hx2;
753: DMDAGetElements(da,&ne,&nc,&e);
754: for (i=0; i<ne; i++) {
755: idx[0] = e[3*i];
756: idx[1] = e[3*i+1];
757: idx[2] = e[3*i+2];
758: MatSetValuesLocal(jac,3,idx,3,idx,(PetscScalar*)values,ADD_VALUES);
759: }
760: DMDARestoreElements(da,&ne,&nc,&e);
761: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
762: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
763: return(0);
764: }
768: PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u)
769: {
770: Vec uOldLocal, uLocal;
771: PetscScalar *cOld;
772: PetscScalar *c;
773: PetscInt i,ne,nc;
774: const PetscInt *e;
776:
778: VecSet(u,0.0);
779: DMGetLocalVector(da, &uOldLocal);
780: DMGetLocalVector(da, &uLocal);
781: VecSet(uLocal,0.0);
782: DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
783: DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
784: VecGetArray(uOldLocal, &cOld);
785: VecGetArray(uLocal, &c);
787: /* access the list of elements on this processor and loop over them */
788: DMDAGetElements(da,&ne,&nc,&e);
789: for (i=0; i<ne; i++) {
791: /* this is nonsense, but copy each nodal value*/
792: c[e[3*i]] = cOld[e[3*i]];
793: c[e[3*i+1]] = cOld[e[3*i+1]];
794: c[e[3*i+2]] = cOld[e[3*i+2]];
795: }
796: DMDARestoreElements(da,&ne,&nc,&e);
797: VecRestoreArray(uOldLocal, &cOld);
798: VecRestoreArray(uLocal, &c);
799: DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
800: DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
801: DMRestoreLocalVector(da, &uOldLocal);
802: DMRestoreLocalVector(da, &uLocal);
803: return(0);
804: }