[讨论]

纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之a

2022-11-25 09:50:03 370 控制程序

0 !********************************************************** !******************************************************** !纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之a !******************************************************** ! SINOMACH ! ! 纳维斯托克斯两维精准估算解决非压缩时间依赖方程于一任意2D领域之a ! ! 周勇 20221125 !********************************************************** !********************************************************** !c $$$ Run, Output: !********************************************************** !********************************************************** Program main !******************************************************************72 ! !c NS2DE_TEST tests the NS2DE library. ! Implicit None Call timestamp() Write (, '(a)') '' Write (, '(a)') 'NS2DE_PRB' Write (, '(a)') ' FORTRAN90,95,&2018 version' Write (, '(a)') ' Test the NS2DE library.' ! ! Taylor Vortex Flow. ! Call uvp_taylor_test() Call uvp_taylor_test2() Call rhs_taylor_test() Call resid_taylor_test() Call gnuplot_taylor_test() Call parameter_taylor_test() ! ! Spiral Flow. ! Call uvp_spiral_test() Call uvp_spiral_test2() Call rhs_spiral_test() Call resid_spiral_test() Call gnuplot_spiral_test() Call parameter_spiral_test() ! ! Lucas Flow. ! Call uvp_lucas_test() Call uvp_lucas_test2() Call rhs_lucas_test() Call resid_lucas_test() Call gnuplot_lucas_test() ! ! Terminate. ! Write (, '(a)') '' Write (, '(a)') 'NS2DE_TEST' Write (, '(a)') ' Normal end of execution.' Call timestamp() Return End Program main Subroutine uvp_taylor_test() !*******************************************************************72 ! !c UVP_TAYLOR_TEST samples the solution at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_TAYLOR_TEST' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the square centered at (1.5,1.5) with "radius" 1.0,' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.5D+00 r8_hi = +2.5D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call uvp_taylor(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_taylor_test Subroutine uvp_taylor_test2() !******************************************************************72 ! !c UVP_TAYLOR_TEST2 samples the solution on the boundary at the initial time. ! Implicit None Integer n Parameter (n=400) Integer i Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 t = 0.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_TAYLOR_TEST2' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' along the boundary' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the square centered at (1.5,1.5) with "radius" 1.0,' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.5D+00 r8_hi = 2.5D+00 Call r8vec_linspace(100, r8_lo, r8_hi, x(1:100)) Do i = 1, 100 y(i) = r8_lo End Do Do i = 101, 200 x(i) = r8_hi End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(101:200)) Call r8vec_linspace(100, r8_hi, r8_lo, x(201:300)) Do i = 201, 300 y(i) = r8_hi End Do Do i = 301, 400 x(i) = r8_lo End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(301:400)) Call uvp_taylor(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_taylor_test2 Subroutine rhs_taylor_test() !*******************************************************************72 ! !c RHS_TAYLOR_TEST samples the right hand side at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision f(n) Double Precision g(n) Double Precision h(n) Double Precision nu Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RHS_TAYLOR_TEST' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Sample the Navier-Stokes right hand sides' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the square centered at (1.5,1.5) with "radius" 1.0,' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.5D+00 r8_hi = +2.5D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call rhs_taylor(nu, rho, n, x, y, t, f, g, h) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' F: ', r8vec_min(n, f), r8vec_max(n, f) Write (, '(a,2x,g14.6,2x,g14.6)') ' G: ', r8vec_min(n, g), r8vec_max(n, g) Write (, '(a,2x,g14.6,2x,g14.6)') ' H: ', r8vec_min(n, h), r8vec_max(n, h) Return End Subroutine rhs_taylor_test Subroutine resid_taylor_test() !*****************************************************************72 ! !c RESID_TAYLOR_TEST samples the residual at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision pr(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_amax Double Precision r8vec_amin Double Precision rho Integer seed Double Precision t Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RESID_TAYLOR_TEST' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Sample the Navier-Stokes residuals' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the square centered at (1.5,1.5) with "radius" 1.0,' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.5D+00 r8_hi = +2.5D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call resid_taylor(nu, rho, n, x, y, t, ur, vr, pr) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' Ur: ', r8vec_amin(n, ur),& r8vec_amax(n, ur) Write (, '(a,2x,g14.6,2x,g14.6)') ' Vr: ', r8vec_amin(n, vr),& r8vec_amax(n, vr) Write (, '(a,2x,g14.6,2x,g14.6)') ' Pr: ', r8vec_amin(n, pr),& r8vec_amax(n, pr) Return End Subroutine resid_taylor_test Subroutine gnuplot_taylor_test() !******************************************************************72 ! !c GNUPLOT_TAYLOR_TEST generates a field on a regular grid and plots it. ! Implicit None Integer x_num Parameter (x_num=21) Integer y_num Parameter (y_num=21) Character (255) header Integer n Double Precision nu Double Precision p(x_num, y_num) Double Precision rho Double Precision s Integer seed Double Precision t Double Precision u(x_num, y_num) Double Precision v(x_num, y_num) Double Precision x(x_num, y_num) Double Precision x_hi Double Precision x_lo Double Precision y(x_num, y_num) Double Precision y_hi Double Precision y_lo Write (, '(a)') '' Write (, '(a)') 'GNUPLOT_TAYLOR_TEST:' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Generate a Taylor vortex velocity field on a regular grid.' Write (, '(a)') ' Store in GNUPLOT data and command files.' x_lo = 0.5D+00 x_hi = 2.5D+00 y_lo = 0.5D+00 y_hi = 2.5D+00 Call grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y) nu = 1.0D+00 rho = 1.0D+00 n = x_numy_num t = 0.0D+00 Call uvp_taylor(nu, rho, n, x, y, t, u, v, p) header = 'taylor' s = 0.10D+00 Call ns2de_gnuplot(header, n, x, y, u, v, s) Return End Subroutine gnuplot_taylor_test Subroutine parameter_taylor_test() !********************************************************************72 ! !c PARAMETER_TAYLOR_TEST monitors solution norms for various values of NU, RHO. ! Implicit None Integer n Parameter (n=1000) Integer i Integer j Integer k Double Precision nu Double Precision p(n) Double Precision p_norm Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_norm_l2 Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision u_norm Double Precision v(n) Double Precision v_norm Double Precision x(n) Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'PARAMETER_TAYLOR_TEST' Write (, '(a)') ' Taylor Vortex Flow' Write (, '(a)') ' Monitor solution norms over time for' Write (, '(a)') ' varous values of NU, RHO.' r8_lo = 0.5D+00 r8_hi = +2.5D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) ! ! Vary RHO. ! Write (, '(a)') '' Write (, '(a)') ' RHO affects the pressure scaling.' Write (, '(a)') '' Write (, '(a)') ' RHO NU T' //& ' \|\|U\|\| \|\|V\|\| \|\|P\|\|' Write (, '(a)') '' nu = 1.0D+00 rho = 1.0D+00 Do j = 1, 3 `Do k = 0, 5 t = real(k, kind=8)/5.0D+00 Call uvp_taylor(nu, rho, n, x, y, t, u, v, p) u_norm = r8vec_norm_l2(n, u)/dble(n) v_norm = r8vec_norm_l2(n, v)/dble(n) p_norm = r8vec_norm_l2(n, p)/dble(n) Write (, '(2x,g10.4,2x,g10.4,2x,f8.4,2x,g10.4,'//'2x,g10.4,2x,g10.4)')& rho, nu, t, u_norm, v_norm, p_norm End Do Write (, '(a)') '' rho = rho/100.0D+00` End Do ! ! Vary NU. ! Write (, '(a)') '' Write (, '(a)') ' NU affects the time scaling.' Write (, '(a)') '' Write (, '(a)') ' RHO NU T' //& ' \|\|U\|\| \|\|V\|\| \|\|P\|\|' Write (, '(a)') '' nu = 1.0D+00 rho = 1.0D+00 Do i = 1, 4 `Do k = 0, 5 t = real(k, kind=8)/5.0D+00 Call uvp_taylor(nu, rho, n, x, y, t, u, v, p) u_norm = r8vec_norm_l2(n, u)/dble(n) v_norm = r8vec_norm_l2(n, v)/dble(n) p_norm = r8vec_norm_l2(n, p)/dble(n) Write (, '(2x,g10.4,2x,g10.4,2x,f8.4,2x,g10.4,'//'2x,g10.4,2x,g10.4)')& rho, nu, t, u_norm, v_norm, p_norm End Do Write (, '(a)') '' nu = nu/10.0D+00` End Do Return End Subroutine parameter_taylor_test Subroutine uvp_spiral_test() !*******************************************************************72 ! !c UVP_SPIRAL_TEST samples the Spiral solution at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_SPIRAL_TEST' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call uvp_spiral(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_spiral_test Subroutine uvp_spiral_test2() !******************************************************************72 ! !c UVP_SPIRAL_TEST2 samples the Spiral solution on the boundary at the initial time. ! Implicit None Integer n Parameter (n=400) Integer i Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 t = 0.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_SPIRAL_TEST2' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' along the boundary' Write (, '(a)') ' at the initial time T = 0, using a region that is' Write (, '(a)') ' the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 Call r8vec_linspace(100, r8_lo, r8_hi, x(1:100)) Do i = 1, 100 y(i) = r8_lo End Do Do i = 101, 200 x(i) = r8_hi End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(101:200)) Call r8vec_linspace(100, r8_hi, r8_lo, x(201:300)) Do i = 201, 300 y(i) = r8_hi End Do Do i = 301, 400 x(i) = r8_lo End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(301:400)) Call uvp_spiral(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_spiral_test2 Subroutine rhs_spiral_test() !********************************************************************72 ! !c RHS_SPIRAL_TEST samples the right hand side at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision f(n) Double Precision g(n) Double Precision h(n) Double Precision nu Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RHS_SPIRAL_TEST' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Sample the Navier-Stokes right hand sides' Write (, '(a)') ' at the initial time T = 0, in the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call rhs_spiral(nu, rho, n, x, y, t, f, g, h) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' F: ', r8vec_min(n, f), r8vec_max(n, f) Write (, '(a,2x,g14.6,2x,g14.6)') ' G: ', r8vec_min(n, g), r8vec_max(n, g) Write (, '(a,2x,g14.6,2x,g14.6)') ' H: ', r8vec_min(n, h), r8vec_max(n, h) Return End Subroutine rhs_spiral_test Subroutine resid_spiral_test() !********************************************************************72 ! !c RESID_SPIRAL_TEST samples the residual at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision pr(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_amax Double Precision r8vec_amin Double Precision rho Integer seed Double Precision t Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RESID_SPIRAL_TEST' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Sample the Navier-Stokes residuals' Write (, '(a)') ' at the initial time T = 0, in the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call resid_spiral(nu, rho, n, x, y, t, ur, vr, pr) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' Ur: ', r8vec_amin(n, ur),& r8vec_amax(n, ur) Write (, '(a,2x,g14.6,2x,g14.6)') ' Vr: ', r8vec_amin(n, vr),& r8vec_amax(n, vr) Write (, '(a,2x,g14.6,2x,g14.6)') ' Pr: ', r8vec_amin(n, pr),& r8vec_amax(n, pr) Return End Subroutine resid_spiral_test Subroutine gnuplot_spiral_test() !********************************************************************72 ! !c GNUPLOT_SPIRAL_TEST generates a field on a regular grid and plots it. ! Implicit None Integer x_num Parameter (x_num=21) Integer y_num Parameter (y_num=21) Character (255) header Integer n Double Precision nu Double Precision p(x_num, y_num) Double Precision rho Double Precision s Integer seed Double Precision t Double Precision u(x_num, y_num) Double Precision v(x_num, y_num) Double Precision x(x_num, y_num) Double Precision x_hi Double Precision x_lo Double Precision y(x_num, y_num) Double Precision y_hi Double Precision y_lo Write (, '(a)') '' Write (, '(a)') 'GNUPLOT_SPIRAL_TEST:' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Generate a velocity field on a regular grid.' Write (, '(a)') ' Store in GNUPLOT data and command files.' x_lo = 0.0D+00 x_hi = 1.0D+00 y_lo = 0.0D+00 y_hi = 1.0D+00 Call grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y) nu = 1.0D+00 rho = 1.0D+00 n = x_numy_num t = 0.0D+00 Call uvp_spiral(nu, rho, n, x, y, t, u, v, p) header = 'spiral' s = 5.00D+00 Call ns2de_gnuplot(header, n, x, y, u, v, s) Return End Subroutine gnuplot_spiral_test Subroutine parameter_spiral_test() !********************************************************************72 ! !c PARAMETER_SPIRAL_TEST monitors solution norms for various values of NU, RHO. ! Implicit None Integer n Parameter (n=1000) Integer i Integer j Integer k Double Precision nu Double Precision p(n) Double Precision p_norm Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_norm_l2 Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision u_norm Double Precision v(n) Double Precision v_norm Double Precision x(n) Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'PARAMETER_SPIRAL_TEST' Write (, '(a)') ' Spiral Flow' Write (, '(a)') ' Monitor solution norms over time for' Write (, '(a)') ' varous values of NU, RHO.' r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) ! ! Vary RHO. ! Write (, '(a)') '' Write (, '(a)') ' RHO affects the pressure scaling.' Write (, '(a)') '' Write (, '(a)') ' RHO NU T' // & ' \|\|U\|\| \|\|V\|\| \|\|P\|\|' Write (, '(a)') '' nu = 1.0D+00 rho = 1.0D+00 Do j = 1, 3 `Do k = 0, 5 t = real(k, kind=8)/5.0D+00 Call uvp_spiral(nu, rho, n, x, y, t, u, v, p) u_norm = r8vec_norm_l2(n, u)/dble(n) v_norm = r8vec_norm_l2(n, v)/dble(n) p_norm = r8vec_norm_l2(n, p)/dble(n) Write (, '(2x,g10.4,2x,g10.4,2x,f8.4,2x,g10.4,'//'2x,g10.4,2x,g10.4)')& rho, nu, t, u_norm, v_norm, p_norm End Do Write (, '(a)') '' rho = rho/100.0D+00` End Do ! ! Vary NU. ! Write (, '(a)') '' Write (, '(a)') ' NU affects the time scaling.' Write (, '(a)') '' Write (, '(a)') ' RHO NU T' //& ' \|\|U\|\| \|\|V\|\| \|\|P\|\|' Write (, '(a)') '' nu = 1.0D+00 rho = 1.0D+00 Do i = 1, 4 `Do k = 0, 5 t = real(k, kind=8)/5.0D+00 Call uvp_spiral(nu, rho, n, x, y, t, u, v, p) u_norm = r8vec_norm_l2(n, u)/dble(n) v_norm = r8vec_norm_l2(n, v)/dble(n) p_norm = r8vec_norm_l2(n, p)/dble(n) Write (, '(2x,g10.4,2x,g10.4,2x,f8.4,2x,g10.4,'//'2x,g10.4,2x,g10.4)')& rho, nu, t, u_norm, v_norm, p_norm End Do Write (, '(a)') '' nu = nu/10.0D+00` End Do Return End Subroutine parameter_spiral_test Subroutine uvp_lucas_test() !********************************************************************72 ! !c UVP_LUCAS_TEST samples the Lucas Bystricky solution at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_LUCAS_TEST' Write (, '(a)') ' Lucas Bystricky Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' at the initial time T = 0, over the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = +1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call uvp_lucas(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_lucas_test Subroutine uvp_lucas_test2() !*******************************************************************72 ! !c UVP_LUCAS_TEST2 samples the Lucas Bystricky solution on the boundary & !at the initial time. ! Implicit None Integer n Parameter (n=400) Integer i Double Precision nu Double Precision p(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D+00 t = 0.0D+00 Write (, '(a)') '' Write (, '(a)') 'UVP_LUCAS_TEST2' Write (, '(a)') ' Lucas Bystricky Flow' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' along the boundary' Write (, '(a)') ' at the initial time T = 0, over the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 Call r8vec_linspace(100, r8_lo, r8_hi, x(1:100)) Do i = 1, 100 y(i) = r8_lo End Do Do i = 101, 200 x(i) = r8_hi End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(101:200)) Call r8vec_linspace(100, r8_hi, r8_lo, x(201:300)) Do i = 201, 300 y(i) = r8_hi End Do Do i = 301, 400 x(i) = r8_lo End Do Call r8vec_linspace(100, r8_lo, r8_hi, y(301:400)) Call uvp_lucas(nu, rho, n, x, y, t, u, v, p) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' U: ', r8vec_min(n, u), r8vec_max(n, u) Write (, '(a,2x,g14.6,2x,g14.6)') ' V: ', r8vec_min(n, v), r8vec_max(n, v) Write (, '(a,2x,g14.6,2x,g14.6)') ' P: ', r8vec_min(n, p), r8vec_max(n, p) Return End Subroutine uvp_lucas_test2 Subroutine rhs_lucas_test() !******************************************************************72 ! !c RHS_LUCAS_TEST samples the Lucas Bystricky right hand side. ! Implicit None Integer n Parameter (n=1000) Double Precision f(n) Double Precision g(n) Double Precision h(n) Double Precision nu Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_max Double Precision r8vec_min Double Precision rho Integer seed Double Precision t Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RHS_LUCAS_TEST' Write (, '(a)') ' Lucas Bystricky Flow' Write (, '(a)') ' Sample the Navier-Stokes right hand sides' Write (, '(a)') ' at the initial time T = 0, in the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call rhs_lucas(nu, rho, n, x, y, t, f, g, h) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' F: ', r8vec_min(n, f), r8vec_max(n, f) Write (, '(a,2x,g14.6,2x,g14.6)') ' G: ', r8vec_min(n, g), r8vec_max(n, g) Write (, '(a,2x,g14.6,2x,g14.6)') ' H: ', r8vec_min(n, h), r8vec_max(n, h) Return End Subroutine rhs_lucas_test Subroutine resid_lucas_test() !********************************************************************72 ! !c RESID_LUCAS_TEST samples the Lucas Bystricky residual at the initial time. ! Implicit None Integer n Parameter (n=1000) Double Precision nu Double Precision pr(n) Double Precision r8_hi Double Precision r8_lo Double Precision r8vec_amax Double Precision r8vec_amin Double Precision rho Integer seed Double Precision t Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision y(n) nu = 1.0D+00 rho = 1.0D00 Write (, '(a)') '' Write (, '(a)') 'RESID_LUCAS_TEST' Write (, '(a)') ' Lucas Bystricky Flow' Write (, '(a)') ' Sample the Navier-Stokes residuals' Write (, '(a)') ' at the initial time T = 0, in the unit square.' Write (, '(a,g14.6)') ' Kinematic viscosity NU = ', nu Write (, '(a,g14.6)') ' Fluid density RHO = ', rho r8_lo = 0.0D+00 r8_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, x) Call r8vec_uniform_ab(n, r8_lo, r8_hi, seed, y) t = 0.0D+00 Call resid_lucas(nu, rho, n, x, y, t, ur, vr, pr) Write (, '(a)') '' Write (, '(a)') ' Minimum Maximum' Write (, '(a)') '' Write (, '(a,2x,g14.6,2x,g14.6)') ' Ur: ', r8vec_amin(n, ur),& r8vec_amax(n, ur) Write (, '(a,2x,g14.6,2x,g14.6)') ' Vr: ', r8vec_amin(n, vr),& r8vec_amax(n, vr) Write (, '(a,2x,g14.6,2x,g14.6)') ' Pr: ', r8vec_amin(n, pr),& r8vec_amax(n, pr) Return End Subroutine resid_lucas_test 0
2022-11-25 09:50:03　　评论淘帖0 举报相关推荐 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之c 286 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之图形表达 572 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之b 346 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之a 734 • 纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域 463 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b 714 • 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之c 666 • 精准估算压缩性稳态斯托克斯方程于2D单位正方形 217 • 物理学最难的方程之一纳维-斯托克斯方程，解答奖金高达100万美金 27911 • 精准估算压缩性稳态斯托克斯方程于2D单位正方形程序图形示意 262