FreeEOS Code
Brought to you by:
airwin
Menu
▾
▴
Tree [9a5978] master / src / History
Read Me
The input arguments to the generic free_eos routine include abundance
and option suite information that together specify the free-energy
model that is calculated by the free_eos library. And the results of
the EOS calculation depend on certain first and second partial
derivatives of that free-energy model. I have derived those partial
derivatives based on the principles of calculus, e.g., the chain rule,
implicit differentiation, etc., and, in particular those derivatives
are *not* based on numerical difference approximations. Thus, if the
derivation *and* implementation of that derivation of all relevant
partial derivatives is correct for a given free-energy model, then
subject to small significance loss limits the iterative Newton-Raphson
solution of the EOS should converge quadratically to the exact
solution corresponding to that free-energy model, and that solution
should exhibit exact thermodynamic consistency.
Thus, when such quadratic convergence of the Newton-Raphson iteration
is not observed or the EOS results are not thermodynamically
consistent, then such symptoms typically imply there is an error in
the derivation or implementation of one of the many partial
derivatives of the free-energy model, and the rest of this file
documents how I use numerical difference tests to find which of the
derivatives is the culprit.
The numerical difference tests I have implemented are the
exchange_test, fermi_dirac_test, and free_eos_test executables which
respectively test partial derivatives of exchange- and
free-electron-related Fermi-Dirac integrals and overall thermodynamic
functions that have been implemented as part of the free_eos library.
(These partial differientation test executables are documented in more
detail in ../test_FreeEOS/README.simple_test and
../test_FreeEOS/README.free_eos_test.)
Assuming for the given free-energy model that free_eos_test shows
there is a derivative error in one or more thermodynamic functions,
then here are further tests that can be applied to help find *which*
derivative that those thermodynamic function derivatives depends on is
the source of the error.
* Test the derivatives of the exchange- and free-electron-related Fermi-Dirac
integrals *that are adopted for that free-energy model* using exchange_test and
fermi_dirac test.
* If the test failure is for kif /= 0, then try to replicate the issue for kif = 0
* Try to replicate the issue for simplified abundances, e.g., pure H, pure He,
or pure Fe.
* Locally modify free_eos_detailed to drop various parts of the free-energy model to
determine which component of that free-energy has a derivative error.
* Locally hack free_eos_detailed so that free_eos_test effectively ends up
testing other partial derivative results with different independent variables!
I classify all such test code as a hack because it closely depends
on knowledge of the pattern of free_eos calls in free_eos_test (a
first call evaluating thermodynamic functions and their derivatives
wrt match_variable and tl and 32 subsequent calls evaluating centred
numerical differences in match_variable and tl). The hack is to save
the initial match_variable and tl from the first of those calls, and
use the match_variable and tl values from the further 32 calls to
derive the logarithmic changes in two arbitrary independent
variables that are being tested with centred numerical differences
(so 4 calls per set of centred differences in match_variable and tl)
for 8 different orders of magnitude in those differences. Note, for
those 32 calls, match_variable and tl are restored to their initial
values to provide clean numerical difference results for whatever
two (specified by itemp[12] in free_eos_detailed) independent
variables are chosen for these hacked partial derivative tests.
Currently free_eos_detailed and the subroutines it calls illustrate
this testing technique for three different mutually exclusive sets
of partial derivatives where one of the logical parameters
debug_dv_aux, debug_aux_dv, or debug_jacobian (all of which are
normally .false.) is locally modified to be .true. in order that the
hacked test code associated with that logical variable is executed.
In these three cases that have been implemented, the partial
derivatives being tested are calculated by eos_jacobian. Therefore,
to insure a derivative test that is for independent variables that
are a fairly good approximation to what would have been the final
result for the given physical conditions and free-energy model, I
have chosen to call eos_cold_start before eos_jacobian to provide
fairly realistic independent variables for these tests.
Further notes on each of these testing hacks that are labelled by
the associated logical variable.
+ debug_dv_aux = .true.: this hack uses numerical differences to
test the partial derivatives of components of dv (designated by
jtemp[1-6] in free_eos_detailed which corresponds to the species H+,
He+, He++, C+, H2, and H2+ although other choices are possible)
wrt to components of aux_old, fl, or tl (designated by itemp[12]
in free_eos_detailed). dv(aux_old, fl, tl) is the vector of
non-ideal equilibrium constants which are needed in eos_calc to
help calculate number densities of all species; aux_old is the
vector of old auxiliary variables fl is the EFF logarithmic
degeneracy parameter, and tl is ln T; and aux_new(dv(aux_old, fl,
tl), fl, tl) = aux_new(aux_old, fl, tl) is the vector of new
auxiliary variables which are each defined as a uniquely weighted
sum over the number densities of species that are calculated by
eos_calc. (See the "solution.pdf" paper for further details
concerning the use of auxiliary variables to determine the EOS
solution corresponding to non-ideal free-energy minimization.)
+ debug_aux_dv = .true.: this hack uses numerical differences to
test the partial derivatives of components of aux_new(dv, fl,
tl) (designated by jtemp[1-6] in free_eos_detailed with sets of
choices corresponding to all aux components possible with EOS1)
wrt to components of dv (designated by itemp[12] in
free_eos_detailed where sets of itemp[12] values correspond to
the species H+, He+, He++, C+, H2, and H2+ although other
choices are possible). See discussion of debug_dv_aux =
.true. for more discussion of the partial derivatives being
debugged for this case.
+ debug_jacobian = .true.: this hack uses numerical differences to
test the partial derivatives of components of faux (the RHS
vector for the system of linary equations solved in a
Newton-Raphson step where the components of this vector are
designated by sets of jtemp[1-6] variables in free_eos_detailed) wrt
to components of aux_old, fl, or tl (designated by itemp[12] in
free_eos_detailed). Except for its last component which is a
special case, faux depends on the difference between
aux_new(aux_old, fl, tl) and aux_old so ultimately these partial
derivatives depend on the partials of aux_new(aux_old, fl, tl) -
aux_old wrt components of aux_old (whose negative values are
stored in the jacobian matrix) and the partial derivatives of
aux_new(aux_old, fl, tl) wrt fl and tl.
N.B. eos_jacobian normally calculates the partial derivatives of
aux_new(dv, fl, tl) wrt fl and tl so an extra transformation
(still not completely debugged, see FIXME notice in
free_eos_detailed) needs to be applied to determine the desired
partials of aux_new(aux_old, fl, tl) wrt fl or tl. Note the
equivalent calculations in eos_tqft are implemented in a
completely different (ifnr = 0) way but are done correctly so
for normal use (i.e., none of the debug_* options are set to
.true.) the derivatives tested by free_eos_test are calculated
correctly.