Actual source code: ex31.c
2: static char help[] = "A Chebyshev spectral method for the compressible Blasius boundary layer equations.\n\n";
4: /*
5: Include "petscsnes.h" so that we can use SNES solvers. Note that this
6: file automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: petscksp.h - linear solvers
12: Include "petscdt.h" so that we can have support for use of Quadrature formulas
13: */
14: /*F
15: This examples solves the compressible Blasius boundary layer equations
16: 2(\rho\muf'')' + ff'' = 0
17: (\rho\muh')' + Prfh' + Pr(\gamma-1)Ma^{2}\rho\muf''^{2} = 0
18: following Howarth-Dorodnitsyn transformation with boundary conditions
19: f(0) = f'(0) = 0, f'(\infty) = 1, h(\infty) = 1, h = \theta(0). Where \theta = T/T_{\infty}
20: Note: density (\rho) and viscosity (\mu) are treated as constants in this example
21: F*/
22: #include <petscsnes.h>
23: #include <petscdt.h>
25: /*
26: User-defined routines
27: */
29: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
31: typedef struct {
32: PetscReal Ma, Pr, h_0;
33: PetscInt N;
34: PetscReal dx_deta;
35: PetscReal *x;
36: PetscReal gamma;
37: } Blasius;
39: int main(int argc, char **argv)
40: {
41: SNES snes; /* nonlinear solver context */
42: Vec x, r; /* solution, residual vectors */
43: PetscMPIInt size;
44: Blasius *blasius;
45: PetscReal L, *weight; /* L is size of the domain */
48: PetscInitialize(&argc, &argv, (char *)0, help);
49: MPI_Comm_size(PETSC_COMM_WORLD, &size);
52: // Read command-line arguments
53: PetscCalloc1(1, &blasius);
54: blasius->Ma = 2; /* Mach number */
55: blasius->Pr = 0.7; /* Prandtl number */
56: blasius->h_0 = 2.; /* relative temperature at the wall */
57: blasius->N = 10; /* Number of Chebyshev terms */
58: blasius->gamma = 1.4; /* specific heat ratio */
59: L = 5;
60: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Compressible Blasius boundary layer equations", "");
61: PetscOptionsReal("-mach", "Mach number at freestream", "", blasius->Ma, &blasius->Ma, NULL);
62: PetscOptionsReal("-prandtl", "Prandtl number", "", blasius->Pr, &blasius->Pr, NULL);
63: PetscOptionsReal("-h_0", "Relative enthalpy at wall", "", blasius->h_0, &blasius->h_0, NULL);
64: PetscOptionsReal("-gamma", "Ratio of specific heats", "", blasius->gamma, &blasius->gamma, NULL);
65: PetscOptionsInt("-N", "Number of Chebyshev terms for f", "", blasius->N, &blasius->N, NULL);
66: PetscOptionsReal("-L", "Extent of the domain", "", L, &L, NULL);
67: PetscOptionsEnd();
68: blasius->dx_deta = 2 / L; /* this helps to map [-1,1] to [0,L] */
69: PetscMalloc2(blasius->N - 3, &blasius->x, blasius->N - 3, &weight);
70: PetscDTGaussQuadrature(blasius->N - 3, -1., 1., blasius->x, weight);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create nonlinear solver context
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: SNESCreate(PETSC_COMM_WORLD, &snes);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create matrix and vector data structures; set corresponding routines
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: /*
81: Create vectors for solution and nonlinear function
82: */
83: VecCreate(PETSC_COMM_WORLD, &x);
84: VecSetSizes(x, PETSC_DECIDE, 2 * blasius->N - 1);
85: VecSetFromOptions(x);
86: VecDuplicate(x, &r);
88: /*
89: Set function evaluation routine and vector.
90: */
91: SNESSetFunction(snes, r, FormFunction, blasius);
92: {
93: KSP ksp;
94: PC pc;
95: SNESGetKSP(snes, &ksp);
96: KSPSetType(ksp, KSPPREONLY);
97: KSPGetPC(ksp, &pc);
98: PCSetType(pc, PCLU);
99: }
100: /*
101: Set SNES/KSP/KSP/PC runtime options, e.g.,
102: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
103: These options will override those specified above as long as
104: SNESSetFromOptions() is called _after_ any other customization
105: routines.
106: */
107: SNESSetFromOptions(snes);
109: SNESSolve(snes, NULL, x);
110: //VecView(x,PETSC_VIEWER_STDOUT_WORLD);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Free work space. All PETSc objects should be destroyed when they
114: are no longer needed.
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscFree2(blasius->x, weight);
118: PetscFree(blasius);
119: VecDestroy(&x);
120: VecDestroy(&r);
121: SNESDestroy(&snes);
122: PetscFinalize();
123: return 0;
124: }
126: /*
127: Helper function to evaluate Chebyshev polynomials with a set of coefficients
128: with all their derivatives represented as a recurrence table
129: */
130: static void ChebyshevEval(PetscInt N, const PetscScalar *Tf, PetscReal x, PetscReal dx_deta, PetscScalar *f)
131: {
132: PetscScalar table[4][3] = {
133: {1, x, 2 * x * x - 1},
134: {0, 1, 4 * x },
135: {0, 0, 4 },
136: {0, 0, 0 } /* Chebyshev polynomials T_0, T_1, T_2 of the first kind in (-1,1) */
137: };
138: for (int i = 0; i < 4; i++) { f[i] = table[i][0] * Tf[0] + table[i][1] * Tf[1] + table[i][2] * Tf[2]; /* i-th derivative of f */ }
139: for (int i = 3; i < N; i++) {
140: table[0][i % 3] = 2 * x * table[0][(i - 1) % 3] - table[0][(i - 2) % 3]; /* T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x) */
141: /* Differentiate Chebyshev polynomials with the recurrence relation */
142: for (int j = 1; j < 4; j++) { table[j][i % 3] = i * (2 * table[j - 1][(i - 1) % 3] + table[j][(i - 2) % 3] / (i - 2)); /* T'_{n}(x)/n = 2T_{n-1}(x) + T'_{n-2}(x)/n-2 */ }
143: for (int j = 0; j < 4; j++) f[j] += table[j][i % 3] * Tf[i];
144: }
145: for (int i = 1; i < 4; i++) {
146: for (int j = 0; j < i; j++) f[i] *= dx_deta; /* Here happens the physics of the problem */
147: }
148: }
150: /*
151: FormFunction - Evaluates nonlinear function, F(x).
153: Input Parameters:
154: . snes - the SNES context
155: . X - input vector
156: . ctx - optional user-defined context
158: Output Parameter:
159: . R - function vector
160: */
161: PetscErrorCode FormFunction(SNES snes, Vec X, Vec R, void *ctx)
162: {
163: Blasius *blasius = (Blasius *)ctx;
164: const PetscScalar *Tf, *Th; /* Tf and Th are Chebyshev coefficients */
165: PetscScalar *r, f[4], h[4];
166: PetscInt N = blasius->N;
167: PetscReal Ma = blasius->Ma, Pr = blasius->Pr;
169: /*
170: Get pointers to vector data.
171: - For default PETSc vectors, VecGetArray() returns a pointer to
172: the data array. Otherwise, the routine is implementation dependent.
173: - You MUST call VecRestoreArray() when you no longer need access to
174: the array.
175: */
176: VecGetArrayRead(X, &Tf);
177: Th = Tf + N;
178: VecGetArray(R, &r);
180: /* Compute function */
181: ChebyshevEval(N, Tf, -1., blasius->dx_deta, f);
182: r[0] = f[0];
183: r[1] = f[1];
184: ChebyshevEval(N, Tf, 1., blasius->dx_deta, f);
185: r[2] = f[1] - 1; /* Right end boundary condition */
186: for (int i = 0; i < N - 3; i++) {
187: ChebyshevEval(N, Tf, blasius->x[i], blasius->dx_deta, f);
188: r[3 + i] = 2 * f[3] + f[2] * f[0];
189: ChebyshevEval(N - 1, Th, blasius->x[i], blasius->dx_deta, h);
190: r[N + 2 + i] = h[2] + Pr * f[0] * h[1] + Pr * (blasius->gamma - 1) * PetscSqr(Ma * f[2]);
191: }
192: ChebyshevEval(N - 1, Th, -1., blasius->dx_deta, h);
193: r[N] = h[0] - blasius->h_0; /* Left end boundary condition */
194: ChebyshevEval(N - 1, Th, 1., blasius->dx_deta, h);
195: r[N + 1] = h[0] - 1; /* Left end boundary condition */
197: /* Restore vectors */
198: VecRestoreArrayRead(X, &Tf);
199: VecRestoreArray(R, &r);
200: return 0;
201: }
203: /*TEST
205: test:
206: args: -snes_monitor -pc_type svd
207: requires: !single
209: TEST*/