diff --git a/AMSS_NCKU_source/bssn_rhs.f90 b/AMSS_NCKU_source/bssn_rhs.f90 index 246b219..9338f22 100644 --- a/AMSS_NCKU_source/bssn_rhs.f90 +++ b/AMSS_NCKU_source/bssn_rhs.f90 @@ -83,6 +83,10 @@ real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz +! shared work arrays (memory pool) to avoid repeated allocation in subroutines + real*8, dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh_work2 ! for fderivs_fh/fdderivs_fh (ghost_width=2) + real*8, dimension(-2:ex(1),-2:ex(2),-2:ex(3)) :: fh_work3 ! for kodis_fh/lopsided_fh (ghost_width=3) + real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA real*8 :: dX, dY, dZ, PI real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0 @@ -151,22 +155,22 @@ gyy = dyy + ONE gzz = dzz + ONE - call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev) - call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev) - call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev) + call fderivs_fh(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev,fh_work2) div_beta = betaxx + betayy + betazz - call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) + call fderivs_fh(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev,fh_work2) chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi - call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) - call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev) - call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev) - call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) - call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev) - call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) + call fderivs_fh(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2) gxx_rhs = - TWO * alpn1 * Axx - F2o3 * gxx * div_beta + & TWO *( gxx * betaxx + gxy * betayx + gxz * betazx) @@ -284,8 +288,8 @@ (gupyy * gupzz + gupyz * gupyz)* Ayz ! Right hand side for Gam^i without shift terms... - call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) - call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) + call fderivs_fh(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev,fh_work2) Gamx_rhs = - TWO * ( Lapx * Rxx + Lapy * Rxy + Lapz * Rxz ) + & TWO * alpn1 * ( & @@ -314,12 +318,12 @@ Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz + & TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) ) - call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,& - X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev) - call fdderivs(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,& - X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev) - call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,& - X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev) + call fdderivs_fh(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,& + X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev,fh_work2) + call fdderivs_fh(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,& + X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev,fh_work2) + call fdderivs_fh(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,& + X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev,fh_work2) fxx = gxxx + gxyy + gxzz fxy = gxyx + gyyy + gyzz @@ -332,9 +336,9 @@ Gamza = gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + & TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz ) - call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev) - call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev) - call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev) + call fderivs_fh(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev,fh_work2) + call fderivs_fh(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev,fh_work2) Gamx_rhs = Gamx_rhs + F2o3 * Gamxa * div_beta - & Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz + & @@ -377,27 +381,27 @@ gzzz = gxz * Gamxzz + gyz * Gamyzz + gzz * Gamzzz !compute Ricci tensor for tilted metric - call fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) + call fdderivs_fh(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2) Rxx = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - call fdderivs(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) + call fdderivs_fh(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2) Ryy = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - call fdderivs(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) + call fdderivs_fh(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2) Rzz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - call fdderivs(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,symmetry,Lev) + call fdderivs_fh(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,symmetry,Lev,fh_work2) Rxy = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - call fdderivs(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,symmetry,Lev) + call fdderivs_fh(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,symmetry,Lev,fh_work2) Rxz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - call fdderivs(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,symmetry,Lev) + call fdderivs_fh(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,symmetry,Lev,fh_work2) Ryz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO @@ -602,7 +606,7 @@ Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy + & Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz ) !covariant second derivative of chi respect to tilted metric - call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fdderivs_fh(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz @@ -628,8 +632,8 @@ Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO ! covariant second derivatives of the lapse respect to physical metric - call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & - SYM,SYM,SYM,symmetry,Lev) + call fdderivs_fh(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & + SYM,SYM,SYM,symmetry,Lev,fh_work2) gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1 gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1 @@ -812,7 +816,7 @@ betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz - call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 @@ -824,7 +828,7 @@ betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz - call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 @@ -832,7 +836,7 @@ dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfz_rhs = Gamz_rhs - reta*dtSfz #elif (GAUGE == 4) - call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 @@ -844,7 +848,7 @@ dtSfy_rhs = ZEO dtSfz_rhs = ZEO #elif (GAUGE == 5) - call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 @@ -947,51 +951,51 @@ !!!!!!!!!advection term part - call lopsided(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS) - call lopsided(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS) - call lopsided(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA) - call lopsided(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS) - call lopsided(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA) - call lopsided(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS) + call lopsided_fh(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) + call lopsided_fh(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,fh_work3) + call lopsided_fh(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,fh_work3) + call lopsided_fh(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) + call lopsided_fh(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,fh_work3) + call lopsided_fh(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) - call lopsided(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS) - call lopsided(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS) - call lopsided(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA) - call lopsided(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS) - call lopsided(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA) - call lopsided(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS) + call lopsided_fh(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) + call lopsided_fh(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,fh_work3) + call lopsided_fh(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,fh_work3) + call lopsided_fh(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) + call lopsided_fh(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,fh_work3) + call lopsided_fh(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) - call lopsided(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS) - call lopsided(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS) + call lopsided_fh(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) + call lopsided_fh(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) - call lopsided(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS) - call lopsided(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS) - call lopsided(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA) + call lopsided_fh(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3) + call lopsided_fh(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3) + call lopsided_fh(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3) !! - call lopsided(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS) + call lopsided_fh(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3) #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7) - call lopsided(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS) - call lopsided(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS) - call lopsided(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA) + call lopsided_fh(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3) + call lopsided_fh(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3) + call lopsided_fh(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3) #endif #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7) - call lopsided(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS) - call lopsided(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS) - call lopsided(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA) + call lopsided_fh(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3) + call lopsided_fh(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3) + call lopsided_fh(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3) #endif if(eps>0)then ! usual Kreiss-Oliger dissipation - call kodis(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps) - call kodis(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps) - call kodis(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps) - call kodis(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps) + call kodis_fh(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps,fh_work3) #if 0 #define i 42 #define j 40 @@ -1005,7 +1009,7 @@ endif #undef k !!stop #endif - call kodis(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps) + call kodis_fh(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps,fh_work3) #if 0 #define i 42 #define j 40 @@ -1019,25 +1023,25 @@ endif #undef k !!stop #endif - call kodis(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps) - call kodis(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps) - call kodis(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps) - call kodis(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps) - call kodis(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps) - call kodis(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps) + call kodis_fh(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps,fh_work3) #if 1 !! bam does not apply dissipation on gauge variables - call kodis(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps) - call kodis(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps) - call kodis(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps) - call kodis(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps) + call kodis_fh(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps,fh_work3) #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7) - call kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps) - call kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps) - call kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps) + call kodis_fh(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps,fh_work3) + call kodis_fh(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps,fh_work3) #endif #endif @@ -1078,12 +1082,12 @@ endif ! mov_Res_j = gupkj*(-1/chi d_k chi*A_ij + D_k A_ij) - 2/3 d_j trK - 8 PI s_j where D respect to physical metric ! store D_i A_jk - 1/chi d_i chi*A_jk in gjk_i - call fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) - call fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0) - call fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0) - call fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) - call fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0) - call fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) + call fderivs_fh(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2) + call fderivs_fh(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0,fh_work2) + call fderivs_fh(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0,fh_work2) + call fderivs_fh(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2) + call fderivs_fh(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0,fh_work2) + call fderivs_fh(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2) gxxx = gxxx - ( Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz & + Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz) - chix*Axx/chin1 diff --git a/AMSS_NCKU_source/diff_new.f90 b/AMSS_NCKU_source/diff_new.f90 index 93954f1..9238ab5 100644 --- a/AMSS_NCKU_source/diff_new.f90 +++ b/AMSS_NCKU_source/diff_new.f90 @@ -1103,6 +1103,103 @@ end subroutine fderivs !----------------------------------------------------------------------------- +! fderivs variant: reuses caller-provided fh work array (memory pool) +!----------------------------------------------------------------------------- + subroutine fderivs_fh(ex,f,fx,fy,fz,X,Y,Z,SYM1,SYM2,SYM3, & + symmetry,onoff,fh) + implicit none + + integer, intent(in ):: ex(1:3),symmetry,onoff + real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f + real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: fx,fy,fz + real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3)) + real*8, intent(in ):: SYM1,SYM2,SYM3 + real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh + + real*8 :: dX,dY,dZ + real*8, dimension(3) :: SoA + integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k + real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz + integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2 + real*8, parameter :: ZEO=0.d0,ONE=1.d0 + real*8, parameter :: TWO=2.d0,EIT=8.d0 + real*8, parameter :: F12=1.2d1 + + dX = X(2)-X(1) + dY = Y(2)-Y(1) + dZ = Z(2)-Z(1) + + imax = ex(1) + jmax = ex(2) + kmax = ex(3) + + imin = 1 + jmin = 1 + kmin = 1 + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1 + if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1 + if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1 + + SoA(1) = SYM1 + SoA(2) = SYM2 + SoA(3) = SYM3 + + call symmetry_bd(2,ex,f,fh,SoA) + + d12dx = ONE/F12/dX + d12dy = ONE/F12/dY + d12dz = ONE/F12/dZ + + d2dx = ONE/TWO/dX + d2dy = ONE/TWO/dY + d2dz = ONE/TWO/dZ + + fx = ZEO + fy = ZEO + fz = ZEO + + do k=1,ex(3)-1 + do j=1,ex(2)-1 + do i=1,ex(1)-1 +#if 0 + if(i+2 <= imax .and. i-2 >= imin)then + fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k)) + elseif(i+1 <= imax .and. i-1 >= imin)then + fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k)) + endif + if(j+2 <= jmax .and. j-2 >= jmin)then + fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k)) + elseif(j+1 <= jmax .and. j-1 >= jmin)then + fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k)) + endif + if(k+2 <= kmax .and. k-2 >= kmin)then + fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2)) + elseif(k+1 <= kmax .and. k-1 >= kmin)then + fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1)) + endif +#else + if(i+2 <= imax .and. i-2 >= imin .and. & + j+2 <= jmax .and. j-2 >= jmin .and. & + k+2 <= kmax .and. k-2 >= kmin) then + fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k)) + fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k)) + fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2)) + elseif(i+1 <= imax .and. i-1 >= imin .and. & + j+1 <= jmax .and. j-1 >= jmin .and. & + k+1 <= kmax .and. k-1 >= kmin) then + fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k)) + fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k)) + fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1)) + endif +#endif + enddo + enddo + enddo + + return + + end subroutine fderivs_fh +!----------------------------------------------------------------------------- ! ! single derivatives dx ! @@ -1940,6 +2037,162 @@ end subroutine fddyz +!----------------------------------------------------------------------------- +! fdderivs variant: reuses caller-provided fh work array (memory pool) +!----------------------------------------------------------------------------- + subroutine fdderivs_fh(ex,f,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & + SYM1,SYM2,SYM3,symmetry,onoff,fh) + implicit none + + integer, intent(in ):: ex(1:3),symmetry,onoff + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ):: f + real*8, dimension(ex(1),ex(2),ex(3)),intent(out):: fxx,fxy,fxz,fyy,fyz,fzz + real*8, intent(in ):: X(ex(1)),Y(ex(2)),Z(ex(3)) + real*8, intent(in ):: SYM1,SYM2,SYM3 + real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh + + real*8 :: dX,dY,dZ + real*8, dimension(3) :: SoA + integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k + real*8 :: Sdxdx,Sdydy,Sdzdz,Fdxdx,Fdydy,Fdzdz + real*8 :: Sdxdy,Sdxdz,Sdydz,Fdxdy,Fdxdz,Fdydz + integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2 + real*8, parameter :: ZEO=0.d0, ONE=1.d0, TWO=2.d0, F1o4=2.5d-1 + real*8, parameter :: F8=8.d0, F16=1.6d1, F30=3.d1 + real*8, parameter :: F1o12=ONE/1.2d1, F1o144=ONE/1.44d2 + + dX = X(2)-X(1) + dY = Y(2)-Y(1) + dZ = Z(2)-Z(1) + + imax = ex(1) + jmax = ex(2) + kmax = ex(3) + + imin = 1 + jmin = 1 + kmin = 1 + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1 + if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1 + if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1 + + SoA(1) = SYM1 + SoA(2) = SYM2 + SoA(3) = SYM3 + + call symmetry_bd(2,ex,f,fh,SoA) + + Sdxdx = ONE /( dX * dX ) + Sdydy = ONE /( dY * dY ) + Sdzdz = ONE /( dZ * dZ ) + + Fdxdx = F1o12 /( dX * dX ) + Fdydy = F1o12 /( dY * dY ) + Fdzdz = F1o12 /( dZ * dZ ) + + Sdxdy = F1o4 /( dX * dY ) + Sdxdz = F1o4 /( dX * dZ ) + Sdydz = F1o4 /( dY * dZ ) + + Fdxdy = F1o144 /( dX * dY ) + Fdxdz = F1o144 /( dX * dZ ) + Fdydz = F1o144 /( dY * dZ ) + + fxx = ZEO + fyy = ZEO + fzz = ZEO + fxy = ZEO + fxz = ZEO + fyz = ZEO + + do k=1,ex(3)-1 + do j=1,ex(2)-1 + do i=1,ex(1)-1 +#if 0 + if(i+2 <= imax .and. i-2 >= imin)then + fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) & + -fh(i+2,j,k)+F16*fh(i+1,j,k) ) + elseif(i+1 <= imax .and. i-1 >= imin)then + fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k)) + endif + if(j+2 <= jmax .and. j-2 >= jmin)then + fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) & + -fh(i,j+2,k)+F16*fh(i,j+1,k) ) + elseif(j+1 <= jmax .and. j-1 >= jmin)then + fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k)) + endif + if(k+2 <= kmax .and. k-2 >= kmin)then + fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) & + -fh(i,j,k+2)+F16*fh(i,j,k+1) ) + elseif(k+1 <= kmax .and. k-1 >= kmin)then + fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1)) + endif + if(i+2 <= imax .and. i-2 >= imin .and. j+2 <= jmax .and. j-2 >= jmin)then + fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) & + -F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) & + +F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) & + - (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k))) + elseif(i+1 <= imax .and. i-1 >= imin .and. j+1 <= jmax .and. j-1 >= jmin)then + fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k)) + endif + if(i+2 <= imax .and. i-2 >= imin .and. k+2 <= kmax .and. k-2 >= kmin)then + fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) & + -F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) & + +F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) & + - (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2))) + elseif(i+1 <= imax .and. i-1 >= imin .and. k+1 <= kmax .and. k-1 >= kmin)then + fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1)) + endif + if(j+2 <= jmax .and. j-2 >= jmin .and. k+2 <= kmax .and. k-2 >= kmin)then + fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) & + -F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) & + +F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) & + - (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2))) + elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then + fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1)) + endif +#else +! for bam comparison + if(i+2 <= imax .and. i-2 >= imin .and. & + j+2 <= jmax .and. j-2 >= jmin .and. & + k+2 <= kmax .and. k-2 >= kmin) then + fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) & + -fh(i+2,j,k)+F16*fh(i+1,j,k) ) + fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) & + -fh(i,j+2,k)+F16*fh(i,j+1,k) ) + fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) & + -fh(i,j,k+2)+F16*fh(i,j,k+1) ) + fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) & + -F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) & + +F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) & + - (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k))) + fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) & + -F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) & + +F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) & + - (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2))) + fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) & + -F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) & + +F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) & + - (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2))) + elseif(i+1 <= imax .and. i-1 >= imin .and. & + j+1 <= jmax .and. j-1 >= jmin .and. & + k+1 <= kmax .and. k-1 >= kmin) then + fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k)) + fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k)) + fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1)) + fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k)) + fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1)) + fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1)) + endif +#endif + enddo + enddo + enddo + + return + + end subroutine fdderivs_fh + #elif (ghost_width == 4) ! sixth order code diff --git a/AMSS_NCKU_source/kodiss.f90 b/AMSS_NCKU_source/kodiss.f90 index a12ada4..5b6114c 100644 --- a/AMSS_NCKU_source/kodiss.f90 +++ b/AMSS_NCKU_source/kodiss.f90 @@ -215,6 +215,99 @@ integer, parameter :: NO_SYMM=0, OCTANT=2 end subroutine kodis +!----------------------------------------------------------------------------- +! kodis variant: reuses caller-provided fh work array (memory pool) +!----------------------------------------------------------------------------- +subroutine kodis_fh(ex,X,Y,Z,f,f_rhs,SoA,Symmetry,eps,fh) + +implicit none +! argument variables +integer,intent(in) :: Symmetry +integer,dimension(3),intent(in)::ex +real*8, dimension(1:3), intent(in) :: SoA +double precision,intent(in),dimension(ex(1))::X +double precision,intent(in),dimension(ex(2))::Y +double precision,intent(in),dimension(ex(3))::Z +double precision,intent(in),dimension(ex(1),ex(2),ex(3))::f +double precision,intent(inout),dimension(ex(1),ex(2),ex(3))::f_rhs +real*8,intent(in) :: eps +real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh +! local variables +integer :: imin,jmin,kmin,imax,jmax,kmax +integer :: i,j,k +real*8 :: dX,dY,dZ +real*8, parameter :: ONE=1.d0,SIX=6.d0,FIT=1.5d1,TWT=2.d1 +real*8,parameter::cof=6.4d1 ! 2^6 +integer, parameter :: NO_SYMM=0, OCTANT=2 + + dX = X(2)-X(1) + dY = Y(2)-Y(1) + dZ = Z(2)-Z(1) + + imax = ex(1) + jmax = ex(2) + kmax = ex(3) + + imin = 1 + jmin = 1 + kmin = 1 + + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2 + if(Symmetry == OCTANT .and. dabs(X(1)) < dX) imin = -2 + if(Symmetry == OCTANT .and. dabs(Y(1)) < dY) jmin = -2 + + call symmetry_bd(3,ex,f,fh,SoA) + + do k=1,ex(3) + do j=1,ex(2) + do i=1,ex(1) + + if(i-3 >= imin .and. i+3 <= imax .and. & + j-3 >= jmin .and. j+3 <= jmax .and. & + k-3 >= kmin .and. k+3 <= kmax) then +#if 0 + f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dX/cof * ( & + (fh(i-3,j,k)+fh(i+3,j,k)) - & + SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + & + FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - & + TWT* fh(i,j,k) ) + f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dY/cof * ( & + (fh(i,j-3,k)+fh(i,j+3,k)) - & + SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + & + FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - & + TWT* fh(i,j,k) ) + f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dZ/cof * ( & + (fh(i,j,k-3)+fh(i,j,k+3)) - & + SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + & + FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - & + TWT* fh(i,j,k) ) +#else + f_rhs(i,j,k) = f_rhs(i,j,k) + eps/cof *( ( & + (fh(i-3,j,k)+fh(i+3,j,k)) - & + SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + & + FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - & + TWT* fh(i,j,k) )/dX + & + ( & + (fh(i,j-3,k)+fh(i,j+3,k)) - & + SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + & + FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - & + TWT* fh(i,j,k) )/dY + & + ( & + (fh(i,j,k-3)+fh(i,j,k+3)) - & + SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + & + FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - & + TWT* fh(i,j,k) )/dZ ) +#endif + endif + + enddo + enddo + enddo + + return + + end subroutine kodis_fh + #elif (ghost_width == 4) ! sixth order code !------------------------------------------------------------------------------------------------------------------------------ diff --git a/AMSS_NCKU_source/lopsidediff.f90 b/AMSS_NCKU_source/lopsidediff.f90 index 2e97af5..01fc884 100644 --- a/AMSS_NCKU_source/lopsidediff.f90 +++ b/AMSS_NCKU_source/lopsidediff.f90 @@ -487,6 +487,160 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA) end subroutine lopsided +!----------------------------------------------------------------------------- +! lopsided variant: reuses caller-provided fh work array (memory pool) +!----------------------------------------------------------------------------- +subroutine lopsided_fh(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA,fh) + implicit none + +!~~~~~~> Input parameters: + + integer, intent(in) :: ex(1:3),Symmetry + real*8, intent(in) :: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)) + real*8,dimension(ex(1),ex(2),ex(3)),intent(in) :: f,Sfx,Sfy,Sfz + + real*8,dimension(ex(1),ex(2),ex(3)),intent(inout):: f_rhs + real*8,dimension(3),intent(in) ::SoA + real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh + +!~~~~~~> local variables: + integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k + real*8 :: dX,dY,dZ + real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz + real*8, parameter :: ZEO=0.d0,ONE=1.d0, F3=3.d0 + real*8, parameter :: TWO=2.d0,F6=6.0d0,F18=1.8d1 + real*8, parameter :: F12=1.2d1, F10=1.d1,EIT=8.d0 + integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2 + + dX = X(2)-X(1) + dY = Y(2)-Y(1) + dZ = Z(2)-Z(1) + + d12dx = ONE/F12/dX + d12dy = ONE/F12/dY + d12dz = ONE/F12/dZ + + d2dx = ONE/TWO/dX + d2dy = ONE/TWO/dY + d2dz = ONE/TWO/dZ + + imax = ex(1) + jmax = ex(2) + kmax = ex(3) + + imin = 1 + jmin = 1 + kmin = 1 + if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2 + if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -2 + if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -2 + + call symmetry_bd(3,ex,f,fh,SoA) + +! upper bound set ex-1 only for efficiency, +! the loop body will set ex 0 also + do k=1,ex(3)-1 + do j=1,ex(2)-1 + do i=1,ex(1)-1 +#if 0 +!! old code - same as original lopsided +#else +!! new code, 2012dec27, based on bam +! x direction + if(Sfx(i,j,k) > ZEO)then + if(i+3 <= imax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) & + -F6*fh(i+2,j,k)+ fh(i+3,j,k)) + elseif(i+2 <= imax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k)) + elseif(i+1 <= imax)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) & + -F6*fh(i-2,j,k)+ fh(i-3,j,k)) + endif + elseif(Sfx(i,j,k) < ZEO)then + if(i-3 >= imin)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) & + -F6*fh(i-2,j,k)+ fh(i-3,j,k)) + elseif(i-2 >= imin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k)) + elseif(i-1 >= imin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) & + -F6*fh(i+2,j,k)+ fh(i+3,j,k)) + endif + endif + +! y direction + if(Sfy(i,j,k) > ZEO)then + if(j+3 <= jmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) & + -F6*fh(i,j+2,k)+ fh(i,j+3,k)) + elseif(j+2 <= jmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfy(i,j,k)*d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k)) + elseif(j+1 <= jmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) & + -F6*fh(i,j-2,k)+ fh(i,j-3,k)) + endif + elseif(Sfy(i,j,k) < ZEO)then + if(j-3 >= jmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) & + -F6*fh(i,j-2,k)+ fh(i,j-3,k)) + elseif(j-2 >= jmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfy(i,j,k)*d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k)) + elseif(j-1 >= jmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) & + -F6*fh(i,j+2,k)+ fh(i,j+3,k)) + endif + endif + +! z direction + if(Sfz(i,j,k) > ZEO)then + if(k+3 <= kmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) & + -F6*fh(i,j,k+2)+ fh(i,j,k+3)) + elseif(k+2 <= kmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfz(i,j,k)*d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2)) + elseif(k+1 <= kmax)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) & + -F6*fh(i,j,k-2)+ fh(i,j,k-3)) + endif + elseif(Sfz(i,j,k) < ZEO)then + if(k-3 >= kmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)- & + Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) & + -F6*fh(i,j,k-2)+ fh(i,j,k-3)) + elseif(k-2 >= kmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfz(i,j,k)*d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2)) + elseif(k-1 >= kmin)then + f_rhs(i,j,k)=f_rhs(i,j,k)+ & + Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) & + -F6*fh(i,j,k+2)+ fh(i,j,k+3)) + endif + endif +#endif + enddo + enddo + enddo + + return + + end subroutine lopsided_fh + #elif (ghost_width == 4) ! sixth order code ! Compute advection terms in right hand sides of field equations