[FLASH-BUGS] Possible bug in physical_constants.F90
Vincenzo Antonuccio-Delogu
van at ct.astro.it
Thu Jul 22 08:38:02 CDT 2004
Ffrom: V. Antonuccio-Delogu, INAF-Catania Astrophysical Observatory, ITALY
Code: FLASH
Version: 2.3
Description: I believe that the scalings of the "ideal gas constant" and
Boltzmann" constant in the calls of add_constant_to_db at lines 196 and 199
are wrong. Both constants have dimensions: [energy] [temp]^-1 [mole]^-1 =
[length]^2 [time]^-2 [temp]^-1. Then, the dependence on mass should be
dropped, which means that lines 195-197 should read:
call add_constant_to_db ("Boltzmann", 1.38065812E-16, &
& 2., -2., 0., &
& -1., 0.)
rather than:
call add_constant_to_db ("Boltzmann", 1.38065812E-16, &
& 2., -2., 1., &
& -1., 0.)
and analogously for the ideal gas constant.
Note that in eos.F90, eos1d.F90 they are used to deduce the pressure:
pres = gascon * dens * temp * abarinv
If density has dimensions: [mass] [length]^-3, then, in order that pressure
has the right dimensions ([mass] [length]^-1 [time]^-2) it is necessary that
the dimension of gascon be: [length]^2 [time]^-2 [temp]^-1.
I enclose a version of physical_constants.F90 which includes this correction.
In this new version I also define new units, which I denote as "FLY" (see
lines 169-171), defined in such a way that: length_unit = 1 Mpc, time_unit =
2/(3H_0), mass_unit=5.23e12 M_sol, H_0 = 100 km/sec/Mpc. In these units, the
gas constants has a value: gascon = 3.69 * 10^-7. Using the uncorrected
scaling in the original physical_constants.F90 one gets a value of ~ 10^-59
for gascon in FLY units. After the proposed correction, one recovers the
correct values.
*************************************************************
Vincenzo ANTONUCCIO-DELOGU *
*
VISITING COPENHAGEN ASTROPHYSICAL OBSERVATORY *
July 1st - September 8th, 2004 *
Tlf.: +45 - 35 32 59 05 *
*
INAF - Catania Astrophysical Observatory *
ITALY *
Citta' Universitaria *
Via Santa Sofia 78 *
I-95125 Catania - ITALY *
*
Tlf.: +39-095 7332 318 *
*
Fax: + 39-095-33 05 92 *
*
e-mails: van at ct.astro.it, antonucc at tac.dk *
WWW: http://w3.ct.astro.it/CHEAL/cheal.html *
*
*************************************************************
-------------- next part --------------
!*************************************************************************
!
! Module: physical_constants
!
! Description: These routines maintain a database of physical constants
! which can be queried by string name and optionally converted
! to any chosen system of units. The default system of units
! is CGS.
!
! For example, a program using this module might obtain the
! value of Newton's gravitational constant in units of
! Mpc^3 Gyr^-2 Msun^-1 by calling
!
! call get_constant_from_db ("Newton", G, len_unit="Mpc",
! time_unit="Gyr", mass_unit="Msun")
!
! In this example, the local variable G is set to the result,
! 4.4983E-15 (to five significant figures).
!
! Physical constants are taken from the 1998 Review of Particle
! Properties, Eur. Phys. J. C 3, 1 (1998), which in turn takes
! most of its values from Cohen, E. R. and Taylor, B. N.,
! Rev. Mod. Phys. 59, 1121 (1987).
! Provides: get_constant_from_db (n, v[, units])
! Return the value of the physical constant named n
! in the variable v. Optional unit specifications are
! used to convert the result. If the constant name or
! one or more unit names aren't recognized, a value of
! 0 is returned.
! add_constant_to_db (n, v, len, time, mass, temp, chg)
! Add a physical constant to the database. n is the
! name to assign, and v is the value in CGS units.
! len, time, mass, temp, and chg are the exponents of
! the various base units used in defining the unit
! scaling of the constant. For example, a constant
! with units (in CGS) of cm^3 s^-2 g^-1 would have
! len=3, time=-2, mass=-1, temp=0, and chg=0.
! add_unit_to_db (t, n, v)
! Add a unit of measurement to the database. t is the
! type of unit ("length", "time", "mass", "charge",
! "temperature"), n is the name of the unit, and v is
! its value in terms of the corresponding CGS unit.
! init_constants_db ()
! Initialize the constants and units databases.
! Can be called by the user program, but doesn't have
! to be, as it is automatically called when needed
! (ie. if a "get" or "add" is called before init).
! list_constants_db (lun)
! List the constants and units
! databases to the specified logical I/O unit.
! destroy_constants_db ()
! Deallocate the memory used by the constants and units
! databases, requiring another init call before they
! can be accessed again.
module physical_constants
!*************************************************************************
! Public data types, constants, and routines.
implicit none
private
public :: get_constant_from_db, add_constant_to_db, &
& add_unit_to_db, list_constants_db, init_constants_db, &
& destroy_constants_db
!-------------------------------------------------------------------------
! Constants used to set base types.
integer, parameter :: n_base_units = 5
integer, parameter :: len_base_unit = 1, time_base_unit = 2, &
& mass_base_unit = 3, temp_base_unit = 4, &
& chg_base_unit = 5
character(len=3), dimension(n_base_units) :: base_unit_name = &
& (/ "cm ", "s ", "g ", "K ", "esu" /)
! Derived type declarations for physical constants.
type const_node_type
character(len=80) :: name
real :: cgs_value
real :: unit_exp(n_base_units)
type (const_node_type), pointer :: next
end type
type unit_node_type
character(len=80) :: name
integer :: base_unit
real :: cgs_value
type (unit_node_type), pointer :: next
end type
! The constant and unit databases.
type (const_node_type), pointer, save :: const_db
type (unit_node_type), pointer, save :: unit_db
! Has the initialization routine been called yet?
logical, save :: init_was_called = .false.
!-------------------------------------------------------------------------
contains
!*************************************************************************
! Routine: init_constants_db ()
! Description: Initialize the physical constant and unit databases.
! Called automatically by the "get" routines if it hasn't
! been called previously.
subroutine init_constants_db()
!-------------------------------------------------------------------------
! If we're already initialized, don't do it again.
if (init_was_called) then
write (*,*) 'init_constants_db: already initialized'
return
endif
! Make sure the database pointers are sane.
nullify (const_db)
nullify (unit_db)
! init_constants_db has been called.
init_was_called = .true.
! Physical constants are taken from the 1998 Review of Particle
! Properties, Eur. Phys. J. C 3, 1 (1998), which in turn takes
! most of its values from Cohen, E. R. and Taylor, B. N.,
! Rev. Mod. Phys. 59, 1121 (1987).
! Set up the units of measurement.
call add_unit_to_db ("length", "cm", 1.)
call add_unit_to_db ("time", "s", 1.)
call add_unit_to_db ("temperature", "K", 1.)
call add_unit_to_db ("mass", "g", 1.)
call add_unit_to_db ("charge", "esu", 1.)
call add_unit_to_db ("length", "m", 1.E2)
call add_unit_to_db ("length", "km", 1.E5)
call add_unit_to_db ("length", "pc", 3.0856775807E18)
call add_unit_to_db ("length", "kpc", 3.0856775807E21)
call add_unit_to_db ("length", "Mpc", 3.0856775807E24)
call add_unit_to_db ("length", "Gpc", 3.0856775807E27)
call add_unit_to_db ("length", "Rsun", 6.96E10)
call add_unit_to_db ("length", "AU", 1.49597870662E13)
call add_unit_to_db ("time", "yr", 3.15569252E7)
call add_unit_to_db ("time", "Myr", 3.15569252E13)
call add_unit_to_db ("time", "Gyr", 3.15569252E16)
call add_unit_to_db ("mass", "kg", 1.E3)
call add_unit_to_db ("mass", "Msun", 1.9889225E33)
call add_unit_to_db ("mass", "amu", 1.660540210E-24)
call add_unit_to_db ("charge", "C", 2.99792458E9)
! here are FLY units
call add_unit_to_db ("length", "LFLY", 3.0856775807E24)
call add_unit_to_db ("time", "TFLY", 2.05759E17)
call add_unit_to_db ("mass", "MFLY", 9.8847E45)
! Set up the physical constants.
! Exponent ordering: len, time, mass, temp, chg
call add_constant_to_db ("Newton", 6.6725985E-8, &
& 3., -2., -1., &
& 0., 0.)
call add_constant_to_db ("speed of light", 2.99792458E10, &
& 1., -1., 0., &
& 0., 0.)
call add_constant_to_db ("Planck", 6.626075540E-27, &
& 2., -1., 1., &
& 0., 0.)
call add_constant_to_db ("electron charge", 4.803206815E-10, &
& 0., 0., 0., &
& 0., 1.)
call add_constant_to_db ("electron mass", 9.109389754E-28, &
& 0., 0., 1., &
& 0., 0.)
call add_constant_to_db ("proton mass", 1.672623110E-24, &
& 0., 0., 1., &
& 0., 0.)
call add_constant_to_db ("fine-structure", 7.2973530764E-3, &
& 0., 0., 0., &
& 0., 0.)
call add_constant_to_db ("Avogadro", 6.022136736E23, &
& 0., 0., 0., &
& 0., 0.)
call add_constant_to_db ("Boltzmann", 1.38065812E-16, &
& 2., -2., 0., &
& -1., 0.)
call add_constant_to_db ("ideal gas constant", 8.3145119843E7, &
& 2., -2., 0., &
& -1., 0.)
call add_constant_to_db ("Wien", 2.89775624E-1, &
& 1., 0., 0., &
& 1., 0.)
call add_constant_to_db ("Stefan-Boltzmann", 5.6705119E-5, &
& 0., -3., 1., &
& -4., 0.)
call add_constant_to_db ("pi", 3.141592653589793238E0, &
& 0., 0., 0., 0., 0.)
call add_constant_to_db ("e", 2.718281828459045235E0, &
& 0., 0., 0., 0., 0.)
call add_constant_to_db ("Euler", 5.77215664901532861E-1, &
& 0., 0., 0., 0., 0.)
!------------------------------------------------------------------------
return
end subroutine init_constants_db
!************************************************************************
! Routine: destroy_constants_db ()
! Description: Deallocate the memory used by the databases.
subroutine destroy_constants_db()
!------------------------------------------------------------------------
type (const_node_type), pointer :: cptr_this, cptr_next
type (unit_node_type), pointer :: uptr_this, uptr_next
!------------------------------------------------------------------------
! Deallocate the databases.
cptr_this => const_db
do while (associated(cptr_this))
cptr_next => cptr_this%next
deallocate (cptr_this)
cptr_this => cptr_next
enddo
uptr_this => unit_db
do while (associated(uptr_this))
uptr_next => uptr_this%next
deallocate (uptr_this)
uptr_this => uptr_next
enddo
! Will have to call init_constants_db again before using databases.
init_was_called = .false.
!------------------------------------------------------------------------
return
end subroutine destroy_constants_db
!************************************************************************
! Routine: get_constant_from_db()
! Description: Return the value of the physical constant named n
! in the variable v. Optional unit specifications are
! used to convert the result. If the constant name or
! one or more unit names aren't recognized, a value of
! 0 is returned.
subroutine get_constant_from_db (name, var, len_unit, time_unit, &
& mass_unit, temp_unit, chg_unit)
!------------------------------------------------------------------------
character(len=*) :: name
real :: var
character(len=*), optional :: len_unit, time_unit, &
& mass_unit, temp_unit, &
& chg_unit
type (const_node_type), pointer :: cnode
type (unit_node_type), pointer :: unode
real :: len_scale, time_scale, &
& mass_scale, temp_scale, &
& chg_scale
!------------------------------------------------------------------------
! First check to see if the database has been created.
! If not, call init_constants_db to do it.
if (.not. init_was_called) call init_constants_db
! Next, find the constant in the database. If it's not found,
! exit with a value of zero.
call find_const_by_name (name, cnode)
var = 0.e0
if (.not. associated(cnode)) return
! Establish the unit scaling. If a unit is not found (or is
! illegally specified, e.g., time_unit="Mpc"), exit with a
! value of zero.
if (present(len_unit)) then
call find_unit_by_name (len_unit, len_base_unit, unode)
if (.not. associated(unode)) return
len_scale = unode%cgs_value
else
len_scale = 1.e0
endif
if (present(time_unit)) then
call find_unit_by_name (time_unit, time_base_unit, unode)
if (.not. associated(unode)) return
time_scale = unode%cgs_value
else
time_scale = 1.e0
endif
if (present(mass_unit)) then
call find_unit_by_name (mass_unit, mass_base_unit, unode)
if (.not. associated(unode)) return
mass_scale = unode%cgs_value
else
mass_scale = 1.e0
endif
if (present(temp_unit)) then
call find_unit_by_name (temp_unit, temp_base_unit, unode)
if (.not. associated(unode)) return
temp_scale = unode%cgs_value
else
temp_scale = 1.e0
endif
if (present(chg_unit)) then
call find_unit_by_name (chg_unit, chg_base_unit, unode)
if (.not. associated(unode)) return
chg_scale = unode%cgs_value
else
chg_scale = 1.e0
endif
! Convert the constant to the requested units and return it.
var = cnode%cgs_value / &
& ( (len_scale ** cnode%unit_exp(len_base_unit)) * &
& (time_scale ** cnode%unit_exp(time_base_unit)) * &
& (mass_scale ** cnode%unit_exp(mass_base_unit)) * &
& (temp_scale ** cnode%unit_exp(temp_base_unit)) * &
& (chg_scale ** cnode%unit_exp(chg_base_unit)) )
!------------------------------------------------------------------------
return
end subroutine get_constant_from_db
!************************************************************************
! Routine: add_constant_to_db ()
! Description: Add a physical constant to the constants database. You
! must specify the name of the constant (a string), its
! value in CGS units, and the exponents for its scaling with
! length, time, mass, temperature, and charge.
subroutine add_constant_to_db (name, cgs_value, len_exp, &
& time_exp, mass_exp, &
& temp_exp, chg_exp)
!------------------------------------------------------------------------
character(len=*) :: name
real :: cgs_value, len_exp, &
& time_exp, mass_exp, &
& temp_exp, chg_exp
type (const_node_type), pointer :: node, this
integer :: istat
!------------------------------------------------------------------------
! Check to make sure the databases have been created. If we
! have been called by init_constants_db, then init_was_called should
! already have been set to .true.
if (.not. init_was_called) call init_constants_db
! Check to make sure the constant isn't already in the
! database.
call find_const_by_name (name, node)
if (associated(node)) then
write (*,*) 'add_constant_to_db: already exists: ', name
else
! If it isn't, create a node and add it to the constants
! database.
allocate (node, stat=istat)
if (istat /= 0) then
write (*,*) 'add_constant_to_db: could not allocate'
else
node%name = name
node%cgs_value = cgs_value
node%unit_exp(len_base_unit) = len_exp
node%unit_exp(time_base_unit) = time_exp
node%unit_exp(mass_base_unit) = mass_exp
node%unit_exp(temp_base_unit) = temp_exp
node%unit_exp(chg_base_unit) = chg_exp
! call make_lowercase (node%name)
nullify (node%next)
endif
if (.not. associated(const_db)) then
const_db => node
else
this => const_db
do while (associated(this%next))
this => this%next
enddo
this%next => node
endif
endif
!------------------------------------------------------------------------
return
end subroutine add_constant_to_db
!************************************************************************
! Routine: add_unit_to_db ()
! Description: Add a unit of measurement to the units database. You
! must specify the name of the unit (a string), its base
! unit type ("length", "time", "mass", "temperature", or
! "charge"), and its value in CGS units.
subroutine add_unit_to_db (base_unit_kind, name, cgs_value)
!------------------------------------------------------------------------
character(len=*) :: name, base_unit_kind
real :: cgs_value
type (unit_node_type), pointer :: node, this
integer :: istat
integer :: base_unit
!------------------------------------------------------------------------
! Check to make sure the databases have been created. If we
! have been called by init_constants_db, then init_was_called
! should already have been set to .true.
if (.not. init_was_called) call init_constants_db
! Is the base unit kind reasonable?
base_unit = 0
if (base_unit_kind == "length") base_unit = len_base_unit
if (base_unit_kind == "time") base_unit = time_base_unit
if (base_unit_kind == "mass") base_unit = mass_base_unit
if (base_unit_kind == "temperature") base_unit = temp_base_unit
if (base_unit_kind == "charge") base_unit = chg_base_unit
if (base_unit == 0) then
write (*,*) 'add_unit_to_db: base unit ', base_unit_kind, &
& ' unknown'
return
endif
! Check to make sure the unit isn't already in the database.
call find_unit_by_name (name, base_unit, node)
if (associated(node)) then
write (*,*) 'add_unit_to_db: already exists: ', name
else
! If it isn't, create a node and add it to the unit database.
allocate (node, stat=istat)
if (istat /= 0) then
write (*,*) 'add_unit_to_db: could not allocate'
else
node%name = name
node%cgs_value = cgs_value
node%base_unit = base_unit
! call make_lowercase (node%name)
nullify (node%next)
endif
if (.not. associated(unit_db)) then
unit_db => node
else
this => unit_db
do while (associated(this%next))
this => this%next
enddo
this%next => node
endif
endif
!------------------------------------------------------------------------
return
end subroutine add_unit_to_db
!************************************************************************
! Routine: find_const_by_name ()
! Description: Find a constant in a linked list using its name.
! Return a pointer to the constant's node if it's found,
! otherwise return a null pointer.
subroutine find_const_by_name (name, node)
!------------------------------------------------------------------------
character(len=*) :: name
type (const_node_type), pointer :: node
! character(len=len(name)) :: name_lcase
!------------------------------------------------------------------------
! name_lcase = name
! call make_lowercase (name_lcase)
node => const_db
do while (associated(node))
if (node%name == name) exit
node => node%next
enddo
!------------------------------------------------------------------------
return
end subroutine find_const_by_name
!************************************************************************
! Routine: find_unit_by_name ()
! Description: Find a unit in a linked list using its name.
! Return a pointer to the unit's node if it's found,
! otherwise return a null pointer. In looking for the
! unit, we also check to make sure that the base unit
! type matches.
subroutine find_unit_by_name (name, base_unit, node)
!------------------------------------------------------------------------
character(len=*) :: name
type (unit_node_type), pointer :: node
integer :: base_unit
! character(len=len(name)) :: name_lcase
!------------------------------------------------------------------------
! name_lcase = name
! call make_lowercase (name_lcase)
node => unit_db
do while (associated(node))
if ( (node%name == name) .and. &
& (node%base_unit == base_unit) ) exit
node => node%next
enddo
!------------------------------------------------------------------------
return
end subroutine find_unit_by_name
!************************************************************************
! Routine: list_constants_db ()
! Description: Echo the physical constants and units of measurement to a
! specified I/O unit.
subroutine list_constants_db (lun)
!------------------------------------------------------------------------
integer :: lun
type (const_node_type), pointer :: cnode
type (unit_node_type), pointer :: unode
character(len=80) :: fmt_str, unit_str, exp_str
integer :: i, j, k
!------------------------------------------------------------------------
! If the databases haven't been initialized, exit.
if (.not. init_was_called) then
write (lun,*) &
& 'list_constants_db: init_constants_db has not been called'
return
endif
! List the units of measurement.
write (lun,*)
write (lun,*) 'Known units of measurement:'
write (lun,*)
fmt_str = "(2X, A20, ' = ', ES20.13, 1X, A3)"
unode => unit_db
do while (associated(unode))
write (lun, fmt_str) unode%name, unode%cgs_value, &
& base_unit_name(unode%base_unit)
unode => unode%next
enddo
! List the physical constants.
write (lun,*)
write (lun,*) 'Known physical constants:'
write (lun,*)
fmt_str = "(2X, A20, ' = ', ES20.13, 1X, A30)"
cnode => const_db
do while (associated(cnode))
write (unit_str(1:), '(80X)')
j = 1
do i = 1, n_base_units
if (abs(cnode%unit_exp(i)) > 1.e-10) then
write (unit_str(j:), "(A)") &
& base_unit_name(i)(1:len_trim(base_unit_name(i)))
j = j + len_trim(base_unit_name(i))
if (cnode%unit_exp(i) /= 1) then
exp_str = " "
if (cnode%unit_exp(i) == int(cnode%unit_exp(i))) then
write (exp_str, "(I6)") int(cnode%unit_exp(i))
else
write (exp_str, "(F6.3)") cnode%unit_exp(i)
endif
k = index (exp_str(1:len_trim(exp_str)), " ", .true.)
write (unit_str(j:), "('^', A)") &
& exp_str(k+1:len_trim(exp_str))
j = j + 1 + 1 + (len_trim(exp_str)-(k+1) + 1)
else
j = j + 1
endif
endif
enddo
write (lun, fmt_str) cnode%name, cnode%cgs_value, unit_str
cnode => cnode%next
enddo
write (lun,*)
!------------------------------------------------------------------------
return
end subroutine list_constants_db
!************************************************************************
! Routine: make_lowercase()
! Description: Convert a string to lowercase.
subroutine make_lowercase (str)
character(len=*):: str
integer :: i
do i = 1, len_trim(str)
if (lge(str(i:i), 'A') .and. lle(str(i:i), 'Z')) &
& str(i:i) = achar( iachar(str(i:i)) + 32 )
enddo
return
end subroutine make_lowercase
!************************************************************************
end module physical_constants
More information about the flash-bugs
mailing list