diff --git a/AMSS_NCKU_source/prolongrestrict_cell.f90 b/AMSS_NCKU_source/prolongrestrict_cell.f90 index 4f6676b..9fefe48 100644 --- a/AMSS_NCKU_source/prolongrestrict_cell.f90 +++ b/AMSS_NCKU_source/prolongrestrict_cell.f90 @@ -1958,7 +1958,9 @@ real*8 :: tmp_yz(extc(1), 6) ! 存储整条 X 线上 6 个 Y 轴偏置的 Z 向插值结果 real*8 :: tmp_xyz_line(-2:extc(1)) ! 包含 X 向 6 点模板访问所需下界 real*8 :: v1, v2, v3, v4, v5, v6 - integer :: ic, jc, kc, ix_offset,ix,iy,iz,jc_min,jc_max,ic_min,ic_max + integer :: ic, jc, kc, ix_offset,ix,iy,iz,jc_min,jc_max,ic_min,ic_max,kc_min,kc_max + integer :: i_lo, i_hi, j_lo, j_hi, k_lo, k_hi + logical :: need_full_symmetry real*8 :: res_line real*8 :: tmp_z_slab(-2:extc(1), -2:extc(2)) ! 包含 Y/X 向模板访问所需下界 if(wei.ne.3)then @@ -2063,20 +2065,35 @@ endif enddo - maxcx = maxval(cix(imino:imaxo)) - maxcy = maxval(ciy(jmino:jmaxo)) - maxcz = maxval(ciz(kmino:kmaxo)) + ic_min = minval(cix(imino:imaxo)) + ic_max = maxval(cix(imino:imaxo)) + jc_min = minval(ciy(jmino:jmaxo)) + jc_max = maxval(ciy(jmino:jmaxo)) + kc_min = minval(ciz(kmino:kmaxo)) + kc_max = maxval(ciz(kmino:kmaxo)) + + maxcx = ic_max + maxcy = jc_max + maxcz = kc_max if(maxcx+3 > extc(1) .or. maxcy+3 > extc(2) .or. maxcz+3 > extc(3))then write(*,*)"error in prolong" return endif - call symmetry_bd(3,extc,func,funcc,SoA) + i_lo = ic_min - 2 + i_hi = ic_max + 3 + j_lo = jc_min - 2 + j_hi = jc_max + 3 + k_lo = kc_min - 2 + k_hi = kc_max + 3 + need_full_symmetry = (i_lo < 1) .or. (j_lo < 1) .or. (k_lo < 1) + if(need_full_symmetry)then + call symmetry_bd(3,extc,func,funcc,SoA) + else + funcc(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi) = func(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi) + endif + ! 对每个 k(pz, kc 固定)预计算 Z 向插值的 2D 切片 -jc_min = minval(ciy(jmino:jmaxo)) -jc_max = maxval(ciy(jmino:jmaxo)) -ic_min = minval(cix(imino:imaxo)) -ic_max = maxval(cix(imino:imaxo)) do k = kmino, kmaxo pz = piz(k); kc = ciz(k) @@ -2357,6 +2374,8 @@ end do real*8 :: tmp_x_line(-1:extf(1)) integer :: fi, fj, fk, ii, jj, kk integer :: fi_min, fi_max, ii_lo, ii_hi + integer :: fj_min, fj_max, fk_min, fk_max, jj_lo, jj_hi, kk_lo, kk_hi + logical :: need_full_symmetry if(wei.ne.3)then write(*,*)"prolongrestrict.f90::restrict3: this routine only surport 3 dimension" @@ -2436,19 +2455,34 @@ end do stop endif - call symmetry_bd(2,extf,funf,funff,SoA) - ! 仅计算 X 向最终写回所需的窗口: ! func(i,j,k) 只访问 tmp_x_line(fi-2:fi+3) fi_min = 2*(imino + lbc(1) - 1) - 1 - lbf(1) + 1 fi_max = 2*(imaxo + lbc(1) - 1) - 1 - lbf(1) + 1 + fj_min = 2*(jmino + lbc(2) - 1) - 1 - lbf(2) + 1 + fj_max = 2*(jmaxo + lbc(2) - 1) - 1 - lbf(2) + 1 + fk_min = 2*(kmino + lbc(3) - 1) - 1 - lbf(3) + 1 + fk_max = 2*(kmaxo + lbc(3) - 1) - 1 - lbf(3) + 1 ii_lo = fi_min - 2 ii_hi = fi_max + 3 - if(ii_lo < -1 .or. ii_hi > extf(1))then - write(*,*)"restrict3: invalid ii window",ii_lo,ii_hi - write(*,*)"imino,imaxo,lbc(1),lbf(1),extf(1) = ",imino,imaxo,lbc(1),lbf(1),extf(1) + jj_lo = fj_min - 2 + jj_hi = fj_max + 3 + kk_lo = fk_min - 2 + kk_hi = fk_max + 3 + if(ii_lo < -1 .or. ii_hi > extf(1) .or. & + jj_lo < -1 .or. jj_hi > extf(2) .or. & + kk_lo < -1 .or. kk_hi > extf(3))then + write(*,*)"restrict3: invalid stencil window" + write(*,*)"ii=",ii_lo,ii_hi," jj=",jj_lo,jj_hi," kk=",kk_lo,kk_hi + write(*,*)"extf=",extf stop endif + need_full_symmetry = (ii_lo < 1) .or. (jj_lo < 1) .or. (kk_lo < 1) + if(need_full_symmetry)then + call symmetry_bd(2,extf,funf,funff,SoA) + else + funff(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi) = funf(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi) + endif !~~~~~~> restriction start... do k = kmino, kmaxo