Actual source code: ex31.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Concepts: KSP^semi-implicit
5: Processors: n
6: T*/
8: /*
9: This is intended to be a prototypical example of the semi-implicit algorithm for
10: a PDE. We have three phases:
12: 1) An explicit predictor step
14: u^{k+1/3} = P(u^k)
16: 2) An implicit corrector step
18: \Delta u^{k+2/3} = F(u^{k+1/3})
20: 3) An explicit update step
22: u^{k+1} = C(u^{k+2/3})
24: We will solve on the unit square with Dirichlet boundary conditions
26: u = f(x,y) for x = 0, x = 1, y = 0, y = 1
28: Although we are using a DMDA, and thus have a structured mesh, we will discretize
29: the problem with finite elements, splitting each cell of the DMDA into two
30: triangles.
32: This uses multigrid to solve the linear system
33: */
35: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
37: #include <petscdmda.h>
38: #include <petscksp.h>
39: #include <petscpcmg.h>
40: #include <petscdmmg.h>
50: typedef struct {
51: Vec rho; /* The mass solution \rho */
52: Vec rho_u; /* The x-momentum solution \rho u */
53: Vec rho_v; /* The y-momentum solution \rho v */
54: Vec rho_e; /* The energy solution \rho e_t */
55: Vec p; /* The pressure solution P */
56: Vec t; /* The temperature solution T */
57: Vec u; /* The x-velocity solution u */
58: Vec v; /* The y-velocity solution v */
59: } SolutionContext;
61: typedef struct {
62: SolutionContext sol_n; /* The solution at time t^n */
63: SolutionContext sol_phi; /* The element-averaged solution at time t^{n+\phi} */
64: SolutionContext sol_np1; /* The solution at time t^{n+1} */
65: Vec mu; /* The dynamic viscosity \mu(T) at time n */
66: Vec kappa; /* The thermal conductivity \kappa(T) at time n */
67: PetscScalar phi; /* The time weighting parameter */
68: PetscScalar dt; /* The timestep \Delta t */
69: } UserContext;
73: int main(int argc,char **argv)
74: {
75: DMMG *dmmg;
76: DM da;
77: UserContext user;
79: PetscInt l;
81: PetscInitialize(&argc,&argv,(char *)0,help);
83: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
84: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
85: DMMGSetDM(dmmg,(DM)da);
86: DMDestroy(&da);
87: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
88: DMMGSetUser(dmmg,l,&user);
89: }
91: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PCICE", "DMMG");
92: user.phi = 0.5;
93: PetscOptionsScalar("-phi", "The time weighting parameter", "ex31.c", user.phi, &user.phi, PETSC_NULL);
94: user.dt = 0.1;
95: PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, PETSC_NULL);
96: PetscOptionsEnd();
98: CreateStructures(DMMGGetFine(dmmg));
99: ComputeInitialGuess(DMMGGetFine(dmmg));
100: ComputePredictor(DMMGGetFine(dmmg));
102: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
103: DMMGSetInitialGuess(dmmg, DMMGInitialGuessCurrent);
104: DMMGSolve(dmmg);
106: ComputeCorrector(DMMGGetFine(dmmg), DMMGGetx(dmmg), DMMGGetr(dmmg));
108: DestroyStructures(DMMGGetFine(dmmg));
109: DMMGDestroy(dmmg);
110: PetscFinalize();
111: return 0;
112: }
116: PetscErrorCode CreateStructures(DMMG dmmg)
117: {
118: DM da = dmmg->dm;
119: UserContext *user = (UserContext *) dmmg->user;
120: const PetscInt *necon;
121: PetscInt ne,nc;
122: PetscErrorCode ierr;
125: DMDAGetElements(da,&ne,&nc,&necon);
126: DMDARestoreElements(da,&ne,&nc,&necon);
127: DMCreateGlobalVector(da, &user->sol_n.rho);
128: DMCreateGlobalVector(da, &user->sol_n.rho_u);
129: DMCreateGlobalVector(da, &user->sol_n.rho_v);
130: DMCreateGlobalVector(da, &user->sol_n.rho_e);
131: DMCreateGlobalVector(da, &user->sol_n.p);
132: DMCreateGlobalVector(da, &user->sol_n.u);
133: DMCreateGlobalVector(da, &user->sol_n.v);
134: VecCreate(PETSC_COMM_WORLD, &user->sol_phi.rho_u);
135: VecSetSizes(user->sol_phi.rho_u, ne, PETSC_DECIDE);
136: VecSetType(user->sol_phi.rho_u,VECMPI);
137: VecDuplicate(user->sol_phi.rho_u, &user->sol_phi.rho_v);
138: VecDuplicate(user->sol_phi.rho_u, &user->sol_phi.u);
139: VecDuplicate(user->sol_phi.rho_u, &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: }
154: PetscErrorCode DestroyStructures(DMMG dmmg)
155: {
156: UserContext *user = (UserContext *) dmmg->user;
160: VecDestroy(&user->sol_n.rho);
161: VecDestroy(&user->sol_n.rho_u);
162: VecDestroy(&user->sol_n.rho_v);
163: VecDestroy(&user->sol_n.rho_e);
164: VecDestroy(&user->sol_n.p);
165: VecDestroy(&user->sol_n.u);
166: VecDestroy(&user->sol_n.v);
167: VecDestroy(&user->sol_phi.rho_u);
168: VecDestroy(&user->sol_phi.rho_v);
169: VecDestroy(&user->sol_phi.u);
170: VecDestroy(&user->sol_phi.v);
171: VecDestroy(&user->sol_np1.rho);
172: VecDestroy(&user->sol_np1.rho_u);
173: VecDestroy(&user->sol_np1.rho_v);
174: VecDestroy(&user->sol_np1.rho_e);
175: VecDestroy(&user->sol_np1.p);
176: VecDestroy(&user->sol_np1.u);
177: VecDestroy(&user->sol_np1.v);
178: VecDestroy(&user->mu);
179: VecDestroy(&user->kappa);
180: return(0);
181: }
185: PetscErrorCode ComputeInitialGuess(DMMG dmmg)
186: {
188: return(0);
189: }
193: /* Average the velocity (u,v) at time t^n over each element for time n+\phi */
194: PetscErrorCode CalculateElementVelocity(DM da, UserContext *user)
195: {
196: PetscScalar *u_n, *v_n;
197: PetscScalar *u_phi, *v_phi;
198: const PetscInt *necon;
199: PetscInt j, e, ne, nc;
200: PetscErrorCode ierr;
203: DMDAGetElements(da, &ne, &nc, &necon);
204: VecGetArray(user->sol_n.u, &u_n);
205: VecGetArray(user->sol_n.v, &v_n);
206: PetscMalloc(ne*sizeof(PetscScalar),&u_phi);
207: PetscMalloc(ne*sizeof(PetscScalar),&v_phi);
208: for(e = 0; e < ne; e++) {
209: u_phi[e] = 0.0;
210: v_phi[e] = 0.0;
211: for(j = 0; j < 3; j++) {
212: u_phi[e] += u_n[necon[3*e+j]];
213: v_phi[e] += v_n[necon[3*e+j]];
214: }
215: u_phi[e] /= 3.0;
216: v_phi[e] /= 3.0;
217: }
218: PetscFree(u_phi);
219: PetscFree(v_phi);
220: DMDARestoreElements(da, &ne,&nc, &necon);
221: VecRestoreArray(user->sol_n.u, &u_n);
222: VecRestoreArray(user->sol_n.v, &v_n);
223: return(0);
224: }
228: /* This is equation 32,
230: 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
232: which is really simple for linear elements
234: 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
236: where
238: U^{n+\phi}_E = {\rho \rho u \rho v}^{n+\phi}_E
240: and the x and y components of the convective fluxes F are
242: f^n = {\rho u \rho u^2 \rho uv}^n g^n = {\rho v \rho uv \rho v^2}^n
243: */
244: PetscErrorCode TaylorGalerkinStepI(DM da, UserContext *user)
245: {
246: PetscScalar phi_dt = user->phi*user->dt;
247: PetscScalar *u_n, *v_n;
248: PetscScalar *rho_n, *rho_u_n, *rho_v_n;
249: PetscScalar *rho_phi, *rho_u_phi, *rho_v_phi;
250: PetscScalar Fx_x, Fy_y;
251: PetscScalar psi_x[3], psi_y[3];
252: PetscInt idx[3];
253: PetscReal hx, hy;
254: const PetscInt *necon;
255: PetscInt j, e, ne, nc, mx, my;
256: PetscErrorCode ierr;
259: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
260: hx = 1.0 / (PetscReal)(mx-1);
261: hy = 1.0 / (PetscReal)(my-1);
262: VecSet(user->sol_phi.rho,0.0);
263: VecSet(user->sol_phi.rho_u,0.0);
264: VecSet(user->sol_phi.rho_v,0.0);
265: VecGetArray(user->sol_n.u, &u_n);
266: VecGetArray(user->sol_n.v, &v_n);
267: VecGetArray(user->sol_n.rho, &rho_n);
268: VecGetArray(user->sol_n.rho_u, &rho_u_n);
269: VecGetArray(user->sol_n.rho_v, &rho_v_n);
270: VecGetArray(user->sol_phi.rho, &rho_phi);
271: VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
272: VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
273: DMDAGetElements(da, &ne, &nc, &necon);
274: for(e = 0; e < ne; e++) {
275: /* Average the existing fields over the element */
276: for(j = 0; j < 3; j++) {
277: idx[j] = necon[3*e+j];
278: rho_phi[e] += rho_n[idx[j]];
279: rho_u_phi[e] += rho_u_n[idx[j]];
280: rho_v_phi[e] += rho_v_n[idx[j]];
281: }
282: rho_phi[e] /= 3.0;
283: rho_u_phi[e] /= 3.0;
284: rho_v_phi[e] /= 3.0;
285: /* Get basis function deriatives (we need the orientation of the element here) */
286: if (idx[1] > idx[0]) {
287: psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
288: psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
289: } else {
290: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
291: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
292: }
293: /* Determine the convective fluxes for \rho^{n+\phi} */
294: Fx_x = 0.0; Fy_y = 0.0;
295: for(j = 0; j < 3; j++) {
296: Fx_x += psi_x[j]*rho_u_n[idx[j]];
297: Fy_y += psi_y[j]*rho_v_n[idx[j]];
298: }
299: rho_phi[e] -= phi_dt*(Fx_x + Fy_y);
300: /* Determine the convective fluxes for (\rho u)^{n+\phi} */
301: Fx_x = 0.0; Fy_y = 0.0;
302: for(j = 0; j < 3; j++) {
303: Fx_x += psi_x[j]*rho_u_n[idx[j]]*u_n[idx[j]];
304: Fy_y += psi_y[j]*rho_v_n[idx[j]]*u_n[idx[j]];
305: }
306: rho_u_phi[e] -= phi_dt*(Fx_x + Fy_y);
307: /* Determine the convective fluxes for (\rho v)^{n+\phi} */
308: Fx_x = 0.0; Fy_y = 0.0;
309: for(j = 0; j < 3; j++) {
310: Fx_x += psi_x[j]*rho_u_n[idx[j]]*v_n[idx[j]];
311: Fy_y += psi_y[j]*rho_v_n[idx[j]]*v_n[idx[j]];
312: }
313: rho_v_phi[e] -= phi_dt*(Fx_x + Fy_y);
314: }
315: DMDARestoreElements(da, &ne, &nc, &necon);
316: VecRestoreArray(user->sol_n.u, &u_n);
317: VecRestoreArray(user->sol_n.v, &v_n);
318: VecRestoreArray(user->sol_n.rho, &rho_n);
319: VecRestoreArray(user->sol_n.rho_u, &rho_u_n);
320: VecRestoreArray(user->sol_n.rho_v, &rho_v_n);
321: VecRestoreArray(user->sol_phi.rho, &rho_phi);
322: VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
323: VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
324: return(0);
325: }
329: /*
330: The element stiffness matrix for the identity in linear elements is
332: 1 /2 1 1\
333: - |1 2 1|
334: 12 \1 1 2/
336: no matter what the shape of the triangle. */
337: PetscErrorCode TaylorGalerkinStepIIMomentum(DM da, UserContext *user)
338: {
339: MPI_Comm comm;
340: KSP ksp;
341: Mat mat;
342: Vec rhs_u, rhs_v;
343: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
344: 0.08333333333, 0.16666666667, 0.08333333333,
345: 0.08333333333, 0.08333333333, 0.16666666667};
346: PetscScalar *u_n, *v_n, *mu_n;
347: PetscScalar *u_phi, *v_phi;
348: PetscScalar *rho_u_phi, *rho_v_phi;
349: PetscInt idx[3];
350: PetscScalar values_u[3];
351: PetscScalar values_v[3];
352: PetscScalar psi_x[3], psi_y[3];
353: PetscScalar mu, tau_xx, tau_xy, tau_yy;
354: PetscReal hx, hy, area;
355: const PetscInt *necon;
356: PetscInt j, k, e, ne, nc, mx, my;
357: PetscErrorCode ierr;
360: PetscObjectGetComm((PetscObject) da, &comm);
361: DMGetMatrix(da, MATAIJ, &mat);
362: DMGetGlobalVector(da, &rhs_u);
363: DMGetGlobalVector(da, &rhs_v);
364: KSPCreate(comm, &ksp);
365: KSPSetFromOptions(ksp);
367: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
368: hx = 1.0 / (PetscReal)(mx-1);
369: hy = 1.0 / (PetscReal)(my-1);
370: area = 0.5*hx*hy;
371: VecGetArray(user->sol_n.u, &u_n);
372: VecGetArray(user->sol_n.v, &v_n);
373: VecGetArray(user->mu, &mu_n);
374: VecGetArray(user->sol_phi.u, &u_phi);
375: VecGetArray(user->sol_phi.v, &v_phi);
376: VecGetArray(user->sol_phi.rho_u, &rho_u_phi);
377: VecGetArray(user->sol_phi.rho_v, &rho_v_phi);
378: DMDAGetElements(da, &ne, &nc, &necon);
379: for(e = 0; e < ne; e++) {
380: for(j = 0; j < 3; j++) {
381: idx[j] = necon[3*e+j];
382: values_u[j] = 0.0;
383: values_v[j] = 0.0;
384: }
385: /* Get basis function deriatives (we need the orientation of the element here) */
386: if (idx[1] > idx[0]) {
387: psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
388: psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
389: } else {
390: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
391: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
392: }
393: /* <\nabla\psi, F^{n+\phi}_e>: Divergence of the element-averaged convective fluxes */
394: for(j = 0; j < 3; j++) {
395: values_u[j] += psi_x[j]*rho_u_phi[e]*u_phi[e] + psi_y[j]*rho_u_phi[e]*v_phi[e];
396: values_v[j] += psi_x[j]*rho_v_phi[e]*u_phi[e] + psi_y[j]*rho_v_phi[e]*v_phi[e];
397: }
398: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
399: for(j = 0; j < 3; j++) {
400: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
401: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
402: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
403: mu = 0.0;
404: tau_xx = 0.0;
405: tau_xy = 0.0;
406: tau_yy = 0.0;
407: for(k = 0; k < 3; k++) {
408: mu += mu_n[idx[k]];
409: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
410: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
411: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
412: }
413: mu /= 3.0;
414: tau_xx *= (2.0/3.0)*mu;
415: tau_xy *= mu;
416: tau_yy *= (2.0/3.0)*mu;
417: values_u[j] -= area*(psi_x[j]*tau_xx + psi_y[j]*tau_xy);
418: values_v[j] -= area*(psi_x[j]*tau_xy + psi_y[j]*tau_yy);
419: }
420: /* Accumulate to global structures */
421: VecSetValuesLocal(rhs_u, 3, idx, values_u, ADD_VALUES);
422: VecSetValuesLocal(rhs_v, 3, idx, values_v, ADD_VALUES);
423: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
424: }
425: DMDARestoreElements(da, &ne, &nc, &necon);
426: VecRestoreArray(user->sol_n.u, &u_n);
427: VecRestoreArray(user->sol_n.v, &v_n);
428: VecRestoreArray(user->mu, &mu_n);
429: VecRestoreArray(user->sol_phi.u, &u_phi);
430: VecRestoreArray(user->sol_phi.v, &v_phi);
431: VecRestoreArray(user->sol_phi.rho_u, &rho_u_phi);
432: VecRestoreArray(user->sol_phi.rho_v, &rho_v_phi);
434: VecAssemblyBegin(rhs_u);
435: VecAssemblyBegin(rhs_v);
436: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
437: VecAssemblyEnd(rhs_u);
438: VecAssemblyEnd(rhs_v);
439: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
440: VecScale(rhs_u,user->dt);
441: VecScale(rhs_v,user->dt);
443: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
444: KSPSolve(ksp, rhs_u, user->sol_np1.rho_u);
445: KSPSolve(ksp, rhs_v, user->sol_np1.rho_v);
446: KSPDestroy(&ksp);
447: MatDestroy(&mat);
448: DMRestoreGlobalVector(da, &rhs_u);
449: DMRestoreGlobalVector(da, &rhs_v);
450: return(0);
451: }
455: /* Notice that this requires the previous momentum solution.
457: The element stiffness matrix for the identity in linear elements is
459: 1 /2 1 1\
460: - |1 2 1|
461: 12 \1 1 2/
463: no matter what the shape of the triangle. */
464: PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
465: {
466: MPI_Comm comm;
467: Mat mat;
468: Vec rhs_m, rhs_e;
469: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
470: 0.08333333333, 0.16666666667, 0.08333333333,
471: 0.08333333333, 0.08333333333, 0.16666666667};
472: PetscScalar *u_n, *v_n, *p_n, *t_n, *mu_n, *kappa_n;
473: PetscScalar *rho_n, *rho_u_n, *rho_v_n, *rho_e_n;
474: PetscScalar *u_phi, *v_phi;
475: PetscScalar *rho_u_np1, *rho_v_np1;
476: PetscInt idx[3];
477: PetscScalar psi_x[3], psi_y[3];
478: PetscScalar values_m[3];
479: PetscScalar values_e[3];
480: PetscScalar phi = user->phi;
481: PetscScalar mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
482: PetscReal hx, hy, area;
483: KSP ksp;
484: const PetscInt *necon;
485: PetscInt j, k, e, ne, nc, mx, my;
486: PetscErrorCode ierr;
489: PetscObjectGetComm((PetscObject) da, &comm);
490: DMGetMatrix(da, MATAIJ, &mat);
491: DMGetGlobalVector(da, &rhs_m);
492: DMGetGlobalVector(da, &rhs_e);
493: KSPCreate(comm, &ksp);
494: KSPSetFromOptions(ksp);
496: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
497: hx = 1.0 / (PetscReal)(mx-1);
498: hy = 1.0 / (PetscReal)(my-1);
499: area = 0.5*hx*hy;
500: VecGetArray(user->sol_n.u, &u_n);
501: VecGetArray(user->sol_n.v, &v_n);
502: VecGetArray(user->sol_n.p, &p_n);
503: VecGetArray(user->sol_n.t, &t_n);
504: VecGetArray(user->mu, &mu_n);
505: VecGetArray(user->kappa, &kappa_n);
506: VecGetArray(user->sol_n.rho, &rho_n);
507: VecGetArray(user->sol_n.rho_u, &rho_u_n);
508: VecGetArray(user->sol_n.rho_v, &rho_v_n);
509: VecGetArray(user->sol_n.rho_e, &rho_e_n);
510: VecGetArray(user->sol_phi.u, &u_phi);
511: VecGetArray(user->sol_phi.v, &v_phi);
512: VecGetArray(user->sol_np1.rho_u, &rho_u_np1);
513: VecGetArray(user->sol_np1.rho_v, &rho_v_np1);
514: DMDAGetElements(da, &ne, &nc, &necon);
515: for(e = 0; e < ne; e++) {
516: for(j = 0; j < 3; j++) {
517: idx[j] = necon[3*e+j];
518: values_m[j] = 0.0;
519: values_e[j] = 0.0;
520: }
521: /* Get basis function deriatives (we need the orientation of the element here) */
522: if (idx[1] > idx[0]) {
523: psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
524: psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
525: } else {
526: psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
527: psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
528: }
529: /* <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
530: for(j = 0; j < 3; j++) {
531: 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;
532: values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
533: }
534: /* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
535: for(j = 0; j < 3; j++) {
536: /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
537: /* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
538: /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
539: /* q_x = -\kappa(T) {\partial T\over\partial x} */
540: /* q_y = -\kappa(T) {\partial T\over\partial y} */
542: /* above code commeted out - causing ininitialized variables. */
543: q_x =0; q_y =0;
545: mu = 0.0;
546: kappa = 0.0;
547: tau_xx = 0.0;
548: tau_xy = 0.0;
549: tau_yy = 0.0;
550: for(k = 0; k < 3; k++) {
551: mu += mu_n[idx[k]];
552: kappa += kappa_n[idx[k]];
553: tau_xx += 2.0*psi_x[k]*u_n[idx[k]] - psi_y[k]*v_n[idx[k]];
554: tau_xy += psi_y[k]*u_n[idx[k]] + psi_x[k]*v_n[idx[k]];
555: tau_yy += 2.0*psi_y[k]*v_n[idx[k]] - psi_x[k]*u_n[idx[k]];
556: q_x += psi_x[k]*t_n[idx[k]];
557: q_y += psi_y[k]*t_n[idx[k]];
558: }
559: mu /= 3.0;
560: kappa /= 3.0;
561: tau_xx *= (2.0/3.0)*mu;
562: tau_xy *= mu;
563: tau_yy *= (2.0/3.0)*mu;
564: 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));
565: }
566: /* Accumulate to global structures */
567: VecSetValuesLocal(rhs_m, 3, idx, values_m, ADD_VALUES);
568: VecSetValuesLocal(rhs_e, 3, idx, values_e, ADD_VALUES);
569: MatSetValuesLocal(mat, 3, idx, 3, idx, identity, ADD_VALUES);
570: }
571: DMDARestoreElements(da, &ne, &nc, &necon);
572: VecRestoreArray(user->sol_n.u, &u_n);
573: VecRestoreArray(user->sol_n.v, &v_n);
574: VecRestoreArray(user->sol_n.p, &p_n);
575: VecRestoreArray(user->sol_n.t, &t_n);
576: VecRestoreArray(user->mu, &mu_n);
577: VecRestoreArray(user->kappa, &kappa_n);
578: VecRestoreArray(user->sol_n.rho, &rho_n);
579: VecRestoreArray(user->sol_n.rho_u, &rho_u_n);
580: VecRestoreArray(user->sol_n.rho_v, &rho_v_n);
581: VecRestoreArray(user->sol_n.rho_e, &rho_e_n);
582: VecRestoreArray(user->sol_phi.u, &u_phi);
583: VecRestoreArray(user->sol_phi.v, &v_phi);
584: VecRestoreArray(user->sol_np1.rho_u, &rho_u_np1);
585: VecRestoreArray(user->sol_np1.rho_v, &rho_v_np1);
587: VecAssemblyBegin(rhs_m);
588: VecAssemblyBegin(rhs_e);
589: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
590: VecAssemblyEnd(rhs_m);
591: VecAssemblyEnd(rhs_e);
592: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
593: VecScale(rhs_m, user->dt);
594: VecScale(rhs_e, user->dt);
596: KSPSetOperators(ksp, mat, mat, DIFFERENT_NONZERO_PATTERN);
597: KSPSolve(ksp, rhs_m, user->sol_np1.rho);
598: KSPSolve(ksp, rhs_e, user->sol_np1.rho_e);
599: KSPDestroy(&ksp);
600: MatDestroy(&mat);
601: DMRestoreGlobalVector(da, &rhs_m);
602: DMRestoreGlobalVector(da, &rhs_e);
603: return(0);
604: }
608: PetscErrorCode ComputePredictor(DMMG dmmg)
609: {
610: DM da = dmmg->dm;
611: UserContext *user = (UserContext *) dmmg->user;
612: Vec uOldLocal, uLocal,uOld;
613: PetscScalar *pOld;
614: PetscScalar *p;
616:
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: DMRestoreLocalVector(da, &uOldLocal);
644: DMRestoreLocalVector(da, &uLocal);
645: #endif
646: return(0);
647: }
651: /*
652: We integrate over each cell
654: (i, j+1)----(i+1, j+1)
655: | \ |
656: | \ |
657: | \ |
658: | \ |
659: | \ |
660: | \ |
661: | \ |
662: (i, j)----(i+1, j)
663: */
664: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
665: {
666: DM da = dmmg->dm;
667: UserContext *user = (UserContext *) dmmg->user;
668: PetscScalar phi = user->phi;
669: PetscScalar *array;
670: PetscInt ne,nc,i;
671: const PetscInt *e;
673: Vec blocal;
676: /* access a local vector with room for the ghost points */
677: DMGetLocalVector(da,&blocal);
678: VecGetArray(blocal, (PetscScalar **) &array);
680: /* access the list of elements on this processor and loop over them */
681: DMDAGetElements(da,&ne,&nc,&e);
682: for (i=0; i<ne; i++) {
684: /* this is nonsense, but set each nodal value to phi (will actually do integration over element */
685: array[e[3*i]] = phi;
686: array[e[3*i+1]] = phi;
687: array[e[3*i+2]] = phi;
688: }
689: VecRestoreArray(blocal, (PetscScalar **) &array);
690: DMDARestoreElements(da,&ne,&nc,&e);
692: /* add our partial sums over all processors into b */
693: DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
694: DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);
695: DMRestoreLocalVector(da,&blocal);
696: return(0);
697: }
701: /*
702: We integrate over each cell
704: (i, j+1)----(i+1, j+1)
705: | \ |
706: | \ |
707: | \ |
708: | \ |
709: | \ |
710: | \ |
711: | \ |
712: (i, j)----(i+1, j)
714: However, the element stiffness matrix for the identity in linear elements is
716: 1 /2 1 1\
717: - |1 2 1|
718: 12 \1 1 2/
720: no matter what the shape of the triangle. The Laplacian stiffness matrix is
722: 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)\
723: - |-(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)|
724: 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 /
726: where A is the area of the triangle, and (x_i, y_i) is its i'th vertex.
727: */
728: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J,Mat jac)
729: {
730: DM da = dmmg->dm;
731: UserContext *user = (UserContext *) dmmg->user;
732: /* not being used!
733: PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
734: 0.08333333333, 0.16666666667, 0.08333333333,
735: 0.08333333333, 0.08333333333, 0.16666666667};
736: */
737: PetscScalar values[3][3];
738: PetscInt idx[3];
739: PetscScalar hx, hy, hx2, hy2, area,phi_dt2;
740: PetscInt i,mx,my,xm,ym,xs,ys;
741: PetscInt ne,nc;
742: const PetscInt *e;
746: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
747: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
748: hx = 1.0 / (mx-1);
749: hy = 1.0 / (my-1);
750: area = 0.5*hx*hy;
751: phi_dt2 = user->phi*user->dt*user->dt;
752: hx2 = hx*hx/area*phi_dt2;
753: hy2 = hy*hy/area*phi_dt2;
755: /* initially all elements have identical geometry so all element stiffness are identical */
756: values[0][0] = hx2 + hy2; values[0][1] = -hy2; values[0][2] = -hx2;
757: values[1][0] = -hy2; values[1][1] = hy2; values[1][2] = 0.0;
758: values[2][0] = -hx2; values[2][1] = 0.0; values[2][2] = hx2;
760: DMDAGetElements(da,&ne,&nc,&e);
761: for (i=0; i<ne; i++) {
762: idx[0] = e[3*i];
763: idx[1] = e[3*i+1];
764: idx[2] = e[3*i+2];
765: MatSetValuesLocal(jac,3,idx,3,idx,(PetscScalar*)values,ADD_VALUES);
766: }
767: DMDARestoreElements(da,&ne,&nc,&e);
768: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
769: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
770: return(0);
771: }
775: PetscErrorCode ComputeCorrector(DMMG dmmg, Vec uOld, Vec u)
776: {
777: DM da = dmmg->dm;
778: Vec uOldLocal, uLocal;
779: PetscScalar *cOld;
780: PetscScalar *c;
781: PetscInt i,ne,nc;
782: const PetscInt *e;
784:
786: VecSet(u,0.0);
787: DMGetLocalVector(da, &uOldLocal);
788: DMGetLocalVector(da, &uLocal);
789: VecSet(uLocal,0.0);
790: DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);
791: DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);
792: VecGetArray(uOldLocal, &cOld);
793: VecGetArray(uLocal, &c);
795: /* access the list of elements on this processor and loop over them */
796: DMDAGetElements(da,&ne,&nc,&e);
797: for (i=0; i<ne; i++) {
799: /* this is nonsense, but copy each nodal value*/
800: c[e[3*i]] = cOld[e[3*i]];
801: c[e[3*i+1]] = cOld[e[3*i+1]];
802: c[e[3*i+2]] = cOld[e[3*i+2]];
803: }
804: DMDARestoreElements(da,&ne,&nc,&e);
805: VecRestoreArray(uOldLocal, &cOld);
806: VecRestoreArray(uLocal, &c);
807: DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);
808: DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);
809: DMRestoreLocalVector(da, &uOldLocal);
810: DMRestoreLocalVector(da, &uLocal);
811: return(0);
812: }