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