Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module ideal
21 :
22 : use const_def, only: dp
23 : use math_lib
24 : use auto_diff
25 : use eos_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: get_ideal_eos_results, get_ideal_alfa, get_ideal_for_eosdt
31 :
32 : contains
33 :
34 0 : subroutine get_ideal_alfa( &
35 : rq, logRho, logT, Z, abar, zbar, &
36 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
37 : ierr)
38 : use const_def, only: dp
39 : use eos_blend
40 : type (EoS_General_Info), pointer :: rq
41 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
42 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
43 : integer, intent(out) :: ierr
44 :
45 0 : alfa = 0d0
46 0 : d_alfa_dlogRho = 0d0
47 0 : d_alfa_dlogT = 0d0
48 0 : end subroutine get_ideal_alfa
49 :
50 0 : subroutine get_ideal_for_eosdt(handle, dbg, Z, X, abar, zbar, species, chem_id, net_iso, xa, &
51 0 : rho, logRho, T, logT, remaining_fraction, res, d_dlnd, d_dlnT, d_dxa, skip, ierr)
52 : integer, intent(in) :: handle
53 : logical, intent(in) :: dbg
54 : real(dp), intent(in) :: &
55 : Z, X, abar, zbar, remaining_fraction
56 : integer, intent(in) :: species
57 : integer, pointer :: chem_id(:), net_iso(:)
58 : real(dp), intent(in) :: xa(:)
59 : real(dp), intent(in) :: rho, logRho, T, logT
60 : real(dp), intent(inout), dimension(nv) :: res, d_dlnd, d_dlnT
61 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
62 : logical, intent(out) :: skip
63 : integer, intent(out) :: ierr
64 : type (EoS_General_Info), pointer :: rq
65 :
66 : rq => eos_handles(handle)
67 :
68 : call get_ideal_eos_results(rq, Z, X, abar, zbar, rho, logRho, T, logT, species, chem_id, xa, &
69 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
70 0 : skip = .false.
71 :
72 : ! zero all components
73 0 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
74 0 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
75 0 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
76 :
77 : ! mark this one
78 0 : res(i_frac_ideal) = 1.0d0
79 :
80 0 : end subroutine get_ideal_for_eosdt
81 :
82 0 : subroutine get_ideal_eos_results( &
83 : rq, Z, X, abar, zbar, Rho, logRho, T, logT, &
84 0 : species, chem_id, xa, res, d_dlnd, d_dlnT, d_dxa, ierr)
85 : type (EoS_General_Info), pointer :: rq
86 : real(dp), intent(in) :: Z, X, abar, zbar
87 : real(dp), intent(in) :: Rho, logRho, T, logT
88 : integer, intent(in) :: species
89 : integer, pointer :: chem_id(:)
90 : real(dp), intent(in) :: xa(:)
91 : integer, intent(out) :: ierr
92 : real(dp), intent(out), dimension(nv) :: res, d_dlnd, d_dlnT
93 : real(dp), intent(out), dimension(nv, species) :: d_dxa
94 :
95 : ierr = 0
96 :
97 : call ideal_eos( &
98 : T, Rho, X, abar, zbar, &
99 : species, chem_id, xa, &
100 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
101 :
102 : ! composition derivatives not provided
103 0 : d_dxa = 0
104 :
105 0 : end subroutine get_ideal_eos_results
106 :
107 0 : subroutine ideal_eos( &
108 : temp_in, den_in, Xfrac, abar, zbar, &
109 0 : species, chem_id, xa, &
110 : res, d_dlnd, d_dlnT, d_dxa, ierr)
111 :
112 : use eos_def
113 : use const_def, only: dp
114 : use utils_lib, only: is_bad
115 : use chem_def, only: chem_isos
116 : use ion_offset, only: compute_ion_offset
117 : use skye_ideal
118 : use skye_thermodynamics
119 : use auto_diff
120 :
121 : integer :: j
122 : integer, intent(in) :: species
123 : integer, pointer :: chem_id(:)
124 : real(dp), intent(in) :: xa(:)
125 : real(dp), intent(in) :: temp_in, den_in
126 : real(dp), intent(in) :: Xfrac, abar, zbar
127 : integer, intent(out) :: ierr
128 : real(dp), intent(out), dimension(nv) :: res, d_dlnd, d_dlnT
129 : real(dp), intent(out), dimension(nv, species) :: d_dxa
130 : real(dp), parameter :: mass_fraction_limit = 1d-10
131 :
132 : integer :: relevant_species
133 : type(auto_diff_real_2var_order3) :: temp, den
134 0 : real(dp) :: ACMI(species), A(species), ya(species), select_xa(species), norm
135 : type(auto_diff_real_2var_order3) :: etaele, xnefer, phase, latent_ddlnT, latent_ddlnRho
136 : type(auto_diff_real_2var_order3) :: F_rad, F_ideal_ion, F_ele, F_coul
137 :
138 : ! Set up
139 0 : ierr = 0
140 0 : F_ele = 0d0
141 0 : F_coul = 0d0
142 :
143 : ! No electrons, so extreme negative chemical potential
144 0 : etaele = -1d99
145 0 : xnefer = 1d-20
146 :
147 : ! no latent heat, no phase
148 0 : latent_ddlnT = 0d0
149 0 : latent_ddlnRho = 0d0
150 0 : phase = 0d0
151 :
152 : ! Construct rho,T and partials
153 0 : temp = temp_in
154 0 : temp%d1val1 = 1d0
155 0 : den = den_in
156 0 : den%d1val2 = 1d0
157 :
158 : ! Count and pack relevant species for Coulomb corrections. Relevant means mass fraction above limit.
159 0 : relevant_species = 0
160 0 : norm = 0d0
161 0 : do j=1,species
162 0 : if (xa(j) > mass_fraction_limit) then
163 0 : relevant_species = relevant_species + 1
164 0 : ACMI(relevant_species) = chem_isos% W(chem_id(j))
165 0 : A(relevant_species) = chem_isos% Z_plus_N(chem_id(j))
166 0 : select_xa(relevant_species) = xa(j)
167 0 : norm = norm + xa(j)
168 : end if
169 : end do
170 :
171 : ! Normalize
172 0 : do j=1,relevant_species
173 0 : select_xa(j) = select_xa(j) / norm
174 : end do
175 :
176 : ! Compute number fractions
177 : norm = 0d0
178 0 : do j=1,relevant_species
179 0 : ya(j) = select_xa(j) / A(j)
180 0 : norm = norm + ya(j)
181 : end do
182 0 : do j=1,relevant_species
183 0 : ya(j) = ya(j) / norm
184 : end do
185 :
186 : ! Ideal ion free energy, only depends on abar
187 0 : F_ideal_ion = compute_F_ideal_ion(temp, den, abar, relevant_species, ACMI, ya)
188 :
189 : ! The compute_ion_offset correction should only be applied to ionized matter.
190 : ! Right now, ideal always assumes neutral, so leave this off unless/until we consider
191 : ! adding assumption of full ionizatino for high T regions on ideal gas.
192 : !F_ideal_ion = F_ideal_ion + compute_ion_offset(species, xa, chem_id) ! Offset so ion ground state energy is zero.
193 :
194 : ! Radiation free energy, independent of composition
195 0 : F_rad = compute_F_rad(temp, den)
196 :
197 :
198 : call pack_for_export(F_ideal_ion, F_coul, F_rad, F_ele, temp, den, xnefer, etaele, abar, zbar, &
199 0 : phase, latent_ddlnT, latent_ddlnRho, res, d_dlnd, d_dlnT, ierr)
200 0 : if(ierr/=0) return
201 :
202 0 : res(i_mu) = abar ! ideal assumes neutral matter, whereas pack_for_export assumes ionized matter. So we patch it up here.
203 :
204 0 : end subroutine ideal_eos
205 :
206 : end module ideal
|
