💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Multi-Surface Enclosure (N surfaces)
Core Numerical Engine in Fortran 90 • 36 total downloads
multisurface_enclosure.f90
! =========================================================================
! Source File: multisurface_enclosure.f90
! =========================================================================
program multisurface_enclosure
implicit none
integer :: n, i, j
real(8), allocatable :: A(:), eps(:), BC_val(:)
integer, allocatable :: BC_type(:)
real(8), allocatable :: F(:, :)
real(8), allocatable :: AM(:, :), CM(:), J_vec(:)
real(8), allocatable :: G(:), q_net(:), Q_net_val(:), Eb(:), T_C(:)
real(8) :: sigma, q_sum
logical :: error
sigma = 5.670374d-8
! Read N
read(*, *) n
allocate(A(n), eps(n), BC_type(n), BC_val(n))
allocate(F(n, n))
allocate(AM(n, n), CM(n), J_vec(n))
allocate(G(n), q_net(n), Q_net_val(n), Eb(n), T_C(n))
! Read surface properties
do i = 1, n
read(*, *) A(i), eps(i), BC_type(i), BC_val(i)
end do
! Read view factor matrix
do i = 1, n
read(*, *) F(i, :)
end do
! Construct matrix equations
do i = 1, n
if (BC_type(i) == 1) then
! Temperature imposed (BC_val is in C, convert to K)
T_C(i) = BC_val(i)
Eb(i) = sigma * (T_C(i) + 273.15d0)**4
do j = 1, n
if (i == j) then
AM(i, j) = 1.0d0 - (1.0d0 - eps(i)) * F(i, j)
else
AM(i, j) = - (1.0d0 - eps(i)) * F(i, j)
end if
end do
CM(i) = eps(i) * Eb(i)
else
! Net flux density imposed (BC_val is in W/m2)
q_net(i) = BC_val(i)
T_C(i) = 0.0d0 ! Will be calculated later
Eb(i) = 0.0d0 ! Will be calculated later
do j = 1, n
if (i == j) then
AM(i, j) = 1.0d0 - F(i, j)
else
AM(i, j) = - F(i, j)
end if
end do
CM(i) = q_net(i)
end if
end do
! Solve linear system
call solve_system(n, AM, CM, J_vec, error)
if (error) then
print *, "ERROR: Linear system is singular or poorly conditioned."
stop
end if
! Post processing
q_sum = 0.0d0
do i = 1, n
! Irradiation G(i) = sum( F(i, j) * J(j) )
G(i) = 0.0d0
do j = 1, n
G(i) = G(i) + F(i, j) * J_vec(j)
end do
if (BC_type(i) == 1) then
! Temp was imposed, calculate flux
! q_net = J - G
q_net(i) = J_vec(i) - G(i)
Q_net_val(i) = A(i) * q_net(i)
else
! Flux was imposed, calculate temperature
! q_net = J - G (should match input flux, let's use input flux for stability)
if (eps(i) > 1.0d-6) then
Eb(i) = J_vec(i) + ((1.0d0 - eps(i)) / eps(i)) * q_net(i)
if (Eb(i) > 0.0d0) then
T_C(i) = (Eb(i) / sigma)**0.25d0 - 273.15d0
else
T_C(i) = -273.15d0
end if
else
Eb(i) = 0.0d0
T_C(i) = -273.15d0 ! Absolute zero if perfect reflector
end if
Q_net_val(i) = A(i) * q_net(i)
end if
q_sum = q_sum + Q_net_val(i)
end do
! Output results in key-value format
print *, "N=", n
do i = 1, n
print '(A,I0,A,F20.6)', "SURF_", i, "_AREA=", A(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_EPS=", eps(i)
print '(A,I0,A,I0)', "SURF_", i, "_BCTYPE=", BC_type(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_T=", T_C(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_QNET_DEN=", q_net(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_QNET=", Q_net_val(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_J=", J_vec(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_G=", G(i)
print '(A,I0,A,F20.6)', "SURF_", i, "_EB=", Eb(i)
end do
print '(A,F20.8)', "Q_SUM=", q_sum
contains
subroutine solve_system(n, A_mat, C_vec, X_vec, err)
integer, intent(in) :: n
real(8), intent(inout) :: A_mat(n, n)
real(8), intent(inout) :: C_vec(n)
real(8), intent(out) :: X_vec(n)
logical, intent(out) :: err
integer :: i, j, k, max_row
real(8) :: factor, pivot, temp_val
real(8) :: temp_row(n)
err = .false.
X_vec = 0.0d0
do i = 1, n
! Find pivot with partial pivoting
max_row = i
do k = i + 1, n
if (abs(A_mat(k, i)) > abs(A_mat(max_row, i))) then
max_row = k
end if
end do
! Swap rows if necessary
if (max_row /= i) then
temp_row = A_mat(i, :)
A_mat(i, :) = A_mat(max_row, :)
A_mat(max_row, :) = temp_row
temp_val = C_vec(i)
C_vec(i) = C_vec(max_row)
C_vec(max_row) = temp_val
end if
pivot = A_mat(i, i)
if (abs(pivot) < 1.0d-12) then
err = .true.
return
end if
! Normalize row i
C_vec(i) = C_vec(i) / pivot
A_mat(i, :) = A_mat(i, :) / pivot
! Eliminate other rows
do j = 1, n
if (j /= i) then
factor = A_mat(j, i)
A_mat(j, :) = A_mat(j, :) - factor * A_mat(i, :)
C_vec(j) = C_vec(j) - factor * C_vec(i)
end if
end do
end do
X_vec = C_vec
end subroutine solve_system
end program multisurface_enclosure
Solver Description
Solves radiation heat transfer within a closed enclosure of N gray, diffuse, and opaque surfaces using matrix radiosity.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 multisurface_enclosure.f90 -o multisurface_enclosure
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
multisurface_enclosure < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! N
3
! Area_i Emissivity_i BC_type_i BC_val_i
1.0 0.8 1 500.0
! Form factor matrix row by row
2.0 0.5 1 50.0
! Parameter 4
1.5 0.2 2 0.0
! Parameter 5
0.0 1.0 0.0
! Parameter 6
0.5 0.1 0.4
! Parameter 7
0.0 0.53333333 0.46666667
3
! Area_i Emissivity_i BC_type_i BC_val_i
1.0 0.8 1 500.0
! Form factor matrix row by row
2.0 0.5 1 50.0
! Parameter 4
1.5 0.2 2 0.0
! Parameter 5
0.0 1.0 0.0
! Parameter 6
0.5 0.1 0.4
! Parameter 7
0.0 0.53333333 0.46666667