Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2017-2019 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 eosdt_support
21 : use eos_def
22 : use const_def, only: avo, crad, ln10, arg_not_provided, mp, kerg, dp, qp, one_sixth
23 : use utils_lib, only: is_bad, mesa_error
24 : use math_lib
25 :
26 : implicit none
27 :
28 : integer, parameter :: sz = sz_per_eos_point
29 :
30 : contains
31 :
32 1685204 : subroutine Do_EoS_Interpolations( &
33 : nvlo, nvhi, n, nx, x, ny, y, fin1, i, j, &
34 : x0, xget, x1, y0, yget, y1, &
35 : fval, df_dx, df_dy, ierr)
36 : integer, intent(in) :: nvlo, nvhi, n, nx, ny
37 : real(dp), intent(in) :: x(:) ! (nx)
38 : real(dp), intent(in) :: y(:) ! (ny)
39 : real(dp), intent(in), pointer :: fin1(:) ! =(4,n,nx,ny)
40 : integer, intent(in) :: i, j ! target cell in f
41 : real(dp), intent(in) :: x0, xget, x1 ! x0 <= xget <= x1; x0 = xs(i), x1 = xs(i+1)
42 : real(dp), intent(in) :: y0, yget, y1 ! y0 <= yget <= y1; y0 = ys(j), y1 = ys(j+1)
43 : real(dp), intent(inout), dimension(nv) :: fval, df_dx, df_dy
44 : integer, intent(out) :: ierr
45 :
46 1685204 : real(dp) :: xp, xpi, xp2, xpi2, ax, axbar, bx, bxbar, cx, cxi, hx2, cxd, cxdi, hx, hxi
47 1685204 : real(dp) :: yp, ypi, yp2, ypi2, ay, aybar, by, bybar, cy, cyi, hy2, cyd, cydi, hy, hyi
48 1685204 : real(dp) :: sixth_hx2, sixth_hy2, z36th_hx2_hy2
49 1685204 : real(dp) :: sixth_hx, sixth_hxi_hy2, z36th_hx_hy2
50 1685204 : real(dp) :: sixth_hx2_hyi, sixth_hy, z36th_hx2_hy
51 : integer :: k, ip1, jp1
52 1685204 : real(dp), pointer :: fin(:,:,:,:)
53 :
54 : include 'formats'
55 :
56 1685204 : ierr = 0
57 :
58 : fin(1:sz_per_eos_point,1:n,1:nx,1:ny) => &
59 1685204 : fin1(1:sz_per_eos_point*n*nx*ny)
60 :
61 1685204 : hx=x1-x0
62 1685204 : hxi=1d0/hx
63 1685204 : hx2=hx*hx
64 :
65 1685204 : xp=(xget-x0)*hxi
66 :
67 1685204 : xpi=1d0-xp
68 1685204 : xp2=xp*xp
69 1685204 : xpi2=xpi*xpi
70 :
71 1685204 : ax=xp2*(3d0-2d0*xp)
72 1685204 : axbar=1d0-ax
73 :
74 1685204 : bx=-xp2*xpi
75 1685204 : bxbar=xpi2*xp
76 :
77 1685204 : cx=xp*(xp2-1d0)
78 1685204 : cxi=xpi*(xpi2-1d0)
79 1685204 : cxd=3d0*xp2-1d0
80 1685204 : cxdi=-3d0*xpi2+1d0
81 :
82 1685204 : hy=y1-y0
83 1685204 : hyi=1d0/hy
84 1685204 : hy2=hy*hy
85 :
86 1685204 : yp=(yget-y0)*hyi
87 :
88 1685204 : ypi=1d0-yp
89 1685204 : yp2=yp*yp
90 1685204 : ypi2=ypi*ypi
91 :
92 1685204 : ay=yp2*(3d0-2d0*yp)
93 1685204 : aybar=1d0-ay
94 :
95 1685204 : by=-yp2*ypi
96 1685204 : bybar=ypi2*yp
97 :
98 1685204 : cy=yp*(yp2-1d0)
99 1685204 : cyi=ypi*(ypi2-1d0)
100 1685204 : cyd=3d0*yp2-1d0
101 1685204 : cydi=-3d0*ypi2+1d0
102 :
103 1685204 : sixth_hx2 = one_sixth*hx2
104 1685204 : sixth_hy2 = one_sixth*hy2
105 1685204 : z36th_hx2_hy2 = sixth_hx2*sixth_hy2
106 :
107 1685204 : sixth_hx = one_sixth*hx
108 1685204 : sixth_hxi_hy2 = sixth_hy2*hxi
109 1685204 : z36th_hx_hy2 = sixth_hx*sixth_hy2
110 :
111 1685204 : sixth_hx2_hyi = sixth_hx2*hyi
112 1685204 : sixth_hy = one_sixth*hy
113 1685204 : z36th_hx2_hy = sixth_hx2*sixth_hy
114 :
115 1685204 : ip1 = i+1
116 1685204 : jp1 = j+1
117 :
118 3370408 : !$omp simd
119 : do k = nvlo, nvhi
120 : ! bicubic spline interpolation
121 :
122 : ! f(1,i,j) = f(x(i),y(j))
123 : ! f(2,i,j) = d2f/dx2(x(i),y(j))
124 : ! f(3,i,j) = d2f/dy2(x(i),y(j))
125 : ! f(4,i,j) = d4f/dx2dy2(x(i),y(j))
126 :
127 : fval(k) = &
128 : xpi*( &
129 : ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1)) &
130 : +xp*(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1)) &
131 : +sixth_hx2*( &
132 : cxi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+ &
133 : cx*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1))) &
134 : +sixth_hy2*( &
135 : xpi*(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1))+ &
136 : xp*(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1))) &
137 : +z36th_hx2_hy2*( &
138 : cxi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+ &
139 43815304 : cx*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1)))
140 :
141 : ! derivatives of bicubic splines
142 : df_dx(k) = &
143 : hxi*( &
144 : -(ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1)) &
145 : +(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1))) &
146 : +sixth_hx*( &
147 : cxdi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+ &
148 : cxd*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1))) &
149 : +sixth_hxi_hy2*( &
150 : -(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1)) &
151 : +(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1))) &
152 : +z36th_hx_hy2*( &
153 : cxdi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+ &
154 43815304 : cxd*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1)))
155 :
156 : df_dy(k) = &
157 : hyi*( &
158 : xpi*(-fin(1,k,i,j) +fin(1,k,i,jp1))+ &
159 : xp*(-fin(1,k,ip1,j)+fin(1,k,ip1,jp1))) &
160 : +sixth_hx2_hyi*( &
161 : cxi*(-fin(2,k,i,j) +fin(2,k,i,jp1))+ &
162 : cx*(-fin(2,k,ip1,j)+fin(2,k,ip1,jp1))) &
163 : +sixth_hy*( &
164 : xpi*(cydi*fin(3,k,i,j) +cyd*fin(3,k,i,jp1))+ &
165 : xp*(cydi*fin(3,k,ip1,j)+cyd*fin(3,k,ip1,jp1))) &
166 : +z36th_hx2_hy*( &
167 : cxi*(cydi*fin(4,k,i,j) +cyd*fin(4,k,i,jp1))+ &
168 43815304 : cx*(cydi*fin(4,k,ip1,j)+cyd*fin(4,k,ip1,jp1)))
169 :
170 : end do
171 :
172 1685204 : end subroutine Do_EoS_Interpolations
173 :
174 :
175 14940 : subroutine Do_Blend( &
176 : rq, species, Rho, logRho, T, logT, &
177 : alfa_in, d_alfa_dlogT_in, d_alfa_dlogRho_in, linear_blend, &
178 14940 : res_1, d_dlnd_1, d_dlnT_1, d_dxa_1, &
179 14940 : res_2, d_dlnd_2, d_dlnT_2, d_dxa_2, &
180 14940 : res, dlnd, dlnT, d_dxa)
181 : type (EoS_General_Info), pointer :: rq
182 : integer, intent(in) :: species
183 : real(dp), intent(in) :: Rho, logRho, T, logT, &
184 : alfa_in, d_alfa_dlogT_in, d_alfa_dlogRho_in
185 : logical, intent(in) :: linear_blend
186 : real(dp), intent(in), dimension(nv) :: res_1, d_dlnd_1, d_dlnT_1, res_2, d_dlnd_2, d_dlnT_2
187 : real(dp), intent(in), dimension(nv, species) :: d_dxa_1, d_dxa_2
188 : real(dp), intent(out), dimension(nv) :: res, dlnd, dlnT
189 : real(dp), intent(out), dimension(nv, species) :: d_dxa
190 :
191 14940 : real(dp) :: alfa0, d_alfa0_dlnT, d_alfa0_dlnd
192 14940 : real(dp) :: alfa, beta, A, &
193 14940 : d_alfa_dlnd, d_alfa_dlnT, &
194 14940 : d_beta_dlnd, d_beta_dlnT
195 : integer :: j, k
196 :
197 14940 : if (.not. linear_blend) then
198 :
199 : ! smooth the transitions near alfa = 0.0 and 1.0
200 : ! quintic smoothing function with 1st and 2nd derivs = 0 at ends
201 :
202 14940 : alfa0 = alfa_in
203 14940 : d_alfa0_dlnT = d_alfa_dlogT_in/ln10
204 14940 : d_alfa0_dlnd = d_alfa_dlogRho_in/ln10
205 14940 : alfa = -alfa0*alfa0*alfa0*(-10d0 + alfa0*(15d0 - 6d0*alfa0))
206 14940 : A = 30d0*(alfa0 - 1d0)*(alfa0 - 1d0)*alfa0*alfa0
207 14940 : d_alfa_dlnd = A*d_alfa0_dlnd
208 14940 : d_alfa_dlnT = A*d_alfa0_dlnT
209 :
210 : else
211 :
212 0 : alfa = alfa_in
213 0 : d_alfa_dlnT = d_alfa_dlogT_in/ln10
214 0 : d_alfa_dlnd = d_alfa_dlogRho_in/ln10
215 :
216 : end if
217 :
218 14940 : beta = 1d0 - alfa
219 14940 : d_beta_dlnT = -d_alfa_dlnT
220 14940 : d_beta_dlnd = -d_alfa_dlnd
221 :
222 403380 : do j=1,nv
223 388440 : res(j) = alfa*res_1(j) + beta*res_2(j)
224 : dlnd(j) = &
225 : alfa*d_dlnd_1(j) + beta*d_dlnd_2(j) + &
226 388440 : d_alfa_dlnd*res_1(j) + d_beta_dlnd*res_2(j)
227 : dlnT(j) = &
228 : alfa*d_dlnT_1(j) + beta*d_dlnT_2(j) + &
229 403380 : d_alfa_dlnT*res_1(j) + d_beta_dlnT*res_2(j)
230 : end do
231 :
232 134454 : do k = 1, species
233 3241818 : do j = 1, nv
234 : d_dxa(j,k) = &
235 3226878 : alfa*d_dxa_1(j,k) + beta*d_dxa_2(j,k) ! ignores blend derivatives
236 : end do
237 : end do
238 :
239 14940 : end subroutine Do_Blend
240 :
241 : end module eosdt_support
|
