!************************************************************
!************************************************************
!纳维斯托克斯两维精准估算解决非压缩时间依赖方程(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_num*y_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_num*y_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
更多回帖
