# include # include # include # include # include "navier_stokes_3d_exact.h" /******************************************************************************/ double r8_erf ( double x ) /******************************************************************************/ /* Purpose: R8_ERF evaluates the error function ERF(X). Licensing: This code is distributed under the MIT license. Modified: 24 May 2007 Author: Original FORTRAN77 version by William Cody. C version by John Burkardt. Reference: William Cody, "Rational Chebyshev approximations for the error function", Mathematics of Computation, 1969, pages 631-638. Parameters: Input, double X, the argument of the error function. Output, double R8_ERF, the value of the error function. */ { double a[5] = { 3.16112374387056560, 1.13864154151050156E+02, 3.77485237685302021E+02, 3.20937758913846947E+03, 1.85777706184603153E-01 }; double b[4] = { 2.36012909523441209E+01, 2.44024637934444173E+02, 1.28261652607737228E+03, 2.84423683343917062E+03 }; double c[9] = { 5.64188496988670089E-01, 8.88314979438837594, 6.61191906371416295E+01, 2.98635138197400131E+02, 8.81952221241769090E+02, 1.71204761263407058E+03, 2.05107837782607147E+03, 1.23033935479799725E+03, 2.15311535474403846E-08 }; double d[8] = { 1.57449261107098347E+01, 1.17693950891312499E+02, 5.37181101862009858E+02, 1.62138957456669019E+03, 3.29079923573345963E+03, 4.36261909014324716E+03, 3.43936767414372164E+03, 1.23033935480374942E+03 }; double del; double erfx; int i; double p[6] = { 3.05326634961232344E-01, 3.60344899949804439E-01, 1.25781726111229246E-01, 1.60837851487422766E-02, 6.58749161529837803E-04, 1.63153871373020978E-02 }; double q[5] = { 2.56852019228982242, 1.87295284992346047, 5.27905102951428412E-01, 6.05183413124413191E-02, 2.33520497626869185E-03 }; double sqrpi = 0.56418958354775628695; double thresh = 0.46875; double xbig = 26.543; double xabs; double xden; double xnum; double xsmall = 1.11E-16; double xsq; xabs = fabs ( x ); /* Evaluate ERF(X) for |X| <= 0.46875. */ if ( xabs <= thresh ) { if ( xsmall < xabs ) { xsq = xabs * xabs; } else { xsq = 0.0; } xnum = a[4] * xsq; xden = xsq; for ( i = 0; i < 3; i++ ) { xnum = ( xnum + a[i] ) * xsq; xden = ( xden + b[i] ) * xsq; } erfx = x * ( xnum + a[3] ) / ( xden + b[3] ); } /* Evaluate ERFC(X) for 0.46875 <= |X| <= 4.0. */ else if ( xabs <= 4.0 ) { xnum = c[8] * xabs; xden = xabs; for ( i = 0; i < 7; i++ ) { xnum = ( xnum + c[i] ) * xabs; xden = ( xden + d[i] ) * xabs; } erfx = ( xnum + c[7] ) / ( xden + d[7] ); xsq = ( double ) ( ( int ) ( ( xabs * 16.0 ) / 16.0 ) ); del = ( xabs - xsq ) * ( xabs + xsq ); erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx; erfx = ( 0.5 - erfx ) + 0.5; if ( x < 0.0 ) { erfx = - erfx; } } /* Evaluate ERFC(X) for 4.0 < |X|. */ else { if ( xbig <= xabs ) { if ( 0.0 < x ) { erfx = 1.0; } else { erfx = - 1.0; } } else { xsq = 1.0 / ( xabs * xabs ); xnum = p[5] * xsq; xden = xsq; for ( i = 0; i < 4; i++ ) { xnum = ( xnum + p[i] ) * xsq; xden = ( xden + q[i] ) * xsq; } erfx = xsq * ( xnum + p[4] ) / ( xden + q[4] ); erfx = ( sqrpi - erfx ) / xabs; xsq = ( double ) ( ( int ) ( ( xabs * 16.0 ) / 16.0 ) ); del = ( xabs - xsq ) * ( xabs + xsq ); erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx; erfx = ( 0.5 - erfx ) + 0.5; if ( x < 0.0 ) { erfx = - erfx; } } } return erfx; } /******************************************************************************/ double r8vec_amax ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_AMAX returns the maximum absolute value in an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 31 May 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double A[N], the array. Output, double AMAX, the value of the entry of largest magnitude. */ { double amax; int i; amax = 0.0; for ( i = 0; i < n; i++ ) { if ( amax < fabs ( a[i] ) ) { amax = fabs ( a[i] ); } } return amax; } /******************************************************************************/ double r8vec_amin ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_AMIN returns the minimum absolute value in an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 31 May 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double A[N], the array. Output, double R8VEC_AMIN, the value of the entry of smallest magnitude. */ { double amin; int i; amin = HUGE_VAL; for ( i = 0; i < n; i++ ) { if ( fabs ( a[i] ) < amin ) { amin = fabs ( a[i] ); } } return amin; } /******************************************************************************/ double r8vec_max ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MAX returns the value of the maximum element in a R8VEC. Licensing: This code is distributed under the MIT license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], a pointer to the first entry of the array. Output, double R8VEC_MAX, the value of the maximum element. This is set to 0.0 if N <= 0. */ { int i; double value; if ( n <= 0 ) { value = 0.0; return value; } value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( value < r8vec[i] ) { value = r8vec[i]; } } return value; } /******************************************************************************/ double r8vec_min ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MIN returns the value of the minimum element in a R8VEC. Licensing: This code is distributed under the MIT license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], the array to be checked. Output, double R8VEC_MIN, the value of the minimum element. */ { int i; double value; value = r8vec[0]; for ( i = 1; i < n; i++ ) { if ( r8vec[i] < value ) { value = r8vec[i]; } } return value; } /******************************************************************************/ double *r8vec_uniform_ab_new ( int n, double a, double b, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_AB_NEW returns a scaled pseudorandom R8VEC. Discussion: Each dimension ranges from A to B. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input, double A, B, the lower and upper limits of the pseudorandom values. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_AB_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_AB_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = a + ( b - a ) * ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ double *r8vec_zeros_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_ZEROS_NEW returns a zero R8VEC. Licensing: This code is distributed under the MIT license. Modified: 30 July 2015 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double R8VEC_ZERO_NEW[N], the vector of zeros. */ { int i; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { r[i] = 0.0; } return r; } /******************************************************************************/ void resid_burgers ( double nu, int n, double x[], double y[], double z[], double t[], double ur[], double vr[], double wr[], double pr[] ) /******************************************************************************/ /* Purpose: RESID_BURGERS evaluates the Burgers residual. Licensing: This code is distributed under the MIT license. Modified: 30 July 2015 Author: John Burkardt Reference: Martin Bazant, Henry Moffatt, Exact solutions of the Navier-Stokes equations having steady vortex structures, Journal of Fluid Mechanics, Volume 541, pages 55-64, 2005. Johannes Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics, Volume 1, pages 171-199, 1948. Parameters: Input, double NU, the kinematic viscosity. Input, int N, the number of points at which the solution is to be evaluated. Input, double X[N], Y[N], Z[N], the coordinates of the points. Input, double T[N], the time coordinates. Output, double UR[N], VR[N], WR[N], PR[N], the residuals. */ { int i; double px; double py; double pz; const double r8_pi = 3.141592653589793; double u; double ut; double ux; double uxx; double uy; double uyy; double uz; double uzz; double v; double vt; double vx; double vxx; double vy; double vyy; double vz; double vzz; double w; double wt; double wx; double wxx; double wy; double wyy; double wz; double wzz; for ( i = 0; i < n; i++ ) { /* Form the functions and derivatives. */ u = 2.0 * x[i]; ux = 2.0; uxx = 0.0; uy = 0.0; uyy = 0.0; uz = 0.0; uzz = 0.0; ut = 0.0; v = - 2.0 * y[i]; vx = 0.0; vxx = 0.0; vy = - 2.0; vyy = 0.0; vz = 0.0; vzz = 0.0; vt = 0.0; w = r8_erf ( y[i] / sqrt ( nu ) ); wx = 0.0; wxx = 0.0; wy = 2.0 * sqrt ( 1.0 / nu / r8_pi ) * exp ( - y[i] * y[i] / nu ); wyy = - 4.0 * sqrt ( 1.0 / nu / r8_pi ) * y[i] * exp ( - y[i] * y[i] / nu ) / nu; wz = 0.0; wzz = 0.0; wt = 0.0; /* p = - 2.0 * ( x[i] * x[i] + y[i] * y[i] ); */ px = - 4.0 * x[i]; py = - 4.0 * y[i]; pz = 0.0; /* Evaluate the residuals. */ ur[i] = ut + u * ux + v * uy + w * uz + px - nu * ( uxx + uyy + uzz ); vr[i] = vt + u * vx + v * vy + w * vz + py - nu * ( vxx + vyy + vzz ); wr[i] = wt + u * wx + v * wy + w * wz + pz - nu * ( wxx + wyy + wzz ); pr[i] = ux + vy + wz; } return; } /******************************************************************************/ void resid_ethier ( double a, double d, int n, double x[], double y[], double z[], double t[], double ur[], double vr[], double wr[], double pr[] ) /******************************************************************************/ /* Purpose: RESID_ETHIER evaluates the Ethier residual. Licensing: This code is distributed under the MIT license. Modified: 17 January 2015 Author: John Burkardt Reference: C Ross Ethier, David Steinman, Exact fully 3D Navier-Stokes solutions for benchmarking, International Journal for Numerical Methods in Fluids, Volume 19, Number 5, March 1994, pages 369-375. Parameters: Input, double A, D, the parameters. Sample values are A = PI/4 and D = PI/2. Input, int N, the number of points at which the solution is to be evaluated. Input, double X[N], Y[N], Z[N], the coordinates of the points. Input, double T[N], the time coordinates. Output, double UR[N], VR[N], WR[N], PR[N], the residuals. */ { double cxy; double cyz; double czx; double e2t; double e2x; double e2y; double e2z; double e4t; double ex; double exy; double ey; double eyz; double ez; double ezx; int i; double px; double py; double pz; double sxy; double syz; double szx; double u; double ut; double ux; double uxx; double uy; double uyy; double uz; double uzz; double v; double vt; double vx; double vxx; double vy; double vyy; double vz; double vzz; double w; double wt; double wx; double wxx; double wy; double wyy; double wz; double wzz; /* Make some temporaries. */ for ( i = 0; i < n; i++ ) { ex = exp ( a * x[i] ); ey = exp ( a * y[i] ); ez = exp ( a * z[i] ); e2x = exp ( 2.0 * a * x[i] ); e2y = exp ( 2.0 * a * y[i] ); e2z = exp ( 2.0 * a * z[i] ); e2t = exp ( - d * d * t[i] ); e4t = exp ( - 2.0 * d * d * t[i] ); exy = exp ( a * ( x[i] + y[i] ) ); eyz = exp ( a * ( y[i] + z[i] ) ); ezx = exp ( a * ( z[i] + x[i] ) ); sxy = sin ( a * x[i] + d * y[i] ); syz = sin ( a * y[i] + d * z[i] ); szx = sin ( a * z[i] + d * x[i] ); cxy = cos ( a * x[i] + d * y[i] ); cyz = cos ( a * y[i] + d * z[i] ); czx = cos ( a * z[i] + d * x[i] ); /* Form the functions and derivatives. */ u = - a * ( ex * syz + ez * cxy ) * e2t; ux = - a * ( a * ex * syz - a * ez * sxy ) * e2t; uxx = - a * ( a * a * ex * syz - a * a * ez * cxy ) * e2t; uy = - a * ( a * ex * cyz - d * ez * sxy ) * e2t; uyy = - a * ( - a * a * ex * syz - d * d * ez * cxy ) * e2t; uz = - a * ( d * ex * cyz + a * ez * cxy ) * e2t; uzz = - a * ( - d * d * ex * syz + a * a * ez * cxy ) * e2t; ut = + d * d * a * ( ex * syz + ez * cxy ) * e2t; v = - a * ( ey * szx + ex * cyz ) * e2t; vx = - a * ( d * ey * czx + a * ex * cyz ) * e2t; vxx = - a * ( - d * d * ey * szx + a * a * ex * cyz ) * e2t; vy = - a * ( a * ey * szx - a * ex * syz ) * e2t; vyy = - a * ( a * a * ey * szx - a * a * ex * cyz ) * e2t; vz = - a * ( a * ey * czx - d * ex * syz ) * e2t; vzz = - a * ( - a * a * ey * szx - d * d * ex * cyz ) * e2t; vt = + d * d * a * ( ey * szx + ex * cyz ) * e2t; w = - a * ( ez * sxy + ey * czx ) * e2t; wx = - a * ( a * ez * cxy - d * ey * szx ) * e2t; wxx = - a * ( - a * a * ez * sxy - d * d * ey * czx ) * e2t; wy = - a * ( d * ez * cxy + a * ey * czx ) * e2t; wyy = - a * ( - d * d * ez * sxy + a * a * ey * czx ) * e2t; wz = - a * ( a * ez * sxy - a * ey * szx ) * e2t; wzz = - a * ( a * a * ez * sxy - a * a * ey * czx ) * e2t; wt = + d * d * a * ( ez * sxy + ey * czx ) * e2t; /* p = - 0.5 * a * a * e4t * ( + e2x + 2.0 * sxy * czx * eyz + e2y + 2.0 * syz * cxy * ezx + e2z + 2.0 * szx * cyz * exy ); */ px = - 0.5 * a * a * e4t * ( + 2.0 * a * e2x + 2.0 * a * cxy * czx * eyz - 2.0 * d * sxy * szx * eyz - 2.0 * a * syz * sxy * ezx + 2.0 * a * syz * cxy * ezx + 2.0 * d * czx * cyz * exy + 2.0 * a * szx * cyz * exy ); py = - 0.5 * a * a * e4t * ( + 2.0 * d * cxy * czx * eyz + 2.0 * a * sxy * czx * eyz + 2.0 * a * e2y + 2.0 * a * cyz * cxy * ezx - 2.0 * d * syz * sxy * ezx - 2.0 * a * szx * syz * exy + 2.0 * a * szx * cyz * exy ); pz = - 0.5 * a * a * e4t * ( - 2.0 * a * sxy * szx * eyz + 2.0 * a * sxy * czx * eyz + 2.0 * d * cyz * cxy * ezx + 2.0 * a * syz * cxy * ezx + 2.0 * a * e2z + 2.0 * a * czx * cyz * exy - 2.0 * d * szx * syz * exy ); /* Evaluate the residuals. */ ur[i] = ut + u * ux + v * uy + w * uz + px - ( uxx + uyy + uzz ); vr[i] = vt + u * vx + v * vy + w * vz + py - ( vxx + vyy + vzz ); wr[i] = wt + u * wx + v * wy + w * wz + pz - ( wxx + wyy + wzz ); pr[i] = ux + vy + wz; } return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ void uvwp_burgers ( double nu, int n, double x[], double y[], double z[], double t[], double u[], double v[], double w[], double p[] ) /******************************************************************************/ /* Purpose: UVWP_BURGERS evaluates the Burgers solution. Licensing: This code is distributed under the MIT license. Modified: 30 July 2015 Author: John Burkardt Reference: Martin Bazant, Henry Moffatt, Exact solutions of the Navier-Stokes equations having steady vortex structures, Journal of Fluid Mechanics, Volume 541, pages 55-64, 2005. Johannes Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics, Volume 1, pages 171-199, 1948. Parameters: Input, double NU, the kinematic viscosity. Input, int N, the number of points at which the solution is to be evaluated. Input, double X[N], Y[N], Z[N], the coordinates of the points. Input, double T[N], the time coordinates. Output, double U[N], V[N], W[N], P[N], the solution values. */ { int i; for ( i = 0; i < n; i++ ) { /* Form the functions and derivatives. */ u[i] = 2.0 * x[i]; v[i] = - 2.0 * y[i]; w[i] = r8_erf ( y[i] / sqrt ( nu ) ); p[i] = - 2.0 * ( x[i] * x[i] + y[i] * y[i] ); } return; } /******************************************************************************/ void uvwp_ethier ( double a, double d, int n, double x[], double y[], double z[], double t[], double u[], double v[], double w[], double p[] ) /******************************************************************************/ /* Purpose: UVWP_ETHIER evaluates the Ethier solution. Discussion: The reference asserts that the given velocity and pressure fields are exact solutions for the 3D incompressible time-dependent Navier Stokes equations over any region. To define a typical problem, one chooses a bounded spatial region and a starting time, and then imposes boundary and initial conditions by referencing the exact solution appropriately. In the reference paper, a calculation is made for the cube centered at (0,0,0) with a "radius" of 1 unit, and over the time interval from t = 0 to t = 0.1, with parameters a = PI/4 and d = PI/2, and with Dirichlet boundary conditions on all faces of the cube. Licensing: This code is distributed under the MIT license. Modified: 17 January 2015 Author: John Burkardt Reference: C Ross Ethier, David Steinman, Exact fully 3D Navier-Stokes solutions for benchmarking, International Journal for Numerical Methods in Fluids, Volume 19, Number 5, March 1994, pages 369-375. Parameters: Input, double A, D, the parameters. Sample values are A = PI/4 and D = PI/2. Input, int N, the number of points at which the solution is to be evaluated. Input, double X[N], Y[N], Z[N], the coordinates of the points. Input, double T[N], the time coordinates. Output, double U[N], V[N], W[N], P[N], the velocity components and pressure at each of the points. */ { double cxy; double cyz; double czx; double e2t; double ex; double exy; double ey; double eyz; double ez; double ezx; int i; double sxy; double syz; double szx; for ( i = 0; i < n; i++ ) { ex = exp ( a * x[i] ); ey = exp ( a * y[i] ); ez = exp ( a * z[i] ); e2t = exp ( - d * d * t[i] ); exy = exp ( a * ( x[i] + y[i] ) ); eyz = exp ( a * ( y[i] + z[i] ) ); ezx = exp ( a * ( z[i] + x[i] ) ); sxy = sin ( a * x[i] + d * y[i] ); syz = sin ( a * y[i] + d * z[i] ); szx = sin ( a * z[i] + d * x[i] ); cxy = cos ( a * x[i] + d * y[i] ); cyz = cos ( a * y[i] + d * z[i] ); czx = cos ( a * z[i] + d * x[i] ); u[i] = - a * ( ex * syz + ez * cxy ) * e2t; v[i] = - a * ( ey * szx + ex * cyz ) * e2t; w[i] = - a * ( ez * sxy + ey * czx ) * e2t; p[i] = 0.5 * a * a * e2t * e2t * ( + ex * ex + 2.0 * sxy * czx * eyz + ey * ey + 2.0 * syz * cxy * ezx + ez * ez + 2.0 * szx * cyz * exy ); } return; }