[讨论]

精准估算压缩性稳态斯托克斯方程于2D单位正方形

2022-11-29 16:58:36 234 控制

0 !************************************************************ !********************************************************** ! 精准估算压缩性稳态斯托克斯方程于2D单位正方形 !********************************************************** ! SINOMACH ! ! 斯托克斯两维精准估算解答 ! 压缩性稳态斯托克斯方程于2D单位正方形 ! 斯托克斯两维精准估算解答费压缩性稳态斯托克斯方程于2D单位正方形 ! ! 周勇 ! 20221129 !********************************************************** !********************************************************** !c $$$$$Run, Output: !********************************************************** Program main !******************************************************************72 ! !c S2DE_TEST tests the S2DE library. ! Implicit None Call timestamp() Write (, '(a)') '' Write (, '(a)') 'S2DE_PRB' Write (, '(a)') ' FORTRAN95,2013&18 version' Write (, '(a)') ' Test the S2DE library.' Call uvp_stokes1_test() Call resid_stokes1_test() Call gnuplot_stokes1_test() Call uvp_stokes2_test() Call resid_stokes2_test() Call gnuplot_stokes2_test() Call uvp_stokes3_test() Call resid_stokes3_test() Call gnuplot_stokes3_test() ! ! Terminate. ! Write (, '(a)') '' Write (, '(a)') 'S2DE_TEST' Write (, '(a)') ' Normal end of execution.' Call timestamp() Return End Program main Subroutine uvp_stokes1_test() !********************************************************************72 ! !c UVP_STOKES1_TEST samples the solution #1. ! Implicit None Integer n Parameter (n=1000) Double Precision p(n) Double Precision r8vec_max Double Precision r8vec_min Integer seed Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'UVP_STOKES1_TEST' Write (, '(a)') ' Exact Stokes solution #1:' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' using a region that is the unit square.' xy_lo = 0.0D+00 xy_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call uvp_stokes1(n, x, y, 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_stokes1_test Subroutine resid_stokes1_test() !********************************************************************72 ! !c RESID_STOKES1_TEST samples the residual for solution #1. ! Implicit None Integer n Parameter (n=1000) Double Precision pr(n) Double Precision r8vec_amax Double Precision r8vec_amin Integer seed Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'RESID_STOKES1_TEST' Write (, '(a)') ' Exact Stokes solution #1:' Write (, '(a)') ' Sample the Stokes residuals' Write (, '(a)') ' using a region that is the unit square.' xy_lo = 0.0D+00 xy_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call resid_stokes1(n, x, y, 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_stokes1_test Subroutine gnuplot_stokes1_test() !********************************************************************72 ! !c GNUPLOT_STOKES1_TEST plots solution #1 on a regular grid. ! Implicit None Integer x_num Parameter (x_num=21) Integer y_num Parameter (y_num=21) Character (255) header Integer n Double Precision p(x_num, y_num) Double Precision s Integer seed 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_STOKES1_TEST:' Write (, '(a)') ' Exact Stokes solution #1:' Write (, '(a)') ' Generate a Stokes 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) n = x_numy_num Call uvp_stokes1(n, x, y, u, v, p) header = 'stokes1' s = 4.0D+00 Call stokes_gnuplot(header, n, x, y, u, v, s) Return End Subroutine gnuplot_stokes1_test Subroutine uvp_stokes2_test() !********************************************************************72 ! !c UVP_STOKES2_TEST samples the solution #2. ! Implicit None Integer n Parameter (n=1000) Double Precision p(n) Double Precision r8vec_max Double Precision r8vec_min Integer seed Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'UVP_STOKES2_TEST' Write (, '(a)') ' Exact Stokes solution #2:' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' using a region that is the unit square.' xy_lo = 0.0D+00 xy_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call uvp_stokes2(n, x, y, 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_stokes2_test Subroutine resid_stokes2_test() !********************************************************************72 ! !c RESID_STOKES2_TEST samples the residual for solution #1. ! Implicit None Integer n Parameter (n=1000) Double Precision pr(n) Double Precision r8vec_amax Double Precision r8vec_amin Integer seed Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'RESID_STOKES2_TEST' Write (, '(a)') ' Exact Stokes solution #2:' Write (, '(a)') ' Sample the Stokes residuals' Write (, '(a)') ' using a region that is the unit square.' xy_lo = 0.0D+00 xy_hi = 1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call resid_stokes2(n, x, y, 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_stokes2_test Subroutine gnuplot_stokes2_test() !********************************************************************72 ! !c GNUPLOT_STOKES2_TEST plots solution #2 on a regular grid. ! Implicit None Integer x_num Parameter (x_num=21) Integer y_num Parameter (y_num=21) Character (255) header Integer n Double Precision p(x_num, y_num) Double Precision s Integer seed 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_STOKES2_TEST:' Write (, '(a)') ' Exact Stokes solution #2:' Write (, '(a)') ' Generate a Stokes 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) n = x_numy_num Call uvp_stokes2(n, x, y, u, v, p) header = 'stokes2' s = 0.05D+00 Call stokes_gnuplot(header, n, x, y, u, v, s) Return End Subroutine gnuplot_stokes2_test Subroutine uvp_stokes3_test() !********************************************************************72 ! !c UVP_STOKES3_TEST samples the solution #3. ! Implicit None Integer n Parameter (n=1000) Double Precision p(n) Double Precision r8vec_max Double Precision r8vec_min Integer seed Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'UVP_STOKES3_TEST' Write (, '(a)') ' Exact Stokes solution #3:' Write (, '(a)') ' Estimate the range of velocity and pressure' Write (, '(a)') ' using a region that is [-1,+1]x[-1,+1].' xy_lo = -1.0D+00 xy_hi = +1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call uvp_stokes3(n, x, y, 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_stokes3_test Subroutine resid_stokes3_test() !********************************************************************72 ! !c RESID_STOKES3_TEST samples the residual for solution #3. ! Implicit None Integer n Parameter (n=1000) Double Precision pr(n) Double Precision r8vec_amax Double Precision r8vec_amin Integer seed Double Precision ur(n) Double Precision vr(n) Double Precision x(n) Double Precision xy_hi Double Precision xy_lo Double Precision y(n) Write (, '(a)') '' Write (, '(a)') 'RESID_STOKES3_TEST' Write (, '(a)') ' Exact Stokes solution #3:' Write (, '(a)') ' Sample the Stokes residuals' Write (, '(a)') ' using a region that is [-1,+1]x[-1,+1].' xy_lo = -1.0D+00 xy_hi = +1.0D+00 seed = 123456789 Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, x) Call r8vec_uniform_ab(n, xy_lo, xy_hi, seed, y) Call resid_stokes3(n, x, y, 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_stokes3_test Subroutine gnuplot_stokes3_test() !********************************************************************72 ! !c GNUPLOT_STOKES3_TEST plots solution #3 on a regular grid. ! Implicit None Integer x_num Parameter (x_num=21) Integer y_num Parameter (y_num=21) Character (255) header Integer n Double Precision p(x_num, y_num) Double Precision s Integer seed 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_STOKES3_TEST:' Write (, '(a)') ' Exact Stokes solution #3:' Write (, '(a)') ' Generate a Stokes velocity field on [-1,+1]x[-1,+1].' Write (, '(a)') ' Store in GNUPLOT data and command files.' x_lo = -1.0D+00 x_hi = +1.0D+00 y_lo = -1.0D+00 y_hi = +1.0D+00 Call grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y) n = x_numy_num Call uvp_stokes3(n, x, y, u, v, p) header = 'stokes3' s = 0.05D+00 Call stokes_gnuplot(header, n, x, y, u, v, s) Return End Subroutine gnuplot_stokes3_test !*************************** Subroutine get_unit(iunit) !*****************************************************************72 ! !c GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! Implicit None Integer i Integer iunit Logical value iunit = 0 Do i = 1, 99 If (i/=5 .And. i/=6 .And. i/=9) Then Inquire (Unit=i, Opened=value, Err=10) If (.Not. value) Then iunit = i Return End If End If 10 Continue End Do Return End Subroutine get_unit Subroutine grid_2d(x_num, x_lo, x_hi, y_num, y_lo, y_hi, x, y) !******************************************************************72 ! !c GRID_2D returns a regular 2D grid. ! ! Parameters: ! ! Input, integer X_NUM, the number of X values to use. ! ! Input, double precision X_LO, X_HI, the range of X values. ! ! Input, integer Y_NUM, the number of Y values to use. ! ! Input, double precision Y_LO, Y_HI, the range of Y values. ! ! Output, double precision X(X_NUM,Y_NUM), Y(X_NUM,Y_NUM), ! the coordinates of the grid. ! Implicit None Integer x_num Integer y_num Integer i Integer j Double Precision x(x_num, y_num) Double Precision x_hi Double Precision x_lo Double Precision xi Double Precision y(x_num, y_num) Double Precision y_hi Double Precision y_lo Double Precision yj If (x_num==1) Then Do j = 1, y_num Do i = 1, x_num x(i, j) = (x_lo+x_hi)/2.0D+00 End Do End Do Else Do i = 1, x_num xi = (dble(x_num-i)x_lo+dble(i-1)x_hi)/dble(x_num-1) Do j = 1, y_num x(i, j) = xi End Do End Do End If If (y_num==1) Then Do j = 1, y_num Do i = 1, x_num y(i, j) = (y_lo+y_hi)/2.0D+00 End Do End Do Else Do j = 1, y_num yj = (dble(y_num-j)y_lo+dble(j-1)y_hi)/dble(y_num-1) Do i = 1, x_num y(i, j) = yj End Do End Do End If Return End Subroutine grid_2d Function r8vec_amax(n, a) !******************************************************************72 ! !c R8VEC_AMAX returns the maximum absolute value in an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, double precision A(N), the array. ! ! Output, double precision R8VEC_AMAX, the value of the entry ! of largest magnitude. ! Implicit None Integer n Double Precision a(n) Integer i Double Precision r8vec_amax Double Precision value value = 0.0D+00 Do i = 1, n value = max(value, abs(a(i))) End Do r8vec_amax = value Return End Function r8vec_amax Function r8vec_amin(n, a) !*****************************************************************72 ! !c R8VEC_AMIN returns the minimum absolute value in an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, double precisionA(N), the array. ! ! Output, double precision R8VEC_AMIN, the value of the entry ! of smallest magnitude. ! Implicit None Integer n Double Precision a(n) Integer i Double Precision r8_huge Parameter (r8_huge=1.79769313486231571D+308) Double Precision r8vec_amin Double Precision value value = r8_huge Do i = 1, n value = min(value, abs(a(i))) End Do r8vec_amin = value Return End Function r8vec_amin Function r8vec_max(n, a) !*****************************************************************72 ! !c R8VEC_MAX returns the maximum value in an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, double precision A(N), the array. ! ! Output, double precision R8VEC_MAX, the value of the largest entry. ! Implicit None Integer n Double Precision a(n) Integer i Double Precision r8vec_max Double Precision value value = a(1) Do i = 2, n value = max(value, a(i)) End Do r8vec_max = value Return End Function r8vec_max Function r8vec_min(n, a) !*****************************************************************72 ! !c R8VEC_MIN returns the minimum value in an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, double precision A(N), the array. ! ! Output, double precision R8VEC_MIN, the value of the smallest entry. ! Implicit None Integer n Double Precision a(n) Integer i Double Precision r8vec_min Double Precision value value = a(1) Do i = 2, n value = min(value, a(i)) End Do r8vec_min = value Return End Function r8vec_min Function r8vec_norm_l2(n, a) !******************************************************************72 ! !c R8VEC_NORM_L2 returns the L2 norm of an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8 values. ! ! The vector L2 norm is defined as: ! ! R8VEC_NORM_L2 = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ). ! ! Parameters: ! ! Input, integer N, the number of entries in A. ! ! Input, double precision A(N), the vector whose L2 norm is desired. ! ! Output, double precision R8VEC_NORM_L2, the L2 norm of A. ! Implicit None Integer n Double Precision a(n) Integer i Double Precision r8vec_norm_l2 Double Precision value value = 0.0D+00 Do i = 1, n value = value + a(i)a(i) End Do value = sqrt(value) r8vec_norm_l2 = value Return End Function r8vec_norm_l2 Subroutine r8vec_uniform_ab(n, a, b, seed, r) !********************************************************************72 ! !c R8VEC_UNIFORM_AB returns a scaled pseudorandom R8VEC. ! ! Discussion: ! ! Each dimension ranges from A to B. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, double precision A, B, the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which should NOT be 0. ! On output, SEED has been updated. ! ! Output, double precision R(N), the vector of pseudorandom values. ! Implicit None Integer n Double Precision a Double Precision b Integer i Integer i4_huge Parameter (i4_huge=2147483647) Integer k Integer seed Double Precision r(n) If (seed==0) Then Write (, '(a)') ' ' Write (, '(a)') 'R8VEC_UNIFORM_AB - Fatal error!' Write (, '(a)') ' Input value of SEED = 0.' Stop 1 End If Do i = 1, n k = seed/127773 seed = 16807(seed-k127773) - k2836 If (seed<0) Then seed = seed + i4_huge End If r(i) = a + (b-a)dble(seed)4.656612875D-10 End Do Return End Subroutine r8vec_uniform_ab Subroutine resid_stokes1(n, x, y, ur, vr, pr) !*******************************************************************72 ! !c RESID_STOKES1 returns residuals of the exact Stokes solution #1. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision UR(N), VR(N), PR(N), the residuals in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision p Double Precision pr(n) Double Precision px Double Precision py Double Precision u Double Precision ur(n) Double Precision ux Double Precision uxx Double Precision uy Double Precision uyy Double Precision v Double Precision vr(n) Double Precision vx Double Precision vxx Double Precision vy Double Precision vyy Double Precision x(n) Double Precision y(n) ! ! Get the right hand sides. ! Call rhs_stokes1(n, x, y, f, g, h) ! ! Form the functions and derivatives. ! Do i = 1, n u = -2.0D+00x(i)*2(x(i)-1.0D+00)*2y(i)(y(i)-1.0D+00)& (2.0D+00y(i)-1.0D+00) ux = -2.0D+00(4.0D+00x(i)3-6.0D+00x(i)*2+2.0D+00x(i))y(i)& (y(i)-1.0D+00)(2.0D+00y(i)-1.0D+00) uxx = -2.0D+00(12.0D+00x(i)*2-12.0D+00x(i)+2.0D+00)& (2.0D+00y(i)*3-3.0D+00y(i)*2+y(i)) uy = -2.0D+00x(i)*2(x(i)-1.0D+00)*2& (6.0D+00y(i)2-3.0D+00y(i)+1.0D+00) uyy = -2.0D+00(x(i)4-2.0D+00x(i)3+x(i)2)(12.0D+00y(i)-6.0D+00) v = 2.0D+00x(i)(x(i)-1.0D+00)(2.0D+00x(i)-1.0D+00)y(i)2& (y(i)-1.0D+00)*2 vx = 2.0D+00(6.0D+00x(i)2-6.0D+00x(i)+1.0D+00)y(i)2& (y(i)-1.0D+00)*2 vxx = 2.0D+00(12.0D+00x(i)-6.0D+00)y(i)*2(y(i)-1.0D+00)*2 vy = 2.0D+00x(i)(x(i)-1.0D+00)(2.0D+00x(i)-1.0D+00)& (4.0D+00y(i)3-6.0D+00y(i)*2+2.0D+00y(i)) vyy = 2.0D+00x(i)(x(i)-1.0D+00)(2.0D+00x(i)-1.0D+00)& (12.0D+00y(i)*2-12.0D+00y(i)+2.0D+00) p = 0.0D+00 px = 0.0D+00 py = 0.0D+00 ur(i) = px - (uxx+uyy) - f(i) vr(i) = py - (vxx+vyy) - g(i) pr(i) = ux + vy - h(i) End Do Return End Subroutine resid_stokes1 Subroutine resid_stokes2(n, x, y, ur, vr, pr) !********************************************************************72 ! !c RESID_STOKES2 returns residuals of the exact Stokes solution #2. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision UR(N), VR(N), PR(N), the residuals in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision p Double Precision pr(n) Double Precision px Double Precision py Double Precision r8_pi Parameter (r8_pi=3.141592653589793D+00) Double Precision u Double Precision ur(n) Double Precision ux Double Precision uxx Double Precision uy Double Precision uyy Double Precision v Double Precision vr(n) Double Precision vx Double Precision vxx Double Precision vy Double Precision vyy Double Precision x(n) Double Precision y(n) ! ! Get the right hand sides. ! Call rhs_stokes2(n, x, y, f, g, h) Do i = 1, n u = 2.0D+00sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) ux = 4.0D+00r8_picos(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) uxx = -8.0D+00r8_pi2sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) uy = -4.0D+00r8_pisin(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) uyy = -8.0D+00r8_pi2sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) v = -2.0D+00cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vx = 4.0D+00r8_pisin(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vxx = 8.0D+00r8_pi2cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vy = -4.0D+00r8_picos(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) vyy = 8.0D+00r8_pi2cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) p = x(i)2 + y(i)2 px = 2.0D+00x(i) py = 2.0D+00y(i) ur(i) = px - (uxx+uyy) - f(i) vr(i) = py - (vxx+vyy) - g(i) pr(i) = ux + vy - h(i) End Do Return End Subroutine resid_stokes2 Subroutine resid_stokes3(n, x, y, ur, vr, pr) !*******************************************************************72 ! !c RESID_STOKES3 returns residuals of the exact Stokes solution #3. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision UR(N), VR(N), PR(N), the residuals in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision p Double Precision pr(n) Double Precision px Double Precision py Double Precision u Double Precision ur(n) Double Precision ux Double Precision uxx Double Precision uy Double Precision uyy Double Precision v Double Precision vr(n) Double Precision vx Double Precision vxx Double Precision vy Double Precision vyy Double Precision x(n) Double Precision y(n) ! ! Get the right hand sides. ! Call rhs_stokes3(n, x, y, f, g, h) ! ! Form the functions and derivatives. ! Do i = 1, n u = 20.0D+00x(i)y(i)3 ux = 20.0D+00y(i)*3 uxx = 0.0D+00 uy = 60.0D+00x(i)y(i)2 uyy = 120.0D+00x(i)y(i) v = 5.0D+00(x(i)4-y(i)4) vx = 20.0D+00x(i)3 vxx = 60.0D+00x(i)*2 vy = -20.0D+00y(i)*3 vyy = -60.0D+00y(i)*2 p = 60.0D+00x(i)*2y(i) - 20.0D+00y(i)3 + 10.0D+00 px = 120.0D+00x(i)y(i) py = 60.0D+00x(i)*2 - 60.0D+00y(i)2 ur(i) = px - (uxx+uyy) - f(i) vr(i) = py - (vxx+vyy) - g(i) pr(i) = ux + vy - h(i) End Do Return End Subroutine resid_stokes3 Subroutine rhs_stokes1(n, x, y, f, g, h) !******************************************************************72 ! !c RHS_STOKES1 returns the right hand sides of the exact Stokes solution #1. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision F(N), G(N), H(N), the right hand sides in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision x(n) Double Precision y(n) Do i = 1, n f(i) = +2.0D+00(12.0D+00x(i)2-12.0D+00x(i)+2.0D+00)& (2.0D+00y(i)*3-3.0D+00y(i)*2+y(i)) + 2.0D+00& (x(i)*4-2.0D+00x(i)3+x(i)2)(12.0D+00y(i)-6.0D+00) g(i) = -2.0D+00(12.0D+00x(i)-6.0D+00)(y(i)4-2.0D+00y(i)3+y(i)2)& - 2.0D+00(2.0D+00x(i)*3-3.0D+00x(i)*2+x(i))& (12.0D+00y(i)2-12.0D+00y(i)+2.0D+00) h(i) = 0.0D+00 End Do Return End Subroutine rhs_stokes1 Subroutine rhs_stokes2(n, x, y, f, g, h) !********************************************************************72 ! !c RHS_STOKES2 returns the right hand sides of the exact Stokes solution #2. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision F(N), G(N), H(N), the right hand sides in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision p Double Precision px Double Precision py Double Precision r8_pi Parameter (r8_pi=3.141592653589793D+00) Double Precision u Double Precision ux Double Precision uxx Double Precision uy Double Precision uyy Double Precision v Double Precision vx Double Precision vxx Double Precision vy Double Precision vyy Double Precision x(n) Double Precision y(n) Do i = 1, n u = 2.0D+00sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) ux = 4.0D+00r8_picos(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) uxx = -8.0D+00r8_pi2sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) uy = -4.0D+00r8_pisin(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) uyy = -8.0D+00r8_pi2sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) v = -2.0D+00cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vx = 4.0D+00r8_pisin(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vxx = 8.0D+00r8_pi2cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) vy = -4.0D+00r8_picos(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) vyy = 8.0D+00r8_pi2cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) p = x(i)2 + y(i)2 px = 2.0D+00x(i) py = 2.0D+00y(i) f(i) = px - (uxx+uyy) g(i) = py - (vxx+vyy) h(i) = ux + vy End Do Return End Subroutine rhs_stokes2 Subroutine rhs_stokes3(n, x, y, f, g, h) !******************************************************************72 ! !c RHS_STOKES3 returns the right hand sides of the exact Stokes solution #3. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision F(N), G(N), H(N), the right hand sides in the U, ! V and P equations. ! Implicit None Integer n Double Precision f(n) Double Precision g(n) Double Precision h(n) Integer i Double Precision x(n) Double Precision y(n) Do i = 1, n f(i) = 0.0D+00 g(i) = 0.0D+00 h(i) = 0.0D+00 End Do Return End Subroutine rhs_stokes3 Subroutine stokes_gnuplot(header, n, x, y, u, v, s) !******************************************************************72 ! !c STOKES_GNUPLOT writes the Stokes velocity field to files for GNUPLOT. ! ! Parameters: ! ! Input, character ( * ) HEADER, a header to be used to name the files. ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the ! evaluation points. ! ! Input, double precision U(N), V(N), the velocity components. ! ! Input, double precision S, a scale factor for the velocity vectors. ! Implicit None Integer n Character (255) command_filename Integer command_unit Character (255) data_filename Integer data_unit Character () header Integer i Character (255) plot_filename Double Precision s Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) ! ! Write the data file. ! data_filename = trim(header) // '_data.txt' Call get_unit(data_unit) Open (Unit=data_unit, File=data_filename, Status='replace') Do i = 1, n Write (data_unit, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)')& x(i), y(i), su(i), sv(i) End Do Close (Unit=data_unit) Write (, '(a)') '' Write (, '(a)') ' Data written to "' // trim(data_filename) // '".' ! ! Write the command file. ! command_filename = trim(header) // '_commands.txt' Call get_unit(command_unit) plot_filename = trim(header) // '.png' Open (Unit=command_unit, File=command_filename, Status='replace') Write (command_unit, '(a)') '# ' // trim(command_filename) Write (command_unit, '(a)') '#' Write (command_unit, '(a)') 'set term png' Write (command_unit, '(a)') 'set output "' // trim(plot_filename) // '"' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') '# Add s and labels.' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') 'set xlabel "<--- X --->"' Write (command_unit, '(a)') 'set ylabel "<--- Y --->"' Write (command_unit, '(a)') 'set "Stokes flow"' Write (command_unit, '(a)') 'unset key' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') '# Add grid lines.' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') 'set grid' Write (command_unit, '(a)') 'set size ratio -1' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') '# Timestamp the plot.' Write (command_unit, '(a)') '#' Write (command_unit, '(a)') 'set timestamp' Write (command_unit, '(a)') 'plot "' // trim(data_filename) // & '" using 1:2:3:4 with vectors ' Write (command_unit, '(a)') ' head filled lt 2 linecolor rgb "blue"' Write (command_unit, '(a)') 'quit' Close (Unit=command_unit) Write (, '(a)') ' Commands written to "' // trim(command_filename) // '".' Return End Subroutine stokes_gnuplot Subroutine timestamp() !********************************************************************72 ! !c TIMESTAMP prints out the current YMDHMS date as a timestamp. ! ! Parameters: ! ! None ! Implicit None Character (8) ampm Integer d Character (8) date Integer h Integer m Integer mm Character (9) month(12) Integer n Integer s Character (10) time Integer y Save month Data month/'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', 'September',& 'October ', 'November ', 'December '/ Call date_and_time(date, time) Read (date, '(i4,i2,i2)') y, m, d Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm If (h<12) Then ampm = 'AM' Else If (h==12) Then If (n==0 .And. s==0) Then ampm = 'Noon' Else ampm = 'PM' End If Else h = h - 12 If (h<12) Then ampm = 'PM' Else If (h==12) Then If (n==0 .And. s==0) Then ampm = 'Midnight' Else ampm = 'AM' End If End If End If Write (, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,& month(m), y, h, ':', n, ':', s, '.', mm, ampm Return End Subroutine timestamp Subroutine uvp_stokes1(n, x, y, u, v, p) !****************************************************************************80 ! !c UVP_STOKES1 evaluates the exact Stokes solution #1. ! ! Discussion: ! ! The solution is defined over the unit square [0,1]x[0,1]. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision U(N), V(N), P(N), the velocity components and ! pressure at each of the points. ! Implicit None Integer n Integer i Double Precision p(n) Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) Do i = 1, n u(i) = -2.0D+00x(i)*2(x(i)-1.0D+00)*2y(i)(y(i)-1.0D+00)& (2.0D+00y(i)-1.0D+00) v(i) = 2.0D+00x(i)(x(i)-1.0D+00)(2.0D+00x(i)-1.0D+00)y(i)*2& (y(i)-1.0D+00)2 p(i) = 0.0D+00 End Do Return End Subroutine uvp_stokes1 Subroutine uvp_stokes2(n, x, y, u, v, p) !**************************************************************************80 ! !c UVP_STOKES2 evaluates the exact Stokes solution #2. ! ! Discussion: ! ! The solution is defined over the unit square [0,1]x[0,1]. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision U(N), V(N), P(N), the velocity components and ! pressure at each of the points. ! Implicit None Integer n Integer i Double Precision p(n) Double Precision r8_pi Parameter (r8_pi=3.141592653589793D+00) Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) Do i = 1, n u(i) = 2.0D+00sin(2.0D+00r8_pix(i))cos(2.0D+00r8_piy(i)) v(i) = -2.0D+00cos(2.0D+00r8_pix(i))sin(2.0D+00r8_piy(i)) p(i) = x(i)2 + y(i)2 End Do Return End Subroutine uvp_stokes2 Subroutine uvp_stokes3(n, x, y, u, v, p) !*******************************************************************72 ! !c UVP_STOKES3 evaluates the exact Stokes solution #3. ! ! Discussion: ! ! The solution is defined over the unit square [-1,+1]x[-1,+1]. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, double precision X(N), Y(N), the coordinates of the points. ! ! Output, double precision U(N), V(N), P(N), the velocity components and ! pressure at each of the points. ! Implicit None Integer n Integer i Double Precision p(n) Double Precision u(n) Double Precision v(n) Double Precision x(n) Double Precision y(n) Do i = 1, n u(i) = 20.0D+00x(i)y(i)3 v(i) = 5.0D+00(x(i)4-y(i)4) p(i) = 60.0D+00x(i)2y(i) - 2.0D+00y(i)3 + 10.0D+00 End Do Return End Subroutine uvp_stokes3 !*********************************************************** !c $$$$$Run, Output: !********************************************************** ******************************* # stokes1_commands.txt # set term png set output "stokes1.png" ............ ************************** 0.00000 0.00000 -0.00000 0.00000 0.500000E-01 0.00000 -0.00000 0.00000 0.100000 0.00000 -0.00000 0.00000 .................... 1.00000 1.00000 -0.00000 0.00000 ******************************* # stokes2_commands.txt # set term png set output "stokes2.png" # ............. ******************************* 0.00000 0.00000 0.00000 -0.00000 0.500000E-01 0.00000 0.309017E-01 -0.00000 0.100000 0.00000 0.587785E-01 -0.00000 ............................. 0.900000 1.00000 -0.587785E-01 0.198152E-16 0.950000 1.00000 -0.309017E-01 0.232942E-16 1.00000 1.00000 -0.244929E-16 0.244929E-16 ******************************* -1.00000 -1.00000 1.00000 0.00000 -0.900000 -1.00000 0.900000 -0.859750E-01 -0.800000 -1.00000 0.800000 -0.147600 ......................... 0.800000 1.00000 0.800000 -0.147600 0.900000 1.00000 0.900000 -0.859750E-01 1.00000 1.00000 1.00000 0.00000 ******************************* ********************************** S2DE_PRB FORTRAN95,2913&18 version Test the S2DE library. UVP_STOKES1_TEST Exact Stokes solution #1: Estimate the range of velocity and pressure using a region that is the unit square. Minimum Maximum U: -0.119905E-01 0.119441E-01 V: -0.119910E-01 0.120032E-01 P: 0.00000 0.00000 RESID_STOKES1_TEST Exact Stokes solution #1: Sample the Stokes residuals using a region that is the unit square. Minimum Maximum Ur: 0.00000 0.00000 Vr: 0.00000 0.246331E-14 Pr: 0.00000 0.142247E-15 GNUPLOT_STOKES1_TEST: Exact Stokes solution #1: Generate a Stokes velocity field on a regular grid. Store in GNUPLOT data and command files. Data written to "stokes1_data.txt". Commands written to "stokes1_commands.txt". UVP_STOKES2_TEST Exact Stokes solution #2: Estimate the range of velocity and pressure using a region that is the unit square. Minimum Maximum U: -1.99485 1.98591 V: -1.99754 1.99662 P: 0.301972E-03 1.95553 RESID_STOKES2_TEST Exact Stokes solution #2: Sample the Stokes residuals using a region that is the unit square. Minimum Maximum Ur: 0.00000 0.00000 Vr: 0.00000 0.00000 Pr: 0.00000 0.00000 GNUPLOT_STOKES2_TEST: Exact Stokes solution #2: Generate a Stokes velocity field on a regular grid. Store in GNUPLOT data and command files. Data written to "stokes2_data.txt". Commands written to "stokes2_commands.txt". UVP_STOKES3_TEST Exact Stokes solution #3: Estimate the range of velocity and pressure using a region that is [-1,+1]x[-1,+1]. Minimum Maximum U: -18.9476 18.3862 V: -4.87497 4.89690 P: -44.5447 66.6661 RESID_STOKES3_TEST Exact Stokes solution #3: Sample the Stokes residuals using a region that is [-1,+1]x[-1,+1]. Minimum Maximum Ur: 0.00000 0.00000 Vr: 0.00000 0.00000 Pr: 0.00000 0.00000 GNUPLOT_STOKES3_TEST: Exact Stokes solution #3: Generate a Stokes velocity field on [-1,+1]x[-1,+1]. Store in GNUPLOT data and command files. Data written to "stokes3_data.txt". Commands written to "stokes3_commands.txt". S2DE_TEST Normal end of execution. 22 April 2022 8:55:27.899 AM ...Program finished with exit code 0 Press ENTER to exit console. ****************************** ****************************** 斯托克斯两维精准估算解答压缩性稳态斯托克斯方程于2D单位正方形, 运算如上,图形表达示意如下P01,P02,和P03. 流体动力力学斯托克斯速度压力场 ******************************** 广州 0
2022-11-29 16:58:36　　评论淘帖0 举报相关推荐 • 精准估算压缩性稳态斯托克斯方程于2D单位正方形程序图形示意 282 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之c 303 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之a 421 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之图形表达 592 • 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之b 393 • 纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域 508 • Lumileds推出了一款具有正方形发光面的LUXEON 3030 2D 2001 • C语言教程之计算正方形的周长 0 • 针对非正方形图像置乱的算法 9 • 基于三轴设备画正方形程序设计 1126