Actual source code: ex31.c

  1: static char help[] = "A Chebyshev spectral method for the compressible Blasius boundary layer equations.\n\n";

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>

24: /*
25:    User-defined routines
26: */

28: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);

30: typedef struct {
31:   PetscReal  Ma, Pr, h_0;
32:   PetscInt   N;
33:   PetscReal  dx_deta;
34:   PetscReal *x;
35:   PetscReal  gamma;
36: } Blasius;

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 */

46:   PetscFunctionBeginUser;
47:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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");

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));

71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72:      Create nonlinear solver context
73:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

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));

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));

108:   PetscCall(SNESSolve(snes, NULL, x));
109:   //PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Free work space.  All PETSc objects should be destroyed when they
113:      are no longer needed.
114:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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: }

125: /*
126:    Helper function to evaluate Chebyshev polynomials with a set of coefficients
127:    with all their derivatives represented as a recurrence table
128: */
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: }

149: /*
150:    FormFunction - Evaluates nonlinear function, F(x).

152:    Input Parameters:
153: .  snes - the SNES context
154: .  X    - input vector
155: .  ctx  - optional user-defined context

157:    Output Parameter:
158: .  R - function vector
159:  */
160: PetscErrorCode FormFunction(SNES snes, Vec X, Vec R, void *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;

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:    */
177:   Th = Tf + N;
178:   PetscCall(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 */