MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_interpolation.t
Go to the documentation of this file.
3 use mod_comm_lib, only: mpistop
4 implicit none
5
6contains
7
8 subroutine interp_linear(x_table,y_table,n_table,x_itp,y_itp,n_itp)
9 ! linear interpolation
10 ! 1D method
11
12 integer :: n_table,n_itp
13 double precision :: x_table(n_table),y_table(n_table)
14 double precision :: x_itp(n_itp),y_itp(n_itp)
15
16 integer :: i,j
17
18 !if (x_table(1)>x_itp(1) .or. x_table(n_table)<x_itp(n_itp)) then
19 ! call mpistop("out of interpolation table")
20 !endif
21
22 do i=1,n_itp
23 if(x_itp(i)<x_table(1)) then
24 y_itp(i)=y_table(1)
25 exit
26 end if
27 if(x_itp(i)>=x_table(n_table)) then
28 y_itp(i)=y_table(n_table)
29 exit
30 end if
31 loopt: do j=1,n_table-1
32 if (x_itp(i)>=x_table(j) .and. x_itp(i)<x_table(j+1)) then
33 y_itp(i)=y_table(j) + (x_itp(i)-x_table(j))*&
34 (y_table(j+1)-y_table(j))/(x_table(j+1)-x_table(j))
35 exit loopt
36 endif
37 enddo loopt
38 enddo
39
40 end subroutine interp_linear
41
42 subroutine interp_cubic_spline(x_table,y_table,n_table,x_itp,y_itp,n_itp)
43 ! interpolation function fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
44 ! first order derivative and second order derivative is continous
45 ! 1D method
46
47 integer :: n_table,n_itp
48 double precision :: x_table(n_table),y_table(n_table)
49 double precision :: x_itp(n_itp),y_itp(n_itp)
50
51 double precision :: at(n_table),bt(n_table),ct(n_table),dt(n_table)
52 double precision :: dxi
53 integer :: i,j
54
55 if (x_table(1)>x_itp(1) .or. x_table(n_table)<x_itp(n_itp)) then
56 call mpistop("out of interpolation table")
57 endif
58
59 call get_cubic_para(x_table,y_table,at,bt,ct,dt,n_table)
60
61 do i=1,n_itp
62 loopt: do j=1,n_table-1
63 if (x_itp(i)>=x_table(j) .and. x_itp(i)<x_table(j+1)) then
64 dxi=x_itp(i)-x_table(j)
65 y_itp(i)=at(j)+bt(j)*dxi+ct(j)*dxi**2+dt(j)*dxi**3
66 exit loopt
67 endif
68 enddo loopt
69 enddo
70
71 if (x_itp(n_itp)==x_table(n_table)) y_itp(n_itp)=y_table(n_table)
72
73 end subroutine interp_cubic_spline
74
75 subroutine get_cubic_para(xi,yi,ai,bi,ci,di,ni)
76 ! get parameters for interpolation
77
78 integer :: ni,i
79 double precision :: xi(ni),yi(ni)
80 double precision :: ai(ni),bi(ni),ci(ni),di(ni)
81
82 double precision :: matrix1(ni,ni),ri1(ni)
83 double precision :: matrix2(ni,ni),ri2(ni)
84 double precision :: mi(ni),hi(ni)
85
86 ! build equations
87 matrix1=0.d0
88 matrix2=0.d0
89
90 do i=1,ni-1
91 hi(i)=xi(i+1)-xi(i)
92 enddo
93
94 matrix1(1,1)=1.0
95 matrix1(ni,ni)=1.0
96 ri1(1)=0.d0
97 ri1(ni)=0.d0
98 do i=2,ni-1
99 matrix1(i-1,i)=hi(i-1)
100 matrix1(i,i)=2*(hi(i-1)+hi(i))
101 matrix1(i+1,i)=hi(i)
102 ri1(i)=(yi(i+1)-yi(i))/hi(i)-(yi(i)-yi(i-1))/hi(i-1)
103 enddo
104 ri1=ri1*6
105
106 ! solve equations
107 do i=1,ni
108 matrix2(i,i)=1.d0
109 enddo
110
111 matrix2(2,1)=matrix1(2,1)/matrix1(1,1)
112 ri2(1)=ri1(1)/matrix1(1,1)
113
114 do i=2,ni-1
115 matrix2(i+1,i)=matrix1(i+1,i)/(matrix1(i,i)-matrix1(i-1,i)*matrix2(i,i-1))
116 enddo
117
118 do i=2,ni
119 ri2(i)=ri1(i)-matrix1(i-1,i)*ri2(i-1)
120 ri2(i)=ri2(i)/(matrix1(i,i)-matrix1(i-1,i)*matrix2(i,i-1))
121 enddo
122
123 mi(ni)=ri2(ni)
124 do i=ni-1,1,-1
125 mi(i)=ri2(i)-matrix2(i+1,i)*mi(i+1)
126 enddo
127
128 ! get parameters for interpolation
129 ai=yi
130 ci=mi(1:ni)
131 do i=1,ni-1
132 bi(i)=(yi(i+1)-yi(i))/hi(i)-hi(i)*mi(i)/2.0-hi(i)*(mi(i+1)-mi(i))/6.0
133 di(i)=(mi(i+1)-mi(i))/(6.0*hi(i))
134 enddo
135
136 end subroutine get_cubic_para
137
138 subroutine get_interp_factor_linear(xc,xp,factor)
139 ! get factor for linear interpolation
140 ! multi-D method
141
142 double precision :: xc(0:1^D&,1:ndim),xp(1:ndim),factor(0:1^D&)
143 double precision :: dxc^D,xd^D
144 integer :: ix^D
145
146 ^d&dxc^d=xc(1^dd&,^d)-xc(0^dd&,^d)\
147 ^d&xd^d=(xp(^d)-xc(0^dd&,^d))/dxc^d\
148 ! interpolation factor
149 {do ix^d=0,1\}
150 factor(ix^d)={abs(1-ix^d-xd^d)*}
151 {enddo\}
152
153 end subroutine get_interp_factor_linear
154
155 subroutine interp_linear_multid(xc,xp,wc,wp,nwc)
156 ! get point values from nearby cells via linear interpolation
157 ! multi-D method
158
159 integer :: nwc
160 double precision :: xc(0:1^D&,1:ndim),xp(1:ndim)
161 double precision :: wc(0:1^D&,1:nwc),wp(1:nwc)
162
163 double precision :: factor(0:1^D&)
164 integer :: ix^D,iw
165
166 call get_interp_factor_linear(xc,xp,factor)
167
168 wp=0.d0
169 do iw=1,nwc
170 {do ix^db=0,1\}
171 wp(iw)=wp(iw)+wc(ix^d,iw)*factor(ix^d)
172 {enddo\}
173 enddo
174
175 end subroutine interp_linear_multid
176
177 subroutine interp_from_grid_linear(igrid,xp,wp)
178 ! get point values from given grid via linear interpolation
179 ! multi-D method
180
181 integer :: igrid
182 double precision :: xp(1:ndim),wp(1:nw)
183
184 double precision :: dxb^D,xb^L
185 double precision :: xc(0:1^D&,1:ndim),wc(0:1^D&,nw)
186 integer :: inblock,ixO^L,j
187 integer :: ixb^D,ixi^D
188
189 ! block/grid boundaries
190 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
191 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
192
193 ! whether or not next point is in this block/grid
194 inblock=0
195 {if (xp(^d)>=xbmin^d .and. xp(^d)<xbmax^d) inblock=inblock+1\}
196 if (inblock/=ndim) then
197 call mpistop('Interpolation error: given point is not in given grid')
198 endif
199
200 ! cell indexes for the point
201 ^d&dxb^d=rnode(rpdx^d_,igrid)\
202 ^d&ixomin^d=ixmlo^d\
203 ^d&ixb^d=floor((xp(^d)-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
204
205 ! nearby cells for interpolation
206 {do ixi^d=0,1\}
207 xc(ixi^d,:)=ps(igrid)%x(ixb^d+ixi^d,:)
208 wc(ixi^d,:)=ps(igrid)%w(ixb^d+ixi^d,:)
209 {enddo\}
210
211 call interp_linear_multid(xc,xp,wc,wp,nw)
212
213 end subroutine interp_from_grid_linear
214
215end module mod_interpolation
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
subroutine interp_cubic_spline(x_table, y_table, n_table, x_itp, y_itp, n_itp)
subroutine get_interp_factor_linear(xc, xp, factor)
subroutine get_cubic_para(xi, yi, ai, bi, ci, di, ni)
subroutine interp_from_grid_linear(igrid, xp, wp)
subroutine interp_linear_multid(xc, xp, wc, wp, nwc)
subroutine interp_linear(x_table, y_table, n_table, x_itp, y_itp, n_itp)