🔬
Solver Purpose & Physical Scope
Simulate gradually varied open-channel water surface profiles (M1, M2, S1, S2 etc.) using Heun integration of the GVF equation.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 473 times
📄 Source File:
gvf_profile.f90
📁 calcul/FluidMechanics /
gvf_profile.f90
program gvf_profile
implicit none
integer :: i, n_steps, channelType, iostat_val, profileClass
double precision :: Q, b, S0, n_mann, z_side, y_start, x_length, g
double precision :: yc, yn, A, P_wet, Rh, Sf, Fr2, dydx
double precision :: x, y, dx, y_new, A2, P2, Rh2, Sf2, Fr2b, dydx2
double precision :: Ac, Pc, Rhc, An, Pn, Rhn
character(len=40) :: className, channelName
read(*,*,iostat=iostat_val) channelType
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid channel type input.'
stop
end if
read(*,*,iostat=iostat_val) Q
read(*,*,iostat=iostat_val) b
read(*,*,iostat=iostat_val) S0
read(*,*,iostat=iostat_val) n_mann
read(*,*,iostat=iostat_val) z_side
read(*,*,iostat=iostat_val) y_start
read(*,*,iostat=iostat_val) x_length
read(*,*,iostat=iostat_val) g
read(*,*,iostat=iostat_val) n_steps
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all GVF inputs.'
stop
end if
if (Q <= 0.0d0 .or. b <= 0.0d0 .or. n_mann <= 0.0d0) then
write(*,*) 'ERROR: Q, b, and Manning n must be positive.'
stop
end if
if (y_start <= 0.0d0) then
write(*,*) 'ERROR: Starting depth must be positive.'
stop
end if
if (g <= 0.0d0) g = 9.81d0
if (n_steps < 10) n_steps = 500
if (x_length <= 0.0d0) x_length = 1000.0d0
select case(channelType)
case(1)
channelName = 'Rectangular'
case(2)
channelName = 'Trapezoidal'
case default
channelName = 'Rectangular'
channelType = 1
end select
! Critical depth by bisection
yc = find_critical_depth(Q, b, z_side, g, channelType)
! Normal depth by bisection
yn = find_normal_depth(Q, b, z_side, S0, n_mann, channelType)
! Classify profile
call classify_profile(S0, y_start, yn, yc, profileClass, className)
write(*,'(A)') '============================================================'
write(*,'(A)') ' GRADUALLY VARIED FLOW (GVF) PROFILE ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Channel Shape = ', trim(channelName)
write(*,'(A,ES12.4,A)') ' Discharge Q = ', Q, ' m3/s'
write(*,'(A,ES12.4,A)') ' Bottom Width b = ', b, ' m'
write(*,'(A,ES12.4)') ' Bed Slope S0 = ', S0
write(*,'(A,ES12.4)') ' Manning n = ', n_mann
write(*,'(A,ES12.4)') ' Side Slope z = ', z_side
write(*,'(A,ES12.4,A)') ' Starting Depth y_start = ', y_start, ' m'
write(*,'(A,ES12.4,A)') ' Reach Length = ', x_length, ' m'
write(*,'(A,I8)') ' Number of Steps = ', n_steps
write(*,*)
write(*,'(A)') '--- CRITICAL / NORMAL DEPTHS --------------------------------'
write(*,'(A,ES12.4,A)') ' Critical Depth yc = ', yc, ' m'
write(*,'(A,ES12.4,A)') ' Normal Depth yn = ', yn, ' m'
write(*,'(A,A)') ' Profile Classification = ', trim(className)
write(*,*)
! Compute specific energy and Froude at critical depth
call channel_geometry(channelType, yc, b, z_side, Ac, Pc, Rhc)
call channel_geometry(channelType, yn, b, z_side, An, Pn, Rhn)
write(*,'(A)') '--- HYDRAULIC DATA ------------------------------------------'
write(*,'(A,ES12.4,A)') ' Critical Area Ac = ', Ac, ' m2'
write(*,'(A,ES12.4,A)') ' Normal Area An = ', An, ' m2'
write(*,'(A,ES12.4,A)') ' Critical Rh = ', Rhc, ' m'
write(*,'(A,ES12.4,A)') ' Normal Rh = ', Rhn, ' m'
write(*,*)
! Euler integration of dy/dx = (S0 - Sf)/(1 - Fr^2)
dx = x_length / dble(n_steps)
x = 0.0d0
y = y_start
write(*,'(A)') '--- WATER SURFACE PROFILE -----------------------------------'
write(*,'(A)') ' x[m] y[m] Sf Fr E[m]'
write(*,'(A)') ' ----------------------------------------------------------------'
do i = 0, n_steps
call channel_geometry(channelType, y, b, z_side, A, P_wet, Rh)
if (A > 0.0d0) then
Sf = (Q * n_mann / (A * Rh**(2.0d0/3.0d0)))**2
Fr2 = Q**2 * (b + 2.0d0*z_side*y) / (g * A**3)
else
Sf = S0
Fr2 = 0.0d0
end if
write(*,'(F12.3,2X,F12.6,2X,ES12.4,2X,F10.5,2X,F12.6)') &
x, y, Sf, sqrt(max(Fr2,0.0d0)), &
y + Q**2/(2.0d0*g*max(A,1.0d-30)**2)
if (i == n_steps) exit
! Euler predictor
if (abs(1.0d0 - Fr2) > 1.0d-8) then
dydx = (S0 - Sf) / (1.0d0 - Fr2)
else
dydx = 0.0d0
end if
! Heun corrector
y_new = y + dydx * dx
if (y_new <= 0.001d0) y_new = 0.001d0
call channel_geometry(channelType, y_new, b, z_side, A2, P2, Rh2)
if (A2 > 0.0d0) then
Sf2 = (Q * n_mann / (A2 * Rh2**(2.0d0/3.0d0)))**2
Fr2b = Q**2 * (b + 2.0d0*z_side*y_new) / (g * A2**3)
else
Sf2 = S0
Fr2b = 0.0d0
end if
if (abs(1.0d0 - Fr2b) > 1.0d-8) then
dydx2 = (S0 - Sf2) / (1.0d0 - Fr2b)
else
dydx2 = 0.0d0
end if
y = y + 0.5d0*(dydx + dydx2)*dx
if (y <= 0.001d0) y = 0.001d0
x = x + dx
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' dy/dx = (S0 - Sf) / (1 - Fr^2).'
write(*,'(A)') ' Sf = [Q n / (A Rh^(2/3))]^2 (Manning).'
write(*,'(A)') ' Fr^2 = Q^2 T / (g A^3); T = b + 2 z y.'
write(*,'(A)') ' Heun (improved Euler) integration scheme.'
contains
subroutine channel_geometry(ctype, depth, bw, zs, Aout, Pout, Rhout)
implicit none
integer, intent(in) :: ctype
double precision, intent(in) :: depth, bw, zs
double precision, intent(out) :: Aout, Pout, Rhout
if (ctype == 2) then
Aout = (bw + zs*depth)*depth
Pout = bw + 2.0d0*depth*sqrt(1.0d0 + zs**2)
else
Aout = bw * depth
Pout = bw + 2.0d0*depth
end if
if (Pout > 0.0d0) then
Rhout = Aout / Pout
else
Rhout = 0.0d0
end if
end subroutine channel_geometry
double precision function find_critical_depth(Qin, bw, zs, gin, ctype)
implicit none
double precision, intent(in) :: Qin, bw, zs, gin
integer, intent(in) :: ctype
double precision :: lo, hi, mid, Amid, Tmid, fval
integer :: it
lo = 0.0001d0
hi = 50.0d0
do it = 1, 120
mid = 0.5d0*(lo+hi)
if (ctype == 2) then
Amid = (bw + zs*mid)*mid
Tmid = bw + 2.0d0*zs*mid
else
Amid = bw*mid
Tmid = bw
end if
fval = Qin**2*Tmid - gin*Amid**3
if (fval > 0.0d0) then
lo = mid
else
hi = mid
end if
end do
find_critical_depth = 0.5d0*(lo+hi)
end function find_critical_depth
double precision function find_normal_depth(Qin, bw, zs, S0in, nin, ctype)
implicit none
double precision, intent(in) :: Qin, bw, zs, S0in, nin
integer, intent(in) :: ctype
double precision :: lo, hi, mid, Amid, Pmid, Rhmid, Qmann, fval
integer :: it
if (S0in <= 0.0d0) then
find_normal_depth = 1.0d30
return
end if
lo = 0.0001d0
hi = 50.0d0
do it = 1, 120
mid = 0.5d0*(lo+hi)
if (ctype == 2) then
Amid = (bw + zs*mid)*mid
Pmid = bw + 2.0d0*mid*sqrt(1.0d0 + zs**2)
else
Amid = bw*mid
Pmid = bw + 2.0d0*mid
end if
Rhmid = Amid / max(Pmid, 1.0d-30)
Qmann = (1.0d0/nin)*Amid*Rhmid**(2.0d0/3.0d0)*sqrt(S0in)
fval = Qmann - Qin
if (fval < 0.0d0) then
lo = mid
else
hi = mid
end if
end do
find_normal_depth = 0.5d0*(lo+hi)
end function find_normal_depth
subroutine classify_profile(S0in, ystart, ynorm, ycrit, pclass, pname)
implicit none
double precision, intent(in) :: S0in, ystart, ynorm, ycrit
integer, intent(out) :: pclass
character(len=40), intent(out) :: pname
if (S0in <= 0.0d0) then
! Adverse or horizontal
if (ystart > ycrit) then
pclass = 10; pname = 'A2 or H2 - Adverse/Horizontal'
else
pclass = 11; pname = 'A3 or H3 - Adverse/Horizontal'
end if
else if (ynorm > ycrit) then
! Mild slope
if (ystart > ynorm) then
pclass = 1; pname = 'M1 - Mild backwater'
else if (ystart > ycrit) then
pclass = 2; pname = 'M2 - Mild drawdown'
else
pclass = 3; pname = 'M3 - Mild supercritical'
end if
else if (abs(ynorm-ycrit) < 1.0d-6) then
pclass = 7; pname = 'C1/C3 - Critical slope'
else
! Steep slope
if (ystart > ycrit) then
pclass = 4; pname = 'S1 - Steep backwater'
else if (ystart > ynorm) then
pclass = 5; pname = 'S2 - Steep drawdown'
else
pclass = 6; pname = 'S3 - Steep supercritical'
end if
end if
end subroutine classify_profile
end program gvf_profile
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 gvf_profile.f90 -o gvf_profile
2. Execution with input.txt redirection:
gvf_profile < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
ct_init\nDischarge Q [m³/s]\nBottom width b [m]\nBed slope S₀\nManning n\nSide slope z (H:V, 0 = rectangular)\nys_init\nxl_init\nGravity g [m/s²]\nManning n