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