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