Skip to content

CFAST Source: fix leakage .out #2593

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
drjfloyd merged 3 commits into from
Oct 19, 2026
Merged

CFAST Source: fix leakage .out #2593

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e68daae
CFAST Source: fix leakage .out
drjfloyd Oct 18, 2026
0c069ab
Merge remote-tracking branch 'github/master'
drjfloyd Oct 18, 2026
229ecb9
CFAST Source: Define mw constants
drjfloyd Oct 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Source/CFAST/cfast_data.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ module cenviro
real(eb), parameter :: cp = 1012.0_eb
real(eb), parameter :: gamma = 1.40_eb
real(eb), parameter :: rgas = (gamma-1.0_eb)/gamma*cp
real(eb), parameter :: mw_n2 = 28.01340e-3_eb, mw_o2 = 31.99880e-3_eb, mw_co2 = 44.00950e-3_eb, mw_h2o = 18.01528e-3_eb, &
mw_c = 12.01070e-3_eb, mw_co = 28.01010e-3_eb, mw_hcn = 27.02534e-3_eb, mw_hcl = 36.46094e-3_eb

end module cenviro

Expand Down
51 changes: 25 additions & 26 deletions Source/CFAST/fire.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ module fire_routines

use cfast_types, only: detector_type, fire_type, room_type, target_type, vent_type

use cenviro, only: cp
use cenviro, only: cp, mw_n2, mw_o2, mw_co2, mw_h2o, mw_c, mw_co, mw_hcn, mw_hcl
use cparams, only: u, l, m, q, ts, fuel, n2, o2, o2f, mxfires, mxrooms, fuel_moles, fuel_q, fuel_n2, fuel_o2, fuel_co2, &
fuel_co, fuel_hcn, fuel_hcl, hcl, fuel_h2o, h2o, fuel_soot, co2, co, hcn, soot, soot_flaming, soot_smolder, ct, &
mx_hsep, t_max, ns_mass, flaming, smoldering, trigger_by_time, trigger_by_temp, trigger_by_flux, idx_tempf_trg, &
Expand Down Expand Up @@ -222,7 +222,7 @@ subroutine interpolate_pyrolysis (ifire,tsec,iroom,mdot_t,area_t,height_t,qdot_t
if (fireptr%n_hcn>0) then
call interp(fireptr%t_hcn,fireptr%y_hcn,fireptr%n_hcn,xxtime,1,y_hcn)
else
y_HCN = fireptr%n_N*0.027028_eb/fireptr%molar_mass
y_HCN = fireptr%n_N*mw_hcn/fireptr%molar_mass
end if
call interp(fireptr%t_trace,fireptr%y_trace,fireptr%n_trace,xxtime,1,y_trace)
call interp(fireptr%t_area,fireptr%area,fireptr%n_area,xxtime,1,area_t)
Expand Down Expand Up @@ -569,11 +569,11 @@ subroutine chemistry (pyrolysis_rate, molar_mass, entrainment_rate, hoc, y_soot,

! now do the chemistry balance with supplied inputs.
nu_soot = molar_mass/0.01201_eb*y_soot
nu_soot_flaming = molar_mass/0.01201_eb*y_soot_flaming
nu_soot_smolder = molar_mass/0.01201_eb*y_soot_smolder
nu_hcn = min(n_N,molar_mass/0.027028_eb*y_hcn)
nu_soot_flaming = molar_mass/mw_c*y_soot_flaming
nu_soot_smolder = molar_mass/mw_c*y_soot_smolder
nu_hcn = min(n_N,molar_mass/mw_hcn*y_hcn)
nu_hcl = n_Cl
nu_co = molar_mass/0.02801_eb*y_co
nu_co = molar_mass/mw_co*y_co
nu_h2o = (n_H - nu_hcl - nu_hcn)/2.0_eb
nu_co2 = n_C - nu_co - nu_hcn - nu_soot
nu_o2 = nu_co2 + (nu_h2o + nu_co - n_O)/2.0_eb
Expand All @@ -582,16 +582,16 @@ subroutine chemistry (pyrolysis_rate, molar_mass, entrainment_rate, hoc, y_soot,
! chemistry balance is molar-based so convert back to mass rates. fuel and o2 are consumed,
! so negative. Others are produced, so positive
net_fuel = -pyrolysis_rate_constrained
net_o2 = -pyrolysis_rate_constrained*nu_o2*0.032_eb/molar_mass
net_co2 = pyrolysis_rate_constrained*nu_co2*0.04401_eb/molar_mass
net_co = pyrolysis_rate_constrained*nu_co*0.02801_eb/molar_mass
net_h2o = pyrolysis_rate_constrained*nu_h2o*0.018016_eb/molar_mass
net_hcl = pyrolysis_rate_constrained*nu_hcl*0.036458_eb/molar_mass
net_hcn = pyrolysis_rate_constrained*nu_hcn*0.027028_eb/molar_mass
net_soot = pyrolysis_rate_constrained*nu_soot*0.01201_eb/molar_mass
net_soot_flaming = pyrolysis_rate_constrained*nu_soot_flaming*0.01201_eb/molar_mass
net_soot_smolder = pyrolysis_rate_constrained*nu_soot_smolder*0.01201_eb/molar_mass
net_n2 = pyrolysis_rate_constrained*nu_n2*0.02802_eb/molar_mass
net_o2 = -pyrolysis_rate_constrained*nu_o2*mw_o2/molar_mass
net_co2 = pyrolysis_rate_constrained*nu_co2*mw_co2/molar_mass
net_co = pyrolysis_rate_constrained*nu_co*mw_co/molar_mass
net_h2o = pyrolysis_rate_constrained*nu_h2o*mw_h2o/molar_mass
net_hcl = pyrolysis_rate_constrained*nu_hcl*mw_hcl/molar_mass
net_hcn = pyrolysis_rate_constrained*nu_hcn*mw_hcn/molar_mass
net_soot = pyrolysis_rate_constrained*nu_soot*mw_c/molar_mass
net_soot_flaming = pyrolysis_rate_constrained*nu_soot_flaming*mw_c/molar_mass
net_soot_smolder = pyrolysis_rate_constrained*nu_soot_smolder*mw_c/molar_mass
net_n2 = pyrolysis_rate_constrained*nu_n2*mw_n2/molar_mass
net_ct = 0.0_eb

if (-net_o2 > o2_available) then
Expand Down Expand Up @@ -800,13 +800,13 @@ subroutine vent_jet_fire (ito, ifrom, netfuel, entrainment_rate, hrr, species_ma
y_smolder_soot = 0.0_eb
y_co = room2ptr%species_mass(u,fuel_co)/room2ptr%species_mass(u,fuel)
y_hcn = room2ptr%species_mass(u,fuel_hcn)/room2ptr%species_mass(u,fuel)
nu_soot = y_soot*xxmolar_mass/0.01201_eb
nu_co = y_co*xxmolar_mass/0.02801_eb
nu_hcn = y_hcn*xxmolar_mass/0.027028_eb
nu_hcl = room2ptr%species_mass(u,fuel_hcl)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.036458_eb
nu_co2 = room2ptr%species_mass(u,fuel_co2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.04401_eb
nu_h2o = room2ptr%species_mass(u,fuel_h2o)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.018016_eb
nu_o2 = room2ptr%species_mass(u,fuel_o2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.032_eb
nu_soot = y_soot*xxmolar_mass/mw_c
nu_co = y_co*xxmolar_mass/mw_co
nu_hcn = y_hcn*xxmolar_mass/mw_hcn
nu_hcl = room2ptr%species_mass(u,fuel_hcl)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_hcl
nu_co2 = room2ptr%species_mass(u,fuel_co2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_co2
nu_h2o = room2ptr%species_mass(u,fuel_h2o)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_h2o
nu_o2 = room2ptr%species_mass(u,fuel_o2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_o2

n_n = nu_hcn
n_cl = nu_hcl
Expand Down Expand Up @@ -1136,11 +1136,10 @@ subroutine update_species (deltt)
real(eb), intent(in) :: deltt

! molar_masses of the species
real(eb), parameter :: molar_mass(ns_mass) = &
(/0.02802_eb,0.032_eb,0.04401_eb,0.02801_eb,0.027028_eb,0.036458_eb,0.01201_eb,0.018016_eb,0.01201_eb/)
real(eb) :: molar_mass(ns_mass) = (/mw_n2, mw_o2, mw_co2, mw_co, mw_hcn, mw_hcl, mw_c, mw_h2o, mw_c/)

! reciprocal of avagadro's number (so you can't have less than an atom of a species)
real(eb), parameter :: avagad = 1.0_eb/6.022e23_eb
real(eb), parameter :: avagad = 1.0_eb/6.02214076e23_eb

real(eb) :: air_moles(2), v(2)
integer i, layer, lsp
Expand Down
4 changes: 2 additions & 2 deletions Source/CFAST/initialization.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,7 +212,7 @@ end subroutine initialize_ambient
subroutine initialize_leakage

integer iroom, counter
real(eb), parameter :: rho_leak = 352.981915_eb/298.15_eb
real(eb), parameter :: rho_leak = default_pressure/rgas/298.15
real(eb) :: area, area2, leak_area
type(room_type), pointer :: roomptr
type(vent_type), pointer :: ventptr
Expand Down Expand Up @@ -744,7 +744,7 @@ subroutine initialize_species
! tables in the handbook of physics and chemistry.
xt = interior_ambient_temperature
xtemp = 23.2_eb - 3.816e3_eb/(xt-46.0_eb)
xh2o = exp(xtemp)/101325.0_eb*(18.016_eb/28.584_eb)
xh2o = exp(xtemp)/default_pressure*(18.016_eb/28.584_eb)
initial_mass_fraction(h2o) = relative_humidity*xh2o

! normalize the compartment atmosphere
Expand Down
32 changes: 6 additions & 26 deletions Source/CFAST/output.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -452,7 +452,7 @@ subroutine results_vent_flows
'Lower Layer Upper Layer Lower Layer',/, &
'Vent From/Bottom To/Top Fraction Inflow Outflow ', &
'Inflow Outflow Inflow Outflow Inflow Outflow',/,153('-'))
5070 format (a4,3x,a12,3x,a12,15x,a)
5070 format (a4,3x,a12,3x,a12,15x,a104)

end subroutine results_vent_flows

Expand Down Expand Up @@ -1161,7 +1161,7 @@ subroutine flwout (outbuf,flow1,flow2,flow3,flow4,flow5,flow6,flow7,flow8)
real(eb), intent(in) :: flow1, flow2, flow3, flow4, flow5, flow6, flow7, flow8
character(len=*), intent(out) :: outbuf

real :: flow(8), x1000,x100,x10,x1,x01
real :: flow(8)

integer :: i

Expand All @@ -1174,34 +1174,14 @@ subroutine flwout (outbuf,flow1,flow2,flow3,flow4,flow5,flow6,flow7,flow8)
flow(6) = real(flow6)
flow(7) = real(flow7)
flow(8) = real(flow8)
x1000 = 1000.0
x100 = 100.0
x10 = 10.0
x1 = 1.0
x01 = 0.1
do i = 1, 8
if (flow(i)>=x1000) then
write (outbuf(13*(i-1)+1:13*i),5000) flow(i)
else if (flow(i)>=x100) then
write (outbuf(13*(i-1)+1:13*i),5010) flow(i)
else if (flow(i)>=x10) then
write (outbuf(13*(i-1)+1:13*i),5020) flow(i)
else if (flow(i)>=x1) then
write (outbuf(13*(i-1)+1:13*i),5030) flow(i)
else if (flow(i)>=x01) then
write (outbuf(13*(i-1)+1:13*i),5040) flow(i)
if (flow(i)<=atol) then
outbuf(13*(i-1)+1:13*i) = ' '
else
write (outbuf(13*(i-1)+1:13*i),5000) flow(i)
end if
if (flow(i)<=atol) outbuf(13*(i-1)+1:13*i) = ' '
if (flow(i).ne.0.0_eb) write (outbuf(13*(i-1)+1:13*i),5050) flow(i)
write (outbuf(13*(i-1)+1:13*i),5050) flow(i)
endif
end do

5000 format (2x,1pg11.3)
5010 format (f6.0,7x)
5020 format (f7.1,6x)
5030 format (f8.2,5x)
5040 format (f9.3,4x)
5050 format (2x,e11.4)
end subroutine flwout

Expand Down
9 changes: 5 additions & 4 deletions Source/CFAST/radiation.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -767,6 +767,7 @@ end subroutine cross_product
! --------------------------- absorb -------------------------------------------

real(eb) function absorb(cmpt, layer)
use cenviro, only: mw_co2, mw_h2o

! calculates absorbance, due to gases (co2 and h2o) and soot, for the specified compartment and layer.

Expand Down Expand Up @@ -820,9 +821,9 @@ real(eb) function absorb(cmpt, layer)

! declare module data

! physical constants [mw in kg/mol; rg in m^3-atm/mol-kelvin]
! physical constants [rg in m^3-atm/mol-kelvin]

real(eb), parameter :: mwco2 = 44.0088e-3_eb, mwh2o = 18.0153e-3_eb, rg = 82.0562e-6_eb, k = 1195.5_eb, rhos = 1800.0_eb
real(eb), parameter :: rg = 82.0562e-6_eb, k = 1195.5_eb, rhos = 1800.0_eb

! log(t) data for co2 [t in k]

Expand Down Expand Up @@ -910,7 +911,7 @@ real(eb) function absorb(cmpt, layer)
ag = 0.0_eb

! absorbance for co2
ng = roomptr%species_mass(layer,co2)/mwco2
ng = roomptr%species_mass(layer,co2)/mw_co2
plg = ng*rtv*l
if (plg>0.0_eb) then
cplg = log10(plg)
Expand All @@ -922,7 +923,7 @@ real(eb) function absorb(cmpt, layer)
end if

! absorbance for h2o
ng = roomptr%species_mass(layer,h2o)/mwh2o
ng = roomptr%species_mass(layer,h2o)/mw_h2o
plg = ng*rtv*l
if (plg>0.0_eb) then
cplg = log10(plg)
Expand Down
17 changes: 9 additions & 8 deletions Source/CFAST/solve.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ module solve_routines
use conduction_routines, only: conduction
use convection_routines, only: convection
use debug_routines, only: output_spreadsheet_residuals
use defaults, only: default_pressure
use exit_routines, only: cfastexit, delete_output_files
use fire_routines, only: fire, vent_jets, integrate_mass, update_species, collect_fire_data_for_smokeview, update_fire_ignition
use isosurface, only: output_isodata
Expand All @@ -25,7 +26,7 @@ module solve_routines

use cfast_types, only: fire_type, room_type, target_type, vent_type

use cenviro, only: odevara, odevarb, odevarc, constvar, cp, rgas, gamma
use cenviro, only: odevara, odevarb, odevarc, constvar, cp, rgas, gamma, mw_co2, mw_h2o
use cparams, only: u, l, m, q, mxrooms, mxtarg, mxnode, mxbranch, mxdiscon, maxeq, ns, check_state, set_state, update_state, &
nwal, ns_mass, vminfrac, n2, o2, co2, co, h2o, w_from_room, w_to_room, w_from_wall, w_to_wall, w_boundary_condition, &
radiation_fix
Expand Down Expand Up @@ -1112,7 +1113,7 @@ subroutine update_data (y_vector,iflag)
! tables in the handbook of physics and chemistry.
xt = exterior_ambient_temperature
xtemp = 23.2_eb - 3.816e3_eb/(xt-46.0_eb)
xh2o = exp(xtemp)/101325.0_eb*(18.016_eb/28.584_eb)
xh2o = exp(xtemp)/default_pressure*(mw_h2o*1000._eb/28.584_eb)
roomptr%species_fraction(u,h2o) = relative_humidity*xh2o
roomptr%species_fraction(l,h2o) = relative_humidity*xh2o

Expand Down Expand Up @@ -1350,19 +1351,19 @@ subroutine update_data (y_vector,iflag)
roomptr => roominfo(iroom)
isof = isof + 1
if (radi_verification_flag) then
if (partial_pressure_co2+partial_pressure_h2o == 101325._eb) then ! Only H2O and CO2 exist
if (partial_pressure_co2+partial_pressure_h2o == default_pressure) then ! Only H2O and CO2 exist
!if (lsp == 1 .or. lsp == 2) then ! N2 and O2 will have to be zero
if (lsp == n2 .or. lsp == o2) then ! N2 and O2 will have to be zero
ppgas = 0._eb
!else if (lsp == 3) then
else if (lsp == co2) then
! partical pressure = species_mass/its molecular weight*ideal gas constant*gas temperature/volume of the medium
roomptr%species_mass(u,lsp) = &
partial_pressure_co2/101325._eb*44.0088e-3_eb/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u)
partial_pressure_co2/default_pressure*mw_co2/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u)
!else if (lsp == 8) then
else if (lsp == h2o) then
roomptr%species_mass(u,lsp) = &
partial_pressure_h2o/101325._eb*18.0153e-3_eb/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u)
partial_pressure_h2o/default_pressure*mw_h2o/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u)
else
ppgas = 0._eb
end if
Expand All @@ -1378,18 +1379,18 @@ subroutine update_data (y_vector,iflag)

isof = isof + 1
if (radi_verification_flag) then
if (partial_pressure_co2+partial_pressure_h2o == 101325._eb) then
if (partial_pressure_co2+partial_pressure_h2o == default_pressure) then
!if (lsp == 1 .or. lsp == 2) then
if (lsp == n2 .or. lsp == o2) then
ppgas = 0._eb
!else if (lsp == 3) then
else if (lsp == co2) then
roomptr%species_mass(l,lsp) = &
partial_pressure_co2/101325._eb*44.0088e-3_eb/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l)
partial_pressure_co2/default_pressure*mw_co2/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l)
!else if (lsp == 8) then
else if (lsp == h2o) then
roomptr%species_mass(l,lsp) = &
partial_pressure_h2o/101325._eb*18.0153e-3_eb/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l)
partial_pressure_h2o/default_pressure*mw_h2o/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l)
else
ppgas = 0._eb
end if
Expand Down
Loading