64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add Lagrange interpolation subroutine and update calls in prolongrestrict modules

Browse Source
This commit is contained in:
ianchb
2026-02-27 12:46:39 +08:00
parent 94d236385d
commit c001939461
3 changed files with 30 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
AMSS_NCKU_source/fmisc.f90
View File

@@ -1177,6 +1177,30 @@ end subroutine d2dump
return
end subroutine polint
subroutine polint0(xa, ya, y, ordn)
! Lagrange interpolation at x=0, O(n) direct formula
implicit none
integer, intent(in) :: ordn
real*8, dimension(ordn), intent(in) :: xa, ya
real*8, intent(out) :: y
integer :: j, k
real*8 :: wj
y = 0.d0
do j = 1, ordn
wj = 1.d0
do k = 1, ordn
if (k .ne. j) then
wj = wj * xa(k) / (xa(k) - xa(j))
endif
enddo
y = y + wj * ya(j)
enddo
return
end subroutine polint0
!------------------------------------------------------------------------------
!
! interpolation in 2 dimensions, follow yx order

7
AMSS_NCKU_source/prolongrestrict_cell.f90
View File

@@ -217,7 +217,6 @@
real*8,dimension(2*ghost_width) :: X,Y,Z
real*8, dimension(2*ghost_width,2*ghost_width) :: tmp2
real*8, dimension(2*ghost_width) :: tmp1
real*8 :: ddy
real*8,dimension(3) :: ccp
#if (ghost_width == 2)
@@ -580,7 +579,7 @@
tmp1(ghost_width-cxI(1)+cxB(1) :ghost_width-cxI(1)+cxT(1) ) = funf(cxB(1):cxT(1),j,k)
endif
call polint(X,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(X,tmp1,funf(i,j,k),2*ghost_width)
! for y direction
elseif(sum(fg).eq.2.and.fg(2) .eq. 0.and. &
@@ -690,7 +689,7 @@
tmp1(ghost_width-cxI(2)+cxB(2) :ghost_width-cxI(2)+cxT(2) ) = funf(i,cxB(2):cxT(2),k)
endif
call polint(Y,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(Y,tmp1,funf(i,j,k),2*ghost_width)
! for z direction
elseif(sum(fg).eq.2.and.fg(3) .eq. 0.and. &
@@ -802,7 +801,7 @@
tmp1(ghost_width-cxI(3)+cxB(3) :ghost_width-cxI(3)+cxT(3) ) = funf(i,j,cxB(3):cxT(3))
endif
call polint(Z,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(Z,tmp1,funf(i,j,k),2*ghost_width)
#else

7
AMSS_NCKU_source/prolongrestrict_vertex.f90
View File

@@ -217,7 +217,6 @@
real*8,dimension(2*ghost_width) :: X,Y,Z
real*8, dimension(2*ghost_width,2*ghost_width) :: tmp2
real*8, dimension(2*ghost_width) :: tmp1
real*8 :: ddy
#if (ghost_width == 2)
real*8, parameter :: C1=-1.d0/16,C2=9.d0/16
@@ -470,7 +469,7 @@
tmp1(cxB(1)+ghost_width-i+1:cxT(1)+ghost_width-i+1) = fh(cxB(1):cxT(1),j,k)
call polint(X,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(X,tmp1,funf(i,j,k),2*ghost_width)
! for y direction
elseif (fg(2) .eq. 0)then
@@ -529,7 +528,7 @@
tmp1(cxB(2)+ghost_width-j+1:cxT(2)+ghost_width-j+1) = fh(i,cxB(2):cxT(2),k)
call polint(Y,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(Y,tmp1,funf(i,j,k),2*ghost_width)
! for z direction
else
@@ -588,7 +587,7 @@
tmp1(cxB(3)+ghost_width-k+1:cxT(3)+ghost_width-k+1) = fh(i,j,cxB(3):cxT(3))
call polint(Z,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
call polint0(Z,tmp1,funf(i,j,k),2*ghost_width)
endif