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

 47:   PetscFunctionBeginUser;
 48:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 49:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 50:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");

 52:   // Read command-line arguments
 53:   PetscCall(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:   PetscCall(PetscOptionsReal("-mach", "Mach number at freestream", "", blasius->Ma, &blasius->Ma, NULL));
 62:   PetscCall(PetscOptionsReal("-prandtl", "Prandtl number", "", blasius->Pr, &blasius->Pr, NULL));
 63:   PetscCall(PetscOptionsReal("-h_0", "Relative enthalpy at wall", "", blasius->h_0, &blasius->h_0, NULL));
 64:   PetscCall(PetscOptionsReal("-gamma", "Ratio of specific heats", "", blasius->gamma, &blasius->gamma, NULL));
 65:   PetscCall(PetscOptionsInt("-N", "Number of Chebyshev terms for f", "", blasius->N, &blasius->N, NULL));
 66:   PetscCall(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:   PetscCall(PetscMalloc2(blasius->N - 3, &blasius->x, blasius->N - 3, &weight));
 70:   PetscCall(PetscDTGaussQuadrature(blasius->N - 3, -1., 1., blasius->x, weight));

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:      Create nonlinear solver context
 74:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 75:   PetscCall(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:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 84:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 2 * blasius->N - 1));
 85:   PetscCall(VecSetFromOptions(x));
 86:   PetscCall(VecDuplicate(x, &r));

 88:   /*
 89:       Set function evaluation routine and vector.
 90:    */
 91:   PetscCall(SNESSetFunction(snes, r, FormFunction, blasius));
 92:   {
 93:     KSP ksp;
 94:     PC  pc;
 95:     PetscCall(SNESGetKSP(snes, &ksp));
 96:     PetscCall(KSPSetType(ksp, KSPPREONLY));
 97:     PetscCall(KSPGetPC(ksp, &pc));
 98:     PetscCall(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:   PetscCall(SNESSetFromOptions(snes));

109:   PetscCall(SNESSolve(snes, NULL, x));
110:   //PetscCall(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:   PetscCall(PetscFree2(blasius->x, weight));
118:   PetscCall(PetscFree(blasius));
119:   PetscCall(VecDestroy(&x));
120:   PetscCall(VecDestroy(&r));
121:   PetscCall(SNESDestroy(&snes));
122:   PetscCall(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:   PetscFunctionBeginUser;
170:   /*
171:    Get pointers to vector data.
172:       - For default PETSc vectors, VecGetArray() returns a pointer to
173:         the data array.  Otherwise, the routine is implementation dependent.
174:       - You MUST call VecRestoreArray() when you no longer need access to
175:         the array.
176:    */
177:   PetscCall(VecGetArrayRead(X, &Tf));
178:   Th = Tf + N;
179:   PetscCall(VecGetArray(R, &r));

181:   /* Compute function */
182:   ChebyshevEval(N, Tf, -1., blasius->dx_deta, f);
183:   r[0] = f[0];
184:   r[1] = f[1];
185:   ChebyshevEval(N, Tf, 1., blasius->dx_deta, f);
186:   r[2] = f[1] - 1; /* Right end boundary condition */
187:   for (int i = 0; i < N - 3; i++) {
188:     ChebyshevEval(N, Tf, blasius->x[i], blasius->dx_deta, f);
189:     r[3 + i] = 2 * f[3] + f[2] * f[0];
190:     ChebyshevEval(N - 1, Th, blasius->x[i], blasius->dx_deta, h);
191:     r[N + 2 + i] = h[2] + Pr * f[0] * h[1] + Pr * (blasius->gamma - 1) * PetscSqr(Ma * f[2]);
192:   }
193:   ChebyshevEval(N - 1, Th, -1., blasius->dx_deta, h);
194:   r[N] = h[0] - blasius->h_0; /* Left end boundary condition */
195:   ChebyshevEval(N - 1, Th, 1., blasius->dx_deta, h);
196:   r[N + 1] = h[0] - 1; /* Left end boundary condition */

198:   /* Restore vectors */
199:   PetscCall(VecRestoreArrayRead(X, &Tf));
200:   PetscCall(VecRestoreArray(R, &r));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: /*TEST

206:    test:
207:       args: -snes_monitor -pc_type svd
208:       requires: !single

210: TEST*/