!************************************************************
!************************************************************
!纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之c
!************************************************************
! SINOMACH
!
! 纳维斯托克斯两维精准估算解决非压缩时间依赖方程于一任意2D领域之c
!
! 周勇 20221125
!**************************************************************
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Subroutine rhs_spiral(nu, rho, n, x, y, t, f, g, h)

!*********************************************************************72
!
!c RHS_SPIRAL evaluates the right hand side of the spiral flow problem.
!
! Discussion:
!
! The right hand side is artificially determined by the requirement
! that the specified values of U, V and P satisfy the discretized
! Navier Stokes and continuity equations.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the fluid density.
!
! Input, integer N, the number of nodes.
!
! Input, double precision X(N), Y(N), the coordinates of nodes.
!
! Input, double precision T, the current time.
!
! Output, double precision F(N), G(N), H(N), the right hand sides.
!
Implicit None

Integer n

Double Precision dpdx(n)
Double Precision dpdy(n)
Double Precision dudt(n)
Double Precision dudx(n)
Double Precision dudxx(n)
Double Precision dudy(n)
Double Precision dudyy(n)
Double Precision dvdt(n)
Double Precision dvdx(n)
Double Precision dvdxx(n)
Double Precision dvdy(n)
Double Precision dvdyy(n)
Double Precision f(n)
Double Precision g(n)
Double Precision h(n)
Integer i
Double Precision nu
Double Precision p(n)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

Do i = 1, n

u(i) = (1.0D+00+nu*t)*2.0D+00*(x(i)**4-2.0D+00*x(i)**3+x(i)**2)*&
(2.0D+00*y(i)**3-3.0D+00*y(i)**2+y(i))

dudt(i) = nu*2.0D+00*(x(i)**4-2.0D+00*x(i)**3+x(i)**2)*&
(2.0D+00*y(i)**3-3.0D+00*y(i)**2+y(i))

dudx(i) = (1.0D+00+nu*t)*2.0D+00*&
(4.0D+00*x(i)**3-6.0D+00*x(i)**2+2.0D+00*x(i))*&
(2.0D+00*y(i)**3-3.0D+00*y(i)**2+y(i))

dudxx(i) = (1.0D+00+nu*t)*2.0D+00*(12.0D+00*x(i)**2-12.0D+00*x(i)+2.0D+00)*&
(2.0D+00*y(i)**3-3.0D+00*y(i)**2+y(i))

dudy(i) = (1.0D+00+nu*t)*2.0D+00*(x(i)**4-2.0D+00*x(i)**3+x(i)**2)*&
(6.0D+00*y(i)**2-6.0D+00*y(i)+1.0D+00)

dudyy(i) = (1.0D+00+nu*t)*2.0D+00*(x(i)**4-2.0D+00*x(i)**3+x(i)**2)*&
(12.0D+00*y(i)-6.0D+00)

v(i) = -(1.0D+00+nu*t)*2.0D+00*(2.0D+00*x(i)**3-3.0D+00*x(i)**2+x(i))*&
(y(i)**4-2.0D+00*y(i)**3+y(i)**2)

dvdt(i) = -nu*2.0D+00*(2.0D+00*x(i)**3-3.0D+00*x(i)**2+x(i))*(y(i)**4-&
2.0D+00*y(i)**3+y(i)**2)

dvdx(i) = -(1.0D+00+nu*t)*2.0D+00*(6.0D+00*x(i)**2-6.0D+00*x(i)+1.0D+00)*&
(y(i)**4-2.0D+00*y(i)**3+y(i)**2)

dvdxx(i) = -(1.0D+00+nu*t)*2.0D+00*(12.0D+00*x(i)-6.0D+00)*&
(y(i)**4-2.0D+00*y(i)**3+y(i)**2)

dvdy(i) = -(1.0D+00+nu*t)*2.0D+00*(2.0D+00*x(i)**3-3.0D+00*x(i)**2+x(i))*&
(4.0D+00*y(i)**3-6.0D+00*y(i)**2+2.0D+00*y(i))

dvdyy(i) = -(1.0D+00+nu*t)*2.0D+00*(2.0D+00*x(i)**3-3.0D+00*x(i)**2+x(i))*&
(12.0D+00*y(i)**2-12.0D+00*y(i)+2.0D+00)

p(i) = rho*y(i)
dpdx(i) = 0.0D+00
dpdy(i) = rho

f(i) = dudt(i) - nu*(dudxx(i)+dudyy(i)) + u(i)*dudx(i) +&
v(i)*dudy(i) + dpdx(i)/rho

g(i) = dvdt(i) - nu*(dvdxx(i)+dvdyy(i)) + u(i)*dvdx(i) +&
v(i)*dvdy(i) + dpdy(i)/rho

h(i) = dudx(i) + dvdy(i)

End Do

Return
End Subroutine rhs_spiral
Subroutine rhs_taylor(nu, rho, n, x, y, t, f, g, h)

!*********************************************************************72
!
!c RHS_TAYLOR returns right hand sides of the Taylor vortex equations.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the density.
!
! Input, integer N, the number of evaluation points.
!
! Input, double precision X(N), Y(N), the coordinates of the points.
!
! Input, double precision T, the time coordinate or coordinates.
!
! 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 nu
Double Precision rho
Double Precision t
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_taylor
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_lucas(nu, rho, n, x, y, t, u, v, p)

!*****************************************************************************80
!
!c UVP_LUCAS evaluates Lucas Bystricky's exact Navier Stokes solution.
!
! Discussion:
!
! There is no time dependence.
!
! The pressure is 0.
!
! The preferred domain is the unit square.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the density.
!
! Input, integer N, the number of evaluation points.
!
! Input, double precision X(N), Y(N), the coordinates of the points.
!
! Input, double precision T, the time coordinate or coordinates.
!
! 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 nu
Double Precision p(n)
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

Do i = 1, n

u(i) = -cos(r8_pi*x(i))/r8_pi
v(i) = -y(i)*sin(r8_pi*x(i))
p(i) = 0.0D+00

End Do

Return
End Subroutine uvp_lucas
Subroutine uvp_spiral(nu, rho, n, x, y, t, u, v, p)

!*********************************************************************72
!
!c UVP_SPIRAL returns velocity and pressure for the spiral flow.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the fluid density.
!
! Input, integer N, the number of nodes.
!
! Input, double precision X(N), Y(N), the coordinates of nodes.
!
! Input, double precision T, the current time.
!
! Output, double precision U(N), V(N), the X and Y velocity.
!
! Output, double precision P(N), the pressure.
!
Implicit None

Integer n

Integer i
Double Precision nu
Double Precision p(n)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

Do i = 1, n

u(i) = (1.0D+00+nu*t)*2.0D+00*x(i)**2*(x(i)-1.0D+00)**2*y(i)*&
(2.0D+00*y(i)-1.0D+00)*(y(i)-1.0D+00)

v(i) = -(1.0D+00+nu*t)*2.0D+00*x(i)*(2.0D+00*x(i)-1.0D+00)*&
(x(i)-1.0D+00)*y(i)**2*(y(i)-1.0D+00)**2

p(i) = rho*y(i)

End Do

Return
End Subroutine uvp_spiral
Subroutine uvp_taylor(nu, rho, n, x, y, t, u, v, p)

!*****************************************************************************80
!
!c UVP_TAYLOR evaluates the Taylor exact Navier Stokes solution.
!
! Discussion:
!
! This flow is known as a Taylor-Green vortex.
!
! The given velocity and pressure fields are exact solutions for the 2D
! incompressible time-dependent Navier Stokes equations over any region.
!
! To define a typical problem, one chooses a bounded spatial region
! and a starting time, and then imposes boundary and initial conditions
! by referencing the exact solution appropriately.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the density.
!
! Input, integer N, the number of evaluation points.
!
! Input, double precision X(N), Y(N), the coordinates of the points.
!
! Input, double precision T, the time coordinate or coordinates.
!
! 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 nu
Double Precision p(n)
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision v(n)
Double Precision x(n)
Double Precision y(n)

Do i = 1, n
u(i) = -cos(r8_pix(i))sin(r8_piy(i))
v(i) = sin(r8_pix(i))cos(r8_piy(i))
p(i) = -0.25D+00rho(cos(2.0D+00r8_pix(i))+cos(2.0D+00r8_piy(i)))
End Do

Do i = 1, n
u(i) = u(i)exp(-2.0D+00r8_pi2nut)
v(i) = v(i)exp(-2.0D+00r8_pi2nut)
p(i) = p(i)exp(-4.0D+00r8_pi**2nut)
End Do

Return
End Subroutine uvp_taylor

set term png
set output "lucas.png"
........................
plot "lucas_data.txt" using 1:2:3:4 with vectors
head filled lt 2 linecolor rgb "blue"
quit

0.00000         0.00000       -0.795775E-01    -0.00000

.....................
1.00000 1.00000 0.795775E-01 -0.306162E-16

spiral_commands.txt

set term png
........................
plot "spiral_data.txt" using 1:2:3:4 with vectors
head filled lt 2 linecolor rgb "blue"
quit

0.00000         0.00000         0.00000        -0.00000

.....................
1.00000 1.00000 0.00000 -0.00000

NS2DE_PRB
FORTRAN90,95,&2018 version
Test the NS2DE library.

UVP_TAYLOR_TEST
Taylor Vortex Flow
Estimate the range of velocity and pressure
at the initial time T = 0, using a region that is
the square centered at (1.5,1.5) with "radius" 1.0,
Kinematic viscosity NU = 1.00000
Fluid density RHO = 1.00000

Minimum       Maximum

U: -0.997425 0.992956
.....................
UVP_SPIRAL_TEST
Spiral Flow
Estimate the range of velocity and pressure
at the initial time T = 0, using a region that is
the unit square.
Kinematic viscosity NU = 1.00000
Fluid density RHO = 1.00000

Minimum       Maximum

U: -0.119441E-01 0.119905E-01

Data written to "spiral_data.txt".
Commands written to "spiral_commands.txt".

PARAMETER_SPIRAL_TEST
Spiral Flow
Monitor solution norms over time for
varous values of NU, RHO.

RHO affects the pressure scaling.

RHO         NU           T     ||U||       ||V||       ||P||

1.000 1.000 0.0000 0.1767E-03 0.1712E-03 0.1798E-01
..................
NU affects the time scaling.

RHO         NU           T     ||U||       ||V||       ||P||

1.000 1.000 0.0000 0.1767E-03 0.1712E-03 0.1798E-01
...............................
UVP_LUCAS_TEST
Lucas Bystricky Flow
Estimate the range of velocity and pressure
at the initial time T = 0, over the unit square.
Kinematic viscosity NU = 1.00000
Fluid density RHO = 1.00000

Minimum       Maximum

U: -0.318305 0.318303
.....................
GNUPLOT_LUCAS_TEST:
Lucas Bystricky Flow
Generate a velocity field on a regular grid.
Store in GNUPLOT data and command files.

Data written to "lucas_data.txt".
Commands written to "lucas_commands.txt".

NS2DE_TEST
Normal end of execution.
26 April 2022 2:00:56.949 AM

...Program finished with exit code 0
Press ENTER to exit console.

纳维斯托克斯两维精准估算解决
非压缩时间依赖方程于一任意2D领域,