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

Optimize prolong3: replace parity branches with coefficient lookup

Browse Source
This commit is contained in:
CGH0S7
2026-02-28 23:59:57 +08:00
parent 2c0a3055d4
commit d4ec69c98a
1 changed files with 40 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

67
AMSS_NCKU_source/prolongrestrict_cell.f90
View File

@@ -1934,18 +1934,21 @@
! when if=1 -> ic=0, this is different to vertex center grid ! when if=1 -> ic=0, this is different to vertex center grid
real*8, dimension(-2:extc(1),-2:extc(2),-2:extc(3)) :: funcc real*8, dimension(-2:extc(1),-2:extc(2),-2:extc(3)) :: funcc
integer,dimension(3) :: cxI integer,dimension(3) :: cxI
integer :: i,j,k,ii,jj,kk integer :: i,j,k,ii,jj,kk,px,py,pz
real*8, dimension(6,6) :: tmp2 real*8, dimension(6,6) :: tmp2
real*8, dimension(6) :: tmp1 real*8, dimension(6) :: tmp1
integer, dimension(extf(1)) :: cix integer, dimension(extf(1)) :: cix
integer, dimension(extf(2)) :: ciy integer, dimension(extf(2)) :: ciy
integer, dimension(extf(3)) :: ciz integer, dimension(extf(3)) :: ciz
logical, dimension(extf(1)) :: evenx integer, dimension(extf(1)) :: pix
logical, dimension(extf(2)) :: eveny integer, dimension(extf(2)) :: piy
logical, dimension(extf(3)) :: evenz integer, dimension(extf(3)) :: piz
real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3 real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3
real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3 real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3
real*8, dimension(6,2), parameter :: WC = reshape((/&
C1,C2,C3,C4,C5,C6,&
C6,C5,C4,C3,C2,C1/), (/6,2/))
integer::imini,imaxi,jmini,jmaxi,kmini,kmaxi integer::imini,imaxi,jmini,jmaxi,kmini,kmaxi
integer::imino,imaxo,jmino,jmaxo,kmino,kmaxo integer::imino,imaxo,jmino,jmaxo,kmino,kmaxo
@@ -2029,17 +2032,29 @@
do i = imino,imaxo do i = imino,imaxo
ii = i + lbf(1) - 1 ii = i + lbf(1) - 1
cix(i) = ii/2 - lbc(1) + 1 cix(i) = ii/2 - lbc(1) + 1
evenx(i) = (ii/2*2 == ii) if(ii/2*2 == ii)then
pix(i) = 1
else
pix(i) = 2
endif
enddo enddo
do j = jmino,jmaxo do j = jmino,jmaxo
jj = j + lbf(2) - 1 jj = j + lbf(2) - 1
ciy(j) = jj/2 - lbc(2) + 1 ciy(j) = jj/2 - lbc(2) + 1
eveny(j) = (jj/2*2 == jj) if(jj/2*2 == jj)then
piy(j) = 1
else
piy(j) = 2
endif
enddo enddo
do k = kmino,kmaxo do k = kmino,kmaxo
kk = k + lbf(3) - 1 kk = k + lbf(3) - 1
ciz(k) = kk/2 - lbc(3) + 1 ciz(k) = kk/2 - lbc(3) + 1
evenz(k) = (kk/2*2 == kk) if(kk/2*2 == kk)then
piz(k) = 1
else
piz(k) = 2
endif
enddo enddo
call symmetry_bd(3,extc,func,funcc,SoA) call symmetry_bd(3,extc,func,funcc,SoA)
@@ -2051,6 +2066,9 @@
cxI(1) = cix(i) cxI(1) = cix(i)
cxI(2) = ciy(j) cxI(2) = ciy(j)
cxI(3) = ciz(k) cxI(3) = ciz(k)
px = pix(i)
py = piy(j)
pz = piz(k)
if(any(cxI+3 > extc)) write(*,*)"error in prolong" if(any(cxI+3 > extc)) write(*,*)"error in prolong"
#if 0 #if 0
@@ -2140,33 +2158,18 @@
endif endif
endif endif
#else #else
if(evenz(k))then tmp2= WC(1,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+& WC(2,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+& WC(3,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+&
C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+& WC(4,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+& WC(5,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+& WC(6,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
else
tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+&
C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
endif
if(eveny(j))then tmp1= WC(1,py)*tmp2(:,1)+WC(2,py)*tmp2(:,2)+WC(3,py)*tmp2(:,3)+&
tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6) WC(4,py)*tmp2(:,4)+WC(5,py)*tmp2(:,5)+WC(6,py)*tmp2(:,6)
else
tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
endif
if(evenx(i))then funf(i,j,k)= WC(1,px)*tmp1(1)+WC(2,px)*tmp1(2)+WC(3,px)*tmp1(3)+&
funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6) WC(4,px)*tmp1(4)+WC(5,px)*tmp1(5)+WC(6,px)*tmp1(6)
else
funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
endif
#endif #endif
enddo enddo
enddo enddo