纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域

纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域

托克斯,3D精准,非压缩方程

纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域

!

!c NS3DE_PRB tests the NS3DE library.

!

Implicit None

Call timestamp()

Write ( *, '(a)') ''

Write (* , '(a)') 'NS3DE_PRB'

Write ( *, '(a)') '  FORTRAN90,95,&2018 version'

Write (* , '(a)') '  Test the NS3DE library.'

Call test01()

Call test02()

!

!  Terminate.

!

Write ( *, '(a)') ''

Write (* , '(a)') 'NS3DE_PRB'

Write (*, '(a)') '  Normal end of execution.'

Call timestamp()

Return

End Program main

Subroutine test01()

!*********************************************************************72

!

!c TEST01 samples the solution at the initial time.

!

Implicit None

Integer n

Parameter (n=1000)

Double Precision a

Double Precision d

Double Precision p(n)

Double Precision r8_pi

Parameter (r8_pi=3.141592653589793D+00)

Double Precision r8vec_max

Double Precision r8vec_min

Integer seed

Double Precision t

Double Precision u(n)

Double Precision v(n)

Double Precision w(n)

Double Precision x(n)

Double Precision xyz_hi

Double Precision xyz_lo

Double Precision y(n)

Double Precision z(n)

a = r8_pi/4.0D+00

d = r8_pi/2.0D+00

Write ( *, '(a)') ''

Write (* , '(a)') 'TEST01'

Write ( *, '(a)') '  Estimate the range of velocity and pressure'

Write (* , '(a)') '  at the initial time T = 0, in a region that is the'

Write ( *, '(a)') '  cube centered at (0,0,0) with "radius" 1.0.'

Write (* , '(a,g14.6)') '  Parameter A = ', a

Write (*, '(a,g14.6)') '  Parameter D = ', d

xyz_lo = -1.0D+00

xyz_hi = +1.0D+00

seed = 123456789

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, x)

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, y)

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, z)

t = 0.0D+00

Call uvwp_ethier(a, d, n, x, y, z, t, u, v, w, p)

Write ( *, '(a)') ''

Write (* , '(a)') '           Minimum       Maximum'

Write ( *, '(a)') ''

Write (* , '(a,g14.6,2x,g14.6)') '  U:  ', r8vec_min(n, u), r8vec_max(n, u)

Write ( *, '(a,g14.6,2x,g14.6)') '  V:  ', r8vec_min(n, v), r8vec_max(n, v)

Write (* , '(a,g14.6,2x,g14.6)') '  W:  ', r8vec_min(n, w), r8vec_max(n, w)

Write (*, '(a,g14.6,2x,g14.6)') '  P:  ', r8vec_min(n, p), r8vec_max(n, p)

Return

End Subroutine test01

Subroutine test02()

!*********************************************************************72

!

!c TEST02 samples the residual at the initial time.

!

Implicit None

Integer n

Parameter (n=1000)

Double Precision a

Double Precision d

Double Precision pr(n)

Double Precision r8_pi

Parameter (r8_pi=3.141592653589793D+00)

Double Precision r8vec_amax

Double Precision r8vec_amin

Integer seed

Double Precision t

Double Precision ur(n)

Double Precision vr(n)

Double Precision wr(n)

Double Precision x(n)

Double Precision xyz_hi

Double Precision xyz_lo

Double Precision y(n)

Double Precision z(n)

a = r8_pi/4.0D+00

d = r8_pi/2.0D+00

Write ( *, '(a)') ''

Write (* , '(a)') 'TEST02'

Write ( *, '(a)') '  Sample the Navier-Stokes residuals'

Write (* , '(a)') '  at the initial time T = 0, using a region that is'

Write ( *, '(a)') '  the cube centered at (0,0,0) with "radius" 1.0,'

Write (* , '(a,g14.6)') '  Parameter A = ', a

Write (*, '(a,g14.6)') '  Parameter D = ', d

xyz_lo = -1.0D+00

xyz_hi = +1.0D+00

seed = 123456789

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, x)

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, y)

Call r8vec_uniform_ab(n, xyz_lo, xyz_hi, seed, z)

t = 0.0D+00

Call resid_ethier(a, d, n, x, y, z, t, ur, vr, wr, pr)

Write ( *, '(a)') ''

Write (* , '(a)') '           Minimum       Maximum'

Write ( *, '(a)') ''

Write (* , '(a,g14.6,2x,g14.6)') '  Ur:  ', r8vec_amin(n, ur), r8vec_amax(n, ur)

Write ( *, '(a,g14.6,2x,g14.6)') '  Vr:  ', r8vec_amin(n, vr), r8vec_amax(n, vr)

Write (* , '(a,g14.6,2x,g14.6)') '  Wr:  ', r8vec_amin(n, wr), r8vec_amax(n, wr)

Write (*, '(a,g14.6,2x,g14.6)') '  Pr:  ', r8vec_amin(n, pr), r8vec_amax(n, pr)

Return

End Subroutine test02

!******************************

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

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-k*127773) - k*2836

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_ethier(a, d, n, x, y, z, t, ur, vr, wr, pr)

!*********************************************************************72

!

!c RESID_ETHIER evaluates the residual of the Ethier exact Navier Stokes solution.

!

!  Parameters:

!

!    Input, double precision A, D, the parameters.  Sample values are A = PI/4

!    and D = PI/2.

!

!    Input, integer N, the number of points at which the solution

!    is to be evaluated.

!

!    Input, double precision X(N), Y(N), Z(N), the coordinates of the points.

!

!    Input, double precision T, the time coordinate or coordinates.

!

!    Output, double precision UR(N), VR(N), WR(N), PR(N), the residuals.

!

Implicit None

Integer n

Double Precision a

Double Precision cxy(n)

Double Precision cyz(n)

Double Precision czx(n)

Double Precision d

Double Precision e2t

Double Precision e2x(n)

Double Precision e2y(n)

Double Precision e2z(n)

Double Precision e4t

Double Precision ex(n)

Double Precision exy(n)

Double Precision ey(n)

Double Precision eyz(n)

Double Precision ez(n)

Double Precision ezx(n)

Integer i

Double Precision p(n)

Double Precision pr(n)

Double Precision px(n)

Double Precision py(n)

Double Precision pz(n)

Double Precision sxy(n)

Double Precision syz(n)

Double Precision szx(n)

Double Precision t

Double Precision u(n)

Double Precision ur(n)

Double Precision ut(n)

Double Precision ux(n)

Double Precision uxx(n)

Double Precision uy(n)

Double Precision uyy(n)

Double Precision uz(n)

Double Precision uzz(n)

Double Precision v(n)

Double Precision vr(n)

Double Precision vt(n)

Double Precision vx(n)

Double Precision vxx(n)

Double Precision vy(n)

Double Precision vyy(n)

Double Precision vz(n)

Double Precision vzz(n)

Double Precision w(n)

Double Precision wr(n)

Double Precision wt(n)

Double Precision wx(n)

Double Precision wxx(n)

Double Precision wy(n)

Double Precision wyy(n)

Double Precision wz(n)

Double Precision wzz(n)

Double Precision x(n)

Double Precision y(n)

Double Precision z(n)

!

!  Make some temporaries.

!

Do i = 1, n

ex(i) = exp(a*x(i))
ey(i) = exp(a*y(i))
ez(i) = exp(a*z(i))

e2x(i) = exp(2.0D+00*a*x(i))
e2y(i) = exp(2.0D+00*a*y(i))
e2z(i) = exp(2.0D+00*a*z(i))

e2t = exp(-d*d*t)
e4t = exp(-2.0*d*d*t)

exy(i) = exp(a*(x(i)+y(i)))
eyz(i) = exp(a*(y(i)+z(i)))
ezx(i) = exp(a*(z(i)+x(i)))

sxy(i) = sin(a*x(i)+d*y(i))
syz(i) = sin(a*y(i)+d*z(i))
szx(i) = sin(a*z(i)+d*x(i))

cxy(i) = cos(a*x(i)+d*y(i))
cyz(i) = cos(a*y(i)+d*z(i))
czx(i) = cos(a*z(i)+d*x(i))


!

!  Form the functions and derivatives.

!

u(i) = -a*(ex(i)*syz(i)+ez(i)*cxy(i)) *e2t

ux(i) = -a* (a*ex(i)*syz(i)-a*ez(i)*sxy(i)) *e2t

uxx(i) = -a* (a*a*ex(i)*syz(i)-a*a*ez(i)*cxy(i)) *e2t

uy(i) = -a* (a*ex(i)*cyz(i)-d*ez(i)*sxy(i)) *e2t

uyy(i) = -a* (-a*a*ex(i)*syz(i)-d*d*ez(i)*cxy(i)) *e2t

uz(i) = -a* (d*ex(i)*cyz(i)+a*ez(i)*cxy(i)) *e2t

uzz(i) = -a* (-d*d*ex(i)*syz(i)+a*a *ez(i)*cxy(i))*e2t

ut(i) = +d*d*a* (ex(i)*syz(i)+ez(i)*cxy(i))*e2t

v(i) = -a*(ey(i)*szx(i)+ex(i)*cyz(i))*e2t
vx(i) = -a*(d*ey(i)*czx(i)+a*ex(i)*cyz(i))*e2t
vxx(i) = -a*(-d*d*ey(i)*szx(i)+a*a*ex(i)*cyz(i))*e2t
vy(i) = -a*(a*ey(i)*szx(i)-a*ex(i)*syz(i))*e2t
vyy(i) = -a*(a*a*ey(i)*szx(i)-a*a*ex(i)*cyz(i))*e2t
vz(i) = -a*(a*ey(i)*czx(i)-d*ex(i)*syz(i))*e2t
vzz(i) = -a*(-a*a*ey(i)*szx(i)-d*d*ex(i)*cyz(i))*e2t
vt(i) = +d*d*a*(ey(i)*szx(i)+ex(i)*cyz(i))*e2t

w(i) = -a*(ez(i)*sxy(i)+ey(i)*czx(i))*e2t
wx(i) = -a*(a*ez(i)*cxy(i)-d*ey(i)*szx(i))*e2t
wxx(i) = -a*(-a*a*ez(i)*sxy(i)-d*d*ey(i)*czx(i))*e2t
wy(i) = -a*(d*ez(i)*cxy(i)+a*ey(i)*czx(i))*e2t
wyy(i) = -a*(-d*d*ez(i)*sxy(i)+a*a*ey(i)*czx(i))*e2t
wz(i) = -a*(a*ez(i)*sxy(i)-a*ey(i)*szx(i))*e2t
wzz(i) = -a*(a*a*ez(i)*sxy(i)-a*a*ey(i)*czx(i))*e2t
wt(i) = +d*d*a*(ez(i)*sxy(i)+ey(i)*czx(i))*e2t

p(i) = -0.5D+00*a*a*e4t*(+e2x(i)+2.0D+00*sxy(i)*czx(i)*eyz(i)+&
e2y(i)+2.0D+00*syz(i)*cxy(i)*ezx(i)+e2z(i)+2.0D+00*szx(i)*cyz(i)*exy(i))

px(i) = -0.5D+00*a*a*e4t*(+2.0D+00*a*e2x(i)+&
2.0D+00*a*cxy(i)*czx(i)*eyz(i)-2.0D+00*d*sxy(i)*szx(i)*eyz(i)-&
2.0D+00*a*syz(i)*sxy(i)*ezx(i)+2.0D+00*a*syz(i)*cxy(i)*ezx(i)+&
2.0D+00*d*czx(i)*cyz(i)*exy(i)+2.0D+00*a*szx(i)*cyz(i)*exy(i))

py(i) = -0.5D+00*a*a*e4t*(+2.0D+00*d*cxy(i)*czx(i)*eyz(i)+&
2.0D+00*a*sxy(i)*czx(i)*eyz(i)+2.0D+00*a*e2y(i)+&
2.0D+00*a*cyz(i)*cxy(i)*ezx(i)-2.0D+00*d*syz(i)*sxy(i)*ezx(i)-&
2.0D+00*a*szx(i)*syz(i)*exy(i)+2.0D+00*a*szx(i)*cyz(i)*exy(i))

pz(i) = -0.5D+00*a*a*e4t*&
(-2.0D+00*a*sxy(i)*szx(i)*eyz(i)+2.0D+00*a*sxy(i)*czx(i)*eyz(i)+&
2.0D+00*d*cyz(i)*cxy(i)*ezx(i)+2.0D+00*a*syz(i)*cxy(i)*ezx(i)+2.0D+00*&
a*e2z(i)+2.0D+00*a*czx(i)*cyz(i)*exy(i)-2.0D+00*d*szx(i)*syz(i)*exy(i))

!

!  Evaluate the residuals.

!

ur(i) = ut(i) + u(i)*ux(i) + v(i)*uy(i) + w(i)*uz(i) + px(i) -&

(uxx(i)+uyy(i)+uzz(i))


vr(i) = vt(i) + u(i)*vx(i) + v(i)*vy(i) + w(i)*vz(i) + py(i) -&
(vxx(i)+vyy(i)+vzz(i))

wr(i) = wt(i) + u(i)*wx(i) + v(i)*wy(i) + w(i)*wz(i) + pz(i) -&
(wxx(i)+wyy(i)+wzz(i))

pr(i) = ux(i) + vy(i) + wz(i)

End Do

Return

End Subroutine resid_ethier

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 uvwp_ethier(a, d, n, x, y, z, t, u, v, w, p)

!*********************************************************************72

!

!c UVWP_ETHIER evaluates the Ethier exact Navier Stokes solution.

!

!  Discussion:

!

!    The reference asserts that the given velocity and pressure fields

!    are exact solutions for the 3D incompressible time-dependent

!    Navier Stokes equations over any region.

!

!    To define a typical problem, one chooses a bounded spatial region

!    and a starting time, and then imposes boundary and initial conditions

!    by referencing the exact solution appropriately.

!

!    In the reference paper, a calculation is made for the cube centered

!    at (0,0,0) with a "radius" of 1 unit, and over the time interval

!    from t = 0 to t = 0.1, with parameters a = PI/4 and d = PI/2,

!    and with Dirichlet boundary conditions on all faces of the cube.

!

!  Parameters:

!

!    Input, double precision A, D, the parameters.  Sample values are A = PI/4

!    and D = PI/2.

!

!    Input, integer ( kind  N, the number of points at which the solution is to

!    be evaluated.

!

!    Input, double precision X(N), Y(N), Z(N), the coordinates of the points.

!

!    Input, double precision T, the time coordinate or coordinates.

!

!    Output, double precision U(N), V(N), W(N), P(N), the velocity components

!    and pressure at each of the points.

!

Implicit None

Integer n

Double Precision a

Double Precision cxy(n)

Double Precision cyz(n)

Double Precision czx(n)

Double Precision d

Double Precision e2t

Double Precision ex(n)

Double Precision exy(n)

Double Precision ey(n)

Double Precision eyz(n)

Double Precision ez(n)

Double Precision ezx(n)

Integer i

Double Precision p(n)

Double Precision sxy(n)

Double Precision syz(n)

Double Precision szx(n)

Double Precision t

Double Precision u(n)

Double Precision v(n)

Double Precision w(n)

Double Precision x(n)

Double Precision y(n)

Double Precision z(n)

Do i = 1, n


ex(i) = exp(a*x(i))
ey(i) = exp(a*y(i))
ez(i) = exp(a*z(i))

e2t = exp(-d*d*t)

exy(i) = exp(a*(x(i)+y(i)))
eyz(i) = exp(a*(y(i)+z(i)))
ezx(i) = exp(a*(z(i)+x(i)))

sxy(i) = sin(a*x(i)+d*y(i))
syz(i) = sin(a*y(i)+d*z(i))
szx(i) = sin(a*z(i)+d*x(i))

cxy(i) = cos(a*x(i)+d*y(i))
cyz(i) = cos(a*y(i)+d*z(i))
czx(i) = cos(a*z(i)+d*x(i))

u(i) = -a*(ex(i)*syz(i)+ez(i)*cxy(i))*e2t
v(i) = -a*(ey(i)*szx(i)+ex(i)*cyz(i))*e2t
w(i) = -a*(ez(i)*sxy(i)+ey(i)*czx(i))*e2t
p(i) = 0.5D+00*a*a*e2t*&
e2t*(+ex(i)*ex(i)+2.0D+00*sxy(i)*czx(i)*eyz(i)+ey(i)*ey(i)+&
2.0D+00*syz(i)*cxy(i)*ezx(i)+ez(i)*ez(i)+2.0D+00*szx(i)*cyz(i)*exy(i))

End Do

Return

End Subroutine uvwp_ethier


NS3DE_PRB

FORTRAN90,95,&2018 version

Test the NS3DE library.

TEST01

Estimate the range of velocity and pressure

at the initial time T = 0, in a region that is the

cube centered at (0,0,0) with "radius" 1.0.

Parameter A =   0.785398

...................


Minimum       Maximum

U:    -2.96676         1.47141

................................

TEST02

Sample the Navier-Stokes residuals

at the initial time T = 0, using a region that is

the cube centered at (0,0,0) with "radius" 1.0,

Parameter A =   0.785398

............


Minimum       Maximum

Ur:     0.00000        0.355271E-14

................................

NS3DE_PRB

Normal end of execution.

26 April     2022   2:12:17.639 AM

...Program finished with exit code 0

Press ENTER to exit console.

纳维斯托克斯三维精准估算非压缩时间依赖方程于一个3D领域.