1 2 static char help[] = "A Chebyshev spectral method for the compressible Blasius boundary layer equations.\n\n"; 3 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> 24 25 /* 26 User-defined routines 27 */ 28 29 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 30 31 typedef struct { 32 PetscReal Ma, Pr, h_0; 33 PetscInt N; 34 PetscReal dx_deta; 35 PetscReal *x; 36 PetscReal gamma; 37 } Blasius; 38 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 */ 46 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"); 51 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)); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create nonlinear solver context 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 76 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)); 87 88 /* 89 Set function evaluation routine and vector. 90 */ 91 PetscCall(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 PetscCall(SNESSetFromOptions(snes)); 108 109 PetscCall(SNESSolve(snes,NULL,x)); 110 //PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Free work space. All PETSc objects should be destroyed when they 114 are no longer needed. 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 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 } 125 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 PetscScalar table[4][3] = { 132 {1, x, 2*x*x - 1}, {0, 1, 4*x}, {0, 0, 4}, {0, 0, 0} /* Chebyshev polynomials T_0, T_1, T_2 of the first kind in (-1,1) */ 133 }; 134 for (int i=0; i<4; i++) { 135 f[i] = table[i][0] * Tf[0] + table[i][1] * Tf[1] + table[i][2] * Tf[2]; /* i-th derivative of f */ 136 } 137 for (int i=3; i<N; i++) { 138 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) */ 139 /* Differentiate Chebyshev polynomials with the recurrence relation */ 140 for (int j=1; j<4; j++) { 141 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 } 143 for (int j=0; j<4; j++) { 144 f[j] += table[j][i%3] * Tf[i]; 145 } 146 } 147 for (int i=1; i<4; i++) { 148 for (int j=0; j<i; j++) f[i] *= dx_deta; /* Here happens the physics of the problem */ 149 } 150 } 151 152 /* ------------------------------------------------------------------- */ 153 /* 154 FormFunction - Evaluates nonlinear function, F(x). 155 156 Input Parameters: 157 . snes - the SNES context 158 . X - input vector 159 . ctx - optional user-defined context 160 161 Output Parameter: 162 . R - function vector 163 */ 164 PetscErrorCode FormFunction(SNES snes,Vec X,Vec R,void *ctx) 165 { 166 Blasius *blasius = (Blasius *)ctx; 167 const PetscScalar *Tf, *Th; /* Tf and Th are Chebyshev coefficients */ 168 PetscScalar *r, f[4], h[4]; 169 PetscInt N = blasius->N; 170 PetscReal Ma = blasius->Ma, Pr = blasius->Pr; 171 172 /* 173 Get pointers to vector data. 174 - For default PETSc vectors, VecGetArray() returns a pointer to 175 the data array. Otherwise, the routine is implementation dependent. 176 - You MUST call VecRestoreArray() when you no longer need access to 177 the array. 178 */ 179 PetscCall(VecGetArrayRead(X,&Tf)); 180 Th = Tf + N; 181 PetscCall(VecGetArray(R,&r)); 182 183 /* Compute function */ 184 ChebyshevEval(N, Tf, -1., blasius->dx_deta, f); 185 r[0] = f[0]; 186 r[1] = f[1]; 187 ChebyshevEval(N, Tf, 1., blasius->dx_deta, f); 188 r[2] = f[1] - 1; /* Right end boundary condition */ 189 for (int i=0; i<N - 3; i++) { 190 ChebyshevEval(N, Tf, blasius->x[i], blasius->dx_deta, f); 191 r[3+i] = 2*f[3] + f[2] * f[0]; 192 ChebyshevEval(N-1, Th, blasius->x[i], blasius->dx_deta, h); 193 r[N+2+i] = h[2] + Pr * f[0] * h[1] + Pr * (blasius->gamma - 1) * PetscSqr(Ma * f[2]); 194 } 195 ChebyshevEval(N-1, Th, -1., blasius->dx_deta, h); 196 r[N] = h[0] - blasius->h_0; /* Left end boundary condition */ 197 ChebyshevEval(N-1, Th, 1., blasius->dx_deta, h); 198 r[N+1] = h[0] - 1; /* Left end boundary condition */ 199 200 /* Restore vectors */ 201 PetscCall(VecRestoreArrayRead(X,&Tf)); 202 PetscCall(VecRestoreArray(R,&r)); 203 return 0; 204 } 205 206 /*TEST 207 208 test: 209 args: -snes_monitor -pc_type svd 210 requires: !single 211 212 TEST*/ 213