<commit_before>module pressure_resistance_flow
!*Brief Description:* This module contains tools that are used to solve systems of equations representing steady pressure, resistance and flow problems in any branching geometry. The subroutines in this module are core subroutines that are used in many problem types and are applicable beyond lung modelling
!
!*LICENSE:*
!TBC
!
!
!*Full Description:*
!
!This module contains tools that are used to solve systems of equations representing steady pressure, resistance and flow problems in any branching geometry. The subroutines in this module are core subroutines that are used in many problem types and are applicable beyond lung modelling
  use solve, only: BICGSTAB_LinSolv,pmgmres_ilu_cr
  use other_consts, only: TOLERANCE
  implicit none
  !Module parameters

  !Module types

  !Module depvar

  !Interfaces
  private
  public evaluate_prq,calculate_ppl
contains
!###################################################################################
!
!*evaluate_PRQ:* Solves for pressure and flow in a rigid or compliant tree structure
  subroutine evaluate_prq(mesh_type,grav_dirn,grav_factor,bc_type,inlet_bc,outlet_bc)
  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EVALUATE_PRQ" :: EVALUATE_PRQ
    use indices
    use capillaryflow,only: cap_flow_ladder
    use arrays,only: dp,num_elems,num_nodes,elem_field,elem_nodes,elem_cnct,node_xyz
    use diagnostics, only: enter_exit
    !local variables
    integer :: mesh_dof,depvar_types
    integer, allocatable :: mesh_from_depvar(:,:,:)
    integer, allocatable :: depvar_at_node(:,:,:)
    integer, allocatable :: depvar_at_elem(:,:,:)
    integer, dimension(0:2,2) :: depvar_totals
    integer, allocatable :: SparseCol(:)
    integer, allocatable :: SparseRow(:)
    integer, allocatable :: update_resistance_entries(:)
    real(dp), allocatable :: SparseVal(:)
    real(dp), allocatable :: RHS(:)
    integer :: num_vars,NonZeros,MatrixSize
    integer :: AllocateStatus

    real(dp), allocatable :: prq_solution(:,:),solver_solution(:)
    real(dp) :: viscosity,density,inlet_bc,outlet_bc,inletbc,outletbc,grav_vect(3),gamma,total_resistance,ERR
    logical, allocatable :: FIX(:)
    logical :: ADD=.FALSE.,CONVERGED=.FALSE.
    character(len=60) :: sub_name,mesh_type,vessel_type,mechanics_type,bc_type
    integer :: grav_dirn,no,depvar,KOUNT,nz,ne,SOLVER_FLAG,ne0,ne1,nj
    real(dp) :: MIN_ERR,N_MIN_ERR,elasticity_parameters(3),mechanics_parameters(2),grav_factor,P1
    real(dp) :: P2,Q01,Rin,Rout,x_cap,y_cap,z_cap,Ppl,LPM_R,Lin,Lout

    sub_name = 'evaluate_prq'
    call enter_exit(sub_name,1)
!!---------DESCRIPTION OF MODEL Types -----------
!mesh_type: can be simple_tree, full_plus_ladder, full_sheet, full_tube The first can be airways, arteries, veins but no special features at the terminal level, the last one has arteries and veins connected by capillary units of some type (lung ladder acinus, lung sheet capillary bed, capillaries are just tubes represented by an element)

!vessel_type:
  !rigid, no elasticity, no parameters required
  !elastic_g0_beta, R=R0*((Ptm/G0)+1.d0)^(1.d0/elasticity_parameters(2)),with an optional maximum pressure beyond which the vessel radius is constant three parameters, g0, elasticity_parameters(2), elasticity_parameters(3)
  !elastic alpha,  R=R0*(alpha*Ptm+1.d0), up to a limit elasticity_parameters(3) two parameters alpha, elasticity_parameters(3)
  !elastic_hooke, two parameters E and h,R=R0+3.0_dp*R0**2*Ptm/(4.0_dp*E*h*R0)

!mechanics type:
  !linear two parmeters, transpulmonary pressure (average) and pleural density (gradient)
  !mechanics, two parameters, pressure and stretch fields

!bc_type:
    !pressure (at inlet and outlets)
    !flow (flow at inlet pressure at outlet).

vessel_type='elastic_alpha'

mechanics_type='linear'

if (vessel_type.eq.'rigid') then
    elasticity_parameters=0.0_dp
elseif (vessel_type.eq.'elastic_g0_beta') then
    elasticity_parameters(1)=6.67e3_dp!G0 (Pa)
    elasticity_parameters(2)=1.0_dp!elasticity_parameters(2)
    elasticity_parameters(3)=32.0_dp*98.07_dp !elasticity_parameters(3) (Pa)
elseif (vessel_type.eq.'elastic_alpha') then
    elasticity_parameters(1)=1.503e-4_dp!alpha (1/Pa)
    elasticity_parameters(2)=32.0_dp*98.07_dp !elasticity_parameters(3) (Pa)
    elasticity_parameters(3)=0.0_dp !Not used
elseif (vessel_type.eq.'elastic_hooke') then
     elasticity_parameters(1)=1.5e6_dp !Pa
     elasticity_parameters(2)=0.1_dp!this is a fraction of the radius so is unitless
     elasticity_parameters(3)=0.0_dp !Not used
else
    print *, 'WARNING: Your chosen vessel type does not seem to be implemented assuming rigid'
    vessel_type='rigid'
    elasticity_parameters=0.0_dp
endif
<commit_msg>Update pressure_resistance_flow.f90<commit_after>module pressure_resistance_flow
!*Brief Description:* This module contains tools that are used to solve systems of equations representing steady pressure, resistance and flow problems in any branching geometry. The subroutines in this module are core subroutines that are used in many problem types and are applicable beyond lung modelling
!
!*LICENSE:*
!TBC
!
!
!*Full Description:*
!
!This module contains tools that are used to solve systems of equations representing steady pressure, resistance and flow problems in any branching geometry. The subroutines in this module are core subroutines that are used in many problem types and are applicable beyond lung modelling
  use solve, only: BICGSTAB_LinSolv,pmgmres_ilu_cr
  use other_consts, only: TOLERANCE
  implicit none
  !Module parameters

  !Module types

  !Module depvar

  !Interfaces
  private
  public evaluate_prq,calculate_ppl
contains
!###################################################################################
!
!*evaluate_PRQ:* Solves for pressure and flow in a rigid or compliant tree structure
  subroutine evaluate_prq(mesh_type,grav_dirn,grav_factor,bc_type,inlet_bc,outlet_bc)
  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EVALUATE_PRQ" :: EVALUATE_PRQ
    use indices
    use capillaryflow,only: cap_flow_ladder
    use arrays,only: dp,num_elems,num_nodes,elem_field,elem_nodes,elem_cnct,node_xyz
    use diagnostics, only: enter_exit
    !local variables
    integer :: mesh_dof,depvar_types
    integer, allocatable :: mesh_from_depvar(:,:,:)
    integer, allocatable :: depvar_at_node(:,:,:)
    integer, allocatable :: depvar_at_elem(:,:,:)
    integer, dimension(0:2,2) :: depvar_totals
    integer, allocatable :: SparseCol(:)
    integer, allocatable :: SparseRow(:)
    integer, allocatable :: update_resistance_entries(:)
    real(dp), allocatable :: SparseVal(:)
    real(dp), allocatable :: RHS(:)
    integer :: num_vars,NonZeros,MatrixSize
    integer :: AllocateStatus

    real(dp), allocatable :: prq_solution(:,:),solver_solution(:)
    real(dp) :: viscosity,density,inlet_bc,outlet_bc,inletbc,outletbc,grav_vect(3),gamma,total_resistance,ERR
    logical, allocatable :: FIX(:)
    logical :: ADD=.FALSE.,CONVERGED=.FALSE.
    character(len=60) :: sub_name,mesh_type,vessel_type,mechanics_type,bc_type
    integer :: grav_dirn,no,depvar,KOUNT,nz,ne,SOLVER_FLAG,ne0,ne1,nj
    real(dp) :: MIN_ERR,N_MIN_ERR,elasticity_parameters(3),mechanics_parameters(2),grav_factor,P1
    real(dp) :: P2,Q01,Rin,Rout,x_cap,y_cap,z_cap,Ppl,LPM_R,Lin,Lout

    sub_name = 'evaluate_prq'
    call enter_exit(sub_name,1)
!!---------DESCRIPTION OF MODEL Types -----------
!mesh_type: can be simple_tree, full_plus_ladder, full_sheet, full_tube The first can be airways, arteries, veins but no special features at the terminal level, the last one has arteries and veins connected by capillary units of some type (lung ladder acinus, lung sheet capillary bed, capillaries are just tubes represented by an element)

!vessel_type:
  !rigid, no elasticity, no parameters required
  !elastic_g0_beta, R=R0*((Ptm/G0)+1.d0)^(1.d0/elasticity_parameters(2)),with an optional maximum pressure beyond which the vessel radius is constant three parameters, g0, elasticity_parameters(2), elasticity_parameters(3)
  !elastic alpha,  R=R0*(alpha*Ptm+1.d0), up to a limit elasticity_parameters(3) two parameters alpha, elasticity_parameters(3)
  !elastic_hooke, two parameters E and h,R=R0+3.0_dp*R0**2*Ptm/(4.0_dp*E*h*R0)

!mechanics type:
  !linear two parmeters, transpulmonary pressure (average) and pleural density (gradient)
  !mechanics, two parameters, pressure and stretch fields

!bc_type:
    !pressure (at inlet and outlets)
    !flow (flow at inlet pressure at outlet).

vessel_type='elastic_g0_beta'

mechanics_type='linear'

if (vessel_type.eq.'rigid') then
    elasticity_parameters=0.0_dp
elseif (vessel_type.eq.'elastic_g0_beta') then
    elasticity_parameters(1)=6.67e3_dp!G0 (Pa)
    elasticity_parameters(2)=1.0_dp!elasticity_parameters(2)
    elasticity_parameters(3)=32.0_dp*98.07_dp !elasticity_parameters(3) (Pa)
elseif (vessel_type.eq.'elastic_alpha') then
    elasticity_parameters(1)=1.503e-4_dp!alpha (1/Pa)
    elasticity_parameters(2)=32.0_dp*98.07_dp !elasticity_parameters(3) (Pa)
    elasticity_parameters(3)=0.0_dp !Not used
elseif (vessel_type.eq.'elastic_hooke') then
     elasticity_parameters(1)=1.5e6_dp !Pa
     elasticity_parameters(2)=0.1_dp!this is a fraction of the radius so is unitless
     elasticity_parameters(3)=0.0_dp !Not used
else
    print *, 'WARNING: Your chosen vessel type does not seem to be implemented assuming rigid'
    vessel_type='rigid'
    elasticity_parameters=0.0_dp
endif

{'fraction_non_alphanumeric': 0.11687077702702703, 'fraction_numerical': 0.01974239864864865, 'mean_word_length': 5.732764747690121, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2960701', 'n_tokens_mistral': 3284, 'n_tokens_neox': 3063, 'n_words': 795}