From 96e96470f1fdd0af568057a84edd159fadee1ff0 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Mon, 4 Oct 2021 10:23:05 +1300 Subject: [PATCH 01/30] initial framework for lymphatics code --- src/bindings/c/CMakeLists.txt | 3 + src/bindings/c/lymphatics.c | 9 + src/bindings/c/lymphatics.f90 | 28 +++ src/bindings/c/lymphatics.h | 8 + src/bindings/interface/lymphatics.i | 10 + src/bindings/python/CMakeLists.txt | 1 + src/lib/CMakeLists.txt | 1 + src/lib/lymphatics.f90 | 330 ++++++++++++++++++++++++++++ 8 files changed, 390 insertions(+) create mode 100644 src/bindings/c/lymphatics.c create mode 100644 src/bindings/c/lymphatics.f90 create mode 100644 src/bindings/c/lymphatics.h create mode 100644 src/bindings/interface/lymphatics.i create mode 100644 src/lib/lymphatics.f90 diff --git a/src/bindings/c/CMakeLists.txt b/src/bindings/c/CMakeLists.txt index e2e2ff32..a723ad17 100644 --- a/src/bindings/c/CMakeLists.txt +++ b/src/bindings/c/CMakeLists.txt @@ -13,6 +13,7 @@ set(C_FORTRAN_LIB_SRCS growtree.f90 indices.f90 imports.f90 + lymphatics.f90 pressure_resistance_flow.f90 species_transport.f90 utils.f90 @@ -32,6 +33,7 @@ set(C_C_LIB_SRCS growtree.c indices.c imports.c + lymphatics.c pressure_resistance_flow.c species_transport.c utils.c @@ -52,6 +54,7 @@ set(C_LIB_HDRS growtree.h indices.h imports.h + lymphatics.h pressure_resistance_flow.h species_transport.h utils.h diff --git a/src/bindings/c/lymphatics.c b/src/bindings/c/lymphatics.c new file mode 100644 index 00000000..01149745 --- /dev/null +++ b/src/bindings/c/lymphatics.c @@ -0,0 +1,9 @@ +#include "lymphatics.h" +#include "string.h" + +void alveolar_capillary_flux_c(int *num_nodes); + +void alveolar_capillary_flux(int num_nodes) +{ + alveolar_capillary_flux_c(&num_nodes); +} diff --git a/src/bindings/c/lymphatics.f90 b/src/bindings/c/lymphatics.f90 new file mode 100644 index 00000000..7740974f --- /dev/null +++ b/src/bindings/c/lymphatics.f90 @@ -0,0 +1,28 @@ +module lymphatics_c + implicit none + + private + +contains +! +!################################################################################### +! +!*alveolar_capillary_flux:* + subroutine alveolar_capillary_flux_c(num_nodes) bind(C, name="alveolar_capillary_flux_c") + use lymphatics,only: alveolar_capillary_flux + implicit none + + integer,intent(in) :: num_nodes + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_alveolar_capillary_flux(num_nodes) +#else + call alveolar_capillary_flux(num_nodes) +#endif + + end subroutine alveolar_capillary_flux_c +! +!################################################################################### +! + +end module lymphatics_c diff --git a/src/bindings/c/lymphatics.h b/src/bindings/c/lymphatics.h new file mode 100644 index 00000000..a0bf39c2 --- /dev/null +++ b/src/bindings/c/lymphatics.h @@ -0,0 +1,8 @@ +#ifndef AETHER_LYMPHATICS_H +#define AETHER_LYMPHATICS_H + +#include "symbol_export.h" + +SHO_PUBLIC void alveolar_capillary_flux(int num_nodes); + +#endif /* AETHER_LYMPHATICS_H */ diff --git a/src/bindings/interface/lymphatics.i b/src/bindings/interface/lymphatics.i new file mode 100644 index 00000000..99363875 --- /dev/null +++ b/src/bindings/interface/lymphatics.i @@ -0,0 +1,10 @@ + +%module(package="aether") lymphatics +%include symbol_export.h +%include lymphatics.h + +%{ +#include "lymphatics.h" +%} + + diff --git a/src/bindings/python/CMakeLists.txt b/src/bindings/python/CMakeLists.txt index 5d72df0f..74fec668 100644 --- a/src/bindings/python/CMakeLists.txt +++ b/src/bindings/python/CMakeLists.txt @@ -29,6 +29,7 @@ set(INTERFACE_SRCS ../interface/growtree.i ../interface/indices.i ../interface/imports.i + ../interface/lymphatics.i ../interface/pressure_resistance_flow.i ../interface/species_transport.i ../interface/ventilation.i diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt index 6ff55d82..b434172f 100644 --- a/src/lib/CMakeLists.txt +++ b/src/lib/CMakeLists.txt @@ -13,6 +13,7 @@ set(LIB_SRCS growtree.f90 indices.f90 imports.f90 + lymphatics.f90 math_utilities.f90 mesh_utilities.f90 other_consts.f90 diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 new file mode 100644 index 00000000..cdcc444a --- /dev/null +++ b/src/lib/lymphatics.f90 @@ -0,0 +1,330 @@ +module lymphatics + !*Brief Description:* This module contains all lymphatic-specific subroutines + ! + !*LICENSE:* + ! + ! + ! + !*Full Description:* + ! + !This module contains code for pulmonary fluid flux within the alveolo-capillary network, + !and lymph transport through lymphatic collecting vessels. + + use arrays + use diagnostics + use indices + use precision ! sets dp for precision + + implicit none + + !Module parameters + + !Module types + + !Module variables + + !Interfaces +! private + public alveolar_capillary_flux +contains + +!!!############################################################################# + + subroutine alveolar_capillary_flux(num_nodes) + !*alveolar_capillary_flux:* calculate fluid flux from blood to interstitium + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ALVEOLAR_CAPILLARY_FLUX" :: ALVEOLAR_CAPILLARY_FLUX + + use other_consts,only: pi + + integer,intent(in) :: num_nodes + + ! Baseline value parameters (eventually will be user-defined?) + integer,parameter :: sex = 0 ! 0 = male, 1 = female + integer,parameter :: au = 32768 ! number of acinar units + real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp + real(dp),parameter :: lymphatic_density = 1.0_dp !to be calculated from CT?? + + ! Capillary parameters + real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp ! #0.809e-7#0.13256e-7*133/20#/20 #0.1278# obtained from literature - Parker range in cm H2O + real(dp),parameter :: open_capillaries = 1.0_dp/6.0_dp + real(dp) :: collecting_volume = 0.0_dp + real(dp) :: collecting_lymphatic_pressure = 5.0_dp + real(dp),parameter :: reflection_coefficient = 0.0_dp + real(dp),parameter :: capillary_SA = 1022.0_dp*0.599998685_dp !#0.970116# ! 1022 standardises capillary to acinar volumes as one interstitium is there for al 1022 capillaries + real(dp),parameter :: transit_time = 1.315728977_dp !#2.813324# + real(dp),parameter :: capillary_pressure = 11.51694179_dp !#23.37539# + + ! Interstitial parameters + real(dp),parameter :: interstitial_resistance = 38000.0_dp !Bertram 2011 (600-12000000) dyn/cm^2 + real(dp),parameter :: diffusion_constant = 0.00000015_dp + + ! Initial lymphatics parameters + real(dp),parameter :: lymphatic_resistance = 1000.0_dp ! Ngo 2019 + real(dp),parameter :: valve_resistance = 0.020_dp ! pressure requirement to overcome the valve resistance and cause flow OR the pressure lost during open the valve + real(dp),parameter :: lymphatic_integrity = 1.0_dp ! a measure of how 'leaky' the lymphatic vessels are and prone to backflow + + ! Osmotic pressure parameters + real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! mol/L in blood plasma + real(dp),parameter :: IGC = 62.3630_dp ! ideal gas constant in mmHg.L.mol-1.K-1 + real(dp),parameter :: T = 310.0_dp ! temperature in Kelvin - based on 37C blood temp + + ! Simulation parameters + real(dp),parameter :: breathing_rate = 15.0_dp + + + ! Local variables + integer :: i,liflowcount,printcount,time_loop + real(dp) :: alveolar_volume,breathing_function,capillary_flow,capillary_vps,capillary_osmotic,capillary_osm_n, & + capillary_volume,capillary_volume_raw,cap_osm_conc,diffusion,excess,flux_a,flux_b,flux_c,gas_diffusion_restriction, & + initial_lymphatic_flow,initial_lymphatic_pressure,initial_lymphatic_surface_area,initial_lymphatic_volume, & + initial_lymph_conc,initial_osm_n,interstitial_capacity,interstitial_capacity_a,interstitial_capacity_b, & + interstitial_osmotic,interstitial_pressure_a,interstitial_pressure_b,interstitial_saturation,interstitial_volume, & + int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,net_flux, & + osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time2,time_dt,time_period,time_variable,total_flux, & + total_hydro_flux,total_osm_flux + + character(len=60) :: sub_name + + ! -------------------------------------------------------------------------- + + sub_name = 'alveolar_capillary_flux' + call enter_exit(sub_name,1) + + ! Calculated values + lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g + capillary_volume_raw = body_mass/0.3474_dp ! Gehr 213 ml and body mass of 74 kg - rough estimate? + capillary_volume = (capillary_volume_raw*open_capillaries)/real(au) + capillary_flow = (0.085_dp/60.0_dp) ! uL/sec per unit - will obtain from perfusion model? + + ! interstitial values + interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(au))*1000.0_dp ! maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) + interstitial_capacity_a = 0.005_dp*interstitial_capacity + interstitial_capacity_b = 0.995_dp*interstitial_capacity + interstitial_volume_a = 0.0_dp + interstitial_volume_b = 0.49_dp*interstitial_capacity + alveolar_volume = 0.0_dp + liflowcount = 0 + + ! initial lymphatic values + initial_lymphatic_volume = 0.0_dp + + ! Osmotic pressures + i = 1.0_dp + capillary_osmotic = i*capillary_molar_conc*IGC*T + capillary_osm_n = 0.0_dp + interstitial_osmotic = 0.0_dp + int_osm_n = 0.0_dp + initial_osm_n = 0.0_dp + total_osm_flux = 0.0_dp + + ! time is the time that each blood particle is in the capillaries + initial_lymphatic_surface_area = capillary_SA*3.648_dp ! #0.00014046 #to be normalised to volume of lung? + + write(*,'('' blood transit time is '',f12.2,'' s'')') transit_time + ! second is how much blood flows through the capillary each second + + capillary_vps = 1.0_dp / transit_time + write(*,'('' volume per second is '',f12.2,'' mL'')') capillary_vps + + time = 0.0_dp + time2 = 1.0_dp + printcount = 0 + total_hydro_flux = 0.0_dp + time_loop = 96 + !times = [] + !int_flux = [] + !interstitial_saturation_as = [] + !initial_lymph_flux = [] + !initial_lymph_sat = [] + !total_flux_list = [] + !total_uptake_list = [] + !collected_volume = [] + !flux_difference = [] + sumflux = 0.0_dp + sumuptake = 0.0_dp + + breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) + + write(*,'('' breathing rate is '',f8.2,'' breaths per minute'')') breathing_rate + write(*,'('' lung mass is '',f8.2)') lung_mass + write(*,'('' cap SA is '',f8.2)') capillary_SA + write(*,'('' int cap is '',f8.2)') interstitial_capacity + write(*,'('' int capA is '',f8.2)') interstitial_capacity_a + write(*,'('' int capB is '',f8.2)') interstitial_capacity_b + write(*,'('' int vol is '',f8.2)') interstitial_volume_b + write(*,'('' lymphSA is '',f8.2)') initial_lymphatic_surface_area + write(*,'('' int cap is '',f8.2)') interstitial_capacity + + write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & + &''init lymph|init lymph| total| alveolar|'')') + write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & + &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') + + do while(time < 32400.0_dp) + time = time + transit_time + do i = 0,95 + interstitial_volume = interstitial_volume_a + interstitial_volume_b + interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% + time_variable = time-transit_time + (transit_time*(time2/time_loop)) + ! calculating flux from capillary into interstitium + interstitial_pressure_a = 1.47_dp * sin(time_variable * breathing_function) + & + ((-3.98_dp * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & + 8.03_dp * (interstitial_volume_a / interstitial_capacity_a) - 6.52_dp) + interstitial_pressure_b = 1.47_dp * sin(time_variable * breathing_function) + & + ((-3.98_dp * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & + 8.03_dp * (interstitial_volume_b / interstitial_capacity_b) - 6.52_dp) + ! pressure determined from saturation equation based of literature (currently linear, but likely not) + if(capillary_pressure > interstitial_pressure_a)then + flux_a = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & + interstitial_pressure_a)) * (transit_time/time_loop) + else + flux_a = 0.0_dp + endif + if(capillary_pressure > interstitial_pressure_b)then + flux_b = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & + interstitial_pressure_b)) * (transit_time / time_loop) + else + flux_b = 0.0_dp + endif + flux_c = flux_a + flux_b + total_hydro_flux = total_hydro_flux + flux_c + + if(interstitial_volume_a + flux_a > interstitial_capacity_a)then + excess = flux_a - (interstitial_capacity_a - interstitial_volume_a) + interstitial_volume_a = interstitial_capacity_a + alveolar_volume = alveolar_volume + 0.5_dp*excess + + if((interstitial_volume_b + 0.5_dp*excess) > interstitial_capacity_b)then + overflow = 0.5_dp * excess - (interstitial_capacity_b - interstitial_volume_b) + alveolar_volume = alveolar_volume + overflow + interstitial_volume_b = interstitial_capacity_b + else + interstitial_volume_b = interstitial_volume_b + 0.5*excess + endif + else + interstitial_volume_a = interstitial_volume_a + flux_a + endif + + interstitial_volume_b = interstitial_volume_b + flux_b + !diffusion = (interstitial_volume_b/(interstitial_capacity-interstitial_capacity_a))*-diffusion_constant + diffusion_constant + diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & + (160_dp*1.1_dp))*(transit_time/time_loop) + interstitial_volume_b = interstitial_volume_b + diffusion + interstitial_volume_a = interstitial_volume_a - diffusion + !alveolar_volume = alveolar_volume - 4.45e-5 + + ! Osmotic + int_osm_conc = int_osm_n/interstitial_volume + interstitial_osmotic = i*int_osm_conc*IGC*T + osm_flux = (reflection_coefficient * capillary_SA * (capillary_osmotic - interstitial_osmotic)) & + * (transit_time / time_loop) + cap_osm_conc = capillary_osm_n / capillary_volume + int_osm_conc = int_osm_n / interstitial_volume + + if(capillary_osmotic > interstitial_osmotic)then + capillary_volume = capillary_volume - osm_flux + interstitial_volume_b = interstitial_volume_b + osm_flux + else + capillary_volume = capillary_volume + osm_flux + interstitial_volume_b = interstitial_volume_b - osm_flux + endif + + net_flux = osm_flux + flux_a + flux_b + if(net_flux > 0.0_dp)then + osm_n_flux = net_flux * cap_osm_conc + int_osm_n = int_osm_n + osm_n_flux + !assuming capillary is constant and does not need updating + else + osm_n_flux = net_flux * int_osm_conc + int_osm_n = int_osm_n - osm_n_flux + endif + + total_osm_flux = total_osm_flux + osm_flux + + !calculating flux from interstitium to initial lymphatics + if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then + lymphatic_conductivity = 0.9_dp * capillary_conductivity + else + lymphatic_conductivity = ((-1625.1_dp * (interstitial_volume_b / interstitial_capacity_b)**5_dp) + & + (3815.1_dp * (interstitial_volume_b / interstitial_capacity_b)**4_dp) + (-3229.5_dp * & + (interstitial_volume_b / interstitial_capacity_b)**3_dp) + (1258_dp * (interstitial_volume_b / & + interstitial_capacity_b)**2_dp) + (-213.23_dp * (interstitial_volume_b / interstitial_capacity_b)) & + + 11.812_dp)* 4.41335e-8_dp !(capillary_conductivity) + endif + !initial_lymphatic_pressure = (math.sin(time_variable * breathing_function + math.pi)) + ((7 * (interstitial_volume_b / interstitial_capacity_b)) - 7) + initial_lymphatic_pressure = ((1.47_dp * sin(time_variable * breathing_function + pi/2_dp)) + & + ((6.82_dp* (interstitial_volume_b / interstitial_capacity_b)**2_dp) + (0.77_dp * (interstitial_volume_b / & + interstitial_capacity_b)) - 6.52_dp)) + ! initial_lymphatic_pressure = math.sin((breathing_function * time_variable) + (math.pi/2)) + ((2 * ((interstitial_volume_b/interstitial_capacity_b) ** 2) +3 * (interstitial_volume_b/interstitial_capacity_b) - 7)+(2*(interstitial_volume_b/interstitial_capacity_b))) ! based on range of values + if(interstitial_volume.le.0.0_dp)then + initial_lymphatic_flow = 0.0_dp + interstitial_volume = 0.0_dp + elseif (interstitial_pressure_b > initial_lymphatic_pressure)then + if((lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b-initial_lymphatic_pressure)) & + * (transit_time/time_loop) > ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time)) then + initial_lymphatic_flow = ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time) + else + initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - & + initial_lymphatic_pressure)) * (transit_time/time_loop) + !write(*,*) (lymphatic_conductivity, initial_lymphatic_surface_area, interstitial_pressure, initial_lymphatic_pressure, transit_time, time_loop) + endif + else + initial_lymphatic_flow = 0.0_dp + endif + + interstitial_volume_b = interstitial_volume_b - initial_lymphatic_flow + initial_lymphatic_volume = initial_lymphatic_volume + initial_lymphatic_flow + + int_osm_conc = int_osm_n/interstitial_volume_b + liflowcount = liflowcount + initial_lymphatic_flow + initial_osm_n = initial_osm_n + (initial_lymphatic_flow*int_osm_conc) + int_osm_n = int_osm_n - (initial_lymphatic_flow*int_osm_conc) + if (initial_lymphatic_volume > 0)then + initial_lymph_conc = initial_osm_n/initial_lymphatic_volume + else + initial_lymph_conc = 0.0_dp + endif + + gas_diffusion_restriction = 0.0000152587890625_dp * exp(13.8629436112_dp*(interstitial_saturation/100.0_dp)) + + time2 = time2 + 1 + enddo + + time2 = 1 + + ! if (interstitial_pressure_b > -0.95 or interstitial_pressure_b < -8.05): + ! raise ValueError('interstitial pressure_b outside of range: ' + str(interstitial_pressure_b)) + ! if initial_lymphatic_pressure > 1.05 or initial_lymphatic_pressure < -8.05: + ! raise ValueError('initial lymphatic pressure outside of range: '+ str(initial_lymphatic_pressure) + ' int pressure: ' + str(interstitial_pressure_b)) + + total_flux = total_hydro_flux ! +total_osm_flux + sumuptake = sumuptake + initial_lymphatic_flow + printcount = printcount + 1 + if (printcount.eq.100)then + write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & + flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + initial_lymphatic_volume,total_flux,alveolar_volume +! write(*,*) 'at ',time, ' s:' +! write(*,*) ' flux per second is ',flux_c/transit_time*1000,' ul' +! write(*,*) ' interstitial volume is ',interstitial_volume,' mL' +! write(*,*) ' interstitial volume a is ',interstitial_volume_a,' mL' +! write(*,*) ' interstitial volume b is ',interstitial_volume_b,' mL' +! write(*,*) ' interstitial saturation is ',100*interstitial_saturation,'%' +! ! write(*,*) ' diffusion is ' + str("%.10f" % diffusion) + ' mL') +! ! write(*,*) ' interstitial pressure is ' + str("%.2f" % interstitial_pressure) + ' mmHg') +! write(*,*) ' initial lymphatic flux is ',initial_lymphatic_flow*1000,' uL' +! write(*,*) ' initial lymphatic volume is ',initial_lymphatic_volume,' mL' +! write(*,*) ' total flux is ',total_flux,' mL' +! write(*,*) ' alveolar volume is ',alveolar_volume + printcount = 0 + endif + enddo ! while + + call enter_exit(sub_name,2) + + end subroutine alveolar_capillary_flux + +!!!############################################################################# + +end module lymphatics From b965956fcafe15920b0de94290b89b126c6f805e Mon Sep 17 00:00:00 2001 From: EdeTede Date: Wed, 13 Oct 2021 13:58:01 +1300 Subject: [PATCH 02/30] test --- src/lib/lymphatics.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index cdcc444a..a44d8180 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -3,7 +3,7 @@ module lymphatics ! !*LICENSE:* ! - ! + !Test Test Test ! !*Full Description:* ! From dbaf8e00bf554fa43843286d0523eba15babebce Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Thu, 14 Oct 2021 10:08:18 +1300 Subject: [PATCH 03/30] lymphatics_transport subroutine and bindings, to evaluate whole system transport. unit_field used to transfer transit time, surface area, and average capillary pressure for lymphatics calculation --- src/lib/capillaryflow.f90 | 14 +++++++--- src/lib/geometry.f90 | 2 ++ src/lib/indices.f90 | 15 ++++++----- src/lib/lymphatics.f90 | 56 ++++++++++++++++++++++++++++++++------- 4 files changed, 69 insertions(+), 18 deletions(-) diff --git a/src/lib/capillaryflow.f90 b/src/lib/capillaryflow.f90 index ed9bf583..bf970f6b 100644 --- a/src/lib/capillaryflow.f90 +++ b/src/lib/capillaryflow.f90 @@ -136,7 +136,7 @@ subroutine evaluate_ladder(ne,NonZeros,MatrixSize,submatrixsize,ngen,& real(dp) :: area_scale,length_scale,alpha_c ! Local variables - integer :: i,iter,j,gen,zone,num_sheet + integer :: i,iter,j,gen,zone,num_sheet,nunit integer, allocatable :: SparseCol(:) integer, allocatable ::SparseRow(:) real(dp) :: area_new,ErrorEstimate,Hart,Hven,Pin_SHEET,Pout_SHEET @@ -332,9 +332,17 @@ subroutine evaluate_ladder(ne,NonZeros,MatrixSize,submatrixsize,ngen,& ! Rtot=9 Pa/mm^3 | Blood_vol=10 mm^3| sheet_area= 11 mm^2 | ave_TT=12 s |ave_H=13 um |Ppl=14 Pa WRITE(20,& '(I6,X,5(F9.2,X),2(F8.5,X),F10.2,X,F8.4,X,F10.4,X,F10.3,X,F8.4,X,F9.4,X)') & - ne,x,y,z,Pin,Pout,Q01_mthrees*1.d9,Qtot*1.d9,Rtot/1000.d0**3.d0,& + ne,x,y,z,Pin,Pout,Q01_mthrees*1.0e+9_dp,Qtot*1.0e+9_dp,Rtot/1000.0_dp**3.0_dp, & TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL,TOTAL_SHEET_H,Ppl - ENDIF + + ! record the unit values for mean pressure, transit time, surface area. + ! ne is the 'linker' element, so nunit is for the parent element + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + unit_field(nu_blood_press,nunit) = (Pin+Pout)/2.0_dp + unit_field(nu_tt,nunit) = TT_TOTAL + unit_field(nu_sa,nunit) = TOTAL_SHEET_SA + + ENDIF ! deallocate (SparseCol, STAT = AllocateStatus) deallocate (SparseRow, STAT = AllocateStatus) diff --git a/src/lib/geometry.f90 b/src/lib/geometry.f90 index b4555cf9..6ba6248f 100644 --- a/src/lib/geometry.f90 +++ b/src/lib/geometry.f90 @@ -414,6 +414,7 @@ subroutine append_units() unit_field=0.0_dp units=0 elem_units_below(1:num_elems) = 0 !initialise the number of terminal units below a branch + elem_field(ne_unit,:) = 0.0_dp nu=0 do ne=1,num_elems @@ -421,6 +422,7 @@ subroutine append_units() nu=nu+1 units(nu)=ne !Set up units array containing terminals elem_units_below(ne)=1 + elem_field(ne_unit,ne) = real(nu) endif enddo diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index ea747393..47007996 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -28,11 +28,11 @@ module indices ne_resist=0,ne_t_resist=0,ne_Vdot=0,ne_Vdot0=0,ne_a_A=0,& ne_dvdt=0,ne_radius_in=0,ne_radius_in0=0,& ne_radius_out=0,ne_radius_out0=0,ne_group=0,ne_Qdot=0, & - ne_vd_bel=0, ne_vol_bel=0 + ne_vd_bel=0, ne_vol_bel=0, ne_unit=0 ! indices for unit_field integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, & nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,& - nu_vd=0,nu_perf=0,nu_blood_press=0 + nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0 !indices for gas exchange field ! indices for gasex_field integer,parameter :: num_gx = 12 @@ -60,12 +60,12 @@ module indices ne_resist,ne_t_resist,ne_Vdot,ne_Vdot0,ne_a_A,& ne_dvdt,ne_radius_in,ne_radius_in0,ne_radius_out,& ne_radius_out0,ne_group,ne_Qdot, & - ne_vd_bel, ne_vol_bel + ne_vd_bel, ne_vol_bel,ne_unit public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, & nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,& nu_conc1,nu_vent,nu_vd,& - nu_perf,nu_blood_press + nu_perf,nu_blood_press,nu_sa,nu_tt public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, & ng_p_cap_o2, ng_p_cap_co2,ng_source_o2,ng_source_co2, & @@ -281,7 +281,7 @@ subroutine perfusion_indices num_nj=1 nj_bv_press=1 !pressure in blood vessel ! indices for elem_field - num_ne=9 + num_ne=10 ne_radius=1 !strained average radius over whole element ne_radius_in=2 !strained radius into an element ne_radius_out=3 !strained radius out of an element @@ -291,10 +291,13 @@ subroutine perfusion_indices ne_Qdot=7 !flow in an element ne_resist=8 !resistance of a blood vessel ne_group=9!Groups vessels into arteries (field=0), capillaries (field=1) and veins(field=2) + ne_unit=10 !store the unit number for terminal element !indices for units - num_nu=2 + num_nu=4 nu_perf=1 nu_blood_press=2 + nu_sa = 3 + nu_tt = 4 call enter_exit(sub_name,2) end subroutine perfusion_indices diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index cdcc444a..3cfe2977 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -26,17 +26,18 @@ module lymphatics !Interfaces ! private public alveolar_capillary_flux + public lymphatic_transport contains !!!############################################################################# - subroutine alveolar_capillary_flux(num_nodes) + subroutine alveolar_capillary_flux(ne) !*alveolar_capillary_flux:* calculate fluid flux from blood to interstitium !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ALVEOLAR_CAPILLARY_FLUX" :: ALVEOLAR_CAPILLARY_FLUX use other_consts,only: pi - integer,intent(in) :: num_nodes + integer,intent(in) :: ne ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 0 ! 0 = male, 1 = female @@ -73,7 +74,7 @@ subroutine alveolar_capillary_flux(num_nodes) ! Local variables - integer :: i,liflowcount,printcount,time_loop + integer :: i,liflowcount,nunit,printcount,time_loop real(dp) :: alveolar_volume,breathing_function,capillary_flow,capillary_vps,capillary_osmotic,capillary_osm_n, & capillary_volume,capillary_volume_raw,cap_osm_conc,diffusion,excess,flux_a,flux_b,flux_c,gas_diffusion_restriction, & initial_lymphatic_flow,initial_lymphatic_pressure,initial_lymphatic_surface_area,initial_lymphatic_volume, & @@ -90,6 +91,14 @@ subroutine alveolar_capillary_flux(num_nodes) sub_name = 'alveolar_capillary_flux' call enter_exit(sub_name,1) + ! get information for the unit fron unit_field + ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + ! unit_field(nu_blood_press,nunit) + ! unit_field(nu_tt,nunit) + ! unit_field(nu_sa,nunit) + ! + ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g capillary_volume_raw = body_mass/0.3474_dp ! Gehr 213 ml and body mass of 74 kg - rough estimate? @@ -300,11 +309,11 @@ subroutine alveolar_capillary_flux(num_nodes) total_flux = total_hydro_flux ! +total_osm_flux sumuptake = sumuptake + initial_lymphatic_flow printcount = printcount + 1 - if (printcount.eq.100)then - write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & - flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume + !if (printcount.eq.100)then + ! write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & + ! flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + ! 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + ! initial_lymphatic_volume,total_flux,alveolar_volume ! write(*,*) 'at ',time, ' s:' ! write(*,*) ' flux per second is ',flux_c/transit_time*1000,' ul' ! write(*,*) ' interstitial volume is ',interstitial_volume,' mL' @@ -318,13 +327,42 @@ subroutine alveolar_capillary_flux(num_nodes) ! write(*,*) ' total flux is ',total_flux,' mL' ! write(*,*) ' alveolar volume is ',alveolar_volume printcount = 0 - endif + !endif enddo ! while + write(*,'(i8,f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') nunit,time, & + flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + initial_lymphatic_volume,total_flux,alveolar_volume call enter_exit(sub_name,2) end subroutine alveolar_capillary_flux +!!!############################################################################# + + subroutine lymphatic_transport() + !*lymphatic_transport:* whole system transport + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LYMPHATIC_TRANSPORT" :: LYMPHATIC_TRANSPORT + + integer :: ne + + character(len=60) :: sub_name + + ! -------------------------------------------------------------------------- + + sub_name = 'lymphatic_transport' + call enter_exit(sub_name,1) + + do ne = 1,num_elems + if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then + call alveolar_capillary_flux(ne) + endif + enddo + + call enter_exit(sub_name,2) + + end subroutine lymphatic_transport + !!!############################################################################# end module lymphatics From ea0357b9cdbc0155c8cb1ec45d5a0da028e6b163 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Thu, 14 Oct 2021 10:11:36 +1300 Subject: [PATCH 04/30] lymphatics_transport subroutine and bindings, to evaluate whole system transport. unit_field used to transfer transit time, surface area, and average capillary pressure for lymphatics calculation --- src/bindings/c/lymphatics.c | 7 +++++++ src/bindings/c/lymphatics.f90 | 16 ++++++++++++++++ src/bindings/c/lymphatics.h | 1 + 3 files changed, 24 insertions(+) diff --git a/src/bindings/c/lymphatics.c b/src/bindings/c/lymphatics.c index 01149745..8a4bb894 100644 --- a/src/bindings/c/lymphatics.c +++ b/src/bindings/c/lymphatics.c @@ -2,8 +2,15 @@ #include "string.h" void alveolar_capillary_flux_c(int *num_nodes); +void lymphatic_transport_c(); + void alveolar_capillary_flux(int num_nodes) { alveolar_capillary_flux_c(&num_nodes); } + +void lymphatic_transport() +{ + lymphatic_transport_c(); +} diff --git a/src/bindings/c/lymphatics.f90 b/src/bindings/c/lymphatics.f90 index 7740974f..ace0f77e 100644 --- a/src/bindings/c/lymphatics.f90 +++ b/src/bindings/c/lymphatics.f90 @@ -24,5 +24,21 @@ end subroutine alveolar_capillary_flux_c ! !################################################################################### ! +!*lymphatic_transport* + subroutine lymphatic_transport_c() bind(C, name="lymphatic_transport_c") + use lymphatics,only: lymphatic_transport + implicit none + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_lymphatic_transport() +#else + call lymphatic_transport() +#endif + + end subroutine lymphatic_transport_c + +! +!################################################################################### +! end module lymphatics_c diff --git a/src/bindings/c/lymphatics.h b/src/bindings/c/lymphatics.h index a0bf39c2..3be674d6 100644 --- a/src/bindings/c/lymphatics.h +++ b/src/bindings/c/lymphatics.h @@ -4,5 +4,6 @@ #include "symbol_export.h" SHO_PUBLIC void alveolar_capillary_flux(int num_nodes); +SHO_PUBLIC void lymphatic_transport(); #endif /* AETHER_LYMPHATICS_H */ From d1ae6239db316fe544c704a943c8360e9b22d7c5 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Sun, 17 Oct 2021 19:14:33 +1300 Subject: [PATCH 05/30] code produces results and redundancies removed --- src/lib/lymphatics.f90 | 116 ++++++++++++++++------------------------- 1 file changed, 44 insertions(+), 72 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 8fce535d..f4ce1d6e 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -46,18 +46,8 @@ subroutine alveolar_capillary_flux(ne) real(dp),parameter :: lymphatic_density = 1.0_dp !to be calculated from CT?? ! Capillary parameters - real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp ! #0.809e-7#0.13256e-7*133/20#/20 #0.1278# obtained from literature - Parker range in cm H2O - real(dp),parameter :: open_capillaries = 1.0_dp/6.0_dp - real(dp) :: collecting_volume = 0.0_dp - real(dp) :: collecting_lymphatic_pressure = 5.0_dp + real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp ! obtained from literature - Parker range in cm H2O real(dp),parameter :: reflection_coefficient = 0.0_dp - real(dp),parameter :: capillary_SA = 1022.0_dp*0.599998685_dp !#0.970116# ! 1022 standardises capillary to acinar volumes as one interstitium is there for al 1022 capillaries - real(dp),parameter :: transit_time = 1.315728977_dp !#2.813324# - real(dp),parameter :: capillary_pressure = 11.51694179_dp !#23.37539# - - ! Interstitial parameters - real(dp),parameter :: interstitial_resistance = 38000.0_dp !Bertram 2011 (600-12000000) dyn/cm^2 - real(dp),parameter :: diffusion_constant = 0.00000015_dp ! Initial lymphatics parameters real(dp),parameter :: lymphatic_resistance = 1000.0_dp ! Ngo 2019 @@ -82,7 +72,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_osmotic,interstitial_pressure_a,interstitial_pressure_b,interstitial_saturation,interstitial_volume, & int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,net_flux, & osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time2,time_dt,time_period,time_variable,total_flux, & - total_hydro_flux,total_osm_flux + total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA character(len=60) :: sub_name @@ -94,16 +84,13 @@ subroutine alveolar_capillary_flux(ne) ! get information for the unit fron unit_field ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) - ! unit_field(nu_blood_press,nunit) - ! unit_field(nu_tt,nunit) - ! unit_field(nu_sa,nunit) + capillary_pressure = unit_field(nu_blood_press,nunit)/133.0_dp + transit_time = unit_field(nu_tt,nunit) + capillary_SA = unit_field(nu_sa,nunit) ! ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g - capillary_volume_raw = body_mass/0.3474_dp ! Gehr 213 ml and body mass of 74 kg - rough estimate? - capillary_volume = (capillary_volume_raw*open_capillaries)/real(au) - capillary_flow = (0.085_dp/60.0_dp) ! uL/sec per unit - will obtain from perfusion model? ! interstitial values interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(au))*1000.0_dp ! maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) @@ -113,6 +100,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume_b = 0.49_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 + initial_lymphatic_surface_area = capillary_SA*3.648_dp ! ! initial lymphatic values initial_lymphatic_volume = 0.0_dp @@ -125,44 +113,24 @@ subroutine alveolar_capillary_flux(ne) int_osm_n = 0.0_dp initial_osm_n = 0.0_dp total_osm_flux = 0.0_dp - - ! time is the time that each blood particle is in the capillaries - initial_lymphatic_surface_area = capillary_SA*3.648_dp ! #0.00014046 #to be normalised to volume of lung? - - write(*,'('' blood transit time is '',f12.2,'' s'')') transit_time - ! second is how much blood flows through the capillary each second - - capillary_vps = 1.0_dp / transit_time - write(*,'('' volume per second is '',f12.2,'' mL'')') capillary_vps - + time = 0.0_dp time2 = 1.0_dp printcount = 0 total_hydro_flux = 0.0_dp time_loop = 96 - !times = [] - !int_flux = [] - !interstitial_saturation_as = [] - !initial_lymph_flux = [] - !initial_lymph_sat = [] - !total_flux_list = [] - !total_uptake_list = [] - !collected_volume = [] - !flux_difference = [] - sumflux = 0.0_dp - sumuptake = 0.0_dp breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) - write(*,'('' breathing rate is '',f8.2,'' breaths per minute'')') breathing_rate - write(*,'('' lung mass is '',f8.2)') lung_mass - write(*,'('' cap SA is '',f8.2)') capillary_SA - write(*,'('' int cap is '',f8.2)') interstitial_capacity - write(*,'('' int capA is '',f8.2)') interstitial_capacity_a - write(*,'('' int capB is '',f8.2)') interstitial_capacity_b - write(*,'('' int vol is '',f8.2)') interstitial_volume_b - write(*,'('' lymphSA is '',f8.2)') initial_lymphatic_surface_area - write(*,'('' int cap is '',f8.2)') interstitial_capacity + !write(*,'('' breathing rate is '',f8.2,'' breaths per minute'')') breathing_rate + !write(*,'('' lung mass is '',f8.2)') lung_mass + !write(*,'('' cap SA is '',f8.2)') capillary_SA + !write(*,'('' int cap is '',f8.2)') interstitial_capacity + !write(*,'('' int capA is '',f8.2)') interstitial_capacity_a + !write(*,'('' int capB is '',f8.2)') interstitial_capacity_b + !write(*,'('' int vol is '',f8.2)') interstitial_volume_b + !write(*,'('' lymphSA is '',f8.2)') initial_lymphatic_surface_area + !write(*,'('' int cap is '',f8.2)') interstitial_capacity write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & &''init lymph|init lymph| total| alveolar|'')') @@ -171,7 +139,10 @@ subroutine alveolar_capillary_flux(ne) do while(time < 32400.0_dp) time = time + transit_time - do i = 0,95 + write(*,'(''code is running at '',f8.2)')time + write(*,'(''time variable '',f8.2)')time_variable + write(*,'(''breathing function '',f8.2)')breathing_function + do while(time2 < 97.0_dp) interstitial_volume = interstitial_volume_a + interstitial_volume_b interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% time_variable = time-transit_time + (transit_time*(time2/time_loop)) @@ -260,11 +231,11 @@ subroutine alveolar_capillary_flux(ne) interstitial_capacity_b)**2_dp) + (-213.23_dp * (interstitial_volume_b / interstitial_capacity_b)) & + 11.812_dp)* 4.41335e-8_dp !(capillary_conductivity) endif - !initial_lymphatic_pressure = (math.sin(time_variable * breathing_function + math.pi)) + ((7 * (interstitial_volume_b / interstitial_capacity_b)) - 7) + initial_lymphatic_pressure = ((1.47_dp * sin(time_variable * breathing_function + pi/2_dp)) + & ((6.82_dp* (interstitial_volume_b / interstitial_capacity_b)**2_dp) + (0.77_dp * (interstitial_volume_b / & interstitial_capacity_b)) - 6.52_dp)) - ! initial_lymphatic_pressure = math.sin((breathing_function * time_variable) + (math.pi/2)) + ((2 * ((interstitial_volume_b/interstitial_capacity_b) ** 2) +3 * (interstitial_volume_b/interstitial_capacity_b) - 7)+(2*(interstitial_volume_b/interstitial_capacity_b))) ! based on range of values + if(interstitial_volume.le.0.0_dp)then initial_lymphatic_flow = 0.0_dp interstitial_volume = 0.0_dp @@ -275,7 +246,6 @@ subroutine alveolar_capillary_flux(ne) else initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - & initial_lymphatic_pressure)) * (transit_time/time_loop) - !write(*,*) (lymphatic_conductivity, initial_lymphatic_surface_area, interstitial_pressure, initial_lymphatic_pressure, transit_time, time_loop) endif else initial_lymphatic_flow = 0.0_dp @@ -294,9 +264,11 @@ subroutine alveolar_capillary_flux(ne) initial_lymph_conc = 0.0_dp endif - gas_diffusion_restriction = 0.0000152587890625_dp * exp(13.8629436112_dp*(interstitial_saturation/100.0_dp)) + !gas_diffusion_restriction = 0.0000152587890625_dp * exp(13.8629436112_dp*(interstitial_saturation/100.0_dp)) time2 = time2 + 1 + write(*,'(''capillary pressure is '',f8.2)')capillary_pressure + write(*,'(''interstitial pressure is '',f8.2)')interstitial_pressure_b enddo time2 = 1 @@ -309,26 +281,26 @@ subroutine alveolar_capillary_flux(ne) total_flux = total_hydro_flux ! +total_osm_flux sumuptake = sumuptake + initial_lymphatic_flow printcount = printcount + 1 - !if (printcount.eq.100)then - ! write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & - ! flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - ! 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - ! initial_lymphatic_volume,total_flux,alveolar_volume -! write(*,*) 'at ',time, ' s:' -! write(*,*) ' flux per second is ',flux_c/transit_time*1000,' ul' -! write(*,*) ' interstitial volume is ',interstitial_volume,' mL' -! write(*,*) ' interstitial volume a is ',interstitial_volume_a,' mL' -! write(*,*) ' interstitial volume b is ',interstitial_volume_b,' mL' -! write(*,*) ' interstitial saturation is ',100*interstitial_saturation,'%' -! ! write(*,*) ' diffusion is ' + str("%.10f" % diffusion) + ' mL') -! ! write(*,*) ' interstitial pressure is ' + str("%.2f" % interstitial_pressure) + ' mmHg') -! write(*,*) ' initial lymphatic flux is ',initial_lymphatic_flow*1000,' uL' -! write(*,*) ' initial lymphatic volume is ',initial_lymphatic_volume,' mL' -! write(*,*) ' total flux is ',total_flux,' mL' -! write(*,*) ' alveolar volume is ',alveolar_volume + if (printcount.eq.100)then + write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & + flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + initial_lymphatic_volume,total_flux,alveolar_volume + write(*,*) 'at ',time, ' s:' + write(*,*) ' flux per second is ',flux_c/transit_time*1000,' ul' + write(*,*) ' interstitial volume is ',interstitial_volume,' mL' + write(*,*) ' interstitial volume a is ',interstitial_volume_a,' mL' + write(*,*) ' interstitial volume b is ',interstitial_volume_b,' mL' + write(*,*) ' interstitial saturation is ',100*interstitial_saturation,'%' + write(*,*) ' diffusion is ',diffusion,' mL' + write(*,*) ' interstitial pressure is ',interstitial_pressure_b,' mmHg' + write(*,*) ' initial lymphatic flux is ',initial_lymphatic_flow*1000,' uL' + write(*,*) ' initial lymphatic volume is ',initial_lymphatic_volume,' mL' + write(*,*) ' total flux is ',total_flux,' mL' + write(*,*) ' alveolar volume is ',alveolar_volume printcount = 0 - !endif - enddo ! while + endif + enddo !while write(*,'(i8,f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') nunit,time, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & From 02391c05eeec43eb6599b810a9bbe374d47b299a Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Tue, 19 Oct 2021 17:27:22 +1300 Subject: [PATCH 06/30] debugging alveolar flux code --- src/lib/lymphatics.f90 | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index f4ce1d6e..a195c5ae 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -139,9 +139,11 @@ subroutine alveolar_capillary_flux(ne) do while(time < 32400.0_dp) time = time + transit_time - write(*,'(''code is running at '',f8.2)')time - write(*,'(''time variable '',f8.2)')time_variable - write(*,'(''breathing function '',f8.2)')breathing_function + !write(*,'(''code is running at '',f8.2)')time + !write(*,'(''time variable '',f8.2)')time_variable + !write(*,'(''breathing function '',f8.2)')breathing_function + write(*,*) 'time=',time,transit_time,time2 + write(*,*) 'ia,b',interstitial_volume_a, interstitial_volume_b do while(time2 < 97.0_dp) interstitial_volume = interstitial_volume_a + interstitial_volume_b interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% @@ -154,6 +156,10 @@ subroutine alveolar_capillary_flux(ne) ((-3.98_dp * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & 8.03_dp * (interstitial_volume_b / interstitial_capacity_b) - 6.52_dp) ! pressure determined from saturation equation based of literature (currently linear, but likely not) + + write(*,*) time_variable,interstitial_pressure_a,interstitial_pressure_b,interstitial_volume_a, & + interstitial_volume_b + if(capillary_pressure > interstitial_pressure_a)then flux_a = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & interstitial_pressure_a)) * (transit_time/time_loop) @@ -172,6 +178,7 @@ subroutine alveolar_capillary_flux(ne) if(interstitial_volume_a + flux_a > interstitial_capacity_a)then excess = flux_a - (interstitial_capacity_a - interstitial_volume_a) interstitial_volume_a = interstitial_capacity_a + write(*,*) 'up1:',interstitial_volume_a alveolar_volume = alveolar_volume + 0.5_dp*excess if((interstitial_volume_b + 0.5_dp*excess) > interstitial_capacity_b)then @@ -183,6 +190,7 @@ subroutine alveolar_capillary_flux(ne) endif else interstitial_volume_a = interstitial_volume_a + flux_a + write(*,*) 'up2:',interstitial_volume_a, flux_a endif interstitial_volume_b = interstitial_volume_b + flux_b @@ -191,6 +199,7 @@ subroutine alveolar_capillary_flux(ne) (160_dp*1.1_dp))*(transit_time/time_loop) interstitial_volume_b = interstitial_volume_b + diffusion interstitial_volume_a = interstitial_volume_a - diffusion + write(*,*) 'up3:',interstitial_volume_a, diffusion !alveolar_volume = alveolar_volume - 4.45e-5 ! Osmotic @@ -267,8 +276,8 @@ subroutine alveolar_capillary_flux(ne) !gas_diffusion_restriction = 0.0000152587890625_dp * exp(13.8629436112_dp*(interstitial_saturation/100.0_dp)) time2 = time2 + 1 - write(*,'(''capillary pressure is '',f8.2)')capillary_pressure - write(*,'(''interstitial pressure is '',f8.2)')interstitial_pressure_b + !write(*,'(''capillary pressure is '',f8.2)')capillary_pressure + !write(*,'(''interstitial pressure is '',f8.2)')interstitial_pressure_b enddo time2 = 1 @@ -281,7 +290,8 @@ subroutine alveolar_capillary_flux(ne) total_flux = total_hydro_flux ! +total_osm_flux sumuptake = sumuptake + initial_lymphatic_flow printcount = printcount + 1 - if (printcount.eq.100)then + !if (printcount.eq.100)then + if (printcount.eq.1)then write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & @@ -300,6 +310,7 @@ subroutine alveolar_capillary_flux(ne) write(*,*) ' alveolar volume is ',alveolar_volume printcount = 0 endif + read(*,*) enddo !while write(*,'(i8,f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') nunit,time, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & From ca29379b29db55262c910df979c4e55b6e576ee8 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Mon, 25 Oct 2021 17:28:31 +1300 Subject: [PATCH 07/30] ability to transfer data from ventilation and perfusion simulations to the lymphatics simulation --- src/bindings/c/arrays.c | 7 ++ src/bindings/c/arrays.f90 | 20 +++++ src/bindings/c/arrays.h | 1 + src/bindings/c/imports.c | 12 +++ src/bindings/c/imports.f90 | 48 +++++++++++ src/bindings/c/imports.h | 2 + src/bindings/c/indices.c | 6 ++ src/bindings/c/indices.f90 | 17 ++++ src/bindings/c/indices.h | 1 + src/lib/arrays.f90 | 34 +++++++- src/lib/exports.f90 | 43 ++++++---- src/lib/imports.f90 | 140 ++++++++++++++++++++++++++++++- src/lib/indices.f90 | 65 +++++++++++++-- src/lib/lymphatics.f90 | 163 ++++++++++++++----------------------- src/lib/ventilation.f90 | 8 ++ 15 files changed, 438 insertions(+), 129 deletions(-) diff --git a/src/bindings/c/arrays.c b/src/bindings/c/arrays.c index e64c8f72..d4818791 100644 --- a/src/bindings/c/arrays.c +++ b/src/bindings/c/arrays.c @@ -1,9 +1,16 @@ #include "arrays.h" +#include "string.h" void set_node_field_value_c(int *row, int *col, double *value); +void update_parameter_c(const char *parameter_name, int *parameter_name_len, double *parameter_value); void set_node_field_value(int row, int col, double value) { set_node_field_value_c(&row, &col, &value); } +void update_parameter(const char *parameter_name, double parameter_value) +{ + int parameter_name_len = (int)strlen(parameter_name); + update_parameter_c(parameter_name, ¶meter_name_len, ¶meter_value); +} diff --git a/src/bindings/c/arrays.f90 b/src/bindings/c/arrays.f90 index 45790b5d..2d7d1349 100644 --- a/src/bindings/c/arrays.f90 +++ b/src/bindings/c/arrays.f90 @@ -30,5 +30,25 @@ subroutine set_node_field_value_c(row, col, value) bind(C, name="set_node_field_ end subroutine set_node_field_value_c + subroutine update_parameter_c(parameter_name, parameter_name_len, parameter_value) bind(C, name="update_parameter_c") + use iso_c_binding, only: c_ptr + use utils_c, only: strncpy + use other_consts, only: max_filename_len + use arrays, only: update_parameter + implicit none + + integer, intent(in) :: parameter_name_len + real(dp), intent(in) :: parameter_value + type(c_ptr),value, intent(in) :: parameter_name + character(len=max_filename_len) :: parameter_name_f + + call strncpy(parameter_name_f, parameter_name, parameter_name_len) +#if defined _WIN32 && defined __INTEL_COMPILER + call so_update_parameter(parameter_name_f, parameter_value) +#else + call update_parameter(parameter_name_f, parameter_value) +#endif + end subroutine update_parameter_c + end module arrays_c diff --git a/src/bindings/c/arrays.h b/src/bindings/c/arrays.h index 1468b841..15495724 100644 --- a/src/bindings/c/arrays.h +++ b/src/bindings/c/arrays.h @@ -5,5 +5,6 @@ #include "symbol_export.h" SHO_PUBLIC void set_node_field_value(int row, int col, double value); +SHO_PUBLIC void update_parameter(const char *parameter_name, double parameter_value); #endif /* AETHER_ARRAYS_H */ diff --git a/src/bindings/c/imports.c b/src/bindings/c/imports.c index 1c7dc2b2..873f51fb 100644 --- a/src/bindings/c/imports.c +++ b/src/bindings/c/imports.c @@ -4,9 +4,21 @@ #include +void import_capillary_c(const char *FLOWFILE,int *FLOWFILE_LEN); +void import_terminal_c(const char *FLOWFILE,int *FLOWFILE_LEN); void import_ventilation_c(const char *FLOWFILE,int *FLOWFILE_LEN); void import_perfusion_c(const char *FLOWFILE,int *FLOWFILE_LEN); +void import_capillary(const char *FLOWFILE) +{ + int filename_len = strlen(FLOWFILE); + import_capillary_c(FLOWFILE, &filename_len); +} +void import_terminal(const char *FLOWFILE) +{ + int filename_len = strlen(FLOWFILE); + import_terminal_c(FLOWFILE, &filename_len); +} void import_ventilation(const char *FLOWFILE) { int filename_len = strlen(FLOWFILE); diff --git a/src/bindings/c/imports.f90 b/src/bindings/c/imports.f90 index 6c45b4b9..34a70530 100644 --- a/src/bindings/c/imports.f90 +++ b/src/bindings/c/imports.f90 @@ -3,6 +3,54 @@ module imports_c private contains +! +!################################################################################### +! + subroutine import_capillary_c(FLOWFILE, filename_len) bind(C, name="import_capillary_c") + + use iso_c_binding, only: c_ptr + use utils_c, only: strncpy + use imports, only: import_capillary + use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN + implicit none + integer,intent(in) :: filename_len + type(c_ptr), value, intent(in) :: FLOWFILE + character(len=MAX_FILENAME_LEN) :: filename_f + + call strncpy(filename_f, FLOWFILE, filename_len) + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_import_capillary(filename_f) +#else + call import_capillary(filename_f) +#endif + + end subroutine import_capillary_c + +! +!################################################################################### +! + subroutine import_terminal_c(FLOWFILE, filename_len) bind(C, name="import_terminal_c") + + use iso_c_binding, only: c_ptr + use utils_c, only: strncpy + use imports, only: import_terminal + use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN + implicit none + integer,intent(in) :: filename_len + type(c_ptr), value, intent(in) :: FLOWFILE + character(len=MAX_FILENAME_LEN) :: filename_f + + call strncpy(filename_f, FLOWFILE, filename_len) + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_import_terminal(filename_f) +#else + call import_terminal(filename_f) +#endif + + end subroutine import_terminal_c + ! !################################################################################### ! diff --git a/src/bindings/c/imports.h b/src/bindings/c/imports.h index 9a5c6124..3232df27 100644 --- a/src/bindings/c/imports.h +++ b/src/bindings/c/imports.h @@ -3,6 +3,8 @@ #include "symbol_export.h" +SHO_PUBLIC void import_capillary(const char *FLOWFILE); +SHO_PUBLIC void import_terminal(const char *FLOWFILE); SHO_PUBLIC void import_ventilation(const char *FLOWFILE); SHO_PUBLIC void import_perfusion(const char *FLOWFILE); diff --git a/src/bindings/c/indices.c b/src/bindings/c/indices.c index 3ab18822..20174baa 100644 --- a/src/bindings/c/indices.c +++ b/src/bindings/c/indices.c @@ -4,6 +4,7 @@ #include void define_problem_type_c(const char *PROBLEMTYPE,int *PROBLEMTYPE_LEN); +void lymphatic_indices_c(); void ventilation_indices_c(); void perfusion_indices_c(); int get_ne_radius_c(); @@ -18,6 +19,11 @@ void define_problem_type(const char *PROBLEMTYPE) } +void lymphatic_indices() +{ + lymphatic_indices_c(); +} + void ventilation_indices() { ventilation_indices_c(); diff --git a/src/bindings/c/indices.f90 b/src/bindings/c/indices.f90 index f53ac8f9..4048cffe 100644 --- a/src/bindings/c/indices.f90 +++ b/src/bindings/c/indices.f90 @@ -43,6 +43,23 @@ end subroutine ventilation_indices_c ! !###################################################################### ! +!> Lymphatic indices + subroutine lymphatic_indices_c() bind(C, name="lymphatic_indices_c") + + use indices, only: lymphatic_indices + implicit none + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_lymphatic_indices() +#else + call lymphatic_indices() +#endif + + end subroutine lymphatic_indices_c + +! +!###################################################################### +! !> Perfusion indices subroutine perfusion_indices_c() bind(C, name="perfusion_indices_c") diff --git a/src/bindings/c/indices.h b/src/bindings/c/indices.h index f1ec8f87..2c9c6907 100644 --- a/src/bindings/c/indices.h +++ b/src/bindings/c/indices.h @@ -3,6 +3,7 @@ #include "symbol_export.h" SHO_PUBLIC void define_problem_type(const char *PROBLEMTYPE); +SHO_PUBLIC void lymphatic_indices(); SHO_PUBLIC void ventilation_indices(); SHO_PUBLIC void perfusion_indices(); SHO_PUBLIC int get_ne_radius(); diff --git a/src/lib/arrays.f90 b/src/lib/arrays.f90 index d9070cd6..0812681c 100644 --- a/src/lib/arrays.f90 +++ b/src/lib/arrays.f90 @@ -108,6 +108,16 @@ module arrays real(dp) :: air_density end type fluid_properties + type default_lymphatic_properties + ! default values for lymphatic model parameters + real(dp) :: lymphatic_density = 1.0_dp !to be calculated from CT?? + real(dp) :: lymphatic_integrity = 1.0_dp ! a measure of how 'leaky' the lymphatic vessels are and prone to backflow + real(dp) :: reflection_coefficient = 0.0_dp + end type default_lymphatic_properties + +!!! arrays that start with default values, updated during simulations/by user + type(default_lymphatic_properties) :: lymphatic_properties + ! temporary, for debugging: real(dp) :: unit_before @@ -120,8 +130,9 @@ module arrays elem_units_below, maxgen,capillary_bf_parameters, zero_tol,loose_tol,gasex_field, & num_lines_2d, lines_2d, line_versn_2d, lines_in_elem, nodes_in_line, elems_2d, & elem_cnct_2d, elem_nodes_2d, elem_versn_2d, elem_lines_2d, elems_at_node_2d, arclength, & - scale_factors_2d, parentlist, fluid_properties, elasticity_vessels, admittance_param, & - elasticity_param, all_admit_param + scale_factors_2d, parentlist, fluid_properties, lymphatic_properties, & + elasticity_vessels, admittance_param, & + elasticity_param, all_admit_param, update_parameter contains subroutine set_node_field_value(row, col, value) @@ -135,5 +146,24 @@ subroutine set_node_field_value(row, col, value) end subroutine set_node_field_value + subroutine update_parameter(parameter_name, parameter_value) + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_UPDATE_PARAMETER" :: UPDATE_PARAMETER + implicit none + real(dp), intent(in) :: parameter_value + character(len=*), intent(in) :: parameter_name + + select case(parameter_name) + +!!! lymphatic_properties + case('lymphatic_density') + lymphatic_properties%lymphatic_density = parameter_value + case('lymphatic_integrity') + lymphatic_properties%lymphatic_integrity = parameter_value + case('reflection_coefficient') + lymphatic_properties%reflection_coefficient = parameter_value + + end select + + end subroutine update_parameter end module arrays diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index 4fde312d..c6dbc1ef 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -503,7 +503,8 @@ subroutine export_terminal_solution(EXNODEFILE, name) !*** Write the field information VALUE_INDEX=1 if(FIRST_NODE)THEN - write(10,'( '' #Fields=5'' )') + write(10,'( '' #Fields=8'' )') + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') do nj=1,3 if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") @@ -512,32 +513,38 @@ subroutine export_terminal_solution(EXNODEFILE, name) write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 VALUE_INDEX=VALUE_INDEX+1 enddo + !element number + write(10,'('' 2) terminal_element, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 !Ventilation (tidal volume/insp time) - write(10,'('' 2) flow, field, rectangular cartesian, #Components=1'')') + write(10,'('' 3) flow, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 - !VALUE_INDEX=VALUE_INDEX+1 - !Volume - !write(10,'('' 3) volume, field, rectangular cartesian, #Components=1'')') - !write(10,'(2X,''1. '')',advance="no") - !write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 - !VALUE_INDEX=VALUE_INDEX+1 - !!Pressure - !write(10,'('' 4) pressure, field, rectangular cartesian, #Components=1'')') - !write(10,'(2X,''1. '')',advance="no") - !write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 !Compliance - write(10,'('' 5) compliance, field, rectangular cartesian, #Components=1'')') + write(10,'('' 4) compliance, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 VALUE_INDEX=VALUE_INDEX+1 !Pleural pressure - write(10,'('' 6) pleural pressure, field, rectangular cartesian, #Components=1'')') + write(10,'('' 5) pleural_pressure, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 VALUE_INDEX=VALUE_INDEX+1 !Tidal volume - write(10,'('' 7) tidal volume, field, rectangular cartesian, #Components=1'')') + write(10,'('' 6) tidal_volume, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + !maximum Pe at a unit + write(10,'('' 7) max_Pe, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + !minimum Pe at a unit + write(10,'('' 8) min_Pe, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 endif !FIRST_NODE @@ -546,12 +553,14 @@ subroutine export_terminal_solution(EXNODEFILE, name) do nj=1,3 write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np)) !Coordinates enddo !njj2 + write(10,'(2X,4(1X,F12.6))') real(ne) !element number write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vent,NOLIST)) !Ventilation - !write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vol,nolist)) !Volume (end expiration) - !write(10,'(2X,4(1X,F12.6))') (unit_field(nu_press,nolist)) !Pressure write(10,'(2X,4(1X,F12.6))') (unit_field(nu_comp,nolist)) !Compliance (end exp) write(10,'(2X,4(1X,F12.6))') (unit_field(nu_pe,nolist)) !Recoil pressure write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vt,nolist)) !Tidal volume + write(10,'(2X,4(1X,F12.6))') (unit_field(nu_Pe_max,nolist)) !maximum elastic recoil + write(10,'(2X,4(1X,F12.6))') (unit_field(nu_Pe_min,nolist)) !minimum elastic recoil + FIRST_NODE=.FALSE. np_last=np enddo !nolist (np) diff --git a/src/lib/imports.f90 b/src/lib/imports.f90 index ecbb9c68..57b30700 100644 --- a/src/lib/imports.f90 +++ b/src/lib/imports.f90 @@ -24,11 +24,62 @@ module imports !Module variables !Interfaces - private + private + public import_capillary public import_ventilation public import_perfusion - + public import_terminal + contains + ! + !########################################################################################### + ! + !>*import_capillary:* This subroutine reads in the results for the micro-circulatory + ! components (capillary bed within 'units') of a perfusion model . + subroutine import_capillary(micro_unit_file) + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_IMPORT_CAPILLARY" :: IMPORT_CAPILLARY + + character(len=MAX_FILENAME_LEN),intent(in) :: micro_unit_file + !local variables + integer :: count_units,ios,ne,nunit + real(dp) :: Pin,Pout,Qtot,temp(7),TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL + + character(len=60) :: sub_name + + sub_name = 'import_capillary' + call enter_exit(sub_name,1) + + open(10, file = micro_unit_file, status='old') + + ! ios is negative if an end of record condition is encountered or if + ! an endfile condition was detected. It is positive if an error was + ! detected. ios is zero otherwise. + ios = 0 + count_units = 0 + do while (ios == 0) + read(10, '(I6,X,5(F9.2,X),2(F8.5,X),F10.2,X,F8.4,X,F10.4,X,F10.3,X,F8.4,X,F9.4,X)', iostat=ios) & + ne,temp(1:3),Pin,Pout,temp(4),Qtot,temp(5),TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL,temp(6:7) + if(ios == 0)then + ! record the unit values for mean pressure, transit time, surface area. + ! ne is the 'linker' element, so nunit is for its parent element + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + unit_field(nu_blood_press,nunit) = (Pin+Pout)/2.0_dp + unit_field(nu_tt,nunit) = TT_TOTAL + unit_field(nu_sa,nunit) = TOTAL_SHEET_SA + count_units = count_units + 1 + endif + enddo + + if(count_units.ne.num_units)then + write(*,'(''WARNING: the number of capillary units ('',i6,'') does not match the geometric model ('',i6,'')'')') & + count_units,num_units + endif + + close(10) + + call enter_exit(sub_name,2) + + end subroutine import_capillary ! !############################################################################## ! @@ -139,4 +190,89 @@ subroutine import_exelemfield(FLOWFILE,field_no) call enter_exit(sub_name,2) end subroutine import_exelemfield +! +!############################################################################## +! +!>*import_terminal:* This subroutine reads in the content of an exnode field file + subroutine import_terminal(EXFILE) + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_IMPORT_TERMINAL" :: IMPORT_TERMINAL + + character(len=MAX_FILENAME_LEN),intent(in) :: EXFILE + !local variables + integer :: count_units,field_label(10),i,ibeg,iend,ierror,ne,nunit,n_fields + real(dp) :: rtemp + character(LEN=132) :: ctemp,label + character(len=300) :: readfile + + character(len=60) :: sub_name + + sub_name = 'import_terminal' + call enter_exit(sub_name,1) + + if(index(EXFILE, ".exnode")> 0) then !full filename is given + readfile = EXFILE + else ! need to append the correct filename extension + readfile = trim(EXFILE)//'.exnode' + endif + + open(10, file=readfile, status='old') + + n_fields = 0 + ierror = 0 + + read_field_labels : do + read(10, fmt="(a)", iostat=ierror) ctemp + if(index(ctemp, "Node:")> 0) exit read_field_labels + if(index(ctemp, ") ")> 0) then + ibeg = index(ctemp, ") ")+1 ! beginning of label + iend = index(ctemp, ",")-1 ! end of label + label = adjustl(ctemp(ibeg:iend)) + n_fields = n_fields + 1 + field_label(n_fields) = 0 + + select case(label) + case('terminal_element') + case('pleural_pressure') + field_label(n_fields) = nu_pe + case('tidal_volume') + field_label(n_fields) = nu_vt + case('max_Pe') + field_label(n_fields) = nu_Pe_max + case('min_Pe') + field_label(n_fields) = nu_Pe_min + end select + endif + + end do read_field_labels + + nunit = 0 + ierror = 0 + + do while (ierror == 0) + if(index(ctemp, "Node:")> 0) then + + do i = 1,3 ! read coordinates; not used + read(10, *, iostat=ierror) ctemp + enddo + read(10, *, iostat=ierror) ctemp ! read element: not used but could be + + nunit = nunit + 1 + + do i = 3,n_fields + read(10, *, iostat=ierror) rtemp + if(field_label(i).gt.0)then + unit_field(field_label(i),nunit) = rtemp + endif + enddo + + read(10, *, iostat=ierror) ctemp + endif + enddo + + close(10) + + call enter_exit(sub_name,2) + + end subroutine import_terminal + end module imports diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index 47007996..cde82141 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -12,8 +12,10 @@ module indices + use arrays use diagnostics use other_consts + use precision implicit none @@ -32,7 +34,8 @@ module indices ! indices for unit_field integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, & nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,& - nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0 + nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0,nu_Pe_max=0, & + nu_Pe_min=0 !indices for gas exchange field ! indices for gasex_field integer,parameter :: num_gx = 12 @@ -65,7 +68,7 @@ module indices public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, & nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,& nu_conc1,nu_vent,nu_vd,& - nu_perf,nu_blood_press,nu_sa,nu_tt + nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, & ng_p_cap_o2, ng_p_cap_co2,ng_source_o2,ng_source_co2, & @@ -76,8 +79,8 @@ module indices !Interfaces private - public define_problem_type,ventilation_indices, perfusion_indices, get_ne_radius, get_nj_conc1, & - growing_indices + public define_problem_type, lymphatic_indices,ventilation_indices, perfusion_indices, & + get_ne_radius, get_nj_conc1, growing_indices contains @@ -101,6 +104,10 @@ subroutine define_problem_type(PROBLEM_TYPE) case ('gas_transfer') print *, 'You are solving a gas transfer model, setting up indices' call exchange_indices + case ('lymphatic_transport') + print *, 'You are solving a lymphatic transport model, setting up indices' + call lymphatic_indices + call update_units case ('perfusion') print *, 'You are solving a static perfusion model, setting up indices' call perfusion_indices @@ -222,7 +229,7 @@ subroutine ventilation_indices ne_vd_bel = 9 ne_vol_bel = 10 ! indices for unit_field - num_nu=10 + num_nu=12 nu_vol=1 nu_comp=2 nu_Vdot0=3 @@ -233,9 +240,34 @@ subroutine ventilation_indices nu_vt=8 nu_air_press=9 nu_vent=10 + nu_Pe_max = 11 + nu_Pe_min = 12 call enter_exit(sub_name,2) end subroutine ventilation_indices +!!!############################################################################# + !> Lymphatic indices + subroutine lymphatic_indices + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LYMPHATIC_INDICES" :: LYMPHATIC_INDICES + + character(len=60) :: sub_name + + sub_name = 'lymphatic_indices' + call enter_exit(sub_name,1) + + ! indices for unit_field + num_nu = 6 + nu_perf = 1 + nu_blood_press = 2 + nu_sa = 3 + nu_tt = 4 + nu_Pe_max = 5 + nu_Pe_min = 6 + + call enter_exit(sub_name,2) + + end subroutine lymphatic_indices + !!!############################################################################# subroutine growing_indices @@ -301,7 +333,30 @@ subroutine perfusion_indices call enter_exit(sub_name,2) end subroutine perfusion_indices + + subroutine update_units + !*update_units:* reallocates unit_field following a change in problem type + + integer :: nu,num_nu_old,nunits + real(dp),allocatable :: unit_temp(:,:) + + if(allocated(unit_field))then + num_nu_old = size(unit_field,1) + allocate(unit_temp(num_nu_old,nunits)) + unit_temp = unit_field + deallocate(unit_field) + allocate(unit_field(num_nu,num_units)) + unit_field = 0.0_dp + do nu = 1,num_nu_old + unit_field(nu,:) = unit_temp(nu,:) + enddo + deallocate(unit_temp) + endif + + end subroutine update_units + + function get_ne_radius() result(res) !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_GET_NE_RADIUS" :: GET_NE_RADIUS diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index a195c5ae..cb9d67b7 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -43,16 +43,9 @@ subroutine alveolar_capillary_flux(ne) integer,parameter :: sex = 0 ! 0 = male, 1 = female integer,parameter :: au = 32768 ! number of acinar units real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp - real(dp),parameter :: lymphatic_density = 1.0_dp !to be calculated from CT?? ! Capillary parameters real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp ! obtained from literature - Parker range in cm H2O - real(dp),parameter :: reflection_coefficient = 0.0_dp - - ! Initial lymphatics parameters - real(dp),parameter :: lymphatic_resistance = 1000.0_dp ! Ngo 2019 - real(dp),parameter :: valve_resistance = 0.020_dp ! pressure requirement to overcome the valve resistance and cause flow OR the pressure lost during open the valve - real(dp),parameter :: lymphatic_integrity = 1.0_dp ! a measure of how 'leaky' the lymphatic vessels are and prone to backflow ! Osmotic pressure parameters real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! mol/L in blood plasma @@ -64,15 +57,16 @@ subroutine alveolar_capillary_flux(ne) ! Local variables - integer :: i,liflowcount,nunit,printcount,time_loop + integer :: i,liflowcount,nunit,n_timesteps,printcount real(dp) :: alveolar_volume,breathing_function,capillary_flow,capillary_vps,capillary_osmotic,capillary_osm_n, & - capillary_volume,capillary_volume_raw,cap_osm_conc,diffusion,excess,flux_a,flux_b,flux_c,gas_diffusion_restriction, & + capillary_volume,capillary_volume_raw,cap_osm_conc,diffusion,dt,excess,flux_a,flux_b,flux_c,gas_diffusion_restriction, & initial_lymphatic_flow,initial_lymphatic_pressure,initial_lymphatic_surface_area,initial_lymphatic_volume, & initial_lymph_conc,initial_osm_n,interstitial_capacity,interstitial_capacity_a,interstitial_capacity_b, & interstitial_osmotic,interstitial_pressure_a,interstitial_pressure_b,interstitial_saturation,interstitial_volume, & - int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,net_flux, & - osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time2,time_dt,time_period,time_variable,total_flux, & - total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA + int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,max_Pe, & + min_Pe,net_flux, & + open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & + time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA character(len=60) :: sub_name @@ -87,27 +81,35 @@ subroutine alveolar_capillary_flux(ne) capillary_pressure = unit_field(nu_blood_press,nunit)/133.0_dp transit_time = unit_field(nu_tt,nunit) capillary_SA = unit_field(nu_sa,nunit) - ! - + max_Pe = unit_field(nu_Pe_max,nunit) + min_Pe = unit_field(nu_Pe_min,nunit) + + write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' cmH2O'')') & + nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/98.0665_dp + ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g + open_capillaries = 1.0_dp/6.0_dp + capillary_volume_raw = body_mass/0.3474_dp ! Gehr 213 ml and body mass of 74 kg - rough estimate? + capillary_volume = (capillary_volume_raw*open_capillaries)/real(num_units) ! interstitial values - interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(au))*1000.0_dp ! maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) + ! interstitial_capacity == maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) + interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp interstitial_capacity_a = 0.005_dp*interstitial_capacity interstitial_capacity_b = 0.995_dp*interstitial_capacity interstitial_volume_a = 0.0_dp interstitial_volume_b = 0.49_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA*3.648_dp ! + initial_lymphatic_surface_area = capillary_SA*3.648_dp ! initial lymphatic values initial_lymphatic_volume = 0.0_dp ! Osmotic pressures - i = 1.0_dp - capillary_osmotic = i*capillary_molar_conc*IGC*T + i = 1 + capillary_osmotic = real(i)*capillary_molar_conc*IGC*T capillary_osm_n = 0.0_dp interstitial_osmotic = 0.0_dp int_osm_n = 0.0_dp @@ -115,39 +117,28 @@ subroutine alveolar_capillary_flux(ne) total_osm_flux = 0.0_dp time = 0.0_dp - time2 = 1.0_dp printcount = 0 total_hydro_flux = 0.0_dp - time_loop = 96 + ! dt or n_timesteps should be controlled by the user + n_timesteps = 96 + dt = transit_time/real(n_timesteps) + breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) - !write(*,'('' breathing rate is '',f8.2,'' breaths per minute'')') breathing_rate - !write(*,'('' lung mass is '',f8.2)') lung_mass - !write(*,'('' cap SA is '',f8.2)') capillary_SA - !write(*,'('' int cap is '',f8.2)') interstitial_capacity - !write(*,'('' int capA is '',f8.2)') interstitial_capacity_a - !write(*,'('' int capB is '',f8.2)') interstitial_capacity_b - !write(*,'('' int vol is '',f8.2)') interstitial_volume_b - !write(*,'('' lymphSA is '',f8.2)') initial_lymphatic_surface_area - !write(*,'('' int cap is '',f8.2)') interstitial_capacity - write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & &''init lymph|init lymph| total| alveolar|'')') write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') - do while(time < 32400.0_dp) - time = time + transit_time - !write(*,'(''code is running at '',f8.2)')time - !write(*,'(''time variable '',f8.2)')time_variable - !write(*,'(''breathing function '',f8.2)')breathing_function - write(*,*) 'time=',time,transit_time,time2 - write(*,*) 'ia,b',interstitial_volume_a, interstitial_volume_b - do while(time2 < 97.0_dp) +! do while(time < 32400.0_dp) + do while(time < 3.0_dp) + time_sum = dt + do while(time_sum < transit_time) interstitial_volume = interstitial_volume_a + interstitial_volume_b interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% - time_variable = time-transit_time + (transit_time*(time2/time_loop)) + time_variable = time + time_sum + ! calculating flux from capillary into interstitium interstitial_pressure_a = 1.47_dp * sin(time_variable * breathing_function) + & ((-3.98_dp * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & @@ -155,20 +146,17 @@ subroutine alveolar_capillary_flux(ne) interstitial_pressure_b = 1.47_dp * sin(time_variable * breathing_function) + & ((-3.98_dp * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & 8.03_dp * (interstitial_volume_b / interstitial_capacity_b) - 6.52_dp) + ! pressure determined from saturation equation based of literature (currently linear, but likely not) - - write(*,*) time_variable,interstitial_pressure_a,interstitial_pressure_b,interstitial_volume_a, & - interstitial_volume_b - - if(capillary_pressure > interstitial_pressure_a)then + if(capillary_pressure > interstitial_pressure_a)then flux_a = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & - interstitial_pressure_a)) * (transit_time/time_loop) + interstitial_pressure_a)) * dt else flux_a = 0.0_dp endif if(capillary_pressure > interstitial_pressure_b)then flux_b = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & - interstitial_pressure_b)) * (transit_time / time_loop) + interstitial_pressure_b)) * dt else flux_b = 0.0_dp endif @@ -178,7 +166,6 @@ subroutine alveolar_capillary_flux(ne) if(interstitial_volume_a + flux_a > interstitial_capacity_a)then excess = flux_a - (interstitial_capacity_a - interstitial_volume_a) interstitial_volume_a = interstitial_capacity_a - write(*,*) 'up1:',interstitial_volume_a alveolar_volume = alveolar_volume + 0.5_dp*excess if((interstitial_volume_b + 0.5_dp*excess) > interstitial_capacity_b)then @@ -186,27 +173,23 @@ subroutine alveolar_capillary_flux(ne) alveolar_volume = alveolar_volume + overflow interstitial_volume_b = interstitial_capacity_b else - interstitial_volume_b = interstitial_volume_b + 0.5*excess + interstitial_volume_b = interstitial_volume_b + 0.5_dp*excess endif else interstitial_volume_a = interstitial_volume_a + flux_a - write(*,*) 'up2:',interstitial_volume_a, flux_a endif interstitial_volume_b = interstitial_volume_b + flux_b - !diffusion = (interstitial_volume_b/(interstitial_capacity-interstitial_capacity_a))*-diffusion_constant + diffusion_constant diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & - (160_dp*1.1_dp))*(transit_time/time_loop) + (160.0_dp*1.1_dp)) * dt interstitial_volume_b = interstitial_volume_b + diffusion interstitial_volume_a = interstitial_volume_a - diffusion - write(*,*) 'up3:',interstitial_volume_a, diffusion - !alveolar_volume = alveolar_volume - 4.45e-5 ! Osmotic int_osm_conc = int_osm_n/interstitial_volume - interstitial_osmotic = i*int_osm_conc*IGC*T - osm_flux = (reflection_coefficient * capillary_SA * (capillary_osmotic - interstitial_osmotic)) & - * (transit_time / time_loop) + interstitial_osmotic = real(i)*int_osm_conc*IGC*T + osm_flux = (lymphatic_properties%reflection_coefficient * capillary_SA * & + (capillary_osmotic - interstitial_osmotic))* dt cap_osm_conc = capillary_osm_n / capillary_volume int_osm_conc = int_osm_n / interstitial_volume @@ -234,15 +217,15 @@ subroutine alveolar_capillary_flux(ne) if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then lymphatic_conductivity = 0.9_dp * capillary_conductivity else - lymphatic_conductivity = ((-1625.1_dp * (interstitial_volume_b / interstitial_capacity_b)**5_dp) + & - (3815.1_dp * (interstitial_volume_b / interstitial_capacity_b)**4_dp) + (-3229.5_dp * & - (interstitial_volume_b / interstitial_capacity_b)**3_dp) + (1258_dp * (interstitial_volume_b / & - interstitial_capacity_b)**2_dp) + (-213.23_dp * (interstitial_volume_b / interstitial_capacity_b)) & + lymphatic_conductivity = ((-1625.1_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + (3815.1_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-3229.5_dp * & + (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (1258.0_dp * (interstitial_volume_b / & + interstitial_capacity_b)**2.0_dp) + (-213.23_dp * (interstitial_volume_b / interstitial_capacity_b)) & + 11.812_dp)* 4.41335e-8_dp !(capillary_conductivity) endif - initial_lymphatic_pressure = ((1.47_dp * sin(time_variable * breathing_function + pi/2_dp)) + & - ((6.82_dp* (interstitial_volume_b / interstitial_capacity_b)**2_dp) + (0.77_dp * (interstitial_volume_b / & + initial_lymphatic_pressure = ((1.47_dp * sin(time_variable * breathing_function + pi/2.0_dp)) + & + ((6.82_dp* (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + (0.77_dp * (interstitial_volume_b / & interstitial_capacity_b)) - 6.52_dp)) if(interstitial_volume.le.0.0_dp)then @@ -250,11 +233,11 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume = 0.0_dp elseif (interstitial_pressure_b > initial_lymphatic_pressure)then if((lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b-initial_lymphatic_pressure)) & - * (transit_time/time_loop) > ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time)) then + * dt > ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time)) then initial_lymphatic_flow = ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time) else initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - & - initial_lymphatic_pressure)) * (transit_time/time_loop) + initial_lymphatic_pressure)) * dt endif else initial_lymphatic_flow = 0.0_dp @@ -267,55 +250,29 @@ subroutine alveolar_capillary_flux(ne) liflowcount = liflowcount + initial_lymphatic_flow initial_osm_n = initial_osm_n + (initial_lymphatic_flow*int_osm_conc) int_osm_n = int_osm_n - (initial_lymphatic_flow*int_osm_conc) - if (initial_lymphatic_volume > 0)then + if (initial_lymphatic_volume > 0.0_dp)then initial_lymph_conc = initial_osm_n/initial_lymphatic_volume else initial_lymph_conc = 0.0_dp endif - !gas_diffusion_restriction = 0.0000152587890625_dp * exp(13.8629436112_dp*(interstitial_saturation/100.0_dp)) + time_sum = time_sum + dt - time2 = time2 + 1 - !write(*,'(''capillary pressure is '',f8.2)')capillary_pressure - !write(*,'(''interstitial pressure is '',f8.2)')interstitial_pressure_b - enddo - - time2 = 1 - - ! if (interstitial_pressure_b > -0.95 or interstitial_pressure_b < -8.05): - ! raise ValueError('interstitial pressure_b outside of range: ' + str(interstitial_pressure_b)) - ! if initial_lymphatic_pressure > 1.05 or initial_lymphatic_pressure < -8.05: - ! raise ValueError('initial lymphatic pressure outside of range: '+ str(initial_lymphatic_pressure) + ' int pressure: ' + str(interstitial_pressure_b)) - - total_flux = total_hydro_flux ! +total_osm_flux - sumuptake = sumuptake + initial_lymphatic_flow - printcount = printcount + 1 - !if (printcount.eq.100)then - if (printcount.eq.1)then - write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time, & + total_flux = total_hydro_flux ! +total_osm_flux + sumuptake = sumuptake + initial_lymphatic_flow + printcount = printcount + 1 + !if (printcount.eq.100)then + write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & initial_lymphatic_volume,total_flux,alveolar_volume - write(*,*) 'at ',time, ' s:' - write(*,*) ' flux per second is ',flux_c/transit_time*1000,' ul' - write(*,*) ' interstitial volume is ',interstitial_volume,' mL' - write(*,*) ' interstitial volume a is ',interstitial_volume_a,' mL' - write(*,*) ' interstitial volume b is ',interstitial_volume_b,' mL' - write(*,*) ' interstitial saturation is ',100*interstitial_saturation,'%' - write(*,*) ' diffusion is ',diffusion,' mL' - write(*,*) ' interstitial pressure is ',interstitial_pressure_b,' mmHg' - write(*,*) ' initial lymphatic flux is ',initial_lymphatic_flow*1000,' uL' - write(*,*) ' initial lymphatic volume is ',initial_lymphatic_volume,' mL' - write(*,*) ' total flux is ',total_flux,' mL' - write(*,*) ' alveolar volume is ',alveolar_volume - printcount = 0 - endif - read(*,*) + !endif + + enddo + + time = time + transit_time + enddo !while - write(*,'(i8,f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') nunit,time, & - flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume call enter_exit(sub_name,2) diff --git a/src/lib/ventilation.f90 b/src/lib/ventilation.f90 index 026c8da8..b5ee8e91 100644 --- a/src/lib/ventilation.f90 +++ b/src/lib/ventilation.f90 @@ -95,6 +95,8 @@ subroutine evaluate_vent sum_tidal = 0.0_dp ! initialise the inspired and expired volumes sum_expid = 0.0_dp last_vol = 0.0_dp + unit_field(nu_Pe_max,:) = -1.0e6_dp + unit_field(nu_Pe_min,:) = 1.0e6_dp !!! set default values for the parameters that control the breathing simulation !!! these should be controlled by user input (showing hard-coded for now) @@ -561,6 +563,12 @@ subroutine tissue_compliance(chest_wall_compliance,undef) !estimate an elastic recoil pressure for the unit unit_field(nu_pe,nunit) = cc/2.0_dp*(3.0_dp*a+b)*(lambda**2.0_dp & -1.0_dp)*exp_term/lambda + ! set minimum and maximum Pe values + unit_field(nu_Pe_max,nunit) = max(unit_field(nu_Pe_max,nunit), & + unit_field(nu_pe,nunit)) + unit_field(nu_Pe_min,nunit) = min(unit_field(nu_Pe_min,nunit), & + unit_field(nu_pe,nunit)) + enddo !nunit call enter_exit(sub_name,2) From aabc7371b59e601f49121d82755e56aa47b12483 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Tue, 26 Oct 2021 19:32:09 +1300 Subject: [PATCH 08/30] updated pressure volume profiles --- src/lib/lymphatics.f90 | 58 +++++++++++++++++++++++++++++------------- 1 file changed, 40 insertions(+), 18 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index cb9d67b7..132f71c7 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -64,9 +64,10 @@ subroutine alveolar_capillary_flux(ne) initial_lymph_conc,initial_osm_n,interstitial_capacity,interstitial_capacity_a,interstitial_capacity_b, & interstitial_osmotic,interstitial_pressure_a,interstitial_pressure_b,interstitial_saturation,interstitial_volume, & int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,max_Pe, & - min_Pe,net_flux, & + min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA + character(len=60) :: sub_name @@ -120,6 +121,16 @@ subroutine alveolar_capillary_flux(ne) printcount = 0 total_hydro_flux = 0.0_dp + intPmax = -1.00_dp + intPmin = -8.00_dp + lymphPmax = 1.00_dp + lymphPmin = -8.00_dp + mx_pe = max_Pe/133.32239_dp + mn_pe = min_Pe/133.32239_dp + fluctuation = ((mx_pe-mn_pe)/2.0_dp) + write(*,'(''max pressure is '',f8.4)')mx_pe + write(*,'(''min pressure is '',f8.4)')mn_pe + write(*,'(''fluctuation '',f8.4)')fluctuation ! dt or n_timesteps should be controlled by the user n_timesteps = 96 dt = transit_time/real(n_timesteps) @@ -131,8 +142,8 @@ subroutine alveolar_capillary_flux(ne) write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') -! do while(time < 32400.0_dp) - do while(time < 3.0_dp) + do while(time < 86400.0_dp) +! do while(time < 3.0_dp) time_sum = dt do while(time_sum < transit_time) interstitial_volume = interstitial_volume_a + interstitial_volume_b @@ -140,12 +151,14 @@ subroutine alveolar_capillary_flux(ne) time_variable = time + time_sum ! calculating flux from capillary into interstitium - interstitial_pressure_a = 1.47_dp * sin(time_variable * breathing_function) + & - ((-3.98_dp * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & - 8.03_dp * (interstitial_volume_a / interstitial_capacity_a) - 6.52_dp) - interstitial_pressure_b = 1.47_dp * sin(time_variable * breathing_function) + & - ((-3.98_dp * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & - 8.03_dp * (interstitial_volume_b / interstitial_capacity_b) - 6.52_dp) + interstitial_pressure_a = fluctuation * sin(time_variable * breathing_function) + & + (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & + ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a) + & + (intPmin + fluctuation)) + interstitial_pressure_b = fluctuation * sin(time_variable * breathing_function) + & + (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & + ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b) + & + (intPmin + fluctuation)) ! pressure determined from saturation equation based of literature (currently linear, but likely not) if(capillary_pressure > interstitial_pressure_a)then @@ -181,7 +194,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume_b = interstitial_volume_b + flux_b diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & - (160.0_dp*1.1_dp)) * dt + (340_dp)) * dt interstitial_volume_b = interstitial_volume_b + diffusion interstitial_volume_a = interstitial_volume_a - diffusion @@ -224,9 +237,9 @@ subroutine alveolar_capillary_flux(ne) + 11.812_dp)* 4.41335e-8_dp !(capillary_conductivity) endif - initial_lymphatic_pressure = ((1.47_dp * sin(time_variable * breathing_function + pi/2.0_dp)) + & - ((6.82_dp* (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + (0.77_dp * (interstitial_volume_b / & - interstitial_capacity_b)) - 6.52_dp)) + initial_lymphatic_pressure = ((fluctuation * sin(time_variable * breathing_function + pi/2.0_dp)) + & + (((lymphPmin-lymphPmax+(fluctuation*2.0_dp))* (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & + (lymphPmin + fluctuation))) if(interstitial_volume.le.0.0_dp)then initial_lymphatic_flow = 0.0_dp @@ -260,18 +273,27 @@ subroutine alveolar_capillary_flux(ne) total_flux = total_hydro_flux ! +total_osm_flux sumuptake = sumuptake + initial_lymphatic_flow - printcount = printcount + 1 + !printcount = printcount + 1 !if (printcount.eq.100)then - write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & - flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume +! write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & +! flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & +! 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & +! initial_lymphatic_volume,total_flux,alveolar_volume !endif enddo time = time + transit_time + printcount = printcount + 1 + if (printcount.eq.100)then + write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & + flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + initial_lymphatic_volume,total_flux,alveolar_volume + printcount = 0 + endif + enddo !while call enter_exit(sub_name,2) From a05c2837a682c01e3a22bcd76c137ac064fbcb5e Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Wed, 27 Oct 2021 14:22:42 +1300 Subject: [PATCH 09/30] testing developments --- src/lib/lymphatics.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index cb9d67b7..6fda4b51 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -84,6 +84,8 @@ subroutine alveolar_capillary_flux(ne) max_Pe = unit_field(nu_Pe_max,nunit) min_Pe = unit_field(nu_Pe_min,nunit) + write(*,*) 'lp',lymphatic_properties%reflection_coefficient + write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' cmH2O'')') & nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/98.0665_dp From 24abbdd52050e75557868d65674459661402082e Mon Sep 17 00:00:00 2001 From: EdeTede Date: Thu, 28 Oct 2021 01:11:26 +1300 Subject: [PATCH 10/30] added run time to arrays and updated lymphatic parameters --- src/lib/arrays.f90 | 3 +++ src/lib/lymphatics.f90 | 17 +++++++++-------- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/lib/arrays.f90 b/src/lib/arrays.f90 index 0812681c..ba43f611 100644 --- a/src/lib/arrays.f90 +++ b/src/lib/arrays.f90 @@ -113,6 +113,7 @@ module arrays real(dp) :: lymphatic_density = 1.0_dp !to be calculated from CT?? real(dp) :: lymphatic_integrity = 1.0_dp ! a measure of how 'leaky' the lymphatic vessels are and prone to backflow real(dp) :: reflection_coefficient = 0.0_dp + real(dp) :: test_time = 86400.0_dp !number of seconds for test to be run end type default_lymphatic_properties !!! arrays that start with default values, updated during simulations/by user @@ -161,6 +162,8 @@ subroutine update_parameter(parameter_name, parameter_value) lymphatic_properties%lymphatic_integrity = parameter_value case('reflection_coefficient') lymphatic_properties%reflection_coefficient = parameter_value + case('test_time') + lymphatic_properties%test_time = parameter_value end select diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 00fffe8e..e981a6f0 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -85,7 +85,8 @@ subroutine alveolar_capillary_flux(ne) max_Pe = unit_field(nu_Pe_max,nunit) min_Pe = unit_field(nu_Pe_min,nunit) - write(*,*) 'lp',lymphatic_properties%reflection_coefficient + write(*,*) 'lp -reflcoef',lymphatic_properties%reflection_coefficient + write(*,*) 'ld',lymphatic_properties%lymphatic_density write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' cmH2O'')') & nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/98.0665_dp @@ -105,7 +106,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume_b = 0.49_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA*3.648_dp + initial_lymphatic_surface_area = capillary_SA*0.64_dp ! initial lymphatic values initial_lymphatic_volume = 0.0_dp @@ -144,7 +145,7 @@ subroutine alveolar_capillary_flux(ne) write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') - do while(time < 86400.0_dp) + do while(time < lymphatic_properties%test_time) ! do while(time < 3.0_dp) time_sum = dt do while(time_sum < transit_time) @@ -232,11 +233,11 @@ subroutine alveolar_capillary_flux(ne) if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then lymphatic_conductivity = 0.9_dp * capillary_conductivity else - lymphatic_conductivity = ((-1625.1_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & - (3815.1_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-3229.5_dp * & - (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (1258.0_dp * (interstitial_volume_b / & - interstitial_capacity_b)**2.0_dp) + (-213.23_dp * (interstitial_volume_b / interstitial_capacity_b)) & - + 11.812_dp)* 4.41335e-8_dp !(capillary_conductivity) + lymphatic_conductivity = ((31.019_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + (-43.674_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-4.609_dp * & + (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (20.561_dp * (interstitial_volume_b / & + interstitial_capacity_b)**2.0_dp) + (1.4381_dp * (interstitial_volume_b / interstitial_capacity_b)) & + + 0.0288_dp)* 4.41335e-8_dp !(capillary_conductivity) endif initial_lymphatic_pressure = ((fluctuation * sin(time_variable * breathing_function + pi/2.0_dp)) + & From 50ffb3dd2edbf5445d0cc83dffe49153d86b2888 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Thu, 28 Oct 2021 19:29:53 +1300 Subject: [PATCH 11/30] updated conductivity parameters --- src/lib/lymphatics.f90 | 54 ++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 28 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index e981a6f0..2d2867c0 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -41,7 +41,7 @@ subroutine alveolar_capillary_flux(ne) ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 0 ! 0 = male, 1 = female - integer,parameter :: au = 32768 ! number of acinar units + integer,parameter :: au = 30676 ! number of acinar units real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp ! Capillary parameters @@ -85,11 +85,11 @@ subroutine alveolar_capillary_flux(ne) max_Pe = unit_field(nu_Pe_max,nunit) min_Pe = unit_field(nu_Pe_min,nunit) - write(*,*) 'lp -reflcoef',lymphatic_properties%reflection_coefficient - write(*,*) 'ld',lymphatic_properties%lymphatic_density + !write(*,*) 'lp -reflcoef',lymphatic_properties%reflection_coefficient + !write(*,*) 'ld',lymphatic_properties%lymphatic_density - write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' cmH2O'')') & - nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/98.0665_dp + write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' mmHg'')') & + nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/133.32239_dp ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g @@ -103,10 +103,10 @@ subroutine alveolar_capillary_flux(ne) interstitial_capacity_a = 0.005_dp*interstitial_capacity interstitial_capacity_b = 0.995_dp*interstitial_capacity interstitial_volume_a = 0.0_dp - interstitial_volume_b = 0.49_dp*interstitial_capacity + interstitial_volume_b = 0.00049_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA*0.64_dp + initial_lymphatic_surface_area = capillary_SA ! initial lymphatic values initial_lymphatic_volume = 0.0_dp @@ -131,8 +131,7 @@ subroutine alveolar_capillary_flux(ne) mx_pe = max_Pe/133.32239_dp mn_pe = min_Pe/133.32239_dp fluctuation = ((mx_pe-mn_pe)/2.0_dp) - write(*,'(''max pressure is '',f8.4)')mx_pe - write(*,'(''min pressure is '',f8.4)')mn_pe + write(*,'(''fluctuation '',f8.4)')fluctuation ! dt or n_timesteps should be controlled by the user n_timesteps = 96 @@ -146,7 +145,6 @@ subroutine alveolar_capillary_flux(ne) &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') do while(time < lymphatic_properties%test_time) -! do while(time < 3.0_dp) time_sum = dt do while(time_sum < transit_time) interstitial_volume = interstitial_volume_a + interstitial_volume_b @@ -162,9 +160,9 @@ subroutine alveolar_capillary_flux(ne) (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b) + & (intPmin + fluctuation)) - + !write(*,'(''Pint: '',f8.4)')interstitial_pressure_b ! pressure determined from saturation equation based of literature (currently linear, but likely not) - if(capillary_pressure > interstitial_pressure_a)then + if(capillary_pressure > interstitial_pressure_a)then flux_a = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - & interstitial_pressure_a)) * dt else @@ -197,7 +195,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume_b = interstitial_volume_b + flux_b diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & - (340_dp)) * dt + (140_dp)) * dt interstitial_volume_b = interstitial_volume_b + diffusion interstitial_volume_a = interstitial_volume_a - diffusion @@ -231,30 +229,30 @@ subroutine alveolar_capillary_flux(ne) !calculating flux from interstitium to initial lymphatics if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then - lymphatic_conductivity = 0.9_dp * capillary_conductivity + lymphatic_conductivity = 1.48_dp * capillary_conductivity else - lymphatic_conductivity = ((31.019_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & - (-43.674_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-4.609_dp * & - (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (20.561_dp * (interstitial_volume_b / & - interstitial_capacity_b)**2.0_dp) + (1.4381_dp * (interstitial_volume_b / interstitial_capacity_b)) & - + 0.0288_dp)* 4.41335e-8_dp !(capillary_conductivity) + !lymphatic_conductivity = ((-595.5_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + !(1140.8_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-716.43_dp * & + !(interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (212.1_dp * (interstitial_volume_b / & + !interstitial_capacity_b)**2.0_dp) + (-20.077_dp * (interstitial_volume_b / interstitial_capacity_b)) & + !- 0.0004_dp)* 4.41335e-8_dp !(capillary_conductivity) + lymphatic_conductivity = ((-57.556_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + (-268.41_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (593.68_dp * & + (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-293.78_dp * (interstitial_volume_b / & + interstitial_capacity_b)**2.0_dp) + (47.955_dp * (interstitial_volume_b / interstitial_capacity_b)) & + - 0.0049_dp)* 4.41335e-8_dp !(capillary_conductivity) endif - initial_lymphatic_pressure = ((fluctuation * sin(time_variable * breathing_function + pi/2.0_dp)) + & - (((lymphPmin-lymphPmax+(fluctuation*2.0_dp))* (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & + initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + & + ((((lymphPmax-lymphPmin-(fluctuation*2.0_dp))* ((interstitial_volume_b / interstitial_capacity_b)**2.0_dp)) + & (lymphPmin + fluctuation))) - + !write(*,'(''Plym: '',f8.4)')initial_lymphatic_pressure if(interstitial_volume.le.0.0_dp)then initial_lymphatic_flow = 0.0_dp interstitial_volume = 0.0_dp elseif (interstitial_pressure_b > initial_lymphatic_pressure)then - if((lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b-initial_lymphatic_pressure)) & - * dt > ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time)) then - initial_lymphatic_flow = ((27.0_dp * capillary_conductivity * capillary_SA) * transit_time) - else - initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - & + initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - & initial_lymphatic_pressure)) * dt - endif else initial_lymphatic_flow = 0.0_dp endif From 3bb91d1450ae32d111690f55c8dffa6f8d38a0a7 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Fri, 29 Oct 2021 16:17:18 +1300 Subject: [PATCH 12/30] another final update on conductivity parameters --- src/lib/lymphatics.f90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 2d2867c0..d00380ae 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -103,7 +103,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_capacity_a = 0.005_dp*interstitial_capacity interstitial_capacity_b = 0.995_dp*interstitial_capacity interstitial_volume_a = 0.0_dp - interstitial_volume_b = 0.00049_dp*interstitial_capacity + interstitial_volume_b = 0.49_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 initial_lymphatic_surface_area = capillary_SA @@ -195,7 +195,7 @@ subroutine alveolar_capillary_flux(ne) interstitial_volume_b = interstitial_volume_b + flux_b diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & - (140_dp)) * dt + (200_dp)) * dt interstitial_volume_b = interstitial_volume_b + diffusion interstitial_volume_a = interstitial_volume_a - diffusion @@ -231,16 +231,16 @@ subroutine alveolar_capillary_flux(ne) if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then lymphatic_conductivity = 1.48_dp * capillary_conductivity else - !lymphatic_conductivity = ((-595.5_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & - !(1140.8_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (-716.43_dp * & - !(interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (212.1_dp * (interstitial_volume_b / & - !interstitial_capacity_b)**2.0_dp) + (-20.077_dp * (interstitial_volume_b / interstitial_capacity_b)) & - !- 0.0004_dp)* 4.41335e-8_dp !(capillary_conductivity) - lymphatic_conductivity = ((-57.556_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & - (-268.41_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (593.68_dp * & - (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-293.78_dp * (interstitial_volume_b / & - interstitial_capacity_b)**2.0_dp) + (47.955_dp * (interstitial_volume_b / interstitial_capacity_b)) & - - 0.0049_dp)* 4.41335e-8_dp !(capillary_conductivity) + !lymphatic_conductivity = ((-57.556_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + ! (-268.41_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (593.68_dp * & + ! (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-293.78_dp * (interstitial_volume_b / & + ! interstitial_capacity_b)**2.0_dp) + (47.955_dp * (interstitial_volume_b / interstitial_capacity_b)) & + ! - 0.0049_dp)* 4.41335e-8_dp !(capillary_conductivity) + lymphatic_conductivity = ((845.87_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & + (-2416.7_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (2388.5_dp * & + (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-922.24_dp * (interstitial_volume_b / & + interstitial_capacity_b)**2.0_dp) + (125.85_dp * (interstitial_volume_b / interstitial_capacity_b)) & + - 0.0067_dp)* 4.41335e-8_dp !(capillary_conductivity) endif initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + & From 5f4bad0e70d92543860ef8cf46b9f8f871583ab1 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Fri, 29 Oct 2021 18:00:25 +1300 Subject: [PATCH 13/30] new bindings to alveolar_flux and lymphatic_transport --- src/bindings/c/lymphatics.c | 13 ++++---- src/bindings/c/lymphatics.f90 | 21 +++++++++---- src/bindings/c/lymphatics.h | 4 +-- src/lib/lymphatics.f90 | 58 ++++++++++++++++++++++++----------- 4 files changed, 64 insertions(+), 32 deletions(-) diff --git a/src/bindings/c/lymphatics.c b/src/bindings/c/lymphatics.c index 8a4bb894..e0d64424 100644 --- a/src/bindings/c/lymphatics.c +++ b/src/bindings/c/lymphatics.c @@ -1,16 +1,17 @@ #include "lymphatics.h" #include "string.h" -void alveolar_capillary_flux_c(int *num_nodes); -void lymphatic_transport_c(); +void alveolar_capillary_flux_c(int *num_nodes, int *write_out); +void lymphatic_transport_c(const char *filename, int *filename_len); -void alveolar_capillary_flux(int num_nodes) +void alveolar_capillary_flux(int num_nodes, int write_out) { - alveolar_capillary_flux_c(&num_nodes); + alveolar_capillary_flux_c(&num_nodes, &write_out); } -void lymphatic_transport() +void lymphatic_transport(const char *filename) { - lymphatic_transport_c(); + int filename_len = (int)strlen(filename); + lymphatic_transport_c(filename, &filename_len); } diff --git a/src/bindings/c/lymphatics.f90 b/src/bindings/c/lymphatics.f90 index ace0f77e..3bcf0054 100644 --- a/src/bindings/c/lymphatics.f90 +++ b/src/bindings/c/lymphatics.f90 @@ -8,16 +8,17 @@ module lymphatics_c !################################################################################### ! !*alveolar_capillary_flux:* - subroutine alveolar_capillary_flux_c(num_nodes) bind(C, name="alveolar_capillary_flux_c") + subroutine alveolar_capillary_flux_c(num_nodes,write_out) bind(C, name="alveolar_capillary_flux_c") use lymphatics,only: alveolar_capillary_flux implicit none integer,intent(in) :: num_nodes + logical,intent(in) :: write_out #if defined _WIN32 && defined __INTEL_COMPILER - call so_alveolar_capillary_flux(num_nodes) + call so_alveolar_capillary_flux(num_nodes,write_out) #else - call alveolar_capillary_flux(num_nodes) + call alveolar_capillary_flux(num_nodes,write_out) #endif end subroutine alveolar_capillary_flux_c @@ -25,14 +26,22 @@ end subroutine alveolar_capillary_flux_c !################################################################################### ! !*lymphatic_transport* - subroutine lymphatic_transport_c() bind(C, name="lymphatic_transport_c") + subroutine lymphatic_transport_c(filename, filename_len) bind(C, name="lymphatic_transport_c") + use iso_c_binding, only: c_ptr + use utils_c, only: strncpy + use other_consts, only: MAX_FILENAME_LEN use lymphatics,only: lymphatic_transport implicit none + integer,intent(in) :: filename_len + type(c_ptr), value, intent(in) :: filename + character(len=MAX_FILENAME_LEN) :: filename_f + + call strncpy(filename_f, filename, filename_len) #if defined _WIN32 && defined __INTEL_COMPILER - call so_lymphatic_transport() + call so_lymphatic_transport(filename_f) #else - call lymphatic_transport() + call lymphatic_transport(filename_f) #endif end subroutine lymphatic_transport_c diff --git a/src/bindings/c/lymphatics.h b/src/bindings/c/lymphatics.h index 3be674d6..83d0e77b 100644 --- a/src/bindings/c/lymphatics.h +++ b/src/bindings/c/lymphatics.h @@ -3,7 +3,7 @@ #include "symbol_export.h" -SHO_PUBLIC void alveolar_capillary_flux(int num_nodes); -SHO_PUBLIC void lymphatic_transport(); +SHO_PUBLIC void alveolar_capillary_flux(int num_nodes, int write_out); +SHO_PUBLIC void lymphatic_transport(const char *filename); #endif /* AETHER_LYMPHATICS_H */ diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index d00380ae..b2b9d7c0 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -13,6 +13,7 @@ module lymphatics use arrays use diagnostics use indices + use other_consts use precision ! sets dp for precision implicit none @@ -31,13 +32,14 @@ module lymphatics !!!############################################################################# - subroutine alveolar_capillary_flux(ne) + subroutine alveolar_capillary_flux(ne,write_out) !*alveolar_capillary_flux:* calculate fluid flux from blood to interstitium !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ALVEOLAR_CAPILLARY_FLUX" :: ALVEOLAR_CAPILLARY_FLUX use other_consts,only: pi integer,intent(in) :: ne + logical,intent(in) :: write_out ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 0 ! 0 = male, 1 = female @@ -88,8 +90,10 @@ subroutine alveolar_capillary_flux(ne) !write(*,*) 'lp -reflcoef',lymphatic_properties%reflection_coefficient !write(*,*) 'ld',lymphatic_properties%lymphatic_density - write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' mmHg'')') & - nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/133.32239_dp + if(write_out)then + write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' mmHg'')') & + nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/133.32239_dp + endif ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g @@ -132,17 +136,19 @@ subroutine alveolar_capillary_flux(ne) mn_pe = min_Pe/133.32239_dp fluctuation = ((mx_pe-mn_pe)/2.0_dp) - write(*,'(''fluctuation '',f8.4)')fluctuation + if(write_out) write(*,'(''fluctuation '',f8.4)')fluctuation ! dt or n_timesteps should be controlled by the user n_timesteps = 96 dt = transit_time/real(n_timesteps) breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) - write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & - &''init lymph|init lymph| total| alveolar|'')') - write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & - &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') + if(write_out) then + write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & + &''init lymph|init lymph| total| alveolar|'')') + write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & + &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') + endif do while(time < lymphatic_properties%test_time) time_sum = dt @@ -286,13 +292,15 @@ subroutine alveolar_capillary_flux(ne) time = time + transit_time - printcount = printcount + 1 - if (printcount.eq.100)then - write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & - flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume - printcount = 0 + if(write_out)then + printcount = printcount + 1 + if (printcount.eq.100)then + write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & + flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + initial_lymphatic_volume,total_flux,alveolar_volume + printcount = 0 + endif endif enddo !while @@ -303,12 +311,14 @@ end subroutine alveolar_capillary_flux !!!############################################################################# - subroutine lymphatic_transport() + subroutine lymphatic_transport(filename) !*lymphatic_transport:* whole system transport !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LYMPHATIC_TRANSPORT" :: LYMPHATIC_TRANSPORT + character(len=MAX_FILENAME_LEN), intent(in) :: filename + ! Local parameters integer :: ne - + character(len=300) :: writefile character(len=60) :: sub_name ! -------------------------------------------------------------------------- @@ -316,12 +326,24 @@ subroutine lymphatic_transport() sub_name = 'lymphatic_transport' call enter_exit(sub_name,1) + if(index(filename, ".oplymph")> 0) then !full filename is given + writefile = filename + else ! need to append the correct filename extension + writefile = trim(filename)//'.oplymph' + endif + + open(10, file=writefile, status='replace') + do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then - call alveolar_capillary_flux(ne) + call alveolar_capillary_flux(ne,.false.) + write(10,'(i8)') ne + write(*,*) ne endif enddo + close(10) + call enter_exit(sub_name,2) end subroutine lymphatic_transport From d65993f6a4a077e978261382e6cc86dcb772fd2a Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Sun, 31 Oct 2021 11:23:05 +1300 Subject: [PATCH 14/30] testing some conditions for stopping the alveolar flux simulation, and writing out whole lung solution --- src/lib/indices.f90 | 7 +++-- src/lib/lymphatics.f90 | 64 +++++++++++++++++++++++++++++++----------- 2 files changed, 52 insertions(+), 19 deletions(-) diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index cde82141..df88ec25 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -35,7 +35,7 @@ module indices integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, & nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,& nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0,nu_Pe_max=0, & - nu_Pe_min=0 + nu_Pe_min=0,nu_flux=0 !indices for gas exchange field ! indices for gasex_field integer,parameter :: num_gx = 12 @@ -67,7 +67,7 @@ module indices public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, & nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,& - nu_conc1,nu_vent,nu_vd,& + nu_conc1,nu_vent,nu_vd,nu_flux, & nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, & @@ -256,13 +256,14 @@ subroutine lymphatic_indices call enter_exit(sub_name,1) ! indices for unit_field - num_nu = 6 + num_nu = 7 nu_perf = 1 nu_blood_press = 2 nu_sa = 3 nu_tt = 4 nu_Pe_max = 5 nu_Pe_min = 6 + nu_flux = 7 call enter_exit(sub_name,2) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index b2b9d7c0..84273e0a 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -69,6 +69,9 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA + real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2) + + logical :: continue character(len=60) :: sub_name @@ -143,14 +146,20 @@ subroutine alveolar_capillary_flux(ne,write_out) breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) - if(write_out) then - write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & - &''init lymph|init lymph| total| alveolar|'')') - write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & - &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') - endif - - do while(time < lymphatic_properties%test_time) +! if(write_out) then +! write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & +! &''init lymph|init lymph| total| alveolar|'')') +! write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & +! &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') +! endf + + continue = .true. + time_av_flux = 0.0_dp + time_av_vol = 0.0_dp + time_av_flow = 0.0_dp + + !do while(time < lymphatic_properties%test_time) + do while (continue) time_sum = dt do while(time_sum < transit_time) interstitial_volume = interstitial_volume_a + interstitial_volume_b @@ -289,22 +298,41 @@ subroutine alveolar_capillary_flux(ne,write_out) !endif enddo - + + time_av_flux(2) = time_av_flux(1) + time_av_flow(2) = time_av_flow(1) + time_av_vol(2) = time_av_vol(1) + time_av_flux(1) = (time_av_flux(1)*time + flux_c)/(time + transit_time) ! time averaged flux/s + time_av_flow(1) = (time_av_flow(1)*time + initial_lymphatic_flow*transit_time)/(time + transit_time) + time_av_vol(1) = (time_av_vol(1)*time + interstitial_volume*transit_time)/(time + transit_time) + time = time + transit_time if(write_out)then printcount = printcount + 1 if (printcount.eq.100)then - write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & + write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume + initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & + time_av_flow(1)*1000.0_dp printcount = 0 endif endif + + if(time.gt.200.0_dp*transit_time)then + if((abs(100.0_dp*(time_av_flux(1)-time_av_flux(2))/time_av_flux(2)).le.0.1_dp).and. & + (abs(100.0_dp*(time_av_flow(1)-time_av_flow(2))/time_av_flow(2)).le.0.1_dp).and. & + (abs(100.0_dp*(time_av_vol(1)-time_av_vol(2))/time_av_vol(2)).le.0.1_dp))then + continue = .false. + endif + endif + if(time .gt.10000.0_dp*transit_time) continue = .false. enddo !while - + + unit_field(nu_flux,nunit) = time_av_flux(1) + call enter_exit(sub_name,2) end subroutine alveolar_capillary_flux @@ -317,7 +345,7 @@ subroutine lymphatic_transport(filename) character(len=MAX_FILENAME_LEN), intent(in) :: filename ! Local parameters - integer :: ne + integer :: ne,np,nunit character(len=300) :: writefile character(len=60) :: sub_name @@ -333,15 +361,19 @@ subroutine lymphatic_transport(filename) endif open(10, file=writefile, status='replace') - + do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then call alveolar_capillary_flux(ne,.false.) - write(10,'(i8)') ne - write(*,*) ne endif enddo + do nunit = 1,num_units + ne = units(nunit) + np = elem_nodes(2,ne) + write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit) + enddo + close(10) call enter_exit(sub_name,2) From e63e9cf063fd6b93dd5c12e5531761d4f0b9bd8f Mon Sep 17 00:00:00 2001 From: EdeTede Date: Wed, 3 Nov 2021 01:11:13 +1300 Subject: [PATCH 15/30] Updated end-point criteria and added outputs --- src/lib/lymphatics.f90 | 46 +++++++++++++++++++++++++++++++----------- 1 file changed, 34 insertions(+), 12 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 84273e0a..b44df5c4 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -69,7 +69,8 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA - real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2) + real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_intsat,nu_time,sat1,sat2,sat3,sat4,sat5, & + nu_lymphflow,nu_av_flux logical :: continue @@ -157,6 +158,11 @@ subroutine alveolar_capillary_flux(ne,write_out) time_av_flux = 0.0_dp time_av_vol = 0.0_dp time_av_flow = 0.0_dp + sat1 = 1.0_dp + sat2 = 2.0_dp + sat3 = 3.0_dp + sat4 = 4.0_dp + sat5 = 5.0_dp !do while(time < lymphatic_properties%test_time) do while (continue) @@ -311,27 +317,42 @@ subroutine alveolar_capillary_flux(ne,write_out) if(write_out)then printcount = printcount + 1 if (printcount.eq.100)then + sat5 = sat4 + sat4 = sat3 + sat3 = sat2 + sat2 = sat1 + sat1 = interstitial_saturation write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & - time_av_flow(1)*1000.0_dp + time_av_flow(1)*1000.0_dp,sat1,sat2,sat3,sat4,sat5,(abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)) printcount = 0 + if(time.gt.200.0_dp*transit_time)then + if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)).le.0.000005)then + continue = .false. + endif + endif endif endif - if(time.gt.200.0_dp*transit_time)then - if((abs(100.0_dp*(time_av_flux(1)-time_av_flux(2))/time_av_flux(2)).le.0.1_dp).and. & - (abs(100.0_dp*(time_av_flow(1)-time_av_flow(2))/time_av_flow(2)).le.0.1_dp).and. & - (abs(100.0_dp*(time_av_vol(1)-time_av_vol(2))/time_av_vol(2)).le.0.1_dp))then - continue = .false. - endif - endif - if(time .gt.10000.0_dp*transit_time) continue = .false. +! if(time.gt.200.0_dp*transit_time)then +! if((abs(100.0_dp*(time_av_flux(1)-time_av_flux(2))/time_av_flux(2)).le.0.01_dp).and. & +! (abs(100.0_dp*(time_av_flow(1)-time_av_flow(2))/time_av_flow(2)).le.0.01_dp).and. & +! (abs(100.0_dp*(time_av_vol(1)-time_av_vol(2))/time_av_vol(2)).le.0.01_dp))then +! continue = .false. +! endif +! endif + + if(time.gt.10000.0_dp*transit_time) continue = .false. enddo !while unit_field(nu_flux,nunit) = time_av_flux(1) + unit_field(nu_intsat,nunit) = interstitial_saturation + !unit_field(nu_time,nunit) = time + !unit_field(nu_av_flux,nunit) = total_flux/time + !unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time call enter_exit(sub_name,2) @@ -348,7 +369,7 @@ subroutine lymphatic_transport(filename) integer :: ne,np,nunit character(len=300) :: writefile character(len=60) :: sub_name - + real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_intsat!,nu_av_flux,nu_lymphflow,nu_time ! -------------------------------------------------------------------------- sub_name = 'lymphatic_transport' @@ -371,7 +392,8 @@ subroutine lymphatic_transport(filename) do nunit = 1,num_units ne = units(nunit) np = elem_nodes(2,ne) - write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit) + write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit)!, & + !unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit),unit_field(nu_time,nunit) enddo close(10) From 4cac6e3b687129f532789f744a96eeb1d351693b Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Wed, 3 Nov 2021 14:54:06 +1300 Subject: [PATCH 16/30] tidying up nu_ indices, adding lymphatic terminal unit export --- src/bindings/c/exports.c | 9 +++++ src/bindings/c/exports.f90 | 25 ++++++++++++++ src/bindings/c/exports.h | 1 + src/lib/exports.f90 | 70 ++++++++++++++++++++++++++++++++++++-- src/lib/indices.f90 | 7 ++-- src/lib/lymphatics.f90 | 4 +-- 6 files changed, 109 insertions(+), 7 deletions(-) diff --git a/src/bindings/c/exports.c b/src/bindings/c/exports.c index 77c46bc5..94b40407 100644 --- a/src/bindings/c/exports.c +++ b/src/bindings/c/exports.c @@ -11,6 +11,7 @@ void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, cons void export_node_field_c(int *nj_field, const char *EXNODEFIELD, int *EXNODEFIELD_LEN, const char *name, int *name_len, const char *field_name, int *field_name_len); void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, const char *name, int *name_len, int *offset_elem, int *offset_node); +void export_terminal_lymphatic_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len); void export_terminal_solution_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len); void export_terminal_perfusion_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len); void export_node_geometry_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len); @@ -54,6 +55,14 @@ void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, export_node_field_c(&nj_field, EXNODEFIELD, &filename_len, name, &name_len, field_name, &field_name_len); } +void export_terminal_lymphatic(const char *EXNODEFILE, const char *name) +{ + int filename_len = strlen(EXNODEFILE); + int name_len = strlen(name); + + export_terminal_lymphatic_c(EXNODEFILE, &filename_len, name, &name_len); +} + void export_terminal_solution(const char *EXNODEFILE, const char *name) { int filename_len = strlen(EXNODEFILE); diff --git a/src/bindings/c/exports.f90 b/src/bindings/c/exports.f90 index aa692f3e..a0b7e80e 100644 --- a/src/bindings/c/exports.f90 +++ b/src/bindings/c/exports.f90 @@ -165,6 +165,31 @@ subroutine export_data_geometry_c(EXDATAFILE, filename_len, name, name_len, offs end subroutine export_data_geometry_c +!!!######################################################################## + + subroutine export_terminal_lymphatic_c(EXNODEFILE, filename_len, name, name_len) bind(C, name="export_terminal_lymphatic_c") + + use iso_c_binding, only: c_ptr + use utils_c, only: strncpy + use exports, only: export_terminal_lymphatic + use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN + implicit none + integer,intent(in) :: filename_len, name_len + type(c_ptr), value, intent(in) :: EXNODEFILE, name + character(len=MAX_FILENAME_LEN) :: filename_f + character(len=MAX_STRING_LEN) :: name_f + + call strncpy(filename_f, EXNODEFILE, filename_len) + call strncpy(name_f, name, name_len) + +#if defined _WIN32 && defined __INTEL_COMPILER + call so_export_terminal_lymphatic(filename_f, name_f) +#else + call export_terminal_lymphatic(filename_f, name_f) +#endif + + end subroutine export_terminal_lymphatic_c + !!!######################################################################## subroutine export_terminal_solution_c(EXNODEFILE, filename_len, name, name_len) bind(C, name="export_terminal_solution_c") diff --git a/src/bindings/c/exports.h b/src/bindings/c/exports.h index 7f8bbdff..03706267 100644 --- a/src/bindings/c/exports.h +++ b/src/bindings/c/exports.h @@ -8,6 +8,7 @@ SHO_PUBLIC void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const SHO_PUBLIC void export_1d_elem_geometry(const char *EXELEMFILE, const char *name); SHO_PUBLIC void export_elem_geometry_2d(const char *EXELEMFILE, const char *name, int offset_elem, int offset_node); SHO_PUBLIC void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, const char *field_name); +SHO_PUBLIC void export_terminal_lymphatic(const char *EXNODEFILE, const char *name); SHO_PUBLIC void export_terminal_solution(const char *EXNODEFILE, const char *name); SHO_PUBLIC void export_terminal_perfusion(const char *EXNODEFILE, const char *name); SHO_PUBLIC void export_terminal_ssgexch(const char *EXNODEFILE, const char *name); diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index c6dbc1ef..a31bfdc1 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -20,8 +20,8 @@ module exports private public export_1d_elem_geometry,export_elem_geometry_2d,export_node_geometry,export_node_geometry_2d,& - export_node_field,export_elem_field,export_terminal_solution,export_terminal_perfusion,& - export_terminal_ssgexch,export_1d_elem_field,export_data_geometry + export_node_field,export_elem_field,export_terminal_lymphatic,export_terminal_solution,& + export_terminal_perfusion,export_terminal_ssgexch,export_1d_elem_field,export_data_geometry contains !!!################################################################ @@ -473,6 +473,72 @@ subroutine export_data_geometry(EXDATAFILE, name, offset) end subroutine export_data_geometry +! +!############################################################################## +! + + subroutine export_terminal_lymphatic(EXNODEFILE, name) + !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_TERMINAL_LYMPHATIC" :: EXPORT_TERMINAL_LYMPHATIC + +!!! Parameters + character(len=MAX_FILENAME_LEN),intent(in) :: EXNODEFILE + character(len=MAX_STRING_LEN),intent(in) :: name + +!!! Local Variables + integer :: len_end,ne,nj,NOLIST,np,np_last,VALUE_INDEX + logical :: FIRST_NODE + + len_end=len_trim(name) + if(num_units.GT.0) THEN + open(10, file=EXNODEFILE, status='replace') + !** write the group name + write(10,'( '' Group name: '',A)') name(:len_end) + FIRST_NODE=.TRUE. + np_last=1 + !*** Exporting Terminal Solution + do nolist=1,num_units + if(nolist.GT.1) np_last = np + ne=units(nolist) + np=elem_nodes(2,ne) + !*** Write the field information + VALUE_INDEX=1 + if(FIRST_NODE)THEN + write(10,'( '' #Fields=3'' )') + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') + do nj=1,3 + if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") + if(nj.eq.2) write(10,'(2X,''y. '')',advance="no") + if(nj.eq.3) write(10,'(2X,''z. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + enddo + + write(10,'('' 2) flux, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + + write(10,'('' 3) int_sat, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + endif !FIRST_NODE + !*** write the node + write(10,'(1X,''Node: '',I12)') np + do nj=1,3 + write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np)) !Coordinates + enddo !njj2 + write(10,'(2X,4(1X,F12.6))') (unit_field(nu_flux,NOLIST)) + write(10,'(2X,4(1X,F12.6))') (unit_field(nu_intsat,nolist)) + + FIRST_NODE=.FALSE. + np_last=np + enddo !nolist (np) + endif !num_nodes + close(10) + + end subroutine export_terminal_lymphatic + ! !############################################################################## ! diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index df88ec25..5071b14d 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -35,7 +35,7 @@ module indices integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, & nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,& nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0,nu_Pe_max=0, & - nu_Pe_min=0,nu_flux=0 + nu_Pe_min=0,nu_flux=0, nu_intsat=0 !indices for gas exchange field ! indices for gasex_field integer,parameter :: num_gx = 12 @@ -67,7 +67,7 @@ module indices public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, & nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,& - nu_conc1,nu_vent,nu_vd,nu_flux, & + nu_conc1,nu_vent,nu_vd,nu_flux,nu_intsat, & nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, & @@ -256,7 +256,7 @@ subroutine lymphatic_indices call enter_exit(sub_name,1) ! indices for unit_field - num_nu = 7 + num_nu = 8 nu_perf = 1 nu_blood_press = 2 nu_sa = 3 @@ -264,6 +264,7 @@ subroutine lymphatic_indices nu_Pe_max = 5 nu_Pe_min = 6 nu_flux = 7 + nu_intsat = 8 call enter_exit(sub_name,2) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index b44df5c4..915dec3a 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -69,7 +69,7 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA - real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_intsat,nu_time,sat1,sat2,sat3,sat4,sat5, & + real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_time,sat1,sat2,sat3,sat4,sat5, & nu_lymphflow,nu_av_flux logical :: continue @@ -369,7 +369,7 @@ subroutine lymphatic_transport(filename) integer :: ne,np,nunit character(len=300) :: writefile character(len=60) :: sub_name - real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_intsat!,nu_av_flux,nu_lymphflow,nu_time + real(dp) :: interstitial_saturation,interstitial_pressure_b !,nu_av_flux,nu_lymphflow,nu_time ! -------------------------------------------------------------------------- sub_name = 'lymphatic_transport' From 1385675d8afc7b32ea819e841ca9e910b62249d6 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Wed, 3 Nov 2021 17:16:29 +1300 Subject: [PATCH 17/30] fixed errors in export_lymphatics and export_terminal_solution, plus moved convergence criteria outside of a write loop in alveolar_capillary_flux --- src/lib/exports.f90 | 6 ++---- src/lib/lymphatics.f90 | 37 +++++++++++++++++-------------------- 2 files changed, 19 insertions(+), 24 deletions(-) diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index a31bfdc1..dbd068f4 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -504,7 +504,6 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) VALUE_INDEX=1 if(FIRST_NODE)THEN write(10,'( '' #Fields=3'' )') - write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') do nj=1,3 if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") @@ -528,8 +527,8 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) do nj=1,3 write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np)) !Coordinates enddo !njj2 - write(10,'(2X,4(1X,F12.6))') (unit_field(nu_flux,NOLIST)) - write(10,'(2X,4(1X,F12.6))') (unit_field(nu_intsat,nolist)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_flux,NOLIST)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_intsat,nolist)) FIRST_NODE=.FALSE. np_last=np @@ -570,7 +569,6 @@ subroutine export_terminal_solution(EXNODEFILE, name) VALUE_INDEX=1 if(FIRST_NODE)THEN write(10,'( '' #Fields=8'' )') - write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') do nj=1,3 if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 915dec3a..f23c7876 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -314,24 +314,24 @@ subroutine alveolar_capillary_flux(ne,write_out) time = time + transit_time - if(write_out)then - printcount = printcount + 1 - if (printcount.eq.100)then - sat5 = sat4 - sat4 = sat3 - sat3 = sat2 - sat2 = sat1 - sat1 = interstitial_saturation + printcount = printcount + 1 + if (printcount.eq.100)then + sat5 = sat4 + sat4 = sat3 + sat3 = sat2 + sat2 = sat1 + sat1 = interstitial_saturation + if(write_out)then write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & time_av_flow(1)*1000.0_dp,sat1,sat2,sat3,sat4,sat5,(abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)) - printcount = 0 - if(time.gt.200.0_dp*transit_time)then - if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)).le.0.000005)then - continue = .false. - endif + endif + printcount = 0 + if(time.gt.200.0_dp*transit_time)then + if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then + continue = .false. endif endif endif @@ -367,6 +367,7 @@ subroutine lymphatic_transport(filename) character(len=MAX_FILENAME_LEN), intent(in) :: filename ! Local parameters integer :: ne,np,nunit + integer :: n_test_run character(len=300) :: writefile character(len=60) :: sub_name real(dp) :: interstitial_saturation,interstitial_pressure_b !,nu_av_flux,nu_lymphflow,nu_time @@ -385,16 +386,12 @@ subroutine lymphatic_transport(filename) do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) call alveolar_capillary_flux(ne,.false.) + np = elem_nodes(2,ne) + write(10,'(i8,5(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit) endif enddo - - do nunit = 1,num_units - ne = units(nunit) - np = elem_nodes(2,ne) - write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit)!, & - !unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit),unit_field(nu_time,nunit) - enddo close(10) From 4259169ed185fd42f075372f208607f4db68b50e Mon Sep 17 00:00:00 2001 From: EdeTede Date: Thu, 4 Nov 2021 16:31:08 +1300 Subject: [PATCH 18/30] updated output parameters --- src/lib/indices.f90 | 8 ++++++-- src/lib/lymphatics.f90 | 44 +++++++++++++++++------------------------- 2 files changed, 24 insertions(+), 28 deletions(-) diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index 5071b14d..892775c8 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -35,7 +35,7 @@ module indices integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, & nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,& nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0,nu_Pe_max=0, & - nu_Pe_min=0,nu_flux=0, nu_intsat=0 + nu_Pe_min=0,nu_flux=0,nu_intsat=0,nu_av_flux=0,nu_lymphflow=0,nu_time=0 !indices for gas exchange field ! indices for gasex_field integer,parameter :: num_gx = 12 @@ -68,7 +68,8 @@ module indices public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, & nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,& nu_conc1,nu_vent,nu_vd,nu_flux,nu_intsat, & - nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min + nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min, & + nu_av_flux,nu_lymphflow,nu_time public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, & ng_p_cap_o2, ng_p_cap_co2,ng_source_o2,ng_source_co2, & @@ -265,6 +266,9 @@ subroutine lymphatic_indices nu_Pe_min = 6 nu_flux = 7 nu_intsat = 8 + nu_av_flux = 9 + nu_lymphflow = 10 + nu_time = 11 call enter_exit(sub_name,2) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 915dec3a..eeaaff86 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -69,8 +69,8 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA - real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_time,sat1,sat2,sat3,sat4,sat5, & - nu_lymphflow,nu_av_flux + real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_time,sat1,sat2,sat3,sat4,sat5!, & + !nu_lymphflow,nu_av_flux logical :: continue @@ -252,11 +252,7 @@ subroutine alveolar_capillary_flux(ne,write_out) if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then lymphatic_conductivity = 1.48_dp * capillary_conductivity else - !lymphatic_conductivity = ((-57.556_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & - ! (-268.41_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (593.68_dp * & - ! (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-293.78_dp * (interstitial_volume_b / & - ! interstitial_capacity_b)**2.0_dp) + (47.955_dp * (interstitial_volume_b / interstitial_capacity_b)) & - ! - 0.0049_dp)* 4.41335e-8_dp !(capillary_conductivity) + lymphatic_conductivity = ((845.87_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & (-2416.7_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (2388.5_dp * & (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-922.24_dp * (interstitial_volume_b / & @@ -295,13 +291,6 @@ subroutine alveolar_capillary_flux(ne,write_out) total_flux = total_hydro_flux ! +total_osm_flux sumuptake = sumuptake + initial_lymphatic_flow - !printcount = printcount + 1 - !if (printcount.eq.100)then -! write(*,'(f12.2, e12.3, f9.3, e12.3, f11.3, f9.3, e12.3, 2(f11.3), e12.3)') time_variable, & -! flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & -! 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & -! initial_lymphatic_volume,total_flux,alveolar_volume - !endif enddo @@ -322,11 +311,11 @@ subroutine alveolar_capillary_flux(ne,write_out) sat3 = sat2 sat2 = sat1 sat1 = interstitial_saturation - write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & - flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & - 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & - time_av_flow(1)*1000.0_dp,sat1,sat2,sat3,sat4,sat5,(abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)) + !write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & + ! flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & + ! 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & + ! initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & + ! time_av_flow(1)*1000.0_dp,sat1,sat2,sat3,sat4,sat5,(abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)) printcount = 0 if(time.gt.200.0_dp*transit_time)then if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)).le.0.000005)then @@ -350,9 +339,12 @@ subroutine alveolar_capillary_flux(ne,write_out) unit_field(nu_flux,nunit) = time_av_flux(1) unit_field(nu_intsat,nunit) = interstitial_saturation - !unit_field(nu_time,nunit) = time - !unit_field(nu_av_flux,nunit) = total_flux/time - !unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time + unit_field(nu_time,nunit) = time + unit_field(nu_av_flux,nunit) = total_flux/time + unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time + write(*,'('' T='',f9.4,'': intsat='',f4.3,'' %; flux='',f8.7,'' ul/s; avFlux='',f8.7,'' ul/s; lyFlo='',f8.7,'' ul/s'')') & + unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& + unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) call enter_exit(sub_name,2) @@ -369,7 +361,7 @@ subroutine lymphatic_transport(filename) integer :: ne,np,nunit character(len=300) :: writefile character(len=60) :: sub_name - real(dp) :: interstitial_saturation,interstitial_pressure_b !,nu_av_flux,nu_lymphflow,nu_time + !real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_av_flux,nu_lymphflow,nu_time ! -------------------------------------------------------------------------- sub_name = 'lymphatic_transport' @@ -385,15 +377,15 @@ subroutine lymphatic_transport(filename) do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then - call alveolar_capillary_flux(ne,.false.) + call alveolar_capillary_flux(ne,.true.) endif enddo do nunit = 1,num_units ne = units(nunit) np = elem_nodes(2,ne) - write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit)!, & - !unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit),unit_field(nu_time,nunit) + write(10,'(i8,4(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit), & + unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit),unit_field(nu_time,nunit) enddo close(10) From 8d8b1bb16d9b5211ed48f6dd4bc755b89c2f4c39 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Thu, 4 Nov 2021 17:57:51 +1300 Subject: [PATCH 19/30] corrected error in export --- src/lib/exports.f90 | 20 +++++++++++++++++++- src/lib/lymphatics.f90 | 28 +++++++++++++++++++++++++--- 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index dbd068f4..9643544b 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -503,7 +503,7 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) !*** Write the field information VALUE_INDEX=1 if(FIRST_NODE)THEN - write(10,'( '' #Fields=3'' )') + write(10,'( '' #Fields=6'' )') write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') do nj=1,3 if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") @@ -521,6 +521,21 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) write(10,'('' 3) int_sat, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + + write(10,'('' 4) transit_time, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + + write(10,'('' 5) cap_pressure, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + + write(10,'('' 6) cap_SA, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 endif !FIRST_NODE !*** write the node write(10,'(1X,''Node: '',I12)') np @@ -529,6 +544,9 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) enddo !njj2 write(10,'(2X,4(1X,e12.6))') (unit_field(nu_flux,NOLIST)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_intsat,nolist)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_tt,nolist)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_blood_press,nolist)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_sa,nolist)) FIRST_NODE=.FALSE. np_last=np diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index f23c7876..e6ed213b 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -366,8 +366,8 @@ subroutine lymphatic_transport(filename) character(len=MAX_FILENAME_LEN), intent(in) :: filename ! Local parameters - integer :: ne,np,nunit - integer :: n_test_run + integer :: i,ne,ne_child,np,nunit + real(dp) :: time_to_run,time_0 character(len=300) :: writefile character(len=60) :: sub_name real(dp) :: interstitial_saturation,interstitial_pressure_b !,nu_av_flux,nu_lymphflow,nu_time @@ -384,6 +384,7 @@ subroutine lymphatic_transport(filename) open(10, file=writefile, status='replace') + call cpu_time(time_0) do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) @@ -392,8 +393,29 @@ subroutine lymphatic_transport(filename) write(10,'(i8,5(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit) endif enddo - + time_to_run = time_0 + call cpu_time(time_0) + time_to_run = time_to_run - time_0 + write(*,*) 'run time=',time_to_run close(10) + + do ne = num_elems,1,-1 + if(elem_field(ne_group,ne).eq.1.0_dp)then + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + elem_field(ne_radius_in0,ne) = unit_field(nu_flux,nunit) + elem_field(ne_radius_out0,ne) = unit_field(nu_intsat,nunit) + else if(elem_field(ne_group,ne).eq.0.0_dp)then ! artery + elem_field(ne_radius_in0,ne) = 0.0_dp + elem_field(ne_radius_out0,ne) = 0.0_dp + do i = 1,elem_cnct(1,0,ne) ! each child branch + ne_child = elem_cnct(1,i,ne) + elem_field(ne_radius_in0,ne) = elem_field(ne_radius_in0,ne) + elem_field(ne_radius_in0,ne_child) + elem_field(ne_radius_out0,ne) = elem_field(ne_radius_out0,ne) + elem_field(ne_radius_out0,ne_child) + enddo + elem_field(ne_radius_in0,ne) = elem_field(ne_radius_in0,ne)/real(elem_cnct(1,0,ne)) + elem_field(ne_radius_out0,ne) = elem_field(ne_radius_out0,ne)/real(elem_cnct(1,0,ne)) + endif + enddo call enter_exit(sub_name,2) From 845987cc660c53ec99ffabe992d5b5d14429a22d Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Thu, 4 Nov 2021 18:23:29 +1300 Subject: [PATCH 20/30] resolving some errors in indices, lymphatics --- src/lib/indices.f90 | 2 +- src/lib/lymphatics.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lib/indices.f90 b/src/lib/indices.f90 index 892775c8..9a980904 100644 --- a/src/lib/indices.f90 +++ b/src/lib/indices.f90 @@ -257,7 +257,7 @@ subroutine lymphatic_indices call enter_exit(sub_name,1) ! indices for unit_field - num_nu = 8 + num_nu = 11 nu_perf = 1 nu_blood_press = 2 nu_sa = 3 diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index e5895c18..adb6c5cf 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -69,7 +69,7 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA - real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),nu_time,sat1,sat2,sat3,sat4,sat5!, & + real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),sat1,sat2,sat3,sat4,sat5!, & !nu_lymphflow,nu_av_flux logical :: continue @@ -342,7 +342,7 @@ subroutine alveolar_capillary_flux(ne,write_out) unit_field(nu_time,nunit) = time unit_field(nu_av_flux,nunit) = total_flux/time unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time - write(*,'('' T='',f9.4,'': intsat='',f4.3,'' %; flux='',f8.7,'' ul/s; avFlux='',f8.7,'' ul/s; lyFlo='',f8.7,'' ul/s'')') & + write(*,'('' T='',e12.3,'': intsat='',e12.3,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.3,'' ul/s'')') & unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) From f045a38920e8248ee8b70f14b4106cbd4b568306 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Tue, 16 Nov 2021 10:39:50 +1300 Subject: [PATCH 21/30] updated parameters --- src/lib/lymphatics.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index adb6c5cf..d8738af6 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -107,14 +107,14 @@ subroutine alveolar_capillary_flux(ne,write_out) ! interstitial values ! interstitial_capacity == maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) - interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp + interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp interstitial_capacity_a = 0.005_dp*interstitial_capacity interstitial_capacity_b = 0.995_dp*interstitial_capacity interstitial_volume_a = 0.0_dp interstitial_volume_b = 0.49_dp*interstitial_capacity alveolar_volume = 0.0_dp liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA + initial_lymphatic_surface_area = capillary_SA ! initial lymphatic values initial_lymphatic_volume = 0.0_dp @@ -250,7 +250,7 @@ subroutine alveolar_capillary_flux(ne,write_out) !calculating flux from interstitium to initial lymphatics if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then - lymphatic_conductivity = 1.48_dp * capillary_conductivity + lymphatic_conductivity = 1.48_dp * 4.41335e-8_dp!capillary_conductivity else lymphatic_conductivity = ((845.87_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & @@ -311,11 +311,10 @@ subroutine alveolar_capillary_flux(ne,write_out) sat2 = sat1 sat1 = interstitial_saturation if(write_out)then - write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), 4(e12.4))') time_variable, & + write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), f12.8,f12.4)') time_variable, & flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, & 100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, & - initial_lymphatic_volume,total_flux,alveolar_volume,time_av_flux(1)*1000.0_dp,time_av_vol(1), & - time_av_flow(1)*1000.0_dp,sat1,sat2,sat3,sat4,sat5,(abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5)-sat1)) + initial_lymphatic_volume,total_flux,alveolar_volume,interstitial_pressure_b endif printcount = 0 if(time.gt.200.0_dp*transit_time)then @@ -342,7 +341,7 @@ subroutine alveolar_capillary_flux(ne,write_out) unit_field(nu_time,nunit) = time unit_field(nu_av_flux,nunit) = total_flux/time unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time - write(*,'('' T='',e12.3,'': intsat='',e12.3,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.3,'' ul/s'')') & + write(*,'('' T='',e12.3,'': intsat='',e12.6,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.6,'' ul/s'')') & unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) @@ -382,7 +381,8 @@ subroutine lymphatic_transport(filename) nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) call alveolar_capillary_flux(ne,.false.) np = elem_nodes(2,ne) - write(10,'(i8,5(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_flux,nunit),unit_field(nu_intsat,nunit) + write(10,'(i8,6(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_av_flux,nunit),unit_field(nu_intsat,nunit), & + unit_field(nu_lymphflow,nunit) endif enddo time_to_run = time_0 From 56a452ba31a695ea88b5dc6b3d74bc1e6498b003 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Tue, 16 Nov 2021 17:34:05 +1300 Subject: [PATCH 22/30] lymphatic flow in export --- src/lib/exports.f90 | 8 +++++++- src/lib/lymphatics.f90 | 1 - 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index 9643544b..e6b5f367 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -503,7 +503,7 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) !*** Write the field information VALUE_INDEX=1 if(FIRST_NODE)THEN - write(10,'( '' #Fields=6'' )') + write(10,'( '' #Fields=7'' )') write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')') do nj=1,3 if(nj.eq.1) write(10,'(2X,''x. '')',advance="no") @@ -536,6 +536,11 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) write(10,'('' 6) cap_SA, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 + VALUE_INDEX=VALUE_INDEX+1 + + write(10,'('' 6) lymf_flow, field, rectangular cartesian, #Components=1'')') + write(10,'(2X,''1. '')',advance="no") + write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 endif !FIRST_NODE !*** write the node write(10,'(1X,''Node: '',I12)') np @@ -547,6 +552,7 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_tt,nolist)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_blood_press,nolist)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_sa,nolist)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_lymphflow,nolist)) FIRST_NODE=.FALSE. np_last=np diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index d8738af6..595ef859 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -331,7 +331,6 @@ subroutine alveolar_capillary_flux(ne,write_out) ! continue = .false. ! endif ! endif - if(time.gt.10000.0_dp*transit_time) continue = .false. enddo !while From 3a9f7fec8d616d5f527eaa9f0d80d31289247b97 Mon Sep 17 00:00:00 2001 From: EdeTede Date: Tue, 18 Jan 2022 14:15:06 +1300 Subject: [PATCH 23/30] documentation --- src/lib/lymphatics.f90 | 47 +++++++++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 595ef859..05b4cff0 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -43,19 +43,21 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 0 ! 0 = male, 1 = female - integer,parameter :: au = 30676 ! number of acinar units - real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp + !sex only determines the weight and therefore size of the lung. Should be updated based on CT + integer,parameter :: au = 30676 ! number of acinar units - currently manually updated but should be imported directly. + real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp !placeholders values + !bodymass used to calculate capillary volume ! Capillary parameters - real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp ! obtained from literature - Parker range in cm H2O + real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp !ml/s/mmHg obtained from Parker (6e-8 cm H2O) ! Osmotic pressure parameters - real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! mol/L in blood plasma + real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! this number is currenlty g/L, likely needs to change for osmolar to work... real(dp),parameter :: IGC = 62.3630_dp ! ideal gas constant in mmHg.L.mol-1.K-1 real(dp),parameter :: T = 310.0_dp ! temperature in Kelvin - based on 37C blood temp ! Simulation parameters - real(dp),parameter :: breathing_rate = 15.0_dp + real(dp),parameter :: breathing_rate = 15.0_dp !constant but should be imported directly from ventilation model ! Local variables @@ -85,7 +87,7 @@ subroutine alveolar_capillary_flux(ne,write_out) ! get information for the unit fron unit_field ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) - capillary_pressure = unit_field(nu_blood_press,nunit)/133.0_dp + capillary_pressure = unit_field(nu_blood_press,nunit)/133.0_dp !converted to mmHg transit_time = unit_field(nu_tt,nunit) capillary_SA = unit_field(nu_sa,nunit) max_Pe = unit_field(nu_Pe_max,nunit) @@ -100,21 +102,23 @@ subroutine alveolar_capillary_flux(ne,write_out) endif ! Calculated values - lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g; gives female lung weight of 639g and male of 840g - open_capillaries = 1.0_dp/6.0_dp - capillary_volume_raw = body_mass/0.3474_dp ! Gehr 213 ml and body mass of 74 kg - rough estimate? + lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g;female lung weight of 639g and male of 840g - should be updated from CT + open_capillaries = 1.0_dp/6.0_dp !based on open capillaries at rest. Should be solved for by perfusion model + capillary_volume_raw = body_mass/0.3474_dp !gives a capillary volume of 213 ml (Gehr 1978) based on a body mass of 74 kg capillary_volume = (capillary_volume_raw*open_capillaries)/real(num_units) ! interstitial values ! interstitial_capacity == maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) - interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp - interstitial_capacity_a = 0.005_dp*interstitial_capacity - interstitial_capacity_b = 0.995_dp*interstitial_capacity + interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp !based on lung mass which should be obtained from CT + interstitial_capacity_a = 0.005_dp*interstitial_capacity !arbitrarily sized - needs further studies on the capillary-lymph interface + interstitial_capacity_b = 0.995_dp*interstitial_capacity !arbitrarily sized - needs further studies on the capillary-lymph interface interstitial_volume_a = 0.0_dp - interstitial_volume_b = 0.49_dp*interstitial_capacity - alveolar_volume = 0.0_dp + interstitial_volume_b = 0.48_dp*interstitial_capacity !assumed to be around 48% saturated at rest + alveolar_volume = 0.0_dp !alveolar volume likely greater at rest, but is lost to respiration - further information needed to put in model liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA + initial_lymphatic_surface_area = capillary_SA !initial lymphatic SA assumed to be the same as capillary SA as no other indication. + !both initial lymphatic SA and initial lymphatic hydraulic conductivity make up the filtration coefficient and currently it is the + !hydraulic conductivity that is adjusted to compensate ! initial lymphatic values initial_lymphatic_volume = 0.0_dp @@ -167,6 +171,9 @@ subroutine alveolar_capillary_flux(ne,write_out) !do while(time < lymphatic_properties%test_time) do while (continue) time_sum = dt + if(transit_time.le.0.001)then + continue =.false. + endif do while(time_sum < transit_time) interstitial_volume = interstitial_volume_a + interstitial_volume_b interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% @@ -177,6 +184,9 @@ subroutine alveolar_capillary_flux(ne,write_out) (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a) + & (intPmin + fluctuation)) + !arbitrarily defined mathematical relationship between interstitial volume and pressure: (same for a and b) + !interstitial pressure changes a lot at low volumes with a small volume change, but at high volumes + !a large volume change is needed to cause a small change in pressure interstitial_pressure_b = fluctuation * sin(time_variable * breathing_function) + & (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + & ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b) + & @@ -250,7 +260,8 @@ subroutine alveolar_capillary_flux(ne,write_out) !calculating flux from interstitium to initial lymphatics if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then - lymphatic_conductivity = 1.48_dp * 4.41335e-8_dp!capillary_conductivity + lymphatic_conductivity = 1.48_dp * 4.41335e-8_dp!all calculated does as a function of capillary_conductivity + !no information on the size of pores or similar for lympatic conductivity so assumed to be similar to capillary. else lymphatic_conductivity = ((845.87_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + & @@ -263,6 +274,8 @@ subroutine alveolar_capillary_flux(ne,write_out) initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + & ((((lymphPmax-lymphPmin-(fluctuation*2.0_dp))* ((interstitial_volume_b / interstitial_capacity_b)**2.0_dp)) + & (lymphPmin + fluctuation))) + !arbitrarily defined mathematical relationship to show that lymphatic pressure does not change much at low volumes with a + !large volume change, but at high volumes only a small volume change is needed to cause a large change in pressure !write(*,'(''Plym: '',f8.4)')initial_lymphatic_pressure if(interstitial_volume.le.0.0_dp)then initial_lymphatic_flow = 0.0_dp @@ -318,7 +331,7 @@ subroutine alveolar_capillary_flux(ne,write_out) endif printcount = 0 if(time.gt.200.0_dp*transit_time)then - if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then + if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then continue = .false. endif endif From 746426d5f211ae74ab1a261efad33d2da95c3cef Mon Sep 17 00:00:00 2001 From: EdeTede Date: Wed, 19 Jan 2022 15:02:57 +1300 Subject: [PATCH 24/30] removed junk - updated doc --- src/lib/.lymphatics.f90.swp | Bin 0 -> 28672 bytes src/lib/lymphatics.f90 | 49 ++++++++---------------------------- 2 files changed, 10 insertions(+), 39 deletions(-) create mode 100644 src/lib/.lymphatics.f90.swp diff --git a/src/lib/.lymphatics.f90.swp b/src/lib/.lymphatics.f90.swp new file mode 100644 index 0000000000000000000000000000000000000000..cde3c7703994c0603fce4c8e18ec7724c8aa65c3 GIT binary patch literal 28672 zcmeI4d5k32UB}zO#)jY+J5E+YL3wRwJ>9*1&9QT=z1Zv7wY~D#9(iVigQc{pyQ^o4 zt*&ZURnHMt#33LBi9}W;0>oCr5d`5#!eJv}5y=jMC^nHvAQ2WsBntmQh(KVhfWzna zUR~WgvoV$>qPSZ6PFKJ4_kQp9&iekQ+kX4md(={Mp~UBhOQkzL_^Z~N+Ydkb%u1`U_r>9f93P#E zPl_$R3k?(+C^S%LpwK{}fkFd?1_})n8YnbSXrR!*_o4=@eyOzZgQe0{ucvjq|4+{U z-?*k!`UvJV$$;0iOaV zzzR45p8A1OX$YgMg4sovNfW^(@y_lw=1Td(Vj~;keuJtf5_>H9Agy)9@L~ z-LZYTrjv@<8Mx!U8C_5;LgQ8YWTp1Dt-$IbV#T$+lZ)-DI&nfJf@R5+eDqJ(-*RY} zB{-O60Gp8X#0<$v$^ zuKWzGXjcBr$)9=ovmk$tm+R^#UAC@W7qVjoj@5PTdgMozt80%RxA*&jpAYauEDbH! zGo>dcT+a`HT>QqVtU7;Wdue|;c34l#lLyZf?eExuRz^4;>bm$Ik0R?l0@|S&V7m22 z1tU9f{Jwc6VQ4B7+LaX7M?h!CiT3Nr zBeG9bo0U)fKA9$qPda0|+BCd7nI4N62Sy(U2tQT1X_wYW(sC1|J8!aJ?8E9wGbV~n zLZ`(`!@Gx9sWqhS!eIgt4QoksaMC<#n~q|39l&GjFhfaKOh=nuny581?KpBepc7cT zV&Itx^qs+A95Tt(`^>$2JtlRrwLRLYzq;nDSKqS2PMWyGuOf;>T?B~uX%ta=)lEs@Yiz}5r+Oh+gp`5zwNKrk1?8R$V-P3DY z7}=CL3P`LY^{6hZHDBhoWC{s&sF~i9rqo#-n8O(wQbE=T>`vCb+2h%&R#qF%aO`S> zRB5*GQsr$@#Qbbqm!ZlkvWH9vV%4m2p1Rw1clc_kt`*vS#hNp}w5s&ts*`GQ4;J+7 zovReRIM0w&#w=)ezTO6X}hdBi9OtmM*od?lYraB^LiOwUkwypJS+cp}D}m;<(*dHceE1 zH?kZLU9mzBmFv)LMw<<_ZUyWX6#JYlt=!&F-8FY#EZAzYQ_4%E_~`h9RlIFGTiX#* z%q|%tsH@#1Xz93QE)Ol&R+fi2wtft3tix83uliz9Oz?~cy$1=MciGAx=x3Q+>1=nu zsqVBrX@R8DDRi?|H)RXdw?ik8g-?nLZ1mjYaoaIQ91n^~>kaKa=`U=d&9_KYSF=1B z>2y=e>UM0+MTd+EcWmeMp(^$cyPGk)oPD*+#!_`{H4g2d(YFVp^joXO)Z;CKo;=y$ zIiB7D><08k#~(+!of6+*&CUDI z-hOWL%!WF>wsz)$_51GMte;vtdurpv4IB4$POZKD%zbO8?(eLgT0gtS!#g^6tetz{ z1{unS#2GEr(%CDsJ(&wn)Yw$JJsIS-BEX?*2{z^y4@@89Pgc_Lkn-4G1+-gqdPuPq z&#tHXkpDG^twDPGGx7uJojCFkSxaS71kMf`AZ_`q zv{lw9DOLjR@%v+pGC!y0f}{P>bjsXJsY`H%BV=Qj?fKds$St<)d8^H|<=b{$z2^}sh8_p_r5 zr=3X1>~$$yP`hKyxI1l!y}%jC8hWCZL~XJhtMpG2T3b82dgk7ZGg-RY?SW$tCd*I< z$}Pv4)Iu8J2>lN*n#`QTQry%uRalbq{|aaMA&~Qb@_zhV{QeU7Tkx~sy+G2qp!gIT zC^S%LpwK{}fkFd?1_})n8YnbSXrRzQp@BjJm(YOCUB$&m@K^uxh}#6IDm(dS_X3^# z%@OxEa{m7V-ci4Y_t|p(Z^Yt{Ip_Ze@9BRLybD|gKFGWH5%>}C7rcKT0vD`4JO`cu9|s=;8{l4W z6kGu=2hZ^yejBvFmy!GTfR_IfsxP7n4HO#qF4TbHlU>?>_Yhyo+y}!A&^Oq;An=vE zd1GQ#F(IX3!xU^7d6!hMVG1^kCWs_2$*1WiX?VegndH$FY#0_3dH4MP&4y9saN8fd zeFfW4*cS4NGo#KxDTsJy?1h>c5vPPoV7Y{9LPK|7zYgWqY@ip&f>6E4?Jd)lx5c(^e4ol(;F-UVin(C|vhvqzVkZt(9^ob0CkgMR-&Cuq&1kwa8HSdpi9!1O|0&Mzmvg3<^FPOY{rNKI|4)G52VHOzI0~)- zkC12q)Ib$n3mzuFlu-qbQQo8AJ>W9%JKO^dK^t5RzQ}#R6X4gu2>4(f)Ib$n3$6su za$oRC@P4oZ9t5|7TR;sQ0WWfw@L%Ao;B(*)!LNXQ@D9)hE%2}0Gdux)4eWymEP~_U zN5NCvEBrKg8~8f+3eSVT1D^pO2akXq&h`<*3G4M_9 zSN6cYU=>^eehVLfCivsa@UH+L0b-W94@@kG!f;KF3I8Oa5?9U0$)`&aHgpDS?^q?5 zzBZ=~v|6oI8p)%oDan)3d(9~7SgMsHLqU3FoM~ch$EqSlj>>JSuT)-M)vu(!S|a#V z#S?NT(-@LNi_4Tmg`wkBCJ5(*WpvPl(>mR(a>^GIs*vpO4Oc9iaH#(x##quBrQ`tz z4I$7cq$-4y-?T&rLmbbgoM=!+ktWYXOgzyY(m^a4DwmfT$6`+Og$*SXEas4QW47Hs zdRV=2=?B({a*rJ(OQ)NZ*7w!!w&U7H;WSU(GU1_;nk=`H&LGL9IQYjd3hHEYV7Z}f zmLnpZK_$=dZ@P9<^H!@c485d)M470nZ)l2%AgYp-9IrD;gq9HLQU(-WFZc6;6=N1Z^siy(1Ab2PHSyo|GF!YzwHn~XP+?zsGrrnSNgV#Ih%xKc_peePzDNf`Mu z@i4`V>nR&u&&<>pPbMh8V!FIGa;mM$qUqCCmQq;gIpgXJ9XW63ExjwrQQ%GDTmTtX z7E99RGDfdpDdr=yoq-rcAf(5>n%E>~D0DDbZZ6ROXfFvhEzbo=lPlV} zrJ&I z?8$RcaJDYTt`w6g(zVDq0$n0fTy}`oM&(p%S;ovEpoZr`ZBw_yA*t~vxS07LdqpJ8 z$)?YA;m9ct-v!g`Y??Lf^g~O@`^_LlCTQ%>m2(r3=4x@1i6HiUlJIPDrN=WZ8W}rX zon%3*uQm}+b$XX$v~WO9*2SF5D1UM<%9 zc-8JM(((D%wiIA*OuPsBo1o%1wUKw%ZW_2D!$D9nVc(SZV#EwjD1@?wD(J5U6Xl7 zsu9=&x98h~LrRfoaaZthU?mkaeutC0S$&P#FqO=qS^UG-?G7<{$8V$N0~6wX8atEA z`~NR-4y|!cmGl3-@%jJfIp?1M;v4W*@KMhAcY)V~&vVBAUGNat1h;`VfUChbIPZTJ z{3ZA#_)YMW;2v-cTn}F44&Wca^Wd}KaqvsvL*U(D8?1vRa0I--UBKtTUxD|7O>ip^ zpZ@E>-*6A`Mero(fVY6BIrIMv*aN%30(XKYxE{O?e3NtkSHQF2Pr>8hVXzH0KoeXC zt^j|}y@32LfH#AG=Dhz^@EhO^coyB=53~+r{S{w@2EMm6AoE!6lid7pH;_9I6t`2k zeN?g*7PnJuQ#ks^&i`^BUffRUf6WnZe~a6xWMV09r;6LD{B>Kh=bot)w^Quob7%JA zcB&_Def2-II`m#Cy)(!uuDG3gxwl_jmc(!WjOi7(QxKRY?n9=plhD(vemiwC;~TWv z&H06;1uhtq^)mY|^Up`%JIu|_L&MP5K0;G4Bns|+mKM}`x%@Hz7=ha=ZuKJj*3$n1 D4*BkJ literal 0 HcmV?d00001 diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 05b4cff0..819da7b1 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -44,12 +44,9 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 0 ! 0 = male, 1 = female !sex only determines the weight and therefore size of the lung. Should be updated based on CT - integer,parameter :: au = 30676 ! number of acinar units - currently manually updated but should be imported directly. - real(dp),parameter :: height = 175.0_dp, weight = 75.0_dp, body_mass = 74.0_dp !placeholders values - !bodymass used to calculate capillary volume ! Capillary parameters - real(dp),parameter :: capillary_conductivity = 4.41335e-8_dp !ml/s/mmHg obtained from Parker (6e-8 cm H2O) + real(dp),parameter :: capillary_conductivity = 4.41335e-8 !ml/s/Pa obtained from Parker (6e-8 cm H2O) ! Osmotic pressure parameters real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! this number is currenlty g/L, likely needs to change for osmolar to work... @@ -71,8 +68,7 @@ subroutine alveolar_capillary_flux(ne,write_out) min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, & open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, & time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA - real(dp) :: time_av_flow(2),time_av_flux(2),time_av_vol(2),sat1,sat2,sat3,sat4,sat5!, & - !nu_lymphflow,nu_av_flux + real(dp) :: sat1,sat2,sat3,sat4,sat5 logical :: continue @@ -87,14 +83,12 @@ subroutine alveolar_capillary_flux(ne,write_out) ! get information for the unit fron unit_field ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) - capillary_pressure = unit_field(nu_blood_press,nunit)/133.0_dp !converted to mmHg + capillary_pressure = unit_field(nu_blood_press,nunit)/133.32239_dp !converted to mmHg transit_time = unit_field(nu_tt,nunit) capillary_SA = unit_field(nu_sa,nunit) max_Pe = unit_field(nu_Pe_max,nunit) min_Pe = unit_field(nu_Pe_min,nunit) - !write(*,*) 'lp -reflcoef',lymphatic_properties%reflection_coefficient - !write(*,*) 'ld',lymphatic_properties%lymphatic_density if(write_out)then write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' mmHg'')') & @@ -103,9 +97,10 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g;female lung weight of 639g and male of 840g - should be updated from CT - open_capillaries = 1.0_dp/6.0_dp !based on open capillaries at rest. Should be solved for by perfusion model - capillary_volume_raw = body_mass/0.3474_dp !gives a capillary volume of 213 ml (Gehr 1978) based on a body mass of 74 kg + open_capillaries = 1.0_dp/6.0_dp !based on open capillaries at rest. Should be solved for by perfusion model? + capillary_volume_raw = 213.00_dp !Gehr 1978 based on having a body mass of 74 kg - should be unique to each person capillary_volume = (capillary_volume_raw*open_capillaries)/real(num_units) + !Only used for the osmotic model at the moment (which isn't operational) can volume be obtained elsewhere? ! interstitial values ! interstitial_capacity == maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002) @@ -151,17 +146,8 @@ subroutine alveolar_capillary_flux(ne,write_out) breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) -! if(write_out) then -! write(*,'(8X,''Time|'',5X,''flux/s| intrstl| a.intrstl| b.intrstl| intrstl|'',X, & -! &''init lymph|init lymph| total| alveolar|'')') -! write(*,'(9X,''(s)|'',7X,''(uL)| vol(mL)| vol(mL)| vol(mL)| sat(%)|'',3X, & -! &''flux(uL)| vol(mL)| flux(mL)| vol(?)|'')') -! endf - continue = .true. - time_av_flux = 0.0_dp - time_av_vol = 0.0_dp - time_av_flow = 0.0_dp + sat1 = 1.0_dp sat2 = 2.0_dp sat3 = 3.0_dp @@ -260,7 +246,7 @@ subroutine alveolar_capillary_flux(ne,write_out) !calculating flux from interstitium to initial lymphatics if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then - lymphatic_conductivity = 1.48_dp * 4.41335e-8_dp!all calculated does as a function of capillary_conductivity + lymphatic_conductivity = 1.48_dp * 4.41335e-8 !all calculated does as a function of capillary_conductivity !no information on the size of pores or similar for lympatic conductivity so assumed to be similar to capillary. else @@ -268,7 +254,7 @@ subroutine alveolar_capillary_flux(ne,write_out) (-2416.7_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (2388.5_dp * & (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-922.24_dp * (interstitial_volume_b / & interstitial_capacity_b)**2.0_dp) + (125.85_dp * (interstitial_volume_b / interstitial_capacity_b)) & - - 0.0067_dp)* 4.41335e-8_dp !(capillary_conductivity) + - 0.0067_dp)* 4.41335e-8 !(capillary_conductivity) endif initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + & @@ -307,13 +293,6 @@ subroutine alveolar_capillary_flux(ne,write_out) enddo - time_av_flux(2) = time_av_flux(1) - time_av_flow(2) = time_av_flow(1) - time_av_vol(2) = time_av_vol(1) - time_av_flux(1) = (time_av_flux(1)*time + flux_c)/(time + transit_time) ! time averaged flux/s - time_av_flow(1) = (time_av_flow(1)*time + initial_lymphatic_flow*transit_time)/(time + transit_time) - time_av_vol(1) = (time_av_vol(1)*time + interstitial_volume*transit_time)/(time + transit_time) - time = time + transit_time printcount = printcount + 1 @@ -336,19 +315,11 @@ subroutine alveolar_capillary_flux(ne,write_out) endif endif endif - -! if(time.gt.200.0_dp*transit_time)then -! if((abs(100.0_dp*(time_av_flux(1)-time_av_flux(2))/time_av_flux(2)).le.0.01_dp).and. & -! (abs(100.0_dp*(time_av_flow(1)-time_av_flow(2))/time_av_flow(2)).le.0.01_dp).and. & -! (abs(100.0_dp*(time_av_vol(1)-time_av_vol(2))/time_av_vol(2)).le.0.01_dp))then -! continue = .false. -! endif -! endif + if(time.gt.10000.0_dp*transit_time) continue = .false. enddo !while - unit_field(nu_flux,nunit) = time_av_flux(1) unit_field(nu_intsat,nunit) = interstitial_saturation unit_field(nu_time,nunit) = time unit_field(nu_av_flux,nunit) = total_flux/time From 1efd9e8c15153aec357568a5f935872a2c164fde Mon Sep 17 00:00:00 2001 From: EdeTede <75711390+EdeTede@users.noreply.github.com> Date: Thu, 20 Jan 2022 17:56:02 +1300 Subject: [PATCH 25/30] Future Directions --- src/lib/lymphatics.f90 | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 05b4cff0..c72030ce 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -428,3 +428,19 @@ end subroutine lymphatic_transport !!!############################################################################# end module lymphatics + + +!FUTURE DIRECTIONS +!input a constant to account for difference between current values and expected values + !Model appeared to be working within the range of the literature but is likely off by a factor of 1000 due to nl to ul conversion error. + !Need to check that outputted units are correct - most things are in ml and mmHg +!Lymphatic network tree + !currently all lymph is returned to the circulation immediately, in reality it moves up a tree of lymphatics against a pressure gradient + !would require excessive modelling perhaps +!impairment of gas diffusion caused by high interstitial saturation + !unclear at what level this would occur +!alveolar flooding changes + !alveolar flooding should be able to move between adjacent compartments + !alveolar fluid should be removed via respiration naturally and therefore should always occur naturally at some low level +!individuality needs to be added in line with the other modules + !currently operates only on preset male/female values From e7de5eaa702a1be5523343c1eb74868c29d4f022 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Thu, 15 Sep 2022 13:00:03 +1200 Subject: [PATCH 26/30] updating some of the outputs from code for including in paper --- src/lib/exports.f90 | 4 ++-- src/lib/lymphatics.f90 | 12 +++++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/lib/exports.f90 b/src/lib/exports.f90 index e6b5f367..a91f343f 100644 --- a/src/lib/exports.f90 +++ b/src/lib/exports.f90 @@ -538,7 +538,7 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 VALUE_INDEX=VALUE_INDEX+1 - write(10,'('' 6) lymf_flow, field, rectangular cartesian, #Components=1'')') + write(10,'('' 7) lymf_flow, field, rectangular cartesian, #Components=1'')') write(10,'(2X,''1. '')',advance="no") write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0 endif !FIRST_NODE @@ -547,7 +547,7 @@ subroutine export_terminal_lymphatic(EXNODEFILE, name) do nj=1,3 write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np)) !Coordinates enddo !njj2 - write(10,'(2X,4(1X,e12.6))') (unit_field(nu_flux,NOLIST)) + write(10,'(2X,4(1X,e12.6))') (unit_field(nu_av_flux,NOLIST)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_intsat,nolist)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_tt,nolist)) write(10,'(2X,4(1X,e12.6))') (unit_field(nu_blood_press,nolist)) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 20bfe705..3126e091 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -47,6 +47,7 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Capillary parameters real(dp),parameter :: capillary_conductivity = 4.41335e-8 !ml/s/Pa obtained from Parker (6e-8 cm H2O) + !real(dp),parameter :: capillary_conductivity = 9.26e-8 !ml/s/Pa obtained from Parker (6e-8 cm H2O) ! Osmotic pressure parameters real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! this number is currenlty g/L, likely needs to change for osmolar to work... @@ -324,9 +325,9 @@ subroutine alveolar_capillary_flux(ne,write_out) unit_field(nu_time,nunit) = time unit_field(nu_av_flux,nunit) = total_flux/time unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time - write(*,'('' T='',e12.3,'': intsat='',e12.6,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.6,'' ul/s'')') & - unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& - unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) + !write(*,'('' T='',e12.3,'': intsat='',e12.6,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.6,'' ul/s'')') & + ! unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& + ! unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) call enter_exit(sub_name,2) @@ -364,8 +365,9 @@ subroutine lymphatic_transport(filename) nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) call alveolar_capillary_flux(ne,.false.) np = elem_nodes(2,ne) - write(10,'(i8,6(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_av_flux,nunit),unit_field(nu_intsat,nunit), & - unit_field(nu_lymphflow,nunit) + write(10,'(i8,11(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_av_flux,nunit),unit_field(nu_intsat,nunit), & + unit_field(nu_lymphflow,nunit),unit_field(nu_blood_press,nunit),unit_field(nu_tt,nunit), & + unit_field(nu_sa,nunit),unit_field(nu_Pe_max,nunit),unit_field(nu_Pe_min,nunit) endif enddo time_to_run = time_0 From d5b43cfc2726f4120cbdac2c6ed3ba7a75b0a317 Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Tue, 20 Sep 2022 14:31:16 +1200 Subject: [PATCH 27/30] nothing --- src/lib/lymphatics.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 3126e091..9a83474f 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -42,7 +42,7 @@ subroutine alveolar_capillary_flux(ne,write_out) logical,intent(in) :: write_out ! Baseline value parameters (eventually will be user-defined?) - integer,parameter :: sex = 0 ! 0 = male, 1 = female + integer,parameter :: sex = 1 ! 0 = male, 1 = female !sex only determines the weight and therefore size of the lung. Should be updated based on CT ! Capillary parameters From 6bee25318a8921f5c5c5ae0dc6fac880edf95bda Mon Sep 17 00:00:00 2001 From: Merryn Tawhai Date: Wed, 17 Apr 2024 08:27:18 +1000 Subject: [PATCH 28/30] remove parentlist from arrays declaration --- src/lib/arrays.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib/arrays.f90 b/src/lib/arrays.f90 index 84094875..fcf5e533 100644 --- a/src/lib/arrays.f90 +++ b/src/lib/arrays.f90 @@ -152,7 +152,7 @@ module arrays elem_units_below, maxgen,capillary_bf_parameters, zero_tol,loose_tol,gasex_field, & num_lines_2d, lines_2d, line_versn_2d, lines_in_elem, nodes_in_line, elems_2d, & elem_cnct_2d, elem_nodes_2d, elem_versn_2d, elem_lines_2d, elems_at_node_2d, arclength, & - scale_factors_2d, parentlist, fluid_properties, lymphatic_properties, elasticity_vessels, admittance_param, & + scale_factors_2d, fluid_properties, lymphatic_properties, elasticity_vessels, admittance_param, & elasticity_param, two_parameter, three_parameter, four_parameter, all_admit_param, update_parameter, & mesh_from_depvar, depvar_at_node, depvar_at_elem, SparseCol, SparseRow, triangle, & update_resistance_entries, vertex_xyz, & From de9948d1a525bb1c042d54ca8ab42de196084854 Mon Sep 17 00:00:00 2001 From: rli787 Date: Thu, 27 Feb 2025 11:23:01 +1300 Subject: [PATCH 29/30] Remove chest wall compliance from tissue compliance. --- src/lib/ventilation.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lib/ventilation.f90 b/src/lib/ventilation.f90 index 9c6961d2..ae040ace 100644 --- a/src/lib/ventilation.f90 +++ b/src/lib/ventilation.f90 @@ -556,8 +556,8 @@ subroutine tissue_compliance(chest_wall_compliance,undef) *(lambda**2+1.0_dp)/lambda**4) unit_field(nu_comp,nunit) = undef/unit_field(nu_comp,nunit) ! V/P ! add the chest wall (proportionately) in parallel - unit_field(nu_comp,nunit) = 1.0_dp/(1.0_dp/unit_field(nu_comp,nunit)& - +1.0_dp/(chest_wall_compliance/dble(num_units))) +! unit_field(nu_comp,nunit) = 1.0_dp/(1.0_dp/unit_field(nu_comp,nunit)& +! +1.0_dp/(chest_wall_compliance/dble(num_units))) !estimate an elastic recoil pressure for the unit unit_field(nu_pe,nunit) = cc/2.0_dp*(3.0_dp*a+b)*(lambda**2.0_dp & -1.0_dp)*exp_term/lambda From 2c506556b2542763d11b2037dfc6481cab39db72 Mon Sep 17 00:00:00 2001 From: Behdad Shaarbaf Ebrahimi Date: Thu, 4 Dec 2025 15:13:16 +1300 Subject: [PATCH 30/30] Fixing the run time statement --- src/lib/lymphatics.f90 | 92 +++++++++++++++++++++--------------------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90 index 9a83474f..960733fb 100644 --- a/src/lib/lymphatics.f90 +++ b/src/lib/lymphatics.f90 @@ -9,42 +9,42 @@ module lymphatics ! !This module contains code for pulmonary fluid flux within the alveolo-capillary network, !and lymph transport through lymphatic collecting vessels. - + use arrays use diagnostics use indices use other_consts use precision ! sets dp for precision - + implicit none - + !Module parameters - + !Module types - + !Module variables - + !Interfaces ! private public alveolar_capillary_flux public lymphatic_transport contains - + !!!############################################################################# - + subroutine alveolar_capillary_flux(ne,write_out) !*alveolar_capillary_flux:* calculate fluid flux from blood to interstitium !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ALVEOLAR_CAPILLARY_FLUX" :: ALVEOLAR_CAPILLARY_FLUX use other_consts,only: pi - + integer,intent(in) :: ne logical,intent(in) :: write_out ! Baseline value parameters (eventually will be user-defined?) integer,parameter :: sex = 1 ! 0 = male, 1 = female !sex only determines the weight and therefore size of the lung. Should be updated based on CT - + ! Capillary parameters real(dp),parameter :: capillary_conductivity = 4.41335e-8 !ml/s/Pa obtained from Parker (6e-8 cm H2O) !real(dp),parameter :: capillary_conductivity = 9.26e-8 !ml/s/Pa obtained from Parker (6e-8 cm H2O) @@ -57,7 +57,7 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Simulation parameters real(dp),parameter :: breathing_rate = 15.0_dp !constant but should be imported directly from ventilation model - + ! Local variables integer :: i,liflowcount,nunit,n_timesteps,printcount real(dp) :: alveolar_volume,breathing_function,capillary_flow,capillary_vps,capillary_osmotic,capillary_osm_n, & @@ -72,19 +72,19 @@ subroutine alveolar_capillary_flux(ne,write_out) real(dp) :: sat1,sat2,sat3,sat4,sat5 logical :: continue - - + + character(len=60) :: sub_name - + ! -------------------------------------------------------------------------- - + sub_name = 'alveolar_capillary_flux' call enter_exit(sub_name,1) - ! get information for the unit fron unit_field + ! get information for the unit fron unit_field ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element - nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) - capillary_pressure = unit_field(nu_blood_press,nunit)/133.32239_dp !converted to mmHg + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + capillary_pressure = unit_field(nu_blood_press,nunit)/133.32239_dp !converted to mmHg transit_time = unit_field(nu_tt,nunit) capillary_SA = unit_field(nu_sa,nunit) max_Pe = unit_field(nu_Pe_max,nunit) @@ -99,7 +99,7 @@ subroutine alveolar_capillary_flux(ne,write_out) ! Calculated values lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp ! g;female lung weight of 639g and male of 840g - should be updated from CT open_capillaries = 1.0_dp/6.0_dp !based on open capillaries at rest. Should be solved for by perfusion model? - capillary_volume_raw = 213.00_dp !Gehr 1978 based on having a body mass of 74 kg - should be unique to each person + capillary_volume_raw = 213.00_dp !Gehr 1978 based on having a body mass of 74 kg - should be unique to each person capillary_volume = (capillary_volume_raw*open_capillaries)/real(num_units) !Only used for the osmotic model at the moment (which isn't operational) can volume be obtained elsewhere? @@ -112,8 +112,8 @@ subroutine alveolar_capillary_flux(ne,write_out) interstitial_volume_b = 0.48_dp*interstitial_capacity !assumed to be around 48% saturated at rest alveolar_volume = 0.0_dp !alveolar volume likely greater at rest, but is lost to respiration - further information needed to put in model liflowcount = 0 - initial_lymphatic_surface_area = capillary_SA !initial lymphatic SA assumed to be the same as capillary SA as no other indication. - !both initial lymphatic SA and initial lymphatic hydraulic conductivity make up the filtration coefficient and currently it is the + initial_lymphatic_surface_area = capillary_SA !initial lymphatic SA assumed to be the same as capillary SA as no other indication. + !both initial lymphatic SA and initial lymphatic hydraulic conductivity make up the filtration coefficient and currently it is the !hydraulic conductivity that is adjusted to compensate ! initial lymphatic values @@ -127,7 +127,7 @@ subroutine alveolar_capillary_flux(ne,write_out) int_osm_n = 0.0_dp initial_osm_n = 0.0_dp total_osm_flux = 0.0_dp - + time = 0.0_dp printcount = 0 total_hydro_flux = 0.0_dp @@ -144,7 +144,7 @@ subroutine alveolar_capillary_flux(ne,write_out) ! dt or n_timesteps should be controlled by the user n_timesteps = 96 dt = transit_time/real(n_timesteps) - + breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate) continue = .true. @@ -154,7 +154,7 @@ subroutine alveolar_capillary_flux(ne,write_out) sat3 = 3.0_dp sat4 = 4.0_dp sat5 = 5.0_dp - + !do while(time < lymphatic_properties%test_time) do while (continue) time_sum = dt @@ -165,7 +165,7 @@ subroutine alveolar_capillary_flux(ne,write_out) interstitial_volume = interstitial_volume_a + interstitial_volume_b interstitial_saturation = interstitial_volume / interstitial_capacity ! saturation as a proportion of 0-100% time_variable = time + time_sum - + ! calculating flux from capillary into interstitium interstitial_pressure_a = fluctuation * sin(time_variable * breathing_function) + & (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + & @@ -210,7 +210,7 @@ subroutine alveolar_capillary_flux(ne,write_out) else interstitial_volume_a = interstitial_volume_a + flux_a endif - + interstitial_volume_b = interstitial_volume_b + flux_b diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ & (200_dp)) * dt @@ -221,7 +221,7 @@ subroutine alveolar_capillary_flux(ne,write_out) int_osm_conc = int_osm_n/interstitial_volume interstitial_osmotic = real(i)*int_osm_conc*IGC*T osm_flux = (lymphatic_properties%reflection_coefficient * capillary_SA * & - (capillary_osmotic - interstitial_osmotic))* dt + (capillary_osmotic - interstitial_osmotic))* dt cap_osm_conc = capillary_osm_n / capillary_volume int_osm_conc = int_osm_n / interstitial_volume @@ -232,7 +232,7 @@ subroutine alveolar_capillary_flux(ne,write_out) capillary_volume = capillary_volume + osm_flux interstitial_volume_b = interstitial_volume_b - osm_flux endif - + net_flux = osm_flux + flux_a + flux_b if(net_flux > 0.0_dp)then osm_n_flux = net_flux * cap_osm_conc @@ -242,9 +242,9 @@ subroutine alveolar_capillary_flux(ne,write_out) osm_n_flux = net_flux * int_osm_conc int_osm_n = int_osm_n - osm_n_flux endif - + total_osm_flux = total_osm_flux + osm_flux - + !calculating flux from interstitium to initial lymphatics if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then lymphatic_conductivity = 1.48_dp * 4.41335e-8 !all calculated does as a function of capillary_conductivity @@ -261,7 +261,7 @@ subroutine alveolar_capillary_flux(ne,write_out) initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + & ((((lymphPmax-lymphPmin-(fluctuation*2.0_dp))* ((interstitial_volume_b / interstitial_capacity_b)**2.0_dp)) + & (lymphPmin + fluctuation))) - !arbitrarily defined mathematical relationship to show that lymphatic pressure does not change much at low volumes with a + !arbitrarily defined mathematical relationship to show that lymphatic pressure does not change much at low volumes with a !large volume change, but at high volumes only a small volume change is needed to cause a large change in pressure !write(*,'(''Plym: '',f8.4)')initial_lymphatic_pressure if(interstitial_volume.le.0.0_dp)then @@ -286,7 +286,7 @@ subroutine alveolar_capillary_flux(ne,write_out) else initial_lymph_conc = 0.0_dp endif - + time_sum = time_sum + dt total_flux = total_hydro_flux ! +total_osm_flux @@ -311,7 +311,7 @@ subroutine alveolar_capillary_flux(ne,write_out) endif printcount = 0 if(time.gt.200.0_dp*transit_time)then - if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then + if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then continue = .false. endif endif @@ -326,13 +326,13 @@ subroutine alveolar_capillary_flux(ne,write_out) unit_field(nu_av_flux,nunit) = total_flux/time unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time !write(*,'('' T='',e12.3,'': intsat='',e12.6,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.6,'' ul/s'')') & - ! unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& + ! unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),& ! unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit) call enter_exit(sub_name,2) - + end subroutine alveolar_capillary_flux - + !!!############################################################################# subroutine lymphatic_transport(filename) @@ -345,9 +345,9 @@ subroutine lymphatic_transport(filename) real(dp) :: time_to_run,time_0 character(len=300) :: writefile character(len=60) :: sub_name - !real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_av_flux,nu_lymphflow,nu_time + !real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_av_flux,nu_lymphflow,nu_time ! -------------------------------------------------------------------------- - + sub_name = 'lymphatic_transport' call enter_exit(sub_name,1) @@ -356,13 +356,13 @@ subroutine lymphatic_transport(filename) else ! need to append the correct filename extension writefile = trim(filename)//'.oplymph' endif - + open(10, file=writefile, status='replace') call cpu_time(time_0) do ne = 1,num_elems if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then - nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) call alveolar_capillary_flux(ne,.false.) np = elem_nodes(2,ne) write(10,'(i8,11(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_av_flux,nunit),unit_field(nu_intsat,nunit), & @@ -371,14 +371,14 @@ subroutine lymphatic_transport(filename) endif enddo time_to_run = time_0 - call cpu_time(time_0) + call cpu_time(time_to_run) time_to_run = time_to_run - time_0 write(*,*) 'run time=',time_to_run close(10) do ne = num_elems,1,-1 - if(elem_field(ne_group,ne).eq.1.0_dp)then - nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) + if(elem_field(ne_group,ne).eq.1.0_dp)then + nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) elem_field(ne_radius_in0,ne) = unit_field(nu_flux,nunit) elem_field(ne_radius_out0,ne) = unit_field(nu_intsat,nunit) else if(elem_field(ne_group,ne).eq.0.0_dp)then ! artery @@ -393,9 +393,9 @@ subroutine lymphatic_transport(filename) elem_field(ne_radius_out0,ne) = elem_field(ne_radius_out0,ne)/real(elem_cnct(1,0,ne)) endif enddo - + call enter_exit(sub_name,2) - + end subroutine lymphatic_transport !!!############################################################################# @@ -405,7 +405,7 @@ end module lymphatics !FUTURE DIRECTIONS !input a constant to account for difference between current values and expected values - !Model appeared to be working within the range of the literature but is likely off by a factor of 1000 due to nl to ul conversion error. + !Model appeared to be working within the range of the literature but is likely off by a factor of 1000 due to nl to ul conversion error. !Need to check that outputted units are correct - most things are in ml and mmHg !Lymphatic network tree !currently all lymph is returned to the circulation immediately, in reality it moves up a tree of lymphatics against a pressure gradient