Line data Source code
1 : module interp_2D_support
2 : use const_def, only: dp, pi
3 : use interp_2d_lib_db
4 : use interp_2d_lib_sg
5 : use math_lib
6 : use utils_lib, only: mesa_error
7 :
8 : implicit none
9 :
10 : integer, parameter :: num_xpts = 31
11 : integer, parameter :: num_ypts = 35
12 : integer, parameter :: sz_per_pt = 4
13 :
14 : real, parameter :: pi_sg = 3.14159
15 :
16 : real(dp), target :: f_db_ary(sz_per_pt*num_xpts*num_ypts)
17 : real(dp), pointer :: f_db(:, :, :), f_db1(:)
18 : real, target :: f_sg_ary(sz_per_pt*num_xpts*num_ypts)
19 : real, pointer :: f_sg(:, :, :), f_sg1(:)
20 : integer :: ilinx ! =1: x grid is "nearly" equally spaced
21 : integer :: iliny ! =1: y grid is "nearly" equally spaced
22 : real, target :: xpts_sg_ary(num_xpts), ypts_sg_ary(num_ypts)
23 : real, pointer :: xpts_sg(:), ypts_sg(:)
24 : real(dp), target :: xpts_db_ary(num_xpts), ypts_db_ary(num_ypts)
25 : real(dp), pointer :: xpts_db(:), ypts_db(:)
26 :
27 : contains
28 :
29 4340 : real(dp) function test_fcn_db(x, y)
30 : real(dp), intent(in) :: x, y
31 4340 : real(dp) :: r, r0
32 4340 : r = dsqrt(x*x + y*y)
33 4340 : r0 = pi/2
34 4340 : test_fcn_db = exp(-r/r0)*sin(x)*cos(y)
35 4340 : end function test_fcn_db
36 :
37 2170 : real function test_fcn_sg(x, y)
38 : real, intent(in) :: x, y
39 2170 : real :: r, r0
40 2170 : r = sqrt(x*x + y*y)
41 2170 : r0 = pi_sg/2
42 2170 : test_fcn_sg = real(test_fcn_db(dble(x), dble(y)))
43 2170 : end function test_fcn_sg
44 :
45 2 : subroutine get_2D_test_values_db
46 :
47 : integer :: i, j
48 2 : real(dp) :: x, xmin, xmax, dx
49 2 : real(dp) :: y, ymin, ymax, dy
50 :
51 2 : xmin = -pi; xmax = pi; dx = (xmax - xmin)/(num_xpts - 1)
52 2 : ymin = -pi; ymax = pi; dy = (ymax - ymin)/(num_ypts - 1)
53 :
54 64 : do i = 1, num_xpts
55 62 : x = xmin + (i - 1)*dx
56 62 : xpts_db(i) = x
57 2234 : do j = 1, num_ypts
58 2170 : y = ymin + (j - 1)*dy
59 2170 : if (i == 1) ypts_db(j) = y
60 2232 : f_db(1, i, j) = test_fcn_db(x, y)
61 : end do
62 : end do
63 :
64 2 : end subroutine get_2D_test_values_db
65 :
66 2 : subroutine get_2D_test_values_sg
67 :
68 : integer :: i, j
69 2 : real :: x, xmin, xmax, dx
70 2 : real :: y, ymin, ymax, dy
71 :
72 2 : xmin = -pi_sg; xmax = pi_sg; dx = (xmax - xmin)/(num_xpts - 1)
73 2 : ymin = -pi_sg; ymax = pi_sg; dy = (ymax - ymin)/(num_ypts - 1)
74 :
75 64 : do i = 1, num_xpts
76 62 : x = xmin + (i - 1)*dx
77 62 : xpts_sg(i) = x
78 2234 : do j = 1, num_ypts
79 2170 : y = ymin + (j - 1)*dy
80 2170 : if (i == 1) ypts_sg(j) = y
81 2232 : f_sg(1, i, j) = test_fcn_sg(x, y)
82 : end do
83 : end do
84 :
85 2 : end subroutine get_2D_test_values_sg
86 :
87 2 : subroutine setup_to_interp_2D_db(bicub_flag)
88 : logical, intent(in) :: bicub_flag
89 : integer :: nf2 ! 2nd dimension of f, nf2.ge.nx
90 : integer :: ibcxmin ! bc flag for x=xmin
91 72 : real(dp) :: bcxmin(num_ypts) ! bc data vs. y at x=xmin
92 : integer :: ibcxmax ! bc flag for x=xmax
93 72 : real(dp) :: bcxmax(num_ypts) ! bc data vs. y at x=xmax
94 : integer :: ibcymin ! bc flag for y=ymin
95 64 : real(dp) :: bcymin(num_xpts) ! bc data vs. x at y=ymin
96 : integer :: ibcymax ! bc flag for y=ymax
97 64 : real(dp) :: bcymax(num_xpts) ! bc data vs. x at y=ymax
98 : integer :: ier ! =0 on exit if there is no error.
99 :
100 2 : nf2 = num_xpts
101 :
102 2 : if (bicub_flag) then
103 : !..just use "not a knot" bc's
104 1 : ibcxmin = 0; bcxmin = 0
105 1 : ibcxmax = 0; bcxmax = 0
106 1 : ibcymin = 0; bcymin = 0
107 1 : ibcymax = 0; bcymax = 0
108 : call interp_mkbicub_db(xpts_db, num_xpts, ypts_db, num_ypts, f_db1, nf2, &
109 1 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ier)
110 : else
111 1 : call interp_mkbipm_db(xpts_db, num_xpts, ypts_db, num_ypts, f_db1, nf2, ier)
112 : end if
113 :
114 2 : if (ier /= 0) then
115 0 : write (*, *) 'error'
116 0 : call mesa_error(__FILE__, __LINE__)
117 : end if
118 :
119 2 : end subroutine setup_to_interp_2D_db
120 :
121 2 : subroutine setup_to_interp_2D_sg(bicub_flag)
122 : logical, intent(in) :: bicub_flag
123 : integer :: nf2 ! 2nd dimension of f, nf2.ge.nx
124 : integer :: ibcxmin ! bc flag for x=xmin
125 72 : real :: bcxmin(num_ypts) ! bc data vs. y at x=xmin
126 : integer :: ibcxmax ! bc flag for x=xmax
127 72 : real :: bcxmax(num_ypts) ! bc data vs. y at x=xmax
128 : integer :: ibcymin ! bc flag for y=ymin
129 64 : real :: bcymin(num_xpts) ! bc data vs. x at y=ymin
130 : integer :: ibcymax ! bc flag for y=ymax
131 64 : real :: bcymax(num_xpts) ! bc data vs. x at y=ymax
132 : integer :: ier ! =0 on exit if there is no error.
133 :
134 2 : nf2 = num_xpts
135 :
136 2 : if (bicub_flag) then
137 : !..just use "not a knot" bc's
138 1 : ibcxmin = 0; bcxmin = 0
139 1 : ibcxmax = 0; bcxmax = 0
140 1 : ibcymin = 0; bcymin = 0
141 1 : ibcymax = 0; bcymax = 0
142 : call interp_mkbicub_sg(xpts_sg, num_xpts, ypts_sg, num_ypts, f_sg1, nf2, &
143 1 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ier)
144 : else
145 1 : call interp_mkbipm_sg(xpts_sg, num_xpts, ypts_sg, num_ypts, f_sg1, nf2, ier)
146 : end if
147 :
148 2 : if (ier /= 0) then
149 0 : write (*, *) 'error'
150 0 : call mesa_error(__FILE__, __LINE__)
151 : end if
152 :
153 2 : end subroutine setup_to_interp_2D_sg
154 :
155 12 : subroutine eval_2D_interp_db(bicub_flag, x, y, z, dz_dx, dz_dy)
156 : logical, intent(in) :: bicub_flag
157 : real(dp), intent(in) :: x, y
158 : real(dp), intent(OUT) :: z, dz_dx, dz_dy
159 :
160 : integer :: ict(6) ! code specifying output desired
161 84 : real(dp) :: fval(6) ! results
162 : integer :: nf2
163 : integer :: ier ! error code =0 ==> no error
164 :
165 12 : nf2 = num_xpts
166 12 : ier = 0
167 12 : fval = 0
168 :
169 12 : if (bicub_flag) then
170 6 : ict(1) = 1 ! want f at (x,y)
171 6 : ict(2) = 1 ! want df_dx at (x,y)
172 6 : ict(3) = 1 ! want df_dy at (x,y)
173 6 : ict(4) = 0 ! skip d2f_dx2
174 6 : ict(5) = 0 ! skip d2f_dy2
175 6 : ict(6) = 0 ! skip d2f_dx_dy
176 6 : call interp_evbicub_db(x, y, xpts_db, num_xpts, ypts_db, num_ypts, ilinx, iliny, f_db1, nf2, ict, fval, ier)
177 6 : z = fval(1)
178 6 : dz_dx = fval(2)
179 6 : dz_dy = fval(3)
180 : else
181 6 : call interp_evbipm_db(x, y, xpts_db, num_xpts, ypts_db, num_ypts, f_db1, nf2, z, ier)
182 6 : dz_dx = 0
183 6 : dz_dy = 0
184 : end if
185 :
186 12 : if (ier /= 0) then
187 0 : write (*, *) 'error'
188 0 : call mesa_error(__FILE__, __LINE__)
189 : end if
190 :
191 12 : end subroutine eval_2D_interp_db
192 :
193 12 : subroutine eval_2D_interp_sg(bicub_flag, x, y, z, dz_dx, dz_dy)
194 : logical, intent(in) :: bicub_flag
195 : real, intent(in) :: x, y
196 : real, intent(OUT) :: z, dz_dx, dz_dy
197 :
198 : integer :: ict(6) ! code specifying output desired
199 84 : real :: fval(6) ! results
200 : integer :: nf2
201 : integer :: ier ! error code =0 ==> no error
202 :
203 12 : nf2 = num_xpts
204 :
205 12 : ier = 0
206 12 : fval = 0
207 12 : if (bicub_flag) then
208 6 : ict(1) = 1 ! want f at (x,y)
209 6 : ict(2) = 1 ! want df_dx at (x,y)
210 6 : ict(3) = 1 ! want df_dy at (x,y)
211 6 : ict(4) = 0 ! skip d2f_dx2
212 6 : ict(5) = 0 ! skip d2f_dy2
213 6 : ict(6) = 0 ! skip d2f_dx_dy
214 6 : call interp_evbicub_sg(x, y, xpts_sg, num_xpts, ypts_sg, num_ypts, ilinx, iliny, f_sg1, nf2, ict, fval, ier)
215 6 : z = fval(1)
216 6 : dz_dx = fval(2)
217 6 : dz_dy = fval(3)
218 : else
219 6 : call interp_evbipm_sg(x, y, xpts_sg, num_xpts, ypts_sg, num_ypts, f_sg1, nf2, z, ier)
220 6 : dz_dx = 0
221 6 : dz_dy = 0
222 : end if
223 :
224 12 : if (ier /= 0) then
225 0 : write (*, *) 'error'
226 0 : call mesa_error(__FILE__, __LINE__)
227 : end if
228 :
229 12 : end subroutine eval_2D_interp_sg
230 :
231 : end module interp_2D_support
|
