#! /usr/bin/env python3 # def navier_stokes_2d_exact_test ( ): #*****************************************************************************80 # ## navier_stokes_2d_exact_test() tests navier_stokes_2d_exact(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # from numpy.random import default_rng import numpy as np import platform print ( '' ) print ( 'navier_stokes_2d_exact_test():' ) print ( ' python version: ' + platform.python_version ( ) ) print ( ' numpy version: ' + np.version.version ) print ( ' Test navier_stokes_2d_exact().' ) rng = default_rng ( ) grid_2d_test ( ) # # cavity # uvp_cavity_test ( rng ) uvp_cavity_test2 ( ) rhs_cavity_test ( rng ) resid_cavity_test ( rng ) velocity_cavity_matplotlib_test ( ) # # exppoly # uvp_exppoly_test ( rng ) uvp_exppoly_test2 ( ) rhs_exppoly_test ( rng ) resid_exppoly_test ( rng ) velocity_exppoly_matplotlib_test ( ) # # exptrig # uvp_exptrig_test ( rng ) uvp_exptrig_test2 ( ) rhs_exptrig_test ( rng ) resid_exptrig_test ( rng ) velocity_exptrig_matplotlib_test ( ) # # gms # uvp_gms_test ( rng ) uvp_gms_test2 ( ) rhs_gms_test ( rng ) resid_gms_test ( rng ) velocity_gms_matplotlib_test ( ) # # Lukas Bystricky # uvp_lukas_test ( rng ) uvp_lukas_test2 ( ) rhs_lukas_test ( rng ) resid_lukas_test ( rng ) velocity_lukas_matplotlib_test ( ) # # noflow # uvp_noflow_test ( rng ) uvp_noflow_test2 ( ) rhs_noflow_test ( rng ) resid_noflow_test ( rng ) pressure_noflow_matplotlib_test ( ) # # Poiseuille # uvp_poiseuille_test ( rng ) uvp_poiseuille_test2 ( ) rhs_poiseuille_test ( rng ) resid_poiseuille_test ( rng ) velocity_poiseuille_matplotlib_test ( ) parameter_poiseuille_test ( rng ) # # Spiral # uvp_spiral_test ( rng ) uvp_spiral_test2 ( ) rhs_spiral_test ( rng ) resid_spiral_test ( rng ) velocity_spiral_matplotlib_test ( ) parameter_spiral_test ( rng ) # # Taylor # uvp_taylor_test ( rng ) uvp_taylor_test2 ( ) rhs_taylor_test ( rng ) resid_taylor_test ( rng ) velocity_taylor_matplotlib_test ( ) parameter_taylor_test ( rng ) # # vortex # uvp_vortex_test ( rng ) uvp_vortex_test2 ( ) rhs_vortex_test ( rng ) resid_vortex_test ( rng ) velocity_vortex_matplotlib_test ( ) parameter_vortex_test ( rng ) # # Terminate. # print ( '' ) print ( 'navier_stokes_2d_exact_test():' ) print ( ' Normal end of execution.' ) return def all_cavity ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## all_cavity() evaluates the variables of the cavity flow. # # Discussion: # # This is NOT the standard driven cavity problem, in which the lid # boundary condition sets the horizontal velocity to 1, and the # vertical velocity to 0. # # f(x) = x^4-2x^3+x^2 # g(y) = y^4 - y^2 # F(x) = 0.2x^5 - 0.5x^4 + x^3/3 # F1(x) = - 4x^6 + 12x^5 - 14x^4 + 8x^3 - 2x^2 # F2(x) = 0.5 * f(x)^2 = 0.5x^8 +2x^7+3x^6+2x^5+0.5*x^4 # G1(y) = -24y^5 + 8y^3 -4y # # u = 8 f(x)g'(y) # v = - 8 f'(x)g(y) # p = 8/RE [F(x)g''(y) + f'(x)g'(y)] + 64 F2(x)*(g(y)g''(y)-g'(y)^2) # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Reference: # # A C Kengni Jotsa, V A Pennati, # A cost-effective FE method for 2D Navier-Stokes equations, # Engineering Applications of Computational Fluid Mechanics, # Volume 9, Number 1, pages 66-93, 2015. # # Tien-Mo Shih, C H Tan, B C Hwang, # Effects of grid staggering on numerical schemes, # International Journal for Numerical Methods of Fluids, # Volume 9, Number 2, pages 193-212, February 1989. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N), # the horizontal velocity values. # # real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N), # the vertical velocity values. # # real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N), # the pressure values. Note that DPDXX and DPDYY are returned as # zero because life is short. # import numpy as np fint = 0.2 * x**5 - 0.5 * x**4 + x**3 / 3.0 f = x**4 - 2.0 * x**3 + x**2 fx = 4.0 * x**3 - 6.0 * x**2 + 2.0 * x fxx = 12.0 * x**2 - 12.0 * x + 2.0 fxxx = 24.0 * x - 12.0 g = y**4 - y**2 gy = 4.0 * y**3 - 2.0 * y gyy = 12.0 * y**2 - 2.0 gyyy = 24.0 * y u = 8.0 * f * gy dudt = np.zeros ( n ) dudx = 8.0 * fx * gy dudxx = 8.0 * fxx * gy dudy = 8.0 * f * gyy dudyy = 8.0 * f * gyyy v = - 8.0 * fx * g dvdt = np.zeros ( n ) dvdx = - 8.0 * fxx * g dvdxx = - 8.0 * fxxx * g dvdy = - 8.0 * fx * gy dvdyy = - 8.0 * fx * gyy p = ( rho / nu ) * 8.0 * ( \ ( 0.2 * x**5 - 0.5 * x**4 + x**3 / 3.0 ) * ( 12.0 * y**2 - 2.0 ) \ + ( 4.0 * x**3 - 6.0 * x**2 + 2.0 * x ) * ( 4.0 * y**3 - 2.0 * y ) ) \ + 32.0 * ( x**8 + 4.0 * x**7 + 6.0 * x**6 + 4.0 * x**5 + x**4 ) \ * ( - 4.0 * y**6 + 2.0 * y**4 - 2.0 * y**2 ) dpdt = np.zeros ( n ) dpdx = ( rho / nu ) * 8.0 * ( \ ( x**4 - 2.0 * x**2 + x**2 ) * ( 12.0 * y**2 - 2.0 ) \ + ( 12.0 * x**2 - 12.0 * x + 2.0 ) * ( 4.0 * y**3 - 2.0 * y ) ) \ + 32.0 * ( 8.0 * x**7 + 28.0 * x**6 + 36.0 * x**5 + 20.0 * x**4 + 4.0 * x**3 ) \ * ( - 4.0 * y**6 + 2.0 * y**4 - 2.0 * y**2 ) dpdxx = np.zeros ( n ) dpdy = ( rho / nu ) * 8.0 * ( \ ( 0.2 * x**5 - 0.5 * x**4 + x**3 / 3.0 ) * ( 24.0 * y ) \ + ( 4.0 * x**3 - 6.0 * x**2 + 2.0 * x ) * ( 12.0 * y**2 - 2.0 ) ) \ + 32.0 * ( x**8 + 4.0 * x**7 + 6.0 * x**6 + 4.0 * x**5 + x**4 ) \ * ( - 24.0 * y**5 + 8.0 * y**3 - 4.0 * y ) dpdyy = np.zeros ( n ) us = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) vs = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) ps = dudx + dvdy return \ u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps def all_exppoly ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## all_exppoly() evaluates the variables of the exppoly flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Reference: # # Xiaoli Li, Jie Shen, # Error analysis of the SAC-MAC scheme for the Navier-Stokes equations, # arXiv:1909.05131v1 [math.NA] 8 Sep 2019 # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N), # the horizontal velocity values. # # real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N), # the vertical velocity values. # # real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N), # the pressure values. Note that DPDXX and DPDYY are returned as # zero because life is short. # import numpy as np u = - np.exp ( t ) * ( x**4 - 2.0 * x**3 + x**2 ) \ * ( 2.0 * y**3 - 3.0 * y**2 + y ) / 256.0 dudt = - np.exp ( t ) * ( x**4 - 2.0 * x**3 + x**2 ) \ * ( 2.0 * y**3 - 3.0 * y**2 + y ) / 256.0 dudx = - np.exp ( t ) * ( 4.0 * x**3 - 6.0 * x**2 + 2.0 * x ) \ * ( 2.0 * y**3 - 3.0 * y**2 + y ) / 256.0 dudxx = - np.exp ( t ) * ( 12.0 * x**2 - 12.0 * x + 2.0 ) \ * ( 2.0 * y**3 - 3.0 * y**2 + y ) / 256.0 dudy = - np.exp ( t ) * ( x**4 - 2.0 * x**3 + x**2 ) \ * ( 6.0 * y**2 - 6.0 * y + 1.0 ) / 256.0 dudyy = - np.exp ( t ) * ( x**4 - 2.0 * x**3 + x**2 ) \ * ( 12.0 * y - 6.0 ) / 256.0 v = np.exp ( t ) * ( 2.0 * x**3 - 3.0 * x**2 + x ) \ * ( y**4 - 2.0 * y**3 + y**2 ) / 256.0 dvdt = np.exp ( t ) * ( 2.0 * x**3 - 3.0 * x**2 + x ) \ * ( y**4 - 2.0 * y**3 + y**2 ) / 256.0 dvdx = np.exp ( t ) * ( 6.0 * x**2 - 6.0 * x + 1.0 ) \ * ( y**4 - 2.0 * y**3 + y**2 ) / 256.0 dvdxx = np.exp ( t ) * ( 12.0 * x - 6.0 * x ) \ * ( y**4 - 2.0 * y**3 + y**2 ) / 256.0 dvdy = np.exp ( t ) * ( 2.0 * x**3 - 3.0 * x**2 + x ) \ * ( 4.0 * y**3 - 6.0 * y**2 + 2.0 * y ) / 256.0 dvdyy = np.exp ( t ) * ( 2.0 * x**3 - 3.0 * x**2 + x ) \ * ( 12.0 * y**2 - 12.0 * y**2 + 2.0 ) / 256.0 p = np.exp ( t ) * ( np.sin ( np.pi * y ) - 2.0 / np.pi ) dpdt = np.exp ( t ) * ( np.sin ( np.pi * y ) - 2.0 / np.pi ) dpdx = np.zeros_like ( x ) dpdxx = np.zeros_like ( x ) dpdy = np.exp ( t ) * np.pi * np.cos ( np.pi * y ) dpdyy = - np.exp ( t ) * np.pi**2 * np.sin ( np.pi * y ) us = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) vs = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) ps = dudx + dvdy return \ u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps def all_exptrig ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## all_exptrig() evaluates the variables of the exptrig flow. # # Discussion: # # The exptrig flow is time dependent, trigonometric in space, exponential # growth in time, zero velocity boundary conditions on unit square at t=1, # this flow forms a spiral. At later times, the exponential growth seems # to make the solution physically absurd and computationally intractable. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Reference: # # Xiaoli Li, Jie Shen, # Error analysis of the SAC-MAC scheme for the Navier-Stokes equations, # arXiv:1909.05131v1 [math.NA] 8 Sep 2019 # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N), # the horizontal velocity values. # # real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N), # the vertical velocity values. # # real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N), # the pressure values. Note that DPDXX and DPDYY are returned as # zero because life is short. # import numpy as np u = np.exp ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) dudt = np.exp ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) dudx = np.pi * np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dudxx = 2.0 * np.pi**2 * np.exp ( t ) * np.cos ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dudy = 2.0 * np.pi * np.exp ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.cos ( 2.0 * np.pi * y ) dudyy = - 4.0 * np.pi**2 * np.exp ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) v = - np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdt = - np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdx = - 2.0 * np.pi * np.exp ( t ) * np.cos ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdxx = 4.0 * np.pi**2 * np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdy = - np.pi * np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dvdyy = - 2.0 * np.pi**2 * np.exp ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.cos ( 2.0 * np.pi * y ) p = rho * np.exp ( t ) * ( np.sin ( np.pi * y ) - 2.0 / np.pi ) dpdt = rho * np.exp ( t ) * ( np.sin ( np.pi * y ) - 2.0 / np.pi ) dpdx = np.zeros_like ( x ) dpdxx = np.zeros_like ( x ) dpdy = np.pi * rho * np.exp ( t ) * np.cos ( np.pi * y ) dpdyy = - np.pi**2 * rho * np.exp ( t ) * np.sin ( np.pi * x ) us = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) vs = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) ps = dudx + dvdy return \ u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps def all_gms ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## all_gms() evaluates the variables of the gms flow. # # Discussion: # # The flow has been modified by a sign change that makes it slightly # more plausible physically. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 August 2020 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N), # the horizontal velocity values. # # real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N), # the vertical velocity values. # # real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N), # the pressure values. # import numpy as np s = ( -1.0 ) ** ( np.ceil ( x ) + np.ceil ( y ) ) u = np.pi * s * np.sin ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) dudt = np.pi * s * np.cos ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) dudx = np.pi**2 * s * np.sin ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dudxx = 2.0 * np.pi**3 * s * np.sin ( t ) * np.cos ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dudy = 2.0 * np.pi**2 * s * np.sin ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.cos ( 2.0 * np.pi * y ) dudyy = - 4.0 * np.pi**3 * s * np.sin ( t ) * ( np.sin ( np.pi * x ) )**2 \ * np.sin ( 2.0 * np.pi * y ) v = - np.pi * s * np.sin ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdt = - np.pi * s * np.cos ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdx = - 2.0 * np.pi**2 * s * np.sin ( t ) * np.cos ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdxx = 4.0 * np.pi**3 * s * np.sin ( t ) * np.sin ( 2.0 * np.pi * x ) \ * ( np.sin ( np.pi * y ) )**2 dvdy = - np.pi**2 * s * np.sin ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.sin ( 2.0 * np.pi * y ) dvdyy = - 2.0 * np.pi**3 * s * np.sin ( t ) * np.sin ( 2.0 * np.pi * x ) \ * np.cos ( 2.0 * np.pi * y ) p = rho * s * np.sin ( t ) * np.cos ( np.pi * x ) \ * np.sin ( np.pi * y ) dpdt = rho * s * np.cos ( t ) * np.cos ( np.pi * x ) \ * np.sin ( np.pi * y ) dpdx = - np.pi * rho * s * np.sin ( t ) * np.sin ( np.pi * x ) \ * np.sin ( np.pi * y ) dpdxx = - np.pi**2 * rho * s * np.sin ( t ) * np.cos ( np.pi * x ) \ * np.sin ( np.pi * y ) dpdy = np.pi * rho * s * np.sin ( t ) * np.cos ( np.pi * x ) \ * np.cos ( np.pi * y ) dpdyy = - np.pi**2 * rho * s * np.sin ( t ) * np.cos ( np.pi * x ) \ * np.sin ( np.pi * y ) us = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) vs = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) ps = dudx + dvdy return \ u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps def all_noflow ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## all_noflow() evaluates the variables of the NoFlow flow. # # Discussion: # # u = 0 # v = 0 # p = arctan ( x * y ) # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N), # the horizontal velocity values. # # real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N), # the vertical velocity values. # # real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N), # the pressure values. # import numpy as np u = np.zeros ( n ) dudt = np.zeros ( n ) dudx = np.zeros ( n ) dudxx = np.zeros ( n ) dudy = np.zeros ( n ) dudyy = np.zeros ( n ) v = np.zeros ( n ) dvdt = np.zeros ( n ) dvdx = np.zeros ( n ) dvdxx = np.zeros ( n ) dvdy = np.zeros ( n ) dvdyy = np.zeros ( n ) p = np.arctan ( x * y ) dpdt = np.zeros ( n ) dpdx = y / ( x**2 * y**2 + 1.0 ) dpdxx = -2.0 * x * y**3 / ( x**2 * y**2 + 1.0 )**2 dpdy = x / ( x**2 * y**2 + 1.0 ) dpdyy = -2.0 * x**3 * y / ( x**2 * y**2 + 1.0 )**2 us = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) vs = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) ps = dudx + dvdy return \ u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps def grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ): #*****************************************************************************80 # ## grid_2d() returns a regular 2D grid. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 January 2015 # # Author: # # John Burkardt # # Input: # # integer X_NUM, the number of X values to use. # # real X_LO, X_HI, the range of X values. # # integer Y_NUM, the number of Y values to use. # # real Y_LO, Y_HI, the range of Y values. # # Output: # # real X(X_NUM*Y_NUM), Y(X_NUM*Y_NUM), the coordinates of the grid. # import numpy as np x = np.zeros ( x_num * y_num ) y = np.zeros ( x_num * y_num ) if ( x_num == 1 ): xi = ( x_lo + x_hi ) / 2.0 k = 0 for j in range ( 0, y_num ): for i in range ( 0, x_num ): x[k] = xi k = k + 1 else: k = 0 for j in range ( 0, y_num ): for i in range ( 0, x_num ): xi = ( float ( x_num - i - 1 ) * x_lo \ + float ( i ) * x_hi ) \ / float ( x_num - 1 ) x[k] = xi k = k + 1 if ( y_num == 1 ): yj = ( y_lo + y_hi ) / 2.0 k = 0 for j in range ( 0, y_num ): for i in range ( 0, x_num ): y[k] = yj k = k + 1 else: k = 0 for j in range ( 0, y_num ): yj = ( float ( y_num - j - 1 ) * y_lo \ + float ( j ) * y_hi ) \ / float ( y_num - 1 ) for i in range ( 0, x_num ): y[k] = yj k = k + 1 return x, y def grid_2d_test ( ): #*****************************************************************************80 # ## grid_2d_test() makes a small 2D grid. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'grid_2d_test():' ) print ( ' grid_2d() generates a regular grid.' ) x_lo = 10.0 x_hi = 20.0 x_num = 5 y_lo = 5.0 y_hi = 6.0 y_num = 3 [ x, y ] = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) print ( '' ) k = 0 for j in range ( 0, y_num ): for i in range ( 0, x_num): print ( ' %2d %2d %2d %14.6f %14.6f' % ( k, i, j, x[k], y[k] ) ) k = k + 1 return def parameter_poiseuille_test ( rng ): #*****************************************************************************80 # ## parameter_poiseuille_test(): Poiseuille solution norms for various values of NU, RHO. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np print ( '' ) print ( 'parameter_poiseuille_test():' ) print ( ' Poiseuille Flow' ) print ( ' Monitor solution norms for various' ) print ( ' values of NU, RHO.' ) n = 1000 x_lo = 0.0 x_hi = 6.0 y_lo = -1.0 y_hi = +1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) # # Vary RHO. # print ( '' ) print ( ' RHO affects the pressure scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for j in range ( 0, 3 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_poiseuille ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) rho = rho / 100.0 # # Vary NU. # print ( '' ) print ( ' NU affects the time scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for i in range ( 0, 4 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_poiseuille ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) nu = nu / 10.0 return def parameter_spiral_test ( rng ): #*****************************************************************************80 # ## parameter_spiral_test(): solution norms over time for various values of NU, RHO. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np print ( '' ) print ( 'parameter_spiral_test():' ) print ( ' Spiral Flow' ) print ( ' Monitor solution norms over time for various' ) print ( ' values of NU, RHO.' ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) # # Vary RHO. # print ( '' ) print ( ' RHO affects the pressure scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for j in range ( 0, 3 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_spiral ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) rho = rho / 100.0 # # Vary NU. # print ( '' ) print ( ' NU affects the time scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for i in range ( 0, 4 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_spiral ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) nu = nu / 10.0 return def parameter_taylor_test ( rng ): #*****************************************************************************80 # ## parameter_taylor_test() monitors Taylor solution norms for various NU, RHO. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np print ( '' ) print ( 'parameter_taylor_test():' ) print ( ' Taylor Flow' ) print ( ' Monitor solution norms over time for various' ) print ( ' values of NU, RHO.' ) n = 1000 x_lo = 0.5 x_hi = 2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) # # Vary RHO. # print ( '' ) print ( ' RHO affects the pressure scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for j in range ( 0, 3 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_taylor ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) rho = rho / 100.0 # # Vary NU. # print ( '' ) print ( ' NU affects the time scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for i in range ( 0, 4 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_taylor ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) nu = nu / 10.0 return def parameter_vortex_test ( rng ): #*****************************************************************************80 # ## parameter_vortex_test() monitors Vortex solution norms for various values of NU, RHO. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np print ( '' ) print ( 'parameter_vortex_test():' ) print ( ' Vortex Flow' ) print ( ' Monitor solution norms over time for various' ) print ( ' values of NU, RHO.' ) n = 1000 x_lo = 0.5 x_hi = 2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) # # Vary RHO. # print ( '' ) print ( ' RHO affects the pressure scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for j in range ( 0, 3 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_vortex ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) rho = rho / 100.0 # # Vary NU. # print ( '' ) print ( ' NU affects the time scaling.' ) print ( '' ) print ( ' RHO NU T ||U|| ||V|| ||P||' ) print ( '' ) nu = 1.0 rho = 1.0 for i in range ( 0, 4 ): for k in range ( 0, 6 ): t = k / 5.0 u, v, p = uvp_vortex ( nu, rho, n, x, y, t ) u_norm = np.linalg.norm ( u ) / n v_norm = np.linalg.norm ( v ) / n p_norm = np.linalg.norm ( p ) / n print ( ' %10.4g %10.4g %8.4g %10.4g %10.4g %10.4g' \ % ( rho, nu, t, u_norm, v_norm, p_norm ) ) print ( '' ) nu = nu / 10.0 return def pressure_matplotlib ( header, n, x, y, p ): #*****************************************************************************80 # ## pressure_matplotlib() plots a pressure field with matplotlib. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # string HEADER, a header to be used to name the files. # # integer N, the number of evaluation points. # # real X(N), Y(N), the coordinates of the evaluation points. # # real P(N), the pressure samples. # from matplotlib import cm import matplotlib.pyplot as plt import numpy as np # # Idiotic contour code requires 2D arrays. # nroot = int ( np.sqrt ( n ) ) x = x.reshape ( [ nroot, nroot ] ) y = y.reshape ( [ nroot, nroot ] ) p = p.reshape ( [ nroot, nroot ] ) # # Form the figure. # levels = 15 ax = plt.figure ( ) plt.contourf ( x, y, p, levels, cmap = cm.coolwarm ) plt.title ( header, fontsize = 16 ) plt.xlabel ( '<--- X --->', fontsize = 16 ) plt.ylabel ( '<--- Y --->', fontsize = 16 ) plt.axis ( 'equal' ) filename = 'pressure_' + header + '_matplotlib.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( block = False ) plt.close ( ) return def pressure_noflow_matplotlib_test ( ): #*****************************************************************************80 # ## pressure_noflow_matplotlib_test() plots a NoFlow pressure field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'pressure_noflow_matplotlib_test():' ) print ( ' Generate a NoFlow pressure field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = -1.0 x_hi = 1.0 x_num = 21 y_lo = -1.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) n = x_num * y_num u, v, p = uvp_noflow ( nu, rho, n, x, y, t ) header = 'noflow' pressure_matplotlib ( header, n, x, y, p ) return def r8vec_print ( n, a, title ): #*****************************************************************************80 # ## r8vec_print() prints an R8VEC. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Input: # # integer N, the dimension of the vector. # # real A(N), the vector to be printed. # # string TITLE, a title. # print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( '%6d: %12g' % ( i, a[i] ) ) return def resid_cavity ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_cavity() evaluates the residuals of the cavity flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real UR(N), VR(N), PR(N), the residuals. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_cavity ( nu, rho, n, x, y, t ) ur = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) - us vr = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) - vs pr = dudx + dvdy - ps return ur, vr, pr def resid_cavity_test ( rng ): #*****************************************************************************80 # ## resid_cavity_test() samples cavity flow residuals at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'resid_cavity_test():' ) print ( ' cavity flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' over the [0,+1]x[0,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) ur, vr, pr = resid_cavity ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_exppoly ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_exppoly() evaluates the residuals of the exppoly flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real UR(N), VR(N), PR(N), the residuals. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exppoly ( nu, rho, n, x, y, t ) ur = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) - us vr = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) - vs pr = dudx + dvdy - ps return ur, vr, pr def resid_exppoly_test ( rng ): #*****************************************************************************80 # ## resid_exppoly_test() samples exppoly residuals at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'resid_exppoly_test():' ) print ( ' exppoly flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) ur, vr, pr = resid_exppoly ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_exptrig ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_exptrig() evaluates the residuals of the exptrig flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real UR(N), VR(N), PR(N), the residuals. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exptrig ( nu, rho, n, x, y, t ) ur = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) - us vr = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) - vs pr = dudx + dvdy - ps return ur, vr, pr def resid_exptrig_test ( rng ): #*****************************************************************************80 # ## resid_exptrig_test() samples exptrig residuals at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'resid_exptrig_test():' ) print ( ' exptrig flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) ur, vr, pr = resid_exptrig ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_gms ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_gms() evaluates the residuals of the GMS flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real UR(N), VR(N), PR(N), the residuals. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_gms ( nu, rho, n, x, y, t ) ur = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) - us vr = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) - vs pr = dudx + dvdy - ps return ur, vr, pr def resid_gms_test ( rng ): #*****************************************************************************80 # ## resid_gms_test() samples GMS residuals at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'resid_gms_test():' ) print ( ' gms flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = 1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) ur, vr, pr = resid_gms ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_lukas ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_lukas() returns residuals of the Lukas Bystricky Flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real UR(N), VR(N), PR(N), the residuals in the U, V and P equations. # import numpy as np ur = np.zeros ( n ) vr = np.zeros ( n ) pr = np.zeros ( n ) # # Get the right hand side functions. # f, g, h = rhs_lukas ( nu, rho, n, x, y, t ); # # Form the functions and derivatives for the left hand side. # u = - np.cos ( np.pi * x ) / np.pi dudt = np.zeros ( n ) dudx = np.sin ( np.pi * x ) dudxx = np.pi * np.cos ( np.pi * x ) dudy = np.zeros ( n ) dudyy = np.zeros ( n ) v = - y * np.sin ( np.pi * x ) dvdt = np.zeros ( n ) dvdx = - np.pi * y * np.cos ( np.pi * x ) dvdxx = + np.pi * np.pi * y * np.sin ( np.pi * x ) dvdy = - np.sin ( np.pi * x ) dvdyy = np.zeros ( n ) p = np.zeros ( n ) dpdx = np.zeros ( n ) dpdy = np.zeros ( n ) # # Evaluate the residuals. # ur = dudt - nu * ( dudxx + dudyy ) + u * dudx + v * dudy + dpdx / rho - f vr = dvdt - nu * ( dvdxx + dvdyy ) + u * dvdx + v * dvdy + dpdy / rho - g pr = dudx + dvdy - h return ur, vr, pr def resid_lukas_test ( rng ): #*****************************************************************************80 # ## resid_lukas_test() samples Lukas Bystricky Flow residuals at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'resid_lukas_test():' ) print ( ' Lukas Bystricky flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 ur, vr, pr = resid_lukas ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_noflow ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_noflow() evaluates the residuals of the NoFlow flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real UR(N), VR(N), PR(N), the residuals. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_noflow ( nu, rho, n, x, y, t ) ur = dudt + u * dudx + v * dudy + dpdx / rho - nu * ( dudxx + dudyy ) - us vr = dvdt + u * dvdx + v * dvdy + dpdy / rho - nu * ( dvdxx + dvdyy ) - vs pr = dudx + dvdy - ps return ur, vr, pr def resid_noflow_test ( rng ): #*****************************************************************************80 # ## resid_noflow_test() samples NoFlow residuals at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'resid_noflow_test():' ) print ( ' NoFlow flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = 1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) ur, vr, pr = resid_noflow ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_poiseuille ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_poiseuille() returns Poiseuille residuals. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real UR(N), VR(N), PR(N), the residuals in the U, V and P equations. # import numpy as np ur = np.zeros ( n ) vr = np.zeros ( n ) pr = np.zeros ( n ) # # Get the right hand side functions. # f, g, h = rhs_poiseuille ( nu, rho, n, x, y, t ); # # Form the functions and derivatives for the left hand side. # u = 1.0 - y ** 2 dudt = 0.0 dudx = 0.0 dudxx = 0.0 dudy = - 2.0 * y dudyy = - 2.0 v = 0.0 dvdt = 0.0 dvdx = 0.0 dvdxx = 0.0 dvdy = 0.0 dvdyy = 0.0 p = - 2.0 * nu * rho * x dpdx = - 2.0 * nu * rho dpdy = 0.0 # # Evaluate the residuals. # ur = dudt - nu * ( dudxx + dudyy ) \ + u * dudx + v * dudy + dpdx / rho - f vr = dvdt - nu * ( dvdxx + dvdyy ) \ + u * dvdx + v * dvdy + dpdy / rho - g pr = dudx + dvdy - h return ur, vr, pr def resid_poiseuille_test ( rng ): #*****************************************************************************80 # ## resid_poiseuille_test() samples the Poiseuille residual. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'resid_poiseuille_test():' ) print ( ' Poiseuille flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' at the initial time T = 0, over a channel region.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 6.0 y_lo = -1.0 y_hi = +1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 ur, vr, pr = resid_poiseuille ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_spiral ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_spiral() returns residuals of the Spiral Flow equations. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Reference: # # Geoffrey Taylor, # On the decay of vortices in a viscous fluid, # Philosophical Magazine, # Volume 46, 1923, pages 671-674. # # Geoffrey Taylor, A E Green, # Mechanism for the production of small eddies from large ones, # Proceedings of the Royal Society of London, # Series A, Volume 158, 1937, pages 499-521. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real UR(N), VR(N), PR(N), the residuals in the U, V and P equations. # import numpy as np ur = np.zeros ( n ) vr = np.zeros ( n ) pr = np.zeros ( n ) # # Get the right hand side functions. # f, g, h = rhs_spiral ( nu, rho, n, x, y, t ); # # Form the functions and derivatives for the left hand side. # u = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudt = nu * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudx = ( 1.0 + nu * t ) * 2.0 \ * ( 4.0 * x ** 3 - 6.0 * x ** 2 + 2.0 * x ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudxx = ( 1.0 + nu * t ) * 2.0 \ * ( 12.0 * x ** 2 - 12.0 * x + 2.0 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudy = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 6.0 * y ** 2 - 6.0 * y + 1.0 ) dudyy = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 12.0 * y - 6.0 ) v = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdt = - nu * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdx = - ( 1.0 + nu * t ) * 2.0 \ * ( 6.0 * x ** 2 - 6.0 * x + 1.0 ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdxx = - ( 1.0 + nu * t ) * 2.0 \ * ( 12.0 * x - 6.0 ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdy = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( 4.0 * y ** 3 - 6.0 * y ** 2 + 2.0 * y ) dvdyy = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( 12.0 * y ** 2 - 12.0 * y + 2.0 ) p = rho * y dpdx = 0.0 dpdy = rho # # Evaluate the residuals. # ur = dudt - nu * ( dudxx + dudyy ) \ + u * dudx + v * dudy + dpdx / rho - f vr = dvdt - nu * ( dvdxx + dvdyy ) \ + u * dvdx + v * dvdy + dpdy / rho - g pr = dudx + dvdy - h return ur, vr, pr def resid_spiral_test ( rng ): #*****************************************************************************80 # ## resid_spiral_test() samples the residuals at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'resid_spiral_test():' ) print ( ' Spiral flow' ) print ( ' Sample the Navier-Stokes residuals' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 ur, vr, pr = resid_spiral ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_taylor ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_taylor() returns Taylor residuals. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Reference: # # Geoffrey Taylor, # On the decay of vortices in a viscous fluid, # Philosophical Magazine, # Volume 46, 1923, pages 671-674. # # Geoffrey Taylor, A E Green, # Mechanism for the production of small eddies from large ones, # Proceedings of the Royal Society of London, # Series A, Volume 158, 1937, pages 499-521. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real UR(N), VR(N), PR(N), the residuals in the U, V and P equations. # import numpy as np # # Get the right hand sides. # f, g, h = rhs_taylor ( nu, rho, n, x, y, t ); # # Make space. # c2x = np.array ( n ) c2y = np.array ( n ) cx = np.array ( n ) cy = np.array ( n ) e2t = np.array ( n ) e4t = np.array ( n ) p = np.array ( n ) px = np.array ( n ) py = np.array ( n ) s2x = np.array ( n ) s2y = np.array ( n ) sx = np.array ( n ) sy = np.array ( n ) u = np.array ( n ) ut = np.array ( n ) ux = np.array ( n ) uxx = np.array ( n ) uy = np.array ( n ) uyy = np.array ( n ) v = np.array ( n ) vt = np.array ( n ) vx = np.array ( n ) vxx = np.array ( n ) vy = np.array ( n ) vyy = np.array ( n ) # # Make some temporaries. # cx = np.cos ( np.pi * x ) cy = np.cos ( np.pi * y ) sx = np.sin ( np.pi * x ) sy = np.sin ( np.pi * y ) e2t = np.exp ( - 2.0 * np.pi * np.pi * nu * t ) c2x = np.cos ( 2.0 * np.pi * x ) c2y = np.cos ( 2.0 * np.pi * y ) s2x = np.sin ( 2.0 * np.pi * x ) s2y = np.sin ( 2.0 * np.pi * y ) e4t = np.exp ( - 4.0 * np.pi * np.pi * nu * t ) # # Form the functions and derivatives. # u = - cx * sy * e2t dudx = np.pi * sx * sy * e2t dudxx = np.pi * np.pi * cx * sy * e2t dudy = - np.pi * cx * cy * e2t dudyy = np.pi * np.pi * cx * sy * e2t dudt = + 2.0 * nu * np.pi * np.pi * cx * sy * e2t v = sx * cy * e2t dvdx = np.pi * cx * cy * e2t dvdxx = - np.pi * np.pi * sx * cy * e2t dvdy = - np.pi * sx * sy * e2t dvdyy = - np.pi * np.pi * sx * cy * e2t dvdt = - 2.0 * nu * np.pi * np.pi * sx * cy * e2t p = - 0.25 * rho * ( c2x + c2y ) * e4t dpdx = + 0.5 * rho * np.pi * s2x * e4t dpdy = + 0.5 * rho * np.pi * s2y * e4t # # Evaluate the residuals. # ur = dudt + u * dudx + v * dudy + ( 1.0 / rho ) * dpdx \ - nu * ( dudxx + dudyy ) - f vr = dvdt + u * dvdx + v * dvdy + ( 1.0 / rho ) * dpdy \ - nu * ( dvdxx + dvdyy ) - g pr = dudx + dvdy - h return ur, vr, pr def resid_taylor_test ( rng ): #*****************************************************************************80 # ## resid_taylor_test() samples the Taylor residual. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'resid_taylor_test():' ) print ( ' Sample the taylor residuals' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 ur, vr, pr = resid_taylor ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def resid_vortex ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## resid_vortex() returns Vortex residuals. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real UR(N), VR(N), PR(N), the residuals in the U, V and P equations. # import numpy as np # # Get the right hand sides. # f, g, h = rhs_vortex ( nu, rho, n, x, y, t ); # # Make space. # c2x = np.array ( n ) c2y = np.array ( n ) cx = np.array ( n ) cy = np.array ( n ) e2t = np.array ( n ) e4t = np.array ( n ) p = np.array ( n ) px = np.array ( n ) py = np.array ( n ) s2x = np.array ( n ) s2y = np.array ( n ) sx = np.array ( n ) sy = np.array ( n ) u = np.array ( n ) ut = np.array ( n ) ux = np.array ( n ) uxx = np.array ( n ) uy = np.array ( n ) uyy = np.array ( n ) v = np.array ( n ) vt = np.array ( n ) vx = np.array ( n ) vxx = np.array ( n ) vy = np.array ( n ) vyy = np.array ( n ) # # Make some temporaries. # cx = np.cos ( np.pi * x ) cy = np.cos ( np.pi * y ) sx = np.sin ( np.pi * x ) sy = np.sin ( np.pi * y ) c2x = np.cos ( 2.0 * np.pi * x ) c2y = np.cos ( 2.0 * np.pi * y ) s2x = np.sin ( 2.0 * np.pi * x ) s2y = np.sin ( 2.0 * np.pi * y ) # # Form the functions and derivatives. # u = - cx * sy dudx = np.pi * sx * sy dudxx = np.pi * np.pi * cx * sy dudy = - np.pi * cx * cy dudyy = np.pi * np.pi * cx * sy dudt = np.zeros ( n ) v = sx * cy dvdx = np.pi * cx * cy dvdxx = - np.pi * np.pi * sx * cy dvdy = - np.pi * sx * sy dvdyy = - np.pi * np.pi * sx * cy dvdt = np.zeros ( n ) p = - 0.25 * rho * ( c2x + c2y ) dpdx = + 0.5 * rho * np.pi * s2x dpdy = + 0.5 * rho * np.pi * s2y # # Evaluate the residuals. # ur = dudt + u * dudx + v * dudy + ( 1.0 / rho ) * dpdx \ - nu * ( dudxx + dudyy ) - f vr = dvdt + u * dvdx + v * dvdy + ( 1.0 / rho ) * dpdy \ - nu * ( dvdxx + dvdyy ) - g pr = dudx + dvdy - h return ur, vr, pr def resid_vortex_test ( rng ): #*****************************************************************************80 # ## resid_vortex_test() samples the Vortex residual. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'resid_vortex_test():' ) print ( ' Sample the vortex residuals' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 ur, vr, pr = resid_vortex ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( np.abs ( ur ) ), np.max ( np.abs ( ur ) ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( np.abs ( vr ) ), np.max ( np.abs ( vr ) ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( np.abs ( pr ) ), np.max ( np.abs ( pr ) ) ) ) return def rhs_cavity ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_cavity() evaluates the source terms of the cavity flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real US(N), VS(N), PS(N), the source terms. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_cavity ( nu, rho, n, x, y, t ) return us, vs, ps def rhs_cavity_test ( rng ): #*****************************************************************************80 # ## rhs_cavity_test() samples the cavity flow source terms at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'rhs_cavity_test():' ) print ( ' cavity flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' over the [0,+1]x[0,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_cavity ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_exppoly ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_exppoly() evaluates the source terms of the exppoly flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real US(N), VS(N), PS(N), the source terms. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exppoly ( nu, rho, n, x, y, t ) return us, vs, ps def rhs_exppoly_test ( rng ): #*****************************************************************************80 # ## rhs_exppoly_test() samples the exppoly source terms at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'rhs_exppoly_test():' ) print ( ' exppoly flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_exppoly ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_exptrig ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_exptrig() evaluates the source terms of the exptrig flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real US(N), VS(N), PS(N), the source terms. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exptrig ( nu, rho, n, x, y, t ) return us, vs, ps def rhs_exptrig_test ( rng ): #*****************************************************************************80 # ## rhs_exptrig_test() samples the exptrig source terms at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'rhs_exptrig_test():' ) print ( ' exptrig flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_exptrig ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_gms ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_gms() evaluates the source terms of the gms flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real US(N), VS(N), PS(N), the source terms. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_gms ( nu, rho, n, x, y, t ) return us, vs, ps def rhs_gms_test ( rng ): #*****************************************************************************80 # ## rhs_gms_test() samples the gms source terms at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'rhs_gms_test():' ) print ( ' gms flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = 1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_gms ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_lukas ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_lukas() returns right hand sides of the Lukas Bystricky equations. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real F(N), G(N), H(N), the right hand sides in the U, V and P equations. # import numpy as np f = np.zeros ( n ) g = np.zeros ( n ) h = np.zeros ( n ) u = - np.cos ( np.pi * x ) / np.pi dudt = np.zeros ( n ) dudx = np.sin ( np.pi * x ) dudxx = np.pi * np.cos ( np.pi * x ) dudy = np.zeros ( n ) dudyy = np.zeros ( n ) v = - y * np.sin ( np.pi * x ) dvdt = np.zeros ( n ) dvdx = - np.pi * y * np.cos ( np.pi * x ) dvdxx = + np.pi * np.pi * y * np.sin ( np.pi * x ) dvdy = - np.sin ( np.pi * x ) dvdyy = np.zeros ( n ) p = np.zeros ( n ) dpdx = np.zeros ( n ) dpdy = np.zeros ( n ) f = dudt - nu * ( dudxx + dudyy ) + u * dudx + v * dudy + dpdx / rho g = dvdt - nu * ( dvdxx + dvdyy ) + u * dvdx + v * dvdy + dpdy / rho h = dudx + dvdy return f, g, h def rhs_lukas_test ( rng ): #*****************************************************************************80 # ## rhs_lukas_test() samples Lukas right hand sides at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'rhs_lukas_test():' ) print ( ' Lukas Flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_lukas ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_noflow ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_noflow() evaluates the source terms of the NoFlow flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real US(N), VS(N), PS(N), the source terms. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_noflow ( nu, rho, n, x, y, t ) return us, vs, ps def rhs_noflow_test ( rng ): #*****************************************************************************80 # ## rhs_noflow_test() samples the NoFlow source terms at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'rhs_noflow_test():' ) print ( ' NoFlow Flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' over the [-1,+1]x[-1,+1] square' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = 1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_noflow ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_poiseuille ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_poiseuille() returns the Poiseuille right hand side. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real F(N), G(N), H(N), the right hand sides in the # U, V and P equations. # import numpy as np f = np.zeros ( n ) g = np.zeros ( n ) h = np.zeros ( n ) return f, g, h def rhs_poiseuille_test ( rng ): #*****************************************************************************80 # ## rhs_poiseuille_test() tests rhs_poiseuille(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'rhs_poiseuille_test():' ) print ( ' Poiseuille Flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' at the initial time T = 0, over a channel region.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 6.0 y_lo = -1.0 y_hi = +1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_poiseuille ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_spiral ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_spiral() returns right hand sides of the Spiral Flow equations. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Reference: # # Geoffrey Taylor, # On the decay of vortices in a viscous fluid, # Philosophical Magazine, # Volume 46, 1923, pages 671-674. # # Geoffrey Taylor, A E Green, # Mechanism for the production of small eddies from large ones, # Proceedings of the Royal Society of London, # Series A, Volume 158, 1937, pages 499-521. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real F(N), G(N), H(N), the right hand sides in the U, V and P equations. # import numpy as np f = np.zeros ( n ) g = np.zeros ( n ) h = np.zeros ( n ) u = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudt = nu * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudx = ( 1.0 + nu * t ) * 2.0 \ * ( 4.0 * x ** 3 - 6.0 * x ** 2 + 2.0 * x ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudxx = ( 1.0 + nu * t ) * 2.0 \ * ( 12.0 * x ** 2 - 12.0 * x + 2.0 ) \ * ( 2.0 * y ** 3 - 3.0 * y ** 2 + y ) dudy = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 6.0 * y ** 2 - 6.0 * y + 1.0 ) dudyy = ( 1.0 + nu * t ) * 2.0 \ * ( x ** 4 - 2.0 * x ** 3 + x ** 2 ) \ * ( 12.0 * y - 6.0 ) v = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdt = - nu * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdx = - ( 1.0 + nu * t ) * 2.0 \ * ( 6.0 * x ** 2 - 6.0 * x + 1.0 ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdxx = - ( 1.0 + nu * t ) * 2.0 \ * ( 12.0 * x - 6.0 ) \ * ( y ** 4 - 2.0 * y ** 3 + y ** 2 ) dvdy = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( 4.0 * y ** 3 - 6.0 * y ** 2 + 2.0 * y ) dvdyy = - ( 1.0 + nu * t ) * 2.0 \ * ( 2.0 * x ** 3 - 3.0 * x ** 2 + x ) \ * ( 12.0 * y ** 2 - 12.0 * y + 2.0 ) p = rho * y dpdx = 0.0 dpdy = rho f = dudt - nu * ( dudxx + dudyy ) \ + u * dudx + v * dudy + dpdx / rho g = dvdt - nu * ( dvdxx + dvdyy ) \ + u * dvdx + v * dvdy + dpdy / rho h = dudx + dvdy return f, g, h def rhs_spiral_test ( rng ): #*****************************************************************************80 # ## rhs_spiral_test() samples the right hand sides at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'rhs_spiral_test():' ) print ( ' Spiral flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = 1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_spiral ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_taylor ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_taylor() returns the Taylor right hand side. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Reference: # # Geoffrey Taylor, # On the decay of vortices in a viscous fluid, # Philosophical Magazine, # Volume 46, 1923, pages 671-674. # # Geoffrey Taylor, A E Green, # Mechanism for the production of small eddies from large ones, # Proceedings of the Royal Society of London, # Series A, Volume 158, 1937, pages 499-521. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real F(N), G(N), H(N), the residuals in the U, V and P equations. # import numpy as np f = np.zeros ( n ) g = np.zeros ( n ) h = np.zeros ( n ) return f, g, h def rhs_taylor_test ( rng ): #*****************************************************************************80 # ## rhs_taylor_test() samples the Taylor right hand side. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'rhs_taylor_test():' ) print ( ' Taylor flow' ) print ( ' Sample the Navier-Stokes right hand sides' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_taylor ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def rhs_vortex ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## rhs_vortex() returns the Vortex right hand side. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T(N), the time coordinate or coordinates. # # Output: # # real F(N), G(N), H(N), the residuals in the U, V and P equations. # import numpy as np f = - 2.0 * nu * ( np.pi ) ** 2 * ( np.cos ( np.pi * x ) * np.sin ( np.pi * y ) ) g = 2.0 * nu * ( np.pi ) ** 2 * ( np.sin ( np.pi * x ) * np.cos ( np.pi * y ) ) h = np.zeros ( n ) return f, g, h def rhs_vortex_test ( rng ): #*****************************************************************************80 # ## rhs_vortex_test() samples the Vortex right hand side. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'rhs_vortex_test():' ) print ( ' Sample the Vortex right hand sides' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 f, g, h = rhs_vortex ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' Ur: %14.6g %14.6g' % ( np.min ( f ), np.max ( f ) ) ) print ( ' Vr: %14.6g %14.6g' % ( np.min ( g ), np.max ( g ) ) ) print ( ' Pr: %14.6g %14.6g' % ( np.min ( h ), np.max ( h ) ) ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return None def uvp_cavity ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_cavity() returns velocity and pressure for the cavity flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Reference: # # Jean-Luc Guermand, Peter Minev, Jie Shen, # An overview of projection methods for incompressible flows, # Computer methods in applied mechanics and engineering, # Volume 105, pages 6011-6045, 2006. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), V(N), the X and Y velocity. # # real P(N), the pressure. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_cavity ( nu, rho, n, x, y, t ) return u, v, p def uvp_cavity_test ( rng ): #*****************************************************************************80 # ## uvp_cavity_test() samples the cavity flow at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_cavity_test():' ) print ( ' cavity flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the [0,+1]x[0,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = +1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) u, v, p = uvp_cavity ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_cavity_test2 ( ): #*****************************************************************************80 # ## uvp_cavity_test2() samples the cavity flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_cavity_test2():' ) print ( ' cavity flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the boundary of the [0,+1]x[0,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_cavity ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_exppoly ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_exppoly() returns velocity and pressure for the exppoly flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Reference: # # Jean-Luc Guermand, Peter Minev, Jie Shen, # An overview of projection methods for incompressible flows, # Computer methods in applied mechanics and engineering, # Volume 105, pages 6011-6045, 2006. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), V(N), the X and Y velocity. # # real P(N), the pressure. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exppoly ( nu, rho, n, x, y, t ) return u, v, p def uvp_exppoly_test ( rng ): #*****************************************************************************80 # ## uvp_exppoly_test() samples the exppoly flow at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_exppoly_test():' ) print ( ' exppoly flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = +1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) u, v, p = uvp_exppoly ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_exppoly_test2 ( ): #*****************************************************************************80 # ## uvp_exppoly_test2() samples the exppoly flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_exppoly_test2():' ) print ( ' exppoly flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the boundary of the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_exppoly ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_exptrig ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_exptrig() returns velocity and pressure for the exptrig flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), V(N), the X and Y velocity. # # real P(N), the pressure. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_exptrig ( nu, rho, n, x, y, t ) return u, v, p def uvp_exptrig_test ( rng ): #*****************************************************************************80 # ## uvp_exptrig_test() samples the exptrig flow at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_exptrig_test():' ) print ( ' exptrig flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = +1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) u, v, p = uvp_exptrig ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_exptrig_test2 ( ): #*****************************************************************************80 # ## uvp_exptrig_test2() samples the exptrig flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_exptrig_test2():' ) print ( ' exptrig flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the boundary of the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_exptrig ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_gms ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_gms() returns velocity and pressure for the gms flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Reference: # # Jean-Luc Guermand, Peter Minev, Jie Shen, # An overview of projection methods for incompressible flows, # Computer methods in applied mechanics and engineering, # Volume 105, pages 6011-6045, 2006. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), V(N), the X and Y velocity. # # real P(N), the pressure. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_gms ( nu, rho, n, x, y, t ) return u, v, p def uvp_gms_test ( rng ): #*****************************************************************************80 # ## uvp_gms_test() samples the gms flow at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_gms_test():' ) print ( ' gms flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = +1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) u, v, p = uvp_gms ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_gms_test2 ( ): #*****************************************************************************80 # ## uvp_gms_test2() samples the gms flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 August 2020 # # Author: # # John Burkardt # import numpy as np r8_lo = -1.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_gms_test2()' ) print ( ' gms flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the boundary of the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_gms ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_lukas ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_lukas() evaluates Lukas Bystricky's exact Navier Stokes solution. # # Discussion: # # There is no time dependence. # # The pressure is 0. # # The preferred domain is the unit square. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T or T(N), the time coordinate or coordinates. # # Output: # # real U(N), V(N), P(N), the velocity components and # pressure at each of the points. # import numpy as np u = np.zeros ( n ) v = np.zeros ( n ) p = np.zeros ( n ) u = - np.cos ( np.pi * x ) / np.pi v = - y * np.sin ( np.pi * x ) return u, v, p def uvp_lukas_test ( rng ): #*****************************************************************************80 # ## uvp_lukas_test() samples the Lukas Bystricky solution at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'uvp_lukas_test():' ) print ( ' Lukas Bystricky flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = +1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 u, v, p = uvp_lukas ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_lukas_test2 ( ): #*****************************************************************************80 # ## uvp_lukas_test2() samples the Lukas Bystricky flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 0.0 print ( '' ) print ( 'uvp_lukas_test2():' ) print ( ' Lukas Bystricky flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' on the boundary' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_lukas ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_noflow ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_noflow() returns velocity and pressure for the NoFlow flow. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the fluid density. # # integer N, the number of nodes. # # real X(N), Y(N), the coordinates of nodes. # # real T, the current time. # # Output: # # real U(N), V(N), the X and Y velocity. # # real P(N), the pressure. # u, dudt, dudx, dudxx, dudy, dudyy, us, \ v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, \ p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps = all_noflow ( nu, rho, n, x, y, t ) return u, v, p def uvp_noflow_test ( rng ): #*****************************************************************************80 # ## uvp_noflow_test() samples the NoFlow flow at a specific time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_noflow_test():' ) print ( ' NoFlow flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = -1.0 x_hi = +1.0 y_lo = -1.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) u, v, p = uvp_noflow ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_noflow_test2 ( ): #*****************************************************************************80 # ## uvp_noflow_test2() samples the NoFlow flow on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 September 2025 # # Author: # # John Burkardt # import numpy as np r8_lo = -1.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'uvp_noflow_test2():' ) print ( ' NoFlow flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' over the boundary of the [-1,+1]x[-1,+1] square,' ) print ( ' at time T = 1.0' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_noflow ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_poiseuille ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_poiseuille() evaluates the Poiseuille solution. # # Discussion: # # This flow is known as a Poiseuille Flow solution. # # The given velocity and pressure fields are exact solutions for the 2D # incompressible time-dependent Navier Stokes equations over the unit square. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T or T(N), the time coordinate or coordinates. # # Output: # # real U(N), V(N), P(N), the velocity components and # pressure at each of the points. # import numpy as np u = np.zeros ( n ) v = np.zeros ( n ) p = np.zeros ( n ) u = 1.0 - y ** 2 # # Can't write it this way or V becomes a scalar! # # v = 0.0; p = -2.0 * rho * nu * x return u, v, p def uvp_poiseuille_test ( rng ): #*****************************************************************************80 # ## uvp_poiseuille_test() samples the Poiseuille solution. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'uvp_poiseuille_test():' ) print ( ' Poiseuille flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' at the initial time T = 0, over a channel region.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = +0.0 x_hi = +6.0 y_lo = -1.0 y_hi = + 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 u, v, p = uvp_poiseuille ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_poiseuille_test2 ( ): #*****************************************************************************80 # ## uvp_poiseuille_test2() samples the Poiseuille solution on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # import numpy as np x_lo = +0.0 x_hi = +6.0 y_lo = -1.0 y_hi = + 1.0 nu = 1.0 rho = 1.0 t = 0.0 print ( '' ) print ( 'uvp_poiseuille_test2():' ) print ( ' Poiseuille flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' on the boundary' ) print ( ' at the initial time T = 0, over a channel region.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( x_lo, x_hi, 100 ) y[0:100] = y_lo x[100:200] = x_hi y[100:200] = np.linspace ( y_lo, y_hi, 100 ) x[200:300] = np.linspace ( x_hi, x_lo, 100 ) y[200:300] = y_hi x[300:400] = x_lo y[300:400] = np.linspace ( y_hi, y_lo, 100 ) u, v, p = uvp_poiseuille ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_spiral ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_spiral() evaluates the Spiral Flow exact Navier Stokes solution. # # Discussion: # # This flow is known as a Spiral Flow solution. # # The given velocity and pressure fields are exact solutions for the 2D # incompressible time-dependent Navier Stokes equations over the unit square. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 January 2015 # # Author: # # John Burkardt # # Reference: # # Maxim Olshanskii, Leo Rebholz, # Application of barycenter refined meshes in linear elasticity # and incompressible fluid dynamics, # ETNA: Electronic Transactions in Numerical Analysis, # Volume 38, pages 258-274, 2011. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T or T(N), the time coordinate or coordinates. # # Output: # # real U(N), V(N), P(N), the velocity components and # pressure at each of the points. # import numpy as np u = np.zeros ( n ) v = np.zeros ( n ) p = np.zeros ( n ) u = ( 1.0 + nu * t ) * 2.0 \ * x ** 2 * ( x - 1.0 ) ** 2 \ * y * ( 2.0 * y - 1.0 ) * ( y - 1.0 ) v = - ( 1.0 + nu * t ) * 2.0 \ * x * ( 2.0 * x - 1.0 ) * ( x - 1.0 ) \ * y ** 2 * ( y - 1.0 ) ** 2 p = rho * y return u, v, p def uvp_spiral_test ( rng ): #*****************************************************************************80 # ## uvp_spiral_test() samples the Spiral Flow solution at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'uvp_spiral_test():' ) print ( ' Spiral flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.0 x_hi = +1.0 y_lo = 0.0 y_hi = 1.0 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 u, v, p = uvp_spiral ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_spiral_test2 ( ): #*****************************************************************************80 # ## uvp_spiral_test2() samples the Spiral Flow on the boundary at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.0 r8_hi = +1.0 nu = 1.0 rho = 1.0 t = 0.0 print ( '' ) print ( 'uvp_spiral_test2():' ) print ( ' Spiral flow' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' on the boundary' ) print ( ' at the initial time T = 0, over the unit square.' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_spiral ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_taylor ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_taylor() evaluates the Taylor exact Navier Stokes solution. # # Discussion: # # This flow is known as a Taylor-Green vortex. # # The given velocity and pressure fields are exact solutions for the 2D # 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. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 17 January 2015 # # Author: # # John Burkardt # # Reference: # # Geoffrey Taylor, # On the decay of vortices in a viscous fluid, # Philosophical Magazine, # Volume 46, 1923, pages 671-674. # # Geoffrey Taylor, A E Green, # Mechanism for the production of small eddies from large ones, # Proceedings of the Royal Society of London, # Series A, Volume 158, 1937, pages 499-521. # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T or T(N), the time coordinate or coordinates. # # Output: # # real U(N), V(N), P(N), the velocity components and # pressure at each of the points. # import numpy as np cx = np.cos ( np.pi * x ) cy = np.cos ( np.pi * y ) c2x = np.cos ( 2.0 * np.pi * x ) c2y = np.cos ( 2.0 * np.pi * y ) sx = np.sin ( np.pi * x ) sy = np.sin ( np.pi * y ) e2t = np.exp ( - 2.0 * np.pi * np.pi * nu * t ) e4t = np.exp ( - 4.0 * np.pi * np.pi * nu * t ) u = np.zeros ( n ) v = np.zeros ( n ) p = np.zeros ( n ) u = - cx * sy * e2t v = sx * cy * e2t p = - 0.25 * rho * ( c2x + c2y ) * e4t return u, v, p def uvp_taylor_test ( rng ): #*****************************************************************************80 # ## uvp_taylor_test() samples the solution at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 January 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'uvp_taylor_test():' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 u, v, p = uvp_taylor ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_taylor_test2 ( ): #*****************************************************************************80 # ## uvp_taylor_test2() samples the solution on the boundary at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.5 r8_hi = +2.5 nu = 1.0 rho = 1.0 t = 0.0 print ( '' ) print ( 'uvp_taylor_test2():' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' on the boundary' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) # # Python is consistent in its willful flouting of sensible conventions. # X[0:100] means X from 0 to 99...! # x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_taylor ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_vortex ( nu, rho, n, x, y, t ): #*****************************************************************************80 # ## uvp_vortex() evaluates the Vortex solution. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # real NU, the kinematic viscosity. # # real RHO, the density. # # integer N, the number of points at which the solution is to # be evaluated. # # real X(N), Y(N), the coordinates of the points. # # real T or T(N), the time coordinate or coordinates. # # Output: # # real U(N), V(N), P(N), the velocity components and # pressure at each of the points. # import numpy as np cx = np.cos ( np.pi * x ) cy = np.cos ( np.pi * y ) c2x = np.cos ( 2.0 * np.pi * x ) c2y = np.cos ( 2.0 * np.pi * y ) sx = np.sin ( np.pi * x ) sy = np.sin ( np.pi * y ) u = np.zeros ( n ) v = np.zeros ( n ) p = np.zeros ( n ) u = - cx * sy v = sx * cy p = - 0.25 * rho * ( c2x + c2y ) return u, v, p def uvp_vortex_test ( rng ): #*****************************************************************************80 # ## uvp_vortex_test() samples the Vortex solution at the initial time. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 July 2015 # # Author: # # John Burkardt # # Input: # # rng(): the current random number generator. # import numpy as np nu = 1.0 rho = 1.0 print ( '' ) print ( 'uvp_vortex_test():' ) print ( ' Sample the vortex flow.' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 1000 x_lo = 0.5 x_hi = +2.5 y_lo = 0.5 y_hi = 2.5 x = x_lo + ( x_hi - x_lo ) * rng.random ( size = n ) y = y_lo + ( y_hi - y_lo ) * rng.random ( size = n ) t = 0.0 u, v, p = uvp_vortex ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def uvp_vortex_test2 ( ): #*****************************************************************************80 # ## uvp_vortex_test2() samples the Vortex solution on the boundary. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # import numpy as np r8_lo = 0.5 r8_hi = +2.5 nu = 1.0 rho = 1.0 t = 0.0 print ( '' ) print ( 'uvp_vortex_test2()' ) print ( ' Sample the Vortex flow.' ) print ( ' Estimate the range of velocity and pressure' ) print ( ' on the boundary' ) print ( ' at the initial time T = 0, using a region that is' ) print ( ' the square centered at (1.5,1.5) with "radius" 1.0,' ) print ( ' Kinematic viscosity NU = %g' % ( nu ) ) print ( ' Fluid density RHO = %g' % ( rho ) ) n = 400 x = np.zeros ( n ) y = np.zeros ( n ) # # Python is consistent in its willful flouting of sensible conventions. # X[0:100] means X from 0 to 99...! # x[0:100] = np.linspace ( r8_lo, r8_hi, 100 ) y[0:100] = r8_lo x[100:200] = r8_hi y[100:200] = np.linspace ( r8_lo, r8_hi, 100 ) x[200:300] = np.linspace ( r8_hi, r8_lo, 100 ) y[200:300] = r8_hi x[300:400] = r8_lo y[300:400] = np.linspace ( r8_lo, r8_hi, 100 ) u, v, p = uvp_vortex ( nu, rho, n, x, y, t ) print ( '' ) print ( ' Minimum Maximum' ) print ( '' ) print ( ' U: %14.6g %14.6g' % ( np.min ( u ), np.max ( u ) ) ) print ( ' V: %14.6g %14.6g' % ( np.min ( v ), np.max ( v ) ) ) print ( ' P: %14.6g %14.6g' % ( np.min ( p ), np.max ( p ) ) ) return def velocity_matplotlib ( header, n, x, y, u, v, p, s ): #*****************************************************************************80 # ## velocity_matplotlib() plots a velocity vector field with matplotlib. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # # Input: # # string HEADER, a header to be used to name the files. # # integer N, the number of evaluation points. # # real X(N), Y(N), the coordinates of the evaluation points. # # real U(N), V(N), P(N), the solution samples. # # real S, a scale factor for the velocity vectors. # import matplotlib.pyplot as plt myplot = plt.figure ( ) ax = plt.gca ( ) ax.quiver ( x, y, u, v ) ax.set_xlabel ( '<--X-->' ) ax.set_ylabel ( '<--Y-->' ) ax.set_title ( header ) ax.axis ( 'equal' ) plt.draw ( ) plot_filename = 'velocity_' + header + '_matplotlib.png' myplot.savefig ( plot_filename ) plt.show ( block = False ) plt.close ( ) return def velocity_cavity_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_cavity_matplotlib_test() plots a cavity velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 September 2025 # # Author: # # John Burkardt # nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'velocity_cavity_matplotlib_test():' ) print ( ' Generate a cavity flow on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 1.0 x_num = 21 y_lo = 0.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) n = x_num * y_num u, v, p = uvp_cavity ( nu, rho, n, x, y, t ) header = 'cavity' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_exppoly_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_exppoly_matplotlib_test() plots the exppoly velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'velocity_exppoly_matplotlib_test():' ) print ( ' Generate the exppoly flow on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 1.0 x_num = 21 y_lo = 0.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) n = x_num * y_num u, v, p = uvp_exppoly ( nu, rho, n, x, y, t ) header = 'exppoly' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_exptrig_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_exptrig_matplotlib_test() plots the exptrig velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 September 2025 # # Author: # # John Burkardt # nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'velocity_exptrig_matplotlib_test():' ) print ( ' Generate the exptrig flow on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 1.0 x_num = 21 y_lo = 0.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) n = x_num * y_num u, v, p = uvp_exptrig ( nu, rho, n, x, y, t ) header = 'exptrig' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_gms_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_gms_matplotlib_test() plots the gms velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2025 # # Author: # # John Burkardt # nu = 1.0 rho = 1.0 t = 1.0 print ( '' ) print ( 'velocity_gms_matplotlib_test():' ) print ( ' Generate the gms flow on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = -1.0 x_hi = 1.0 x_num = 21 y_lo = -1.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) n = x_num * y_num u, v, p = uvp_gms ( nu, rho, n, x, y, t ) header = 'gms' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_lukas_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_lukas_matplotlib_test() plots a Lukas Bystricky velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 March 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'velocity_lukas_matplotlib_test():' ) print ( ' Generate a Lukas Bystricky Flow field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 1.0 x_num = 21 y_lo = 0.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) nu = 1.0 rho = 1.0 n = x_num * y_num t = 0.0 u, v, p = uvp_lukas ( nu, rho, n, x, y, t ) header = 'lukas' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_poiseuille_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_poiseuille_matplotlib_test() plots a Poiseuille velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'velocity_poiseuille_matplotlib_test():' ) print ( ' Generate a Poiseuille velocity field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 6.0 x_num = 61 y_lo = -1.0 y_hi = +1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) nu = 1.0 rho = 1.0 n = x_num * y_num t = 0.0 u, v, p = uvp_poiseuille ( nu, rho, n, x, y, t ) header = 'poiseuille' s = 5.0 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_spiral_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_spiral_matplotlib_test() plots a Spiral velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'velocity_spiral_matplotlib_test():' ) print ( ' Generate a Spiral Flow velocity field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.0 x_hi = 1.0 x_num = 21 y_lo = 0.0 y_hi = 1.0 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) nu = 1.0 rho = 1.0 n = x_num * y_num t = 0.0 u, v, p = uvp_spiral ( nu, rho, n, x, y, t ) header = 'spiral' s = 5.0 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_taylor_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_taylor_matplotlib_test() plots a Taylor Flow field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'velocity_taylor_matplotlib_test():' ) print ( ' Generate a Taylor velocity field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.5 x_hi = 2.5 x_num = 21 y_lo = 0.5 y_hi = 2.5 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) nu = 1.0 rho = 1.0 n = x_num * y_num t = 0.0 u, v, p = uvp_taylor ( nu, rho, n, x, y, t ) header = 'taylor' s = 0.10 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return def velocity_vortex_matplotlib_test ( ): #*****************************************************************************80 # ## velocity_vortex_matplotlib_test() plots a Vortex velocity field. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 July 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'velocity_vortex_matplotlib_test():' ) print ( ' Generate a Vortex velocity field on a regular grid.' ) print ( ' Display it using matplotlib' ) x_lo = 0.5 x_hi = 1.5 x_num = 21 y_lo = 0.5 y_hi = 1.5 y_num = 21 x, y = grid_2d ( x_num, x_lo, x_hi, y_num, y_lo, y_hi ) nu = 1.0 rho = 1.0 n = x_num * y_num t = 0.0 u, v, p = uvp_vortex ( nu, rho, n, x, y, t ) header = 'vortex' s = 0.25 velocity_matplotlib ( header, n, x, y, u, v, p, s ) return if ( __name__ == '__main__' ): timestamp ( ) navier_stokes_2d_exact_test ( ) timestamp ( )