MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_interpolation.t
Go to the documentation of this file.
3  use mod_comm_lib, only: mpistop
4  implicit none
5 
6 contains
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  integer :: ix^D
144  double precision :: dxc^D,xd^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  integer :: ix^D,iw
164  double precision :: factor(0:1^D&)
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  integer :: inblock,ixO^L,j
186  integer :: ixb^D,ixi^D
187  double precision :: xc(0:1^D&,1:ndim),wc(0:1^D&,nw)
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 
215 end module mod_interpolation
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
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)