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

!************************************************************
!************************************************************
! 纳维斯托克斯两维精准估算解决非压缩时间依赖方程(abc部)之b
!************************************************************
! SINOMACH
!
! 纳维斯托克斯两维精准估算解决非压缩时间依赖方程于一任意2D领域之b
!
! 周勇 20221125
!**************************************************************
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Subroutine gnuplot_lucas_test()

!*********************************************************************72
!
!c GNUPLOT_LUCAS_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_LUCAS_TEST:'
Write (, '(a)') ' Lucas Bystricky 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_lucas(nu, rho, n, x, y, t, u, v, p)

header = 'lucas'
s = 0.25D+00
Call ns2de_gnuplot(header, n, x, y, u, v, s)

Return
End Subroutine gnuplot_lucas_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
Subroutine ns2de_gnuplot(header, n, x, y, u, v, s)

!*********************************************************************72
!
!c NS2DE_GNUPLOT writes the Navier-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 "Navier-Stokes velocity field"'
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 ns2de_gnuplot
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
Subroutine r8vec_linspace(n, a, b, x)

!*********************************************************************72
!
!c R8VEC_LINSPACE creates a vector of linearly spaced values.
!
! Discussion:
!
! An R8VEC is a vector of R8's.
!
! 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
!
! In other words, the interval is divided into N-1 even subintervals,
! and the endpoints of intervals are used as the points.
!
! Parameters:
!
! Input, integer N, the number of entries in the vector.
!
! Input, double precision A, B, the first and last entries.
!
! Output, double precision X(N), a vector of linearly spaced data.
!
Implicit None

Integer n

Double Precision a
Double Precision b
Integer i
Double Precision x(n)

If (n==1) Then

x(1) = (a+b)/2.0D+00

Else

Do i = 1, n
  x(i) = (dble(n-i)*a+dble(i-1)*b)/dble(n-1)
End Do

End If

Return
End Subroutine r8vec_linspace
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-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_lucas(nu, rho, n, x, y, t, ur, vr, pr)

!*****************************************************************************80
!
!c RESID_LUCAS returns Lucas Bystricky residuals.
!
! 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 UR(N), VR(N), PR(N), the residuals in the U,
! V and P equations.
!
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 pr(n)
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision ur(n)
Double Precision v(n)
Double Precision vr(n)
Double Precision x(n)
Double Precision y(n)
!
! Get the right hand side functions.
!
Call rhs_lucas(nu, rho, n, x, y, t, f, g, h)
!
! Form the functions and derivatives for the left hand side.
!
Do i = 1, n

u(i) = -cos(r8_pi*x(i))/r8_pi
dudt(i) = 0.0D+00
dudx(i) = sin(r8_pi*x(i))
dudxx(i) = r8_pi*cos(r8_pi*x(i))
dudy(i) = 0.0D+00
dudyy(i) = 0.0D+00

v(i) = -y(i)*sin(r8_pi*x(i))
dvdt(i) = 0.0D+00
dvdx(i) = -r8_pi*y(i)*cos(r8_pi*x(i))
dvdxx(i) = +r8_pi**2*y(i)*sin(r8_pi*x(i))
dvdy(i) = -sin(r8_pi*x(i))
dvdyy(i) = 0.0D+00

p(i) = 0.0D+00
dpdx(i) = 0.0D+00
dpdy(i) = 0.0D+00

!
! Evaluate the residuals.
!
ur(i) = dudt(i) + u(i)*dudx(i) + v(i)*dudy(i) + (1.0D+00/rho)dpdx(i) -&
nu(dudxx(i)+dudyy(i)) - f(i)

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

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

End Do

Return
End Subroutine resid_lucas
Subroutine resid_spiral(nu, rho, n, x, y, t, ur, vr, pr)

!*********************************************************************72
!
!c RESID_SPIRAL evaluates the pointwise residual of the spiral flow problem.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the 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.
!
! Input, double precision U(N), V(N), the X and Y velocity.
!
! Input, double precision P(N), the pressure.
!
! Output, double precision UR(N), the U-equation residual.
!
! Output, double precision VR(N), the V-equation residual.
!
! Output, double precision PR(N), the P-equation residual.
!
Implicit None

Integer (Kind=4) 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 pr(n)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision ur(n)
Double Precision v(n)
Double Precision vr(n)
Double Precision x(n)
Double Precision y(n)
!
! Get the right hand side functions.
!
Call rhs_spiral(nu, rho, n, x, y, t, f, g, h)
!
! Form the functions and derivatives for the left hand side.
!
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

!
! Evaluate the residuals.
!
ur(i) = dudt(i) - nu*(dudxx(i)+dudyy(i)) + u(i)*dudx(i) + v(i)*dudy(i) +&
dpdx(i)/rho - f(i)

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

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

End Do

Return
End Subroutine resid_spiral
Subroutine resid_taylor(nu, rho, n, x, y, t, ur, vr, pr)

!*********************************************************************72
!
!c RESID_TAYLOR returns residuals of the Taylor exact Navier Stokes solution.
!
! 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 UR(N), VR(N), PR(N), the residuals in the U,
! V and P equations.
!
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 pr(n)
Double Precision r8_pi
Parameter (r8_pi=3.141592653589793D+00)
Double Precision rho
Double Precision t
Double Precision u(n)
Double Precision ur(n)
Double Precision v(n)
Double Precision vr(n)
Double Precision x(n)
Double Precision y(n)
!
! Get the right hand sides.
!
Call rhs_taylor(nu, rho, n, x, y, t, f, g, h)
!
! Form the functions and derivatives for the left hand side.
!
Do i = 1, n
u(i) = -cos(r8_pix(i))sin(r8_piy(i))
dudx(i) = r8_pisin(r8_pix(i))sin(r8_piy(i))
dudxx(i) = r8_pir8_picos(r8_pix(i))sin(r8_piy(i))
dudy(i) = -r8_picos(r8_pix(i))cos(r8_piy(i))
dudyy(i) = r8_pir8_picos(r8_pix(i))sin(r8_piy(i))
dudt(i) = +2.0nur8_pir8_picos(r8_pix(i))sin(r8_piy(i))

v(i) = sin(r8_pi*x(i))*cos(r8_pi*y(i))
dvdx(i) = r8_pi*cos(r8_pi*x(i))*cos(r8_pi*y(i))
dvdxx(i) = -r8_pi*r8_pi*sin(r8_pi*x(i))*cos(r8_pi*y(i))
dvdy(i) = -r8_pi*sin(r8_pi*x(i))*sin(r8_pi*y(i))
dvdyy(i) = -r8_pi*r8_pi*sin(r8_pi*x(i))*cos(r8_pi*y(i))
dvdt(i) = -2.0*nu*r8_pi*r8_pi*sin(r8_pi*x(i))*cos(r8_pi*y(i))

p(i) = -0.25D+00*rho*(cos(2.0D+00*r8_pi*x(i))+cos(2.0D+00*r8_pi*y(i)))
dpdx(i) = +0.5D+00*rho*r8_pi*sin(2.0D+00*r8_pi*x(i))
dpdy(i) = +0.5D+00*rho*r8_pi*sin(2.0D+00*r8_pi*y(i))

!
! Time scaling.
!
u(i) = u(i)exp(-2.0r8_pi2nut)
dudx(i) = dudx(i)exp(-2.0r8_pi2nut)
dudxx(i) = dudxx(i)exp(-2.0r8_pi2nut)
dudy(i) = dudy(i)exp(-2.0r8_pi2nut)
dudyy(i) = dudyy(i)exp(-2.0r8_pi2nut)
dudt(i) = dudt(i)exp(-2.0r8_pi2nut)

v(i) = v(i)*exp(-2.0*r8_pi**2*nu*t)
dvdx(i) = dvdx(i)*exp(-2.0*r8_pi**2*nu*t)
dvdxx(i) = dvdxx(i)*exp(-2.0*r8_pi**2*nu*t)
dvdy(i) = dvdy(i)*exp(-2.0*r8_pi**2*nu*t)
dvdyy(i) = dvdyy(i)*exp(-2.0*r8_pi**2*nu*t)
dvdt(i) = dvdt(i)*exp(-2.0*r8_pi**2*nu*t)

p(i) = p(i)*exp(-4.0*r8_pi**2*nu*t)
dpdx(i) = dpdx(i)*exp(-4.0*r8_pi**2*nu*t)
dpdy(i) = dpdy(i)*exp(-4.0*r8_pi**2*nu*t)

!
! Evaluate the residuals.
!
ur(i) = dudt(i) + u(i)*dudx(i) + v(i)*dudy(i) + (1.0D+00/rho)dpdx(i) -&
nu(dudxx(i)+dudyy(i)) - f(i)

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

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

End Do

Return
End Subroutine resid_taylor
Subroutine rhs_lucas(nu, rho, n, x, y, t, f, g, h)

!*****************************************************************************80
!
!c RHS_LUCAS evaluates the right hand side of Lucas Bystricky's problem.
!
! Parameters:
!
! Input, double precision NU, the kinematic viscosity.
!
! Input, double precision RHO, the 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 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
dudt(i) = 0.0D+00
dudx(i) = sin(r8_pi*x(i))
dudxx(i) = r8_pi*cos(r8_pi*x(i))
dudy(i) = 0.0D+00
dudyy(i) = 0.0D+00

v(i) = -y(i)*sin(r8_pi*x(i))
dvdt(i) = 0.0D+00
dvdx(i) = -r8_pi*y(i)*cos(r8_pi*x(i))
dvdxx(i) = +r8_pi**2*y(i)*sin(r8_pi*x(i))
dvdy(i) = -sin(r8_pi*x(i))
dvdyy(i) = 0.0D+00

p(i) = 0.0D+00
dpdx(i) = 0.0D+00
dpdy(i) = 0.0D+00

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_lucas