From 2c0a3055d4ef400ed5975cfc46e915a93c83900a Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Sat, 28 Feb 2026 23:53:30 +0800 Subject: [PATCH] Optimize prolong3: precompute coarse index/parity maps --- AMSS_NCKU_source/prolongrestrict_cell.f90 | 123 ++++++++++++---------- 1 file changed, 68 insertions(+), 55 deletions(-) diff --git a/AMSS_NCKU_source/prolongrestrict_cell.f90 b/AMSS_NCKU_source/prolongrestrict_cell.f90 index 17fc43d..8d2907f 100644 --- a/AMSS_NCKU_source/prolongrestrict_cell.f90 +++ b/AMSS_NCKU_source/prolongrestrict_cell.f90 @@ -1933,10 +1933,16 @@ integer,dimension(3) :: lbc,ubc,lbf,ubf,lbp,ubp,lbpc,ubpc ! 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 - integer,dimension(3) :: cxI - integer :: i,j,k,ii,jj,kk - real*8, dimension(6,6) :: tmp2 - real*8, dimension(6) :: tmp1 + integer,dimension(3) :: cxI + integer :: i,j,k,ii,jj,kk + real*8, dimension(6,6) :: tmp2 + real*8, dimension(6) :: tmp1 + integer, dimension(extf(1)) :: cix + integer, dimension(extf(2)) :: ciy + integer, dimension(extf(3)) :: ciz + logical, dimension(extf(1)) :: evenx + logical, dimension(extf(2)) :: eveny + logical, dimension(extf(3)) :: evenz 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 @@ -2003,10 +2009,10 @@ kmini=lbpc(3)-lbc(3) + 1 kmaxi=ubpc(3)-lbc(3) + 1 - if(imino.lt.1.or.jmino.lt.1.or.kmino.lt.1.or.& - imini.lt.1.or.jmini.lt.1.or.kmini.lt.1.or.& - imaxo.gt.extf(1).or.jmaxo.gt.extf(2).or.kmaxo.gt.extf(3).or.& - imaxi.gt.extc(1)-2.or.jmaxi.gt.extc(2)-2.or.kmaxi.gt.extc(3)-2)then + if(imino.lt.1.or.jmino.lt.1.or.kmino.lt.1.or.& + imini.lt.1.or.jmini.lt.1.or.kmini.lt.1.or.& + imaxo.gt.extf(1).or.jmaxo.gt.extf(2).or.kmaxo.gt.extf(3).or.& + imaxi.gt.extc(1)-2.or.jmaxi.gt.extc(2)-2.or.kmaxi.gt.extc(3)-2)then write(*,*)"error in prolongation for" write(*,*)"from" write(*,*)llbc,uubc @@ -2017,33 +2023,40 @@ write(*,*)"want" write(*,*)llbp,uubp write(*,*)lbp,ubp,lbpc,ubpc - return - endif - - call symmetry_bd(3,extc,func,funcc,SoA) - -!~~~~~~> prolongation start... - do k = kmino,kmaxo - do j = jmino,jmaxo - do i = imino,imaxo - cxI(1) = i - cxI(2) = j - cxI(3) = k -! change to coarse level reference -!|---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*---| -!|=======x===============x===============x===============x=======| - cxI = (cxI+lbf-1)/2 -! change to array index - cxI = cxI - lbc + 1 - - if(any(cxI+3 > extc)) write(*,*)"error in prolong" - ii=i+lbf(1)-1 - jj=j+lbf(2)-1 - kk=k+lbf(3)-1 -#if 0 - if(ii/2*2==ii)then - if(jj/2*2==jj)then - if(kk/2*2==kk)then + return + endif + + do i = imino,imaxo + ii = i + lbf(1) - 1 + cix(i) = ii/2 - lbc(1) + 1 + evenx(i) = (ii/2*2 == ii) + enddo + do j = jmino,jmaxo + jj = j + lbf(2) - 1 + ciy(j) = jj/2 - lbc(2) + 1 + eveny(j) = (jj/2*2 == jj) + enddo + do k = kmino,kmaxo + kk = k + lbf(3) - 1 + ciz(k) = kk/2 - lbc(3) + 1 + evenz(k) = (kk/2*2 == kk) + enddo + + call symmetry_bd(3,extc,func,funcc,SoA) + +!~~~~~~> prolongation start... + do k = kmino,kmaxo + do j = jmino,jmaxo + do i = imino,imaxo + cxI(1) = cix(i) + cxI(2) = ciy(j) + cxI(3) = ciz(k) + + if(any(cxI+3 > extc)) write(*,*)"error in prolong" +#if 0 + if(ii/2*2==ii)then + if(jj/2*2==jj)then + if(kk/2*2==kk)then tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+& C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+& C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+& @@ -2126,11 +2139,11 @@ endif endif endif -#else - if(kk/2*2==kk)then - tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+& - C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+& - C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+& +#else + if(evenz(k))then + tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+& + C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+& + C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3) )+& C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+& C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+& C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3) @@ -2139,21 +2152,21 @@ 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(jj/2*2==jj)then - tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6) - else - tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6) - endif - - if(ii/2*2==ii)then - funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*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 + 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= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*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)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*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 enddo enddo