From 039d263c53816849c0c2db2fa6fccdb7e028b936 Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Thu, 26 Feb 2026 19:37:38 -0500 Subject: [PATCH 01/12] fix TS adjoint --- src/adjoint/outputForward/residuals_d.f90 | 167 ++++++++- src/adjoint/outputReverse/residuals_b.f90 | 172 ++++++++- .../outputReverseFast/residuals_fast_b.f90 | 154 +++++++- src/solver/residuals.F90 | 56 ++- tests/test_ts_adjoint.py | 347 ++++++++++++++++++ 5 files changed, 872 insertions(+), 24 deletions(-) create mode 100644 tests/test_ts_adjoint.py diff --git a/src/adjoint/outputForward/residuals_d.f90 b/src/adjoint/outputForward/residuals_d.f90 index 8e12ed2e6..0590c069f 100644 --- a/src/adjoint/outputForward/residuals_d.f90 +++ b/src/adjoint/outputForward/residuals_d.f90 @@ -508,8 +508,8 @@ end subroutine sourceterms_block ! differentiation of initres_block in forward (tangent) mode (with options i4 dr8 r8): ! variations of useful results: *dw -! with respect to varying inputs: *dw -! rw status of diff variables: *dw:in-out +! with respect to varying inputs: *dw *flowDoms.w *flowDoms.vol +! rw status of diff variables: *dw:in-out *flowDoms.w:in *flowDoms.vol:in ! plus diff mem management of: dw:in subroutine initres_block_d(varstart, varend, nn, sps) ! @@ -534,7 +534,7 @@ subroutine initres_block_d(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: oneoverdt, tmp, tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp ! return immediately of no variables are in the range. @@ -542,8 +542,8 @@ subroutine initres_block_d(varstart, varend, nn, sps) return else ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then @@ -573,6 +573,163 @@ subroutine initres_block_d(varstart, varend, nn, sps) end do end do end if + case (timespectral) +! time spectral computation. + jj = sectionid + if (currentlevel .eq. groundlevel) then +! fine grid: initialize to zero, then accumulate spectral time derivative + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dwd(i, j, k, l) = 0.0_8 + dw(i, j, k, l) = zero + end do + end do + end do + end do +! loop over spectral instances + do mm = 1, ntimeintervalsspectral + ii = 3 * (mm - 1) + do l = varstart, varend + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = 2, kl + do j = 2, jl + do i = 2, il +! tangent of: tmp = dvector * velocities + tmpd = dvector(jj, ll, ii+1) * flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdomsd(nn, 1, mm)%w(i,j,k,ivz) + tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) +! tangent of: dw += tmp * vol_mm * rho_mm + dwd(i,j,k,l) = dwd(i,j,k,l) & + + tmpd * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + + tmp * flowdomsd(nn, 1, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdomsd(nn, 1, mm)%w(i,j,k,irho) + dw(i,j,k,l) = dw(i,j,k,l) & + + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) + end do + end do + end do + else +! scalar variable + do k = 2, kl + do j = 2, jl + do i = 2, il +! tangent of: dw += dscalar * vol_mm * w_mm + dwd(i,j,k,l) = dwd(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdomsd(nn, 1, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdomsd(nn, 1, mm)%w(i,j,k,l) + dw(i,j,k,l) = dw(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) + end do + end do + end do + end if + end do + end do + else +! coarse grid level + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dwd(i, j, k, l) = 0.0_8 + dw(i, j, k, l) = wr(i, j, k, l) + end do + end do + end do + end do +! loop over spectral instances (coarse grid correction) + do mm = 1, ntimeintervalsspectral + ii = 3 * (mm - 1) + do l = varstart, varend + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = 2, kl + do j = 2, jl + do i = 2, il +! tangent of: tmp = dvector * (rho*vel - rho1*vel1) +! w1 is frozen so w1d = 0 + tmpd = dvector(jj, ll, ii+1) & + * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdomsd(nn, 1, mm)%w(i,j,k,ivx)) & + + dvector(jj, ll, ii+2) & + * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdomsd(nn, 1, mm)%w(i,j,k,ivy)) & + + dvector(jj, ll, ii+3) & + * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & + + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdomsd(nn, 1, mm)%w(i,j,k,ivz)) + tmp = dvector(jj, ll, ii+1) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & + + dvector(jj, ll, ii+2) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & + + dvector(jj, ll, ii+3) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) +! tangent of: dw += tmp * vol_mm + dwd(i,j,k,l) = dwd(i,j,k,l) & + + tmpd * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + + tmp * flowdomsd(nn, 1, mm)%vol(i,j,k) + dw(i,j,k,l) = dw(i,j,k,l) & + + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) + end do + end do + end do + else +! scalar variable + do k = 2, kl + do j = 2, jl + do i = 2, il +! tangent of: dw += dscalar * vol_mm * (w_mm - w1_mm) +! w1 is frozen so w1d = 0 + dwd(i,j,k,l) = dwd(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdomsd(nn, 1, mm)%vol(i,j,k) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdomsd(nn, 1, mm)%w(i,j,k,l) + dw(i,j,k,l) = dw(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) + end do + end do + end do + end if + end do + end do + end if end select ! set the residual in the halo cells to zero. this is just ! to avoid possible problems. their values do not matter. diff --git a/src/adjoint/outputReverse/residuals_b.f90 b/src/adjoint/outputReverse/residuals_b.f90 index b16d723d8..c513ba104 100644 --- a/src/adjoint/outputReverse/residuals_b.f90 +++ b/src/adjoint/outputReverse/residuals_b.f90 @@ -531,8 +531,8 @@ end subroutine sourceterms_block ! differentiation of initres_block in reverse (adjoint) mode (with options noisize i4 dr8 r8): ! gradient of useful results: *dw -! with respect to varying inputs: *dw -! rw status of diff variables: *dw:in-out +! with respect to varying inputs: *dw *flowDoms.w *flowDoms.vol +! rw status of diff variables: *dw:in-out *flowDoms.w:out *flowDoms.vol:out ! plus diff mem management of: dw:in subroutine initres_block_b(varstart, varend, nn, sps) ! @@ -557,25 +557,33 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: oneoverdt, tmp, tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp integer :: branch ! return immediately of no variables are in the range. if (varend .ge. varstart) then ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol2b(1) + call pushcontrol3b(1) + else + call pushcontrol3b(0) + end if + case (timespectral) +! time spectral computation. + if (currentlevel .eq. groundlevel) then + call pushcontrol3b(2) else - call pushcontrol2b(0) + call pushcontrol3b(3) end if case default - call pushcontrol2b(2) + call pushcontrol3b(4) end select +! reverse the halo cell zeroing do l=varend,varstart,-1 do j=jl,2,-1 do i=il,2,-1 @@ -602,8 +610,9 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do end do - call popcontrol2b(branch) + call popcontrol3b(branch) if (branch .eq. 0) then +! steady, coarse grid: forward did dw = wr (constant) do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -614,6 +623,151 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do else if (branch .eq. 1) then +! steady, fine grid: forward did dw = zero + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do + end do + end do + end do + else if (branch .eq. 2) then +! time spectral, fine grid +! reverse the spectral time derivative accumulation + jj = sectionid + do mm = ntimeintervalsspectral, 1, -1 + ii = 3 * (mm - 1) + do l = varend, varstart, -1 + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 +! recompute tmp from forward pass + tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) +! adjoint of: dw(i,j,k,l) += tmp * vol_mm * rho_mm + tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & + + tmp * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) +! adjoint of: tmp = dvector * velocities + flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+1) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+2) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & + + dvector(jj, ll, ii+3) * tmpd + end do + end do + end do + else +! scalar variable +! adjoint of: dw(i,j,k,l) += dscalar(jj,sps,mm) * vol_mm * w_mm(l) + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 + flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & + * dwd(i,j,k,l) + end do + end do + end do + end if + end do + end do +! reverse of initialization: dw = zero => dwd = 0 + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do + end do + end do + end do + else if (branch .eq. 3) then +! time spectral, coarse grid +! reverse the coarse grid spectral time derivative accumulation + jj = sectionid + do mm = ntimeintervalsspectral, 1, -1 + ii = 3 * (mm - 1) + do l = varend, varstart, -1 + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 +! recompute tmp from forward pass + tmp = dvector(jj, ll, ii+1) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & + + dvector(jj, ll, ii+2) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & + + dvector(jj, ll, ii+3) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) +! adjoint of: dw += tmp * vol_mm + tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & + + tmp * dwd(i,j,k,l) +! adjoint of tmp = dvector * (rho*vel - rho1*vel1) +! only w is active, w1 is frozen + flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + + (dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz)) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + end do + end do + end do + else +! scalar variable +! adjoint of: dw += dscalar * vol_mm * (w_mm - w1_mm) +! only w is active, w1 is frozen + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 + flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & + + dscalar(jj, sps, mm) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) * dwd(i,j,k,l) + end do + end do + end do + end if + end do + end do +! reverse of initialization: dw = wr (constant) => dwd = 0 do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index edf78f8d1..ce273184c 100644 --- a/src/adjoint/outputReverseFast/residuals_fast_b.f90 +++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90 @@ -480,8 +480,8 @@ end subroutine sourceterms_block ! differentiation of initres_block in reverse (adjoint) mode (with options noisize i4 dr8 r8): ! gradient of useful results: *dw -! with respect to varying inputs: *dw -! rw status of diff variables: *dw:in-out +! with respect to varying inputs: *dw *flowDoms.w +! rw status of diff variables: *dw:in-out *flowDoms.w:out ! plus diff mem management of: dw:in subroutine initres_block_fast_b(varstart, varend, nn, sps) ! @@ -506,24 +506,31 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: oneoverdt, tmp, tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp integer :: branch ! return immediately of no variables are in the range. if (varend .ge. varstart) then ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol2b(1) + call pushcontrol3b(1) + else + call pushcontrol3b(0) + end if + case (timespectral) +! time spectral computation. + if (currentlevel .eq. groundlevel) then + call pushcontrol3b(2) else - call pushcontrol2b(0) + call pushcontrol3b(3) end if case default - call pushcontrol2b(2) + call pushcontrol3b(4) end select do l=varend,varstart,-1 do j=jl,2,-1 @@ -551,8 +558,9 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do end do - call popcontrol2b(branch) + call popcontrol3b(branch) if (branch .eq. 0) then +! steady, coarse grid do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -563,6 +571,134 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do else if (branch .eq. 1) then +! steady, fine grid + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do + end do + end do + end do + else if (branch .eq. 2) then +! time spectral, fine grid (state-only: vol, dscalar, dvector frozen) + jj = sectionid + do mm = ntimeintervalsspectral, 1, -1 + ii = 3 * (mm - 1) + do l = varend, varstart, -1 + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 +! recompute tmp + tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) +! adjoint of: dw += tmp * vol_mm * rho_mm (vol frozen) + tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) + flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) +! adjoint of: tmp = dvector * velocities + flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+1) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+2) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & + + dvector(jj, ll, ii+3) * tmpd + end do + end do + end do + else +! scalar variable (vol and dscalar frozen) + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 + flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * dwd(i,j,k,l) + end do + end do + end do + end if + end do + end do +! reverse of initialization: dw = zero => dwd = 0 + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do + end do + end do + end do + else if (branch .eq. 3) then +! time spectral, coarse grid (state-only: vol, dscalar, dvector, w1 frozen) + jj = sectionid + do mm = ntimeintervalsspectral, 1, -1 + ii = 3 * (mm - 1) + do l = varend, varstart, -1 + if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then +! momentum variable + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 +! recompute tmp + tmp = dvector(jj, ll, ii+1) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & + + dvector(jj, ll, ii+2) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & + + dvector(jj, ll, ii+3) & + * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & + - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & + * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) +! adjoint of: dw += tmp * vol_mm (vol frozen) + tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) +! adjoint of tmp (w1 frozen) + flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & + + (dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz)) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & + + dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & + + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & + + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd + end do + end do + end do + else +! scalar variable (vol, dscalar, w1 frozen) + do k = kl, 2, -1 + do j = jl, 2, -1 + do i = il, 2, -1 + flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & + + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & + * dwd(i,j,k,l) + end do + end do + end do + end if + end do + end do +! reverse of initialization: dw = wr (constant) => dwd = 0 do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 diff --git a/src/solver/residuals.F90 b/src/solver/residuals.F90 index 44d8a39a6..73c976c84 100644 --- a/src/solver/residuals.F90 +++ b/src/solver/residuals.F90 @@ -695,6 +695,7 @@ subroutine initres_block(varStart, varEnd, nn, sps) end select !=========================================================== +#endif case (timeSpectral) @@ -736,8 +737,10 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! Also store in ii the offset needed for vector ! quantities. +#ifndef USE_TAPENADE wsp => flowDoms(nn, currentLevel, mm)%w volsp => flowDoms(nn, currentLevel, mm)%vol +#endif ii = 3 * (mm - 1) ! Loop over the number of variables to be set. @@ -768,17 +771,29 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! Store the matrix vector product with the ! velocity in tmp. +#ifdef USE_TAPENADE + tmp = dvector(jj, ll, ii + 1) * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivx) & + + dvector(jj, ll, ii + 2) * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivy) & + + dvector(jj, ll, ii + 3) * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivz) +#else tmp = dvector(jj, ll, ii + 1) * wsp(i, j, k, ivx) & + dvector(jj, ll, ii + 2) * wsp(i, j, k, ivy) & + dvector(jj, ll, ii + 3) * wsp(i, j, k, ivz) +#endif ! Update the residual. Note the ! multiplication with the density to obtain ! the correct time derivative for the ! momentum variable. +#ifdef USE_TAPENADE + dw(i, j, k, l) = dw(i, j, k, l) & + + tmp * flowDoms(nn, currentLevel, mm)%vol(i, j, k) & + * flowDoms(nn, currentLevel, mm)%w(i, j, k, irho) +#else dw(i, j, k, l) = dw(i, j, k, l) & + tmp * volsp(i, j, k) * wsp(i, j, k, irho) +#endif end do end do @@ -793,9 +808,16 @@ subroutine initres_block(varStart, varEnd, nn, sps) do k = 2, kl do j = 2, jl do i = 2, il +#ifdef USE_TAPENADE + dw(i, j, k, l) = dw(i, j, k, l) & + + dscalar(jj, sps, mm) & + * flowDoms(nn, currentLevel, mm)%vol(i, j, k) & + * flowDoms(nn, currentLevel, mm)%w(i, j, k, l) +#else dw(i, j, k, l) = dw(i, j, k, l) & + dscalar(jj, sps, mm) & * volsp(i, j, k) * wsp(i, j, k, l) +#endif end do end do @@ -836,9 +858,11 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! Furthermore store in ii the offset needed for ! vector quantities. +#ifndef USE_TAPENADE wsp => flowDoms(nn, currentLevel, mm)%w wsp1 => flowDoms(nn, currentLevel, mm)%w1 volsp => flowDoms(nn, currentLevel, mm)%vol +#endif ii = 3 * (mm - 1) ! Loop over the number of variables to be set. @@ -873,6 +897,23 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! Store the matrix vector product with the ! momentum in tmp. +#ifdef USE_TAPENADE + tmp = dvector(jj, ll, ii + 1) & + * (flowDoms(nn, currentLevel, mm)%w(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivx) & + - flowDoms(nn, currentLevel, mm)%w1(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w1(i, j, k, ivx)) & + + dvector(jj, ll, ii + 2) & + * (flowDoms(nn, currentLevel, mm)%w(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivy) & + - flowDoms(nn, currentLevel, mm)%w1(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w1(i, j, k, ivy)) & + + dvector(jj, ll, ii + 3) & + * (flowDoms(nn, currentLevel, mm)%w(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w(i, j, k, ivz) & + - flowDoms(nn, currentLevel, mm)%w1(i, j, k, irho) & + * flowDoms(nn, currentLevel, mm)%w1(i, j, k, ivz)) +#else tmp = dvector(jj, ll, ii + 1) & * (wsp(i, j, k, irho) * wsp(i, j, k, ivx) & - wsp1(i, j, k, irho) * wsp1(i, j, k, ivx)) & @@ -882,13 +923,19 @@ subroutine initres_block(varStart, varEnd, nn, sps) + dvector(jj, ll, ii + 3) & * (wsp(i, j, k, irho) * wsp(i, j, k, ivz) & - wsp1(i, j, k, irho) * wsp1(i, j, k, ivz)) +#endif ! Add tmp to the residual. Multiply it by ! the volume to obtain the finite volume ! formulation of the derivative of the ! momentum. +#ifdef USE_TAPENADE + dw(i, j, k, l) = dw(i, j, k, l) & + + tmp * flowDoms(nn, currentLevel, mm)%vol(i, j, k) +#else dw(i, j, k, l) = dw(i, j, k, l) + tmp * volsp(i, j, k) +#endif end do end do @@ -903,10 +950,18 @@ subroutine initres_block(varStart, varEnd, nn, sps) do k = 2, kl do j = 2, jl do i = 2, il +#ifdef USE_TAPENADE + dw(i, j, k, l) = dw(i, j, k, l) & + + dscalar(jj, sps, mm) & + * flowDoms(nn, currentLevel, mm)%vol(i, j, k) & + * (flowDoms(nn, currentLevel, mm)%w(i, j, k, l) & + - flowDoms(nn, currentLevel, mm)%w1(i, j, k, l)) +#else dw(i, j, k, l) = dw(i, j, k, l) & + dscalar(jj, sps, mm) & * volsp(i, j, k) & * (wsp(i, j, k, l) - wsp1(i, j, k, l)) +#endif end do end do end do @@ -917,7 +972,6 @@ subroutine initres_block(varStart, varEnd, nn, sps) end do timeLoopCoarse end if spectralLevelTest -#endif end select ! Set the residual in the halo cells to zero. This is just diff --git a/tests/test_ts_adjoint.py b/tests/test_ts_adjoint.py new file mode 100644 index 000000000..ee254ef00 --- /dev/null +++ b/tests/test_ts_adjoint.py @@ -0,0 +1,347 @@ +""" +Test script for time spectral adjoint verification. + +Verifies the adjoint of the time spectral residual initialization +using dot product tests (forward vs backward consistency) and +forward mode AD vs finite difference comparison. + +Run with: mpirun -np 2 python test_ts_adjoint.py +""" + +import os +import sys +import copy +import numpy as np + +# Add the test directory to path for reg_default_options +sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "reg_tests")) + +from adflow import ADFLOW + +try: + from idwarp import USMesh + + HAS_IDWARP = True +except ImportError: + HAS_IDWARP = False + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() + +baseDir = os.path.dirname(os.path.abspath(__file__)) +inputDir = os.path.join(baseDir, "../input_files") +outputDir = os.path.join(baseDir, "output_files") + + +def printroot(msg): + if rank == 0: + print(msg, flush=True) + + +def setup_ts_solver(): + """Set up the time spectral solver with mesh deformation.""" + from baseclasses import AeroProblem + + gridFile = os.path.join(inputDir, "naca64A010_euler-L2.cgns") + + if not os.path.isfile(gridFile): + printroot(f"ERROR: Grid file not found: {gridFile}") + printroot("Run: cd input_files && bash get-input-files.sh") + sys.exit(1) + + ntimeintervalsspectral = 3 + + options = { + # Common defaults + "gridFile": gridFile, + "outputDirectory": outputDir, + "writeVolumeSolution": False, + "writeSurfaceSolution": False, + # Solver options + "blocksplitting": True, + "useblockettes": False, + "equationtype": "Euler", + "equationmode": "time spectral", + "mgcycle": "sg", + "l2convergence": 1e-13, + "ncycles": 200000, + "monitorvariables": ["resrho", "cl"], + # NK solver + "usenksolver": True, + "nkswitchtol": 1e-4, + "NKSubSpaceSize": 400, + "applypcsubspacesize": 400, + # ANK solver + "useanksolver": True, + "ankswitchtol": 1e-2, + "anksubspacesize": 50, + # Time spectral + "alphafollowing": False, + "timeintervals": ntimeintervalsspectral, + "useexternaldynamicmesh": True, + "usetsinterpolatedgridvelocity": True, + # Adjoint + "adjointl2convergence": 1e-14, + "adjointmaxiter": 500, + "nkadpc": False, + } + + # Create the solver + CFDSolver = ADFLOW(options=options, debug=True) + + # Setup aeroproblem + ap = AeroProblem( + name="64A010pitchingTS", + alpha=0.0, + mach=0.796, + reynolds=12.56e6, + reynoldsLength=1.0, + T=305.0, + R=287.085, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd", "cmz"], + xRef=0.248, + xRot=0.248, + omegaFourier=112.59, + ) + ap.addDV("alpha") + ap.addDV("mach") + + if HAS_IDWARP: + # Deform the mesh for time spectral + meshOptions = {"gridFile": gridFile} + mesh = USMesh(options=meshOptions) + CFDSolver.setMesh(mesh) + + # Motion history - pitching oscillation + alpha_0 = 1.01 + deltaAlpha = -alpha_0 * np.pi / 180.0 + alpha = np.linspace(0.0, 2.0 * np.pi, ntimeintervalsspectral + 1)[:-1] + alpha = -np.sin(alpha) + alpha *= deltaAlpha + + xRot = 0.25 + + # Get undeformed points + MDGroup = CFDSolver.allWallsGroup + cfdPts0 = CFDSolver.getSurfaceCoordinates(MDGroup, includeZipper=False) + N_pts = cfdPts0.shape[0] + + # Deform mesh for each time instance + for sps in range(ntimeintervalsspectral): + cfdPoints_deformed = np.zeros((N_pts, 3)) + ptch_loc = alpha[sps] + cc = np.cos(ptch_loc) + ss = np.sin(ptch_loc) + for j in range(N_pts): + cfdPoints_deformed[j, 0] = cc * (cfdPts0[j, 0] - xRot) + ss * cfdPts0[j, 1] + xRot + cfdPoints_deformed[j, 1] = -ss * (cfdPts0[j, 0] - xRot) + cc * cfdPts0[j, 1] + cfdPoints_deformed[j, 2] = cfdPts0[j, 2] + + mesh.setSurfaceCoordinates(cfdPoints_deformed) + mesh.warpMesh() + m = mesh.getSolverGrid() + CFDSolver.adflow.warping.setgridforoneinstance(m, sps=sps + 1) + + CFDSolver._updateGeomInfo = True + CFDSolver.updateGeometryInfo() + + return CFDSolver, ap + + +def test_dot_product(CFDSolver, ap): + """ + Dot product test: verifies forward and backward mode consistency. + + For a matrix J = dR/dw: + == + + This is the gold standard test for adjoint correctness. + """ + printroot("\n" + "=" * 60) + printroot("DOT PRODUCT TEST: Forward vs Backward mode consistency") + printroot("=" * 60) + + # Get random perturbation vectors + wDot = CFDSolver.getStatePerturbation(314) + dwBar = CFDSolver.getStatePerturbation(101) + + # Forward mode: resDot = J * wDot + resDot, funcsDot = CFDSolver.computeJacobianVectorProductFwd( + wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False + ) + + # Backward mode: wBar = J^T * dwBar + wBar = CFDSolver.computeJacobianVectorProductBwd( + resBar=dwBar, wDeriv=True, xVDeriv=False, xDvDerivAero=False + ) + + # Compute dot products locally then reduce + dotLocal1 = np.sum(resDot * dwBar) + dotLocal2 = np.sum(wDot * wBar) + + dotGlobal1 = comm.allreduce(dotLocal1, op=MPI.SUM) + dotGlobal2 = comm.allreduce(dotLocal2, op=MPI.SUM) + + rel_error = abs(dotGlobal1 - dotGlobal2) / (abs(dotGlobal1) + 1e-30) + + printroot(f" = {dotGlobal1:.15e}") + printroot(f" = {dotGlobal2:.15e}") + printroot(f" Relative error = {rel_error:.6e}") + + passed = rel_error < 1e-10 + printroot(f" Result: {'PASS' if passed else 'FAIL'}") + + return passed + + +def test_fwd_vs_fd(CFDSolver, ap): + """ + Forward mode AD vs Finite Difference comparison. + + Verifies dR/dw * wDot computed by AD matches the FD approximation: + (R(w + h*wDot) - R(w)) / h + """ + printroot("\n" + "=" * 60) + printroot("FORWARD MODE AD vs FINITE DIFFERENCE") + printroot("=" * 60) + + wDot = CFDSolver.getStatePerturbation(314) + + # AD forward mode + resDot_AD, funcsDot_AD = CFDSolver.computeJacobianVectorProductFwd( + wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False + ) + + # Finite difference + resDot_FD, funcsDot_FD = CFDSolver.computeJacobianVectorProductFwd( + wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False, mode="FD", h=1e-8 + ) + + # Compare residual derivatives + norm_AD = np.sqrt(comm.allreduce(np.sum(resDot_AD**2), op=MPI.SUM)) + norm_FD = np.sqrt(comm.allreduce(np.sum(resDot_FD**2), op=MPI.SUM)) + norm_diff = np.sqrt(comm.allreduce(np.sum((resDot_AD - resDot_FD) ** 2), op=MPI.SUM)) + + rel_error_res = norm_diff / (norm_AD + 1e-30) + printroot(f" Residual derivative:") + printroot(f" ||dR_AD|| = {norm_AD:.6e}") + printroot(f" ||dR_FD|| = {norm_FD:.6e}") + printroot(f" Relative error = {rel_error_res:.6e}") + res_passed = rel_error_res < 1e-3 + printroot(f" Result: {'PASS' if res_passed else 'FAIL'}") + + # Compare function derivatives + func_passed = True + for func in funcsDot_AD: + val_AD = funcsDot_AD[func] + val_FD = funcsDot_FD[func] + if abs(val_AD) > 1e-15: + rel_err = abs(val_AD - val_FD) / abs(val_AD) + else: + rel_err = abs(val_AD - val_FD) + ok = rel_err < 1e-4 + func_passed = func_passed and ok + printroot(f" {func}: AD={val_AD:.10e}, FD={val_FD:.10e}, rel_err={rel_err:.6e} {'PASS' if ok else 'FAIL'}") + + return res_passed and func_passed + + +def test_adjoint_sens(CFDSolver, ap): + """ + Test evalFunctionsSens - the main user-facing adjoint sensitivity method. + This solves the adjoint system and computes dF/dX for all design variables. + """ + printroot("\n" + "=" * 60) + printroot("ADJOINT SENSITIVITY TEST (evalFunctionsSens)") + printroot("=" * 60) + + funcsSens = {} + CFDSolver.evalFunctionsSens(ap, funcsSens) + + # Check adjoint failure + funcsCheck = {} + CFDSolver.checkAdjointFailure(ap, funcsCheck) + adj_converged = not funcsCheck.get("fail", True) + + printroot(f" Adjoint converged: {adj_converged}") + + for func_name, sens_dict in funcsSens.items(): + printroot(f" {func_name}:") + for dv_name, val in sens_dict.items(): + if isinstance(val, np.ndarray): + printroot(f" d/d({dv_name}) = array, norm={np.linalg.norm(val):.6e}") + else: + printroot(f" d/d({dv_name}) = {val:.10e}") + + return adj_converged + + +def main(): + printroot("=" * 60) + printroot("TIME SPECTRAL ADJOINT VERIFICATION TEST") + printroot("=" * 60) + + if not HAS_IDWARP: + printroot("WARNING: IDWarp not available. Skipping mesh deformation.") + printroot("The test will run without mesh deformation (rigid grid).") + + # Setup + printroot("\n--- Setting up solver ---") + CFDSolver, ap = setup_ts_solver() + + # Solve forward + printroot("\n--- Running forward solve ---") + CFDSolver(ap) + + # Check solution + funcs = {} + CFDSolver.checkSolutionFailure(ap, funcs) + if funcs.get("fail", True): + printroot("ERROR: Forward solve failed!") + sys.exit(1) + + # Evaluate functions + funcs = {} + CFDSolver.evalFunctions(ap, funcs) + printroot("\nFunction values:") + for k, v in funcs.items(): + printroot(f" {k} = {v:.10e}") + + # Initialize residual for adjoint tests + CFDSolver.getResidual(ap) + + # Run tests + results = {} + + # Test 1: Dot product test (most important) + results["dot_product"] = test_dot_product(CFDSolver, ap) + + # Test 2: Forward AD vs FD + results["fwd_vs_fd"] = test_fwd_vs_fd(CFDSolver, ap) + + # Test 3: Adjoint sensitivity + results["adjoint_sens"] = test_adjoint_sens(CFDSolver, ap) + + # Summary + printroot("\n" + "=" * 60) + printroot("SUMMARY") + printroot("=" * 60) + all_passed = True + for name, passed in results.items(): + status = "PASS" if passed else "FAIL" + printroot(f" {name}: {status}") + all_passed = all_passed and passed + + printroot(f"\nOverall: {'ALL TESTS PASSED' if all_passed else 'SOME TESTS FAILED'}") + printroot("=" * 60) + + return 0 if all_passed else 1 + + +if __name__ == "__main__": + ret = main() + sys.exit(ret) From a991fe5c1ba04187a2e55154404df352fc9bc306 Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Sun, 8 Mar 2026 13:39:51 -0400 Subject: [PATCH 02/12] more --- tests/test_ts_adjoint.py | 347 --------------------------------------- 1 file changed, 347 deletions(-) delete mode 100644 tests/test_ts_adjoint.py diff --git a/tests/test_ts_adjoint.py b/tests/test_ts_adjoint.py deleted file mode 100644 index ee254ef00..000000000 --- a/tests/test_ts_adjoint.py +++ /dev/null @@ -1,347 +0,0 @@ -""" -Test script for time spectral adjoint verification. - -Verifies the adjoint of the time spectral residual initialization -using dot product tests (forward vs backward consistency) and -forward mode AD vs finite difference comparison. - -Run with: mpirun -np 2 python test_ts_adjoint.py -""" - -import os -import sys -import copy -import numpy as np - -# Add the test directory to path for reg_default_options -sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "reg_tests")) - -from adflow import ADFLOW - -try: - from idwarp import USMesh - - HAS_IDWARP = True -except ImportError: - HAS_IDWARP = False - -from mpi4py import MPI - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() - -baseDir = os.path.dirname(os.path.abspath(__file__)) -inputDir = os.path.join(baseDir, "../input_files") -outputDir = os.path.join(baseDir, "output_files") - - -def printroot(msg): - if rank == 0: - print(msg, flush=True) - - -def setup_ts_solver(): - """Set up the time spectral solver with mesh deformation.""" - from baseclasses import AeroProblem - - gridFile = os.path.join(inputDir, "naca64A010_euler-L2.cgns") - - if not os.path.isfile(gridFile): - printroot(f"ERROR: Grid file not found: {gridFile}") - printroot("Run: cd input_files && bash get-input-files.sh") - sys.exit(1) - - ntimeintervalsspectral = 3 - - options = { - # Common defaults - "gridFile": gridFile, - "outputDirectory": outputDir, - "writeVolumeSolution": False, - "writeSurfaceSolution": False, - # Solver options - "blocksplitting": True, - "useblockettes": False, - "equationtype": "Euler", - "equationmode": "time spectral", - "mgcycle": "sg", - "l2convergence": 1e-13, - "ncycles": 200000, - "monitorvariables": ["resrho", "cl"], - # NK solver - "usenksolver": True, - "nkswitchtol": 1e-4, - "NKSubSpaceSize": 400, - "applypcsubspacesize": 400, - # ANK solver - "useanksolver": True, - "ankswitchtol": 1e-2, - "anksubspacesize": 50, - # Time spectral - "alphafollowing": False, - "timeintervals": ntimeintervalsspectral, - "useexternaldynamicmesh": True, - "usetsinterpolatedgridvelocity": True, - # Adjoint - "adjointl2convergence": 1e-14, - "adjointmaxiter": 500, - "nkadpc": False, - } - - # Create the solver - CFDSolver = ADFLOW(options=options, debug=True) - - # Setup aeroproblem - ap = AeroProblem( - name="64A010pitchingTS", - alpha=0.0, - mach=0.796, - reynolds=12.56e6, - reynoldsLength=1.0, - T=305.0, - R=287.085, - areaRef=1.0, - chordRef=1.0, - evalFuncs=["cl", "cd", "cmz"], - xRef=0.248, - xRot=0.248, - omegaFourier=112.59, - ) - ap.addDV("alpha") - ap.addDV("mach") - - if HAS_IDWARP: - # Deform the mesh for time spectral - meshOptions = {"gridFile": gridFile} - mesh = USMesh(options=meshOptions) - CFDSolver.setMesh(mesh) - - # Motion history - pitching oscillation - alpha_0 = 1.01 - deltaAlpha = -alpha_0 * np.pi / 180.0 - alpha = np.linspace(0.0, 2.0 * np.pi, ntimeintervalsspectral + 1)[:-1] - alpha = -np.sin(alpha) - alpha *= deltaAlpha - - xRot = 0.25 - - # Get undeformed points - MDGroup = CFDSolver.allWallsGroup - cfdPts0 = CFDSolver.getSurfaceCoordinates(MDGroup, includeZipper=False) - N_pts = cfdPts0.shape[0] - - # Deform mesh for each time instance - for sps in range(ntimeintervalsspectral): - cfdPoints_deformed = np.zeros((N_pts, 3)) - ptch_loc = alpha[sps] - cc = np.cos(ptch_loc) - ss = np.sin(ptch_loc) - for j in range(N_pts): - cfdPoints_deformed[j, 0] = cc * (cfdPts0[j, 0] - xRot) + ss * cfdPts0[j, 1] + xRot - cfdPoints_deformed[j, 1] = -ss * (cfdPts0[j, 0] - xRot) + cc * cfdPts0[j, 1] - cfdPoints_deformed[j, 2] = cfdPts0[j, 2] - - mesh.setSurfaceCoordinates(cfdPoints_deformed) - mesh.warpMesh() - m = mesh.getSolverGrid() - CFDSolver.adflow.warping.setgridforoneinstance(m, sps=sps + 1) - - CFDSolver._updateGeomInfo = True - CFDSolver.updateGeometryInfo() - - return CFDSolver, ap - - -def test_dot_product(CFDSolver, ap): - """ - Dot product test: verifies forward and backward mode consistency. - - For a matrix J = dR/dw: - == - - This is the gold standard test for adjoint correctness. - """ - printroot("\n" + "=" * 60) - printroot("DOT PRODUCT TEST: Forward vs Backward mode consistency") - printroot("=" * 60) - - # Get random perturbation vectors - wDot = CFDSolver.getStatePerturbation(314) - dwBar = CFDSolver.getStatePerturbation(101) - - # Forward mode: resDot = J * wDot - resDot, funcsDot = CFDSolver.computeJacobianVectorProductFwd( - wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False - ) - - # Backward mode: wBar = J^T * dwBar - wBar = CFDSolver.computeJacobianVectorProductBwd( - resBar=dwBar, wDeriv=True, xVDeriv=False, xDvDerivAero=False - ) - - # Compute dot products locally then reduce - dotLocal1 = np.sum(resDot * dwBar) - dotLocal2 = np.sum(wDot * wBar) - - dotGlobal1 = comm.allreduce(dotLocal1, op=MPI.SUM) - dotGlobal2 = comm.allreduce(dotLocal2, op=MPI.SUM) - - rel_error = abs(dotGlobal1 - dotGlobal2) / (abs(dotGlobal1) + 1e-30) - - printroot(f" = {dotGlobal1:.15e}") - printroot(f" = {dotGlobal2:.15e}") - printroot(f" Relative error = {rel_error:.6e}") - - passed = rel_error < 1e-10 - printroot(f" Result: {'PASS' if passed else 'FAIL'}") - - return passed - - -def test_fwd_vs_fd(CFDSolver, ap): - """ - Forward mode AD vs Finite Difference comparison. - - Verifies dR/dw * wDot computed by AD matches the FD approximation: - (R(w + h*wDot) - R(w)) / h - """ - printroot("\n" + "=" * 60) - printroot("FORWARD MODE AD vs FINITE DIFFERENCE") - printroot("=" * 60) - - wDot = CFDSolver.getStatePerturbation(314) - - # AD forward mode - resDot_AD, funcsDot_AD = CFDSolver.computeJacobianVectorProductFwd( - wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False - ) - - # Finite difference - resDot_FD, funcsDot_FD = CFDSolver.computeJacobianVectorProductFwd( - wDot=wDot, residualDeriv=True, funcDeriv=True, fDeriv=False, mode="FD", h=1e-8 - ) - - # Compare residual derivatives - norm_AD = np.sqrt(comm.allreduce(np.sum(resDot_AD**2), op=MPI.SUM)) - norm_FD = np.sqrt(comm.allreduce(np.sum(resDot_FD**2), op=MPI.SUM)) - norm_diff = np.sqrt(comm.allreduce(np.sum((resDot_AD - resDot_FD) ** 2), op=MPI.SUM)) - - rel_error_res = norm_diff / (norm_AD + 1e-30) - printroot(f" Residual derivative:") - printroot(f" ||dR_AD|| = {norm_AD:.6e}") - printroot(f" ||dR_FD|| = {norm_FD:.6e}") - printroot(f" Relative error = {rel_error_res:.6e}") - res_passed = rel_error_res < 1e-3 - printroot(f" Result: {'PASS' if res_passed else 'FAIL'}") - - # Compare function derivatives - func_passed = True - for func in funcsDot_AD: - val_AD = funcsDot_AD[func] - val_FD = funcsDot_FD[func] - if abs(val_AD) > 1e-15: - rel_err = abs(val_AD - val_FD) / abs(val_AD) - else: - rel_err = abs(val_AD - val_FD) - ok = rel_err < 1e-4 - func_passed = func_passed and ok - printroot(f" {func}: AD={val_AD:.10e}, FD={val_FD:.10e}, rel_err={rel_err:.6e} {'PASS' if ok else 'FAIL'}") - - return res_passed and func_passed - - -def test_adjoint_sens(CFDSolver, ap): - """ - Test evalFunctionsSens - the main user-facing adjoint sensitivity method. - This solves the adjoint system and computes dF/dX for all design variables. - """ - printroot("\n" + "=" * 60) - printroot("ADJOINT SENSITIVITY TEST (evalFunctionsSens)") - printroot("=" * 60) - - funcsSens = {} - CFDSolver.evalFunctionsSens(ap, funcsSens) - - # Check adjoint failure - funcsCheck = {} - CFDSolver.checkAdjointFailure(ap, funcsCheck) - adj_converged = not funcsCheck.get("fail", True) - - printroot(f" Adjoint converged: {adj_converged}") - - for func_name, sens_dict in funcsSens.items(): - printroot(f" {func_name}:") - for dv_name, val in sens_dict.items(): - if isinstance(val, np.ndarray): - printroot(f" d/d({dv_name}) = array, norm={np.linalg.norm(val):.6e}") - else: - printroot(f" d/d({dv_name}) = {val:.10e}") - - return adj_converged - - -def main(): - printroot("=" * 60) - printroot("TIME SPECTRAL ADJOINT VERIFICATION TEST") - printroot("=" * 60) - - if not HAS_IDWARP: - printroot("WARNING: IDWarp not available. Skipping mesh deformation.") - printroot("The test will run without mesh deformation (rigid grid).") - - # Setup - printroot("\n--- Setting up solver ---") - CFDSolver, ap = setup_ts_solver() - - # Solve forward - printroot("\n--- Running forward solve ---") - CFDSolver(ap) - - # Check solution - funcs = {} - CFDSolver.checkSolutionFailure(ap, funcs) - if funcs.get("fail", True): - printroot("ERROR: Forward solve failed!") - sys.exit(1) - - # Evaluate functions - funcs = {} - CFDSolver.evalFunctions(ap, funcs) - printroot("\nFunction values:") - for k, v in funcs.items(): - printroot(f" {k} = {v:.10e}") - - # Initialize residual for adjoint tests - CFDSolver.getResidual(ap) - - # Run tests - results = {} - - # Test 1: Dot product test (most important) - results["dot_product"] = test_dot_product(CFDSolver, ap) - - # Test 2: Forward AD vs FD - results["fwd_vs_fd"] = test_fwd_vs_fd(CFDSolver, ap) - - # Test 3: Adjoint sensitivity - results["adjoint_sens"] = test_adjoint_sens(CFDSolver, ap) - - # Summary - printroot("\n" + "=" * 60) - printroot("SUMMARY") - printroot("=" * 60) - all_passed = True - for name, passed in results.items(): - status = "PASS" if passed else "FAIL" - printroot(f" {name}: {status}") - all_passed = all_passed and passed - - printroot(f"\nOverall: {'ALL TESTS PASSED' if all_passed else 'SOME TESTS FAILED'}") - printroot("=" * 60) - - return 0 if all_passed else 1 - - -if __name__ == "__main__": - ret = main() - sys.exit(ret) From 3a728b7528780ccb1bb61b8f57f8203aee7fe052 Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Sun, 8 Mar 2026 14:56:49 -0400 Subject: [PATCH 03/12] fix ad --- src/adjoint/Makefile_tapenade | 4 +- src/adjoint/outputForward/residuals_d.f90 | 473 ++++++++++++---- src/adjoint/outputReverse/residuals_b.f90 | 630 ++++++++++++++++------ 3 files changed, 838 insertions(+), 269 deletions(-) diff --git a/src/adjoint/Makefile_tapenade b/src/adjoint/Makefile_tapenade index 963e5bcf5..509054bb0 100644 --- a/src/adjoint/Makefile_tapenade +++ b/src/adjoint/Makefile_tapenade @@ -186,8 +186,8 @@ sa%saResScale(scratch, dw) > \ residuals%sourceTerms_block(w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat, actuatorRegions%volume) > \ (w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat, actuatorRegions%volume) \ \ -residuals%initres_block(dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) > \ - (dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) \ +residuals%initres_block(dw, fw, flowDoms%w, flowDoms%vol) > \ + (dw, fw, flowDoms%w, flowDoms%vol) \ \ fluxes%inviscidCentralFlux(w, timeRef, vol, si, sj, sk, p, dw, sfacei, sfacej, sfacek) > \ (w, timeRef, vol, si, sj, sk, p, dw, sfacei, sfacej, sfacek) \ diff --git a/src/adjoint/outputForward/residuals_d.f90 b/src/adjoint/outputForward/residuals_d.f90 index 0590c069f..6588e8cac 100644 --- a/src/adjoint/outputForward/residuals_d.f90 +++ b/src/adjoint/outputForward/residuals_d.f90 @@ -508,9 +508,12 @@ end subroutine sourceterms_block ! differentiation of initres_block in forward (tangent) mode (with options i4 dr8 r8): ! variations of useful results: *dw -! with respect to varying inputs: *dw *flowDoms.w *flowDoms.vol -! rw status of diff variables: *dw:in-out *flowDoms.w:in *flowDoms.vol:in -! plus diff mem management of: dw:in +! with respect to varying inputs: *(flowdoms.vol) *(flowdoms.w) +! *dw +! rw status of diff variables: *(flowdoms.vol):in *(flowdoms.w):in +! *dw:in-out +! plus diff mem management of: flowdoms.vol:in flowdoms.w:in +! dw:in subroutine initres_block_d(varstart, varend, nn, sps) ! ! initres initializes the given range of the residual. either to @@ -534,16 +537,26 @@ subroutine initres_block_d(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp, tmpd + real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp + real(kind=realtype) :: temp + real(kind=realtype) :: temp0 + real(kind=realtype) :: temp1 + real(kind=realtype) :: temp2 + real(kind=realtype) :: temp3 + real(kind=realtype) :: temp4 + real(kind=realtype) :: temp5 + real(kind=realtype) :: temp6 + real(kind=realtype) :: temp7 ! return immediately of no variables are in the range. if (varend .lt. varstart) then return else ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then @@ -573,11 +586,18 @@ subroutine initres_block_d(varstart, varend, nn, sps) end do end do end if - case (timespectral) -! time spectral computation. + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. jj = sectionid +! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then -! fine grid: initialize to zero, then accumulate spectral time derivative +! finest multigrid level. the residual must be +! initialized to the time derivative. +! initialize it to zero. do l=varstart,varend do k=2,kl do j=2,jl @@ -588,61 +608,88 @@ subroutine initres_block_d(varstart, varend, nn, sps) end do end do end do -! loop over spectral instances - do mm = 1, ntimeintervalsspectral - ii = 3 * (mm - 1) - do l = varstart, varend - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. if (l .eq. ivx) ll = 3*sps - 2 if (l .eq. ivy) ll = 3*sps - 1 if (l .eq. ivz) ll = 3*sps - do k = 2, kl - do j = 2, jl - do i = 2, il -! tangent of: tmp = dvector * velocities - tmpd = dvector(jj, ll, ii+1) * flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdomsd(nn, 1, mm)%w(i,j,k,ivz) - tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) -! tangent of: dw += tmp * vol_mm * rho_mm - dwd(i,j,k,l) = dwd(i,j,k,l) & - + tmpd * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - + tmp * flowdomsd(nn, 1, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdomsd(nn, 1, mm)%w(i,j,k,irho) - dw(i,j,k,l) = dw(i,j,k,l) & - + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + temp = dvector(jj, ll, ii+1) + temp0 = dvector(jj, ll, ii+2) + temp1 = dvector(jj, ll, ii+3) + tmpd = temp*flowdomsd(nn, currentlevel, mm)%w(i, j& +& , k, ivx) + temp0*flowdomsd(nn, currentlevel, mm& +& )%w(i, j, k, ivy) + temp1*flowdomsd(nn, & +& currentlevel, mm)%w(i, j, k, ivz) + tmp = temp*flowdoms(nn, currentlevel, mm)%w(i, j, & +& k, ivx) + temp0*flowdoms(nn, currentlevel, mm)%w& +& (i, j, k, ivy) + temp1*flowdoms(nn, currentlevel& +& , mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + temp1 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho) + temp0 = flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& ) + dwd(i, j, k, l) = dwd(i, j, k, l) + temp1*(temp0*& +& tmpd+tmp*flowdomsd(nn, currentlevel, mm)%vol(i, & +& j, k)) + tmp*temp0*flowdomsd(nn, currentlevel, & +& mm)%w(i, j, k, irho) + dw(i, j, k, l) = dw(i, j, k, l) + tmp*temp0*temp1 end do end do end do else -! scalar variable - do k = 2, kl - do j = 2, jl - do i = 2, il -! tangent of: dw += dscalar * vol_mm * w_mm - dwd(i,j,k,l) = dwd(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdomsd(nn, 1, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdomsd(nn, 1, mm)%w(i,j,k,l) - dw(i,j,k,l) = dw(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) +! scalar variable. loop over the owned cells to +! add the contribution of wsp to the time +! derivative. + do k=2,kl + do j=2,jl + do i=2,il + temp1 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& l) + temp0 = flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& ) + dwd(i, j, k, l) = dwd(i, j, k, l) + dscalar(jj, & +& sps, mm)*(temp1*flowdomsd(nn, currentlevel, mm)%& +& vol(i, j, k)+temp0*flowdomsd(nn, currentlevel, & +& mm)%w(i, j, k, l)) + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*(temp0*temp1) end do end do end do end if - end do - end do + end do varloopfine + end do timeloopfine else -! coarse grid level +! coarse grid level. initialize the owned cells to the +! residual forcing term plus a correction for the +! multigrid treatment of the time derivative term. +! initialization to the residual forcing term. do l=varstart,varend do k=2,kl do j=2,jl @@ -653,82 +700,106 @@ subroutine initres_block_d(varstart, varend, nn, sps) end do end do end do -! loop over spectral instances (coarse grid correction) - do mm = 1, ntimeintervalsspectral - ii = 3 * (mm - 1) - do l = varstart, varend - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. if (l .eq. ivx) ll = 3*sps - 2 if (l .eq. ivy) ll = 3*sps - 1 if (l .eq. ivz) ll = 3*sps - do k = 2, kl - do j = 2, jl - do i = 2, il -! tangent of: tmp = dvector * (rho*vel - rho1*vel1) -! w1 is frozen so w1d = 0 - tmpd = dvector(jj, ll, ii+1) & - * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdomsd(nn, 1, mm)%w(i,j,k,ivx)) & - + dvector(jj, ll, ii+2) & - * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdomsd(nn, 1, mm)%w(i,j,k,ivy)) & - + dvector(jj, ll, ii+3) & - * (flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & - + flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdomsd(nn, 1, mm)%w(i,j,k,ivz)) - tmp = dvector(jj, ll, ii+1) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & - + dvector(jj, ll, ii+2) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & - + dvector(jj, ll, ii+3) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) -! tangent of: dw += tmp * vol_mm - dwd(i,j,k,l) = dwd(i,j,k,l) & - + tmpd * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - + tmp * flowdomsd(nn, 1, mm)%vol(i,j,k) - dw(i,j,k,l) = dw(i,j,k,l) & - + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) +! add the contribution of wps to the correction +! of the time derivative. the difference between +! the current time derivative and the one when +! entering this grid level must be added, because +! the residual forcing term only takes the spatial +! part of the coarse grid residual into account. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! momentum in tmp. + temp1 = dvector(jj, ll, ii+1) + temp0 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivx) + temp = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho) + temp2 = dvector(jj, ll, ii+2) + temp3 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivy) + temp4 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho) + temp5 = dvector(jj, ll, ii+3) + temp6 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivz) + temp7 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho) + tmpd = temp1*(temp0*flowdomsd(nn, currentlevel, mm& +& )%w(i, j, k, irho)+temp*flowdomsd(nn, & +& currentlevel, mm)%w(i, j, k, ivx)) + temp2*(& +& temp3*flowdomsd(nn, currentlevel, mm)%w(i, j, k& +& , irho)+temp4*flowdomsd(nn, currentlevel, mm)%w(& +& i, j, k, ivy)) + temp5*(temp6*flowdomsd(nn, & +& currentlevel, mm)%w(i, j, k, irho)+temp7*& +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz)) + tmp = temp1*(temp*temp0-flowdoms(nn, currentlevel& +& , mm)%w1(i, j, k, irho)*flowdoms(nn, & +& currentlevel, mm)%w1(i, j, k, ivx)) + temp2*(& +& temp4*temp3-flowdoms(nn, currentlevel, mm)%w1(i& +& , j, k, irho)*flowdoms(nn, currentlevel, mm)%w1(& +& i, j, k, ivy)) + temp5*(temp7*temp6-flowdoms(nn& +& , currentlevel, mm)%w1(i, j, k, irho)*flowdoms(& +& nn, currentlevel, mm)%w1(i, j, k, ivz)) +! add tmp to the residual. multiply it by +! the volume to obtain the finite volume +! formulation of the derivative of the +! momentum. + temp7 = flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& ) + dwd(i, j, k, l) = dwd(i, j, k, l) + temp7*tmpd + & +& tmp*flowdomsd(nn, currentlevel, mm)%vol(i, j, k) + dw(i, j, k, l) = dw(i, j, k, l) + tmp*temp7 end do end do end do else -! scalar variable - do k = 2, kl - do j = 2, jl - do i = 2, il -! tangent of: dw += dscalar * vol_mm * (w_mm - w1_mm) -! w1 is frozen so w1d = 0 - dwd(i,j,k,l) = dwd(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdomsd(nn, 1, mm)%vol(i,j,k) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdomsd(nn, 1, mm)%w(i,j,k,l) - dw(i,j,k,l) = dw(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) +! scalar variable. loop over the owned cells +! to add the contribution of wsp to the correction +! of the time derivative. + do k=2,kl + do j=2,jl + do i=2,il + temp7 = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& l) - flowdoms(nn, currentlevel, mm)%w1(i, j, k, & +& l) + temp6 = flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& ) + dwd(i, j, k, l) = dwd(i, j, k, l) + dscalar(jj, & +& sps, mm)*(temp7*flowdomsd(nn, currentlevel, mm)%& +& vol(i, j, k)+temp6*flowdomsd(nn, currentlevel, & +& mm)%w(i, j, k, l)) + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*(temp6*temp7) end do end do end do end if - end do - end do + end do varloopcoarse + end do timeloopcoarse end if end select ! set the residual in the halo cells to zero. this is just @@ -834,6 +905,174 @@ subroutine initres_block(varstart, varend, nn, sps) end do end do end if + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. + jj = sectionid +! determine the currently active multigrid level. + if (currentlevel .eq. groundlevel) then +! finest multigrid level. the residual must be +! initialized to the time derivative. +! initialize it to zero. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = zero + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + tmp = dvector(jj, ll, ii+1)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx) + dvector(jj, & +& ll, ii+2)*flowdoms(nn, currentlevel, mm)%w(i, j& +& , k, ivy) + dvector(jj, ll, ii+3)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho) + end do + end do + end do + else +! scalar variable. loop over the owned cells to +! add the contribution of wsp to the time +! derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*flowdoms(nn, currentlevel, mm)%w(i, j, k, l) + end do + end do + end do + end if + end do varloopfine + end do timeloopfine + else +! coarse grid level. initialize the owned cells to the +! residual forcing term plus a correction for the +! multigrid treatment of the time derivative term. +! initialization to the residual forcing term. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = wr(i, j, k, l) + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! add the contribution of wps to the correction +! of the time derivative. the difference between +! the current time derivative and the one when +! entering this grid level must be added, because +! the residual forcing term only takes the spatial +! part of the coarse grid residual into account. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! momentum in tmp. + tmp = dvector(jj, ll, ii+1)*(flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx)-flowdoms(nn, & +& currentlevel, mm)%w1(i, j, k, irho)*flowdoms(nn& +& , currentlevel, mm)%w1(i, j, k, ivx)) + dvector(& +& jj, ll, ii+2)*(flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, ivy)-flowdoms(nn, currentlevel, mm)%w1(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w1& +& (i, j, k, ivy)) + dvector(jj, ll, ii+3)*(& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, irho)*& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, ivz)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, irho)& +& *flowdoms(nn, currentlevel, mm)%w1(i, j, k, ivz)& +& ) +! add tmp to the residual. multiply it by +! the volume to obtain the finite volume +! formulation of the derivative of the +! momentum. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k) + end do + end do + end do + else +! scalar variable. loop over the owned cells +! to add the contribution of wsp to the correction +! of the time derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*(flowdoms(nn, currentlevel, mm)%w(i, j, k, l)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, l)) + end do + end do + end do + end if + end do varloopcoarse + end do timeloopcoarse + end if end select ! set the residual in the halo cells to zero. this is just ! to avoid possible problems. their values do not matter. diff --git a/src/adjoint/outputReverse/residuals_b.f90 b/src/adjoint/outputReverse/residuals_b.f90 index c513ba104..8a831baae 100644 --- a/src/adjoint/outputReverse/residuals_b.f90 +++ b/src/adjoint/outputReverse/residuals_b.f90 @@ -530,10 +530,14 @@ subroutine sourceterms_block(nn, res, iregion, plocal) end subroutine sourceterms_block ! differentiation of initres_block in reverse (adjoint) mode (with options noisize i4 dr8 r8): -! gradient of useful results: *dw -! with respect to varying inputs: *dw *flowDoms.w *flowDoms.vol -! rw status of diff variables: *dw:in-out *flowDoms.w:out *flowDoms.vol:out -! plus diff mem management of: dw:in +! gradient of useful results: *(flowdoms.vol) *(flowdoms.w) +! *dw +! with respect to varying inputs: *(flowdoms.vol) *(flowdoms.w) +! *dw +! rw status of diff variables: *(flowdoms.vol):incr *(flowdoms.w):incr +! *dw:in-out +! plus diff mem management of: flowdoms.vol:in flowdoms.w:in +! dw:in subroutine initres_block_b(varstart, varend, nn, sps) ! ! initres initializes the given range of the residual. either to @@ -557,33 +561,189 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp, tmpd + real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp + real(kind=realtype) :: tempd + real(kind=realtype) :: temp + real(kind=realtype) :: tempd0 + real(kind=realtype) :: tempd1 integer :: branch ! return immediately of no variables are in the range. if (varend .ge. varstart) then ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol3b(1) + call pushcontrol3b(3) else - call pushcontrol3b(0) + call pushcontrol3b(2) end if - case (timespectral) -! time spectral computation. + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. + jj = sectionid +! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol3b(2) +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + call pushinteger4(ii) + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) then + call pushinteger4(ll) + ll = 3*sps - 2 + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (l .eq. ivy) then + call pushinteger4(ll) + ll = 3*sps - 1 + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (l .eq. ivz) then + call pushinteger4(ll) + ll = 3*sps + call pushcontrol1b(1) + else + call pushcontrol1b(0) + end if +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + call pushreal8(tmp) + tmp = dvector(jj, ll, ii+1)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx) + dvector(jj, & +& ll, ii+2)*flowdoms(nn, currentlevel, mm)%w(i, j& +& , k, ivy) + dvector(jj, ll, ii+3)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + end do + end do + end do + call pushcontrol1b(1) + else + call pushcontrol1b(0) + end if + end do varloopfine + end do timeloopfine + call pushcontrol3b(1) else - call pushcontrol3b(3) +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + call pushinteger4(ii) + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) then + call pushinteger4(ll) + ll = 3*sps - 2 + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (l .eq. ivy) then + call pushinteger4(ll) + ll = 3*sps - 1 + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (l .eq. ivz) then + call pushinteger4(ll) + ll = 3*sps + call pushcontrol1b(1) + else + call pushcontrol1b(0) + end if +! add the contribution of wps to the correction +! of the time derivative. the difference between +! the current time derivative and the one when +! entering this grid level must be added, because +! the residual forcing term only takes the spatial +! part of the coarse grid residual into account. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! momentum in tmp. + call pushreal8(tmp) + tmp = dvector(jj, ll, ii+1)*(flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx)-flowdoms(nn, & +& currentlevel, mm)%w1(i, j, k, irho)*flowdoms(nn& +& , currentlevel, mm)%w1(i, j, k, ivx)) + dvector(& +& jj, ll, ii+2)*(flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, ivy)-flowdoms(nn, currentlevel, mm)%w1(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w1& +& (i, j, k, ivy)) + dvector(jj, ll, ii+3)*(& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, irho)*& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, ivz)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, irho)& +& *flowdoms(nn, currentlevel, mm)%w1(i, j, k, ivz)& +& ) +! add tmp to the residual. multiply it by +! the volume to obtain the finite volume +! formulation of the derivative of the +! momentum. + end do + end do + end do + call pushcontrol1b(1) + else + call pushcontrol1b(0) + end if + end do varloopcoarse + end do timeloopcoarse + call pushcontrol3b(0) end if case default call pushcontrol3b(4) end select -! reverse the halo cell zeroing do l=varend,varstart,-1 do j=jl,2,-1 do i=il,2,-1 @@ -611,83 +771,155 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do call popcontrol3b(branch) - if (branch .eq. 0) then -! steady, coarse grid: forward did dw = wr (constant) - do l=varend,varstart,-1 - do k=kl,2,-1 - do j=jl,2,-1 - do i=il,2,-1 - dwd(i, j, k, l) = 0.0_8 - end do + if (branch .lt. 2) then + if (branch .eq. 0) then + do mm=ntimeintervalsspectral,1,-1 + do l=varend,varstart,-1 + call popcontrol1b(branch) + if (branch .eq. 0) then + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + tempd1 = dscalar(jj, sps, mm)*dwd(i, j, k, l) + flowdomsd(nn, currentlevel, mm)%vol(i, j, k) = & +& flowdomsd(nn, currentlevel, mm)%vol(i, j, k) + (& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, l)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, l))*& +& tempd1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & +& flowdoms(nn, currentlevel, mm)%vol(i, j, k)*& +& tempd1 + end do + end do + end do + else + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + tmpd = flowdoms(nn, currentlevel, mm)%vol(i, j, k)& +& *dwd(i, j, k, l) + flowdomsd(nn, currentlevel, mm)%vol(i, j, k) = & +& flowdomsd(nn, currentlevel, mm)%vol(i, j, k) + & +& tmp*dwd(i, j, k, l) + call popreal8(tmp) + tempd0 = dvector(jj, ll, ii+1)*tmpd + tempd = dvector(jj, ll, ii+2)*tmpd + tempd1 = dvector(jj, ll, ii+3)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivz)*tempd1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivy)*tempd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivx)*tempd0 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd0 + end do + end do + end do + call popcontrol1b(branch) + if (branch .ne. 0) call popinteger4(ll) + call popcontrol1b(branch) + if (branch .eq. 0) call popinteger4(ll) + call popcontrol1b(branch) + if (branch .eq. 0) call popinteger4(ll) + end if end do + call popinteger4(ii) end do - end do - else if (branch .eq. 1) then -! steady, fine grid: forward did dw = zero - do l=varend,varstart,-1 - do k=kl,2,-1 - do j=jl,2,-1 - do i=il,2,-1 - dwd(i, j, k, l) = 0.0_8 + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do end do end do end do - end do - else if (branch .eq. 2) then -! time spectral, fine grid -! reverse the spectral time derivative accumulation - jj = sectionid - do mm = ntimeintervalsspectral, 1, -1 - ii = 3 * (mm - 1) - do l = varend, varstart, -1 - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 -! recompute tmp from forward pass - tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) -! adjoint of: dw(i,j,k,l) += tmp * vol_mm * rho_mm - tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & - + tmp * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) -! adjoint of: tmp = dvector * velocities - flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+1) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+2) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & - + dvector(jj, ll, ii+3) * tmpd + else + do mm=ntimeintervalsspectral,1,-1 + do l=varend,varstart,-1 + call popcontrol1b(branch) + if (branch .eq. 0) then + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + tempd0 = dscalar(jj, sps, mm)*dwd(i, j, k, l) + flowdomsd(nn, currentlevel, mm)%vol(i, j, k) = & +& flowdomsd(nn, currentlevel, mm)%vol(i, j, k) + & +& flowdoms(nn, currentlevel, mm)%w(i, j, k, l)*& +& tempd0 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & +& flowdoms(nn, currentlevel, mm)%vol(i, j, k)*& +& tempd0 + end do end do end do - end do - else -! scalar variable -! adjoint of: dw(i,j,k,l) += dscalar(jj,sps,mm) * vol_mm * w_mm(l) - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 - flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & - * dwd(i,j,k,l) + else + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + temp = flowdoms(nn, currentlevel, mm)%vol(i, j, k) + tempd = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho)*dwd(i, j, k, l) + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + tmp*temp*dwd(i, j, k, l) + tmpd = temp*tempd + flowdomsd(nn, currentlevel, mm)%vol(i, j, k) = & +& flowdomsd(nn, currentlevel, mm)%vol(i, j, k) + & +& tmp*tempd + call popreal8(tmp) + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & +& + dvector(jj, ll, ii+1)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & +& + dvector(jj, ll, ii+2)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & +& + dvector(jj, ll, ii+3)*tmpd + end do end do end do + call popcontrol1b(branch) + if (branch .ne. 0) call popinteger4(ll) + call popcontrol1b(branch) + if (branch .eq. 0) call popinteger4(ll) + call popcontrol1b(branch) + if (branch .eq. 0) call popinteger4(ll) + end if + end do + call popinteger4(ii) + end do + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do end do - end if + end do end do - end do -! reverse of initialization: dw = zero => dwd = 0 + end if + else if (branch .eq. 2) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -698,76 +930,6 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do else if (branch .eq. 3) then -! time spectral, coarse grid -! reverse the coarse grid spectral time derivative accumulation - jj = sectionid - do mm = ntimeintervalsspectral, 1, -1 - ii = 3 * (mm - 1) - do l = varend, varstart, -1 - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 -! recompute tmp from forward pass - tmp = dvector(jj, ll, ii+1) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & - + dvector(jj, ll, ii+2) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & - + dvector(jj, ll, ii+3) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) -! adjoint of: dw += tmp * vol_mm - tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & - + tmp * dwd(i,j,k,l) -! adjoint of tmp = dvector * (rho*vel - rho1*vel1) -! only w is active, w1 is frozen - flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - + (dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz)) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - end do - end do - end do - else -! scalar variable -! adjoint of: dw += dscalar * vol_mm * (w_mm - w1_mm) -! only w is active, w1 is frozen - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 - flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%vol(i,j,k) = flowdomsd(nn, 1, mm)%vol(i,j,k) & - + dscalar(jj, sps, mm) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,l) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,l)) * dwd(i,j,k,l) - end do - end do - end do - end if - end do - end do -! reverse of initialization: dw = wr (constant) => dwd = 0 do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -841,6 +1003,174 @@ subroutine initres_block(varstart, varend, nn, sps) end do end do end if + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. + jj = sectionid +! determine the currently active multigrid level. + if (currentlevel .eq. groundlevel) then +! finest multigrid level. the residual must be +! initialized to the time derivative. +! initialize it to zero. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = zero + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + tmp = dvector(jj, ll, ii+1)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx) + dvector(jj, & +& ll, ii+2)*flowdoms(nn, currentlevel, mm)%w(i, j& +& , k, ivy) + dvector(jj, ll, ii+3)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho) + end do + end do + end do + else +! scalar variable. loop over the owned cells to +! add the contribution of wsp to the time +! derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*flowdoms(nn, currentlevel, mm)%w(i, j, k, l) + end do + end do + end do + end if + end do varloopfine + end do timeloopfine + else +! coarse grid level. initialize the owned cells to the +! residual forcing term plus a correction for the +! multigrid treatment of the time derivative term. +! initialization to the residual forcing term. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = wr(i, j, k, l) + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! add the contribution of wps to the correction +! of the time derivative. the difference between +! the current time derivative and the one when +! entering this grid level must be added, because +! the residual forcing term only takes the spatial +! part of the coarse grid residual into account. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! momentum in tmp. + tmp = dvector(jj, ll, ii+1)*(flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx)-flowdoms(nn, & +& currentlevel, mm)%w1(i, j, k, irho)*flowdoms(nn& +& , currentlevel, mm)%w1(i, j, k, ivx)) + dvector(& +& jj, ll, ii+2)*(flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, ivy)-flowdoms(nn, currentlevel, mm)%w1(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w1& +& (i, j, k, ivy)) + dvector(jj, ll, ii+3)*(& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, irho)*& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, ivz)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, irho)& +& *flowdoms(nn, currentlevel, mm)%w1(i, j, k, ivz)& +& ) +! add tmp to the residual. multiply it by +! the volume to obtain the finite volume +! formulation of the derivative of the +! momentum. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k) + end do + end do + end do + else +! scalar variable. loop over the owned cells +! to add the contribution of wsp to the correction +! of the time derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*(flowdoms(nn, currentlevel, mm)%w(i, j, k, l)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, l)) + end do + end do + end do + end if + end do varloopcoarse + end do timeloopcoarse + end if end select ! set the residual in the halo cells to zero. this is just ! to avoid possible problems. their values do not matter. From e81c9e8df8730144aa358790f07a14604faa887c Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Sun, 8 Mar 2026 15:20:19 -0400 Subject: [PATCH 04/12] clean --- .../outputReverseFast/residuals_fast_b.f90 | 567 ++++++++++++++---- 1 file changed, 435 insertions(+), 132 deletions(-) diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index ce273184c..fca35d480 100644 --- a/src/adjoint/outputReverseFast/residuals_fast_b.f90 +++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90 @@ -479,10 +479,10 @@ subroutine sourceterms_block(nn, res, iregion, plocal) end subroutine sourceterms_block ! differentiation of initres_block in reverse (adjoint) mode (with options noisize i4 dr8 r8): -! gradient of useful results: *dw -! with respect to varying inputs: *dw *flowDoms.w -! rw status of diff variables: *dw:in-out *flowDoms.w:out -! plus diff mem management of: dw:in +! gradient of useful results: *(flowdoms.w) *dw +! with respect to varying inputs: *(flowdoms.w) *dw +! rw status of diff variables: *(flowdoms.w):incr *dw:in-out +! plus diff mem management of: flowdoms.w:in dw:in subroutine initres_block_fast_b(varstart, varend, nn, sps) ! ! initres initializes the given range of the residual. either to @@ -506,28 +506,157 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! local variables. ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m - real(kind=realtype) :: oneoverdt, tmp, tmpd + real(kind=realtype) :: oneoverdt, tmp + real(kind=realtype) :: tmpd real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1 real(kind=realtype), dimension(:, :, :), pointer :: volsp + real(kind=realtype) :: tempd + real(kind=realtype) :: tempd0 + real(kind=realtype) :: tempd1 integer :: branch ! return immediately of no variables are in the range. if (varend .ge. varstart) then ! determine the equation mode and act accordingly. - select case (equationmode) - case (steady) + select case (equationmode) + case (steady) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol3b(1) + call pushcontrol3b(3) else - call pushcontrol3b(0) + call pushcontrol3b(2) end if - case (timespectral) -! time spectral computation. + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. + jj = sectionid +! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol3b(2) +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) then + ll = 3*sps - 2 +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + end if + if (l .eq. ivy) then + ll = 3*sps - 1 +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + end if + if (l .eq. ivz) then + ll = 3*sps +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + end if +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + tmp = dvector(jj, ll, ii+1)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx) + dvector(jj, & +& ll, ii+2)*flowdoms(nn, currentlevel, mm)%w(i, j& +& , k, ivy) + dvector(jj, ll, ii+3)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + end do + end do + end do +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + end if + end do varloopfine + end do timeloopfine + call pushcontrol3b(1) else - call pushcontrol3b(3) +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) then + ll = 3*sps - 2 +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + end if + if (l .eq. ivy) then + ll = 3*sps - 1 +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + end if + if (l .eq. ivz) then + ll = 3*sps +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + end if +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 1 + else +myIntPtr = myIntPtr + 1 + myIntStack(myIntPtr) = 0 + end if + end do varloopcoarse + end do timeloopcoarse + call pushcontrol3b(0) end if case default call pushcontrol3b(4) @@ -559,76 +688,142 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do call popcontrol3b(branch) - if (branch .eq. 0) then -! steady, coarse grid - do l=varend,varstart,-1 - do k=kl,2,-1 - do j=jl,2,-1 - do i=il,2,-1 - dwd(i, j, k, l) = 0.0_8 - end do + if (branch .lt. 2) then + if (branch .eq. 0) then + do mm=ntimeintervalsspectral,1,-1 + do l=varend,varstart,-1 +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) then + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & +& dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & +& mm)%vol(i, j, k)*dwd(i, j, k, l) + end do + end do + end do + else + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + tmpd = flowdoms(nn, currentlevel, mm)%vol(i, j, k)& +& *dwd(i, j, k, l) + tempd = dvector(jj, ll, ii+1)*tmpd + tempd0 = dvector(jj, ll, ii+2)*tmpd + tempd1 = dvector(jj, ll, ii+3)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivz)*tempd1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivy)*tempd0 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd0 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& ivx)*tempd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & +& + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& +& )*tempd + end do + end do + end do +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .ne. 0) call popinteger4(ll) +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) call popinteger4(ll) +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) call popinteger4(ll) + end if end do end do - end do - else if (branch .eq. 1) then -! steady, fine grid - do l=varend,varstart,-1 - do k=kl,2,-1 - do j=jl,2,-1 - do i=il,2,-1 - dwd(i, j, k, l) = 0.0_8 + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do end do end do end do - end do - else if (branch .eq. 2) then -! time spectral, fine grid (state-only: vol, dscalar, dvector frozen) - jj = sectionid - do mm = ntimeintervalsspectral, 1, -1 - ii = 3 * (mm - 1) - do l = varend, varstart, -1 - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 -! recompute tmp - tmp = dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) -! adjoint of: dw += tmp * vol_mm * rho_mm (vol frozen) - tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * dwd(i,j,k,l) - flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - + tmp * flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) -! adjoint of: tmp = dvector * velocities - flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+1) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+2) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & - + dvector(jj, ll, ii+3) * tmpd + else + do mm=ntimeintervalsspectral,1,-1 + do l=varend,varstart,-1 +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) then + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & +& dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & +& mm)%vol(i, j, k)*dwd(i, j, k, l) + end do end do end do - end do - else -! scalar variable (vol and dscalar frozen) - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 - flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * dwd(i,j,k,l) + else + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + tempd = flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*dwd(i, j, k, l) + tmpd = flowdoms(nn, currentlevel, mm)%w(i, j, k, & +& irho)*tempd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& +& ) + tmp*tempd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & +& + dvector(jj, ll, ii+1)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & +& + dvector(jj, ll, ii+2)*tmpd + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & +& + dvector(jj, ll, ii+3)*tmpd + end do end do end do +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .ne. 0) call popinteger4(ll) +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) call popinteger4(ll) +branch = myIntStack(myIntPtr) + myIntPtr = myIntPtr - 1 + if (branch .eq. 0) call popinteger4(ll) + end if + end do + end do + do l=varend,varstart,-1 + do k=kl,2,-1 + do j=jl,2,-1 + do i=il,2,-1 + dwd(i, j, k, l) = 0.0_8 + end do end do - end if + end do end do - end do -! reverse of initialization: dw = zero => dwd = 0 + end if + else if (branch .eq. 2) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -639,66 +834,6 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do else if (branch .eq. 3) then -! time spectral, coarse grid (state-only: vol, dscalar, dvector, w1 frozen) - jj = sectionid - do mm = ntimeintervalsspectral, 1, -1 - ii = 3 * (mm - 1) - do l = varend, varstart, -1 - if (l .eq. ivx .or. l .eq. ivy .or. l .eq. ivz) then -! momentum variable - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 -! recompute tmp - tmp = dvector(jj, ll, ii+1) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivx)) & - + dvector(jj, ll, ii+2) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivy)) & - + dvector(jj, ll, ii+3) & - * (flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz) & - - flowdoms(nn, currentlevel, mm)%w1(i,j,k,irho) & - * flowdoms(nn, currentlevel, mm)%w1(i,j,k,ivz)) -! adjoint of: dw += tmp * vol_mm (vol frozen) - tmpd = flowdoms(nn, currentlevel, mm)%vol(i,j,k) * dwd(i,j,k,l) -! adjoint of tmp (w1 frozen) - flowdomsd(nn, 1, mm)%w(i,j,k,irho) = flowdomsd(nn, 1, mm)%w(i,j,k,irho) & - + (dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,ivz)) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivx) = flowdomsd(nn, 1, mm)%w(i,j,k,ivx) & - + dvector(jj, ll, ii+1) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivy) = flowdomsd(nn, 1, mm)%w(i,j,k,ivy) & - + dvector(jj, ll, ii+2) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - flowdomsd(nn, 1, mm)%w(i,j,k,ivz) = flowdomsd(nn, 1, mm)%w(i,j,k,ivz) & - + dvector(jj, ll, ii+3) * flowdoms(nn, currentlevel, mm)%w(i,j,k,irho) * tmpd - end do - end do - end do - else -! scalar variable (vol, dscalar, w1 frozen) - do k = kl, 2, -1 - do j = jl, 2, -1 - do i = il, 2, -1 - flowdomsd(nn, 1, mm)%w(i,j,k,l) = flowdomsd(nn, 1, mm)%w(i,j,k,l) & - + dscalar(jj, sps, mm) * flowdoms(nn, currentlevel, mm)%vol(i,j,k) & - * dwd(i,j,k,l) - end do - end do - end do - end if - end do - end do -! reverse of initialization: dw = wr (constant) => dwd = 0 do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -772,6 +907,174 @@ subroutine initres_block(varstart, varend, nn, sps) end do end do end if + case (timespectral) +! time spectral computation. the time derivative of the +! current solution is given by a linear combination of +! all other solutions, i.e. a matrix vector product. +! first store the section to which this block belongs +! in jj. + jj = sectionid +! determine the currently active multigrid level. + if (currentlevel .eq. groundlevel) then +! finest multigrid level. the residual must be +! initialized to the time derivative. +! initialize it to zero. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = zero + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopfine:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the volume. +! also store in ii the offset needed for vector +! quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopfine:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! loop over the owned cell centers to add the +! contribution from wsp. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! velocity in tmp. + tmp = dvector(jj, ll, ii+1)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx) + dvector(jj, & +& ll, ii+2)*flowdoms(nn, currentlevel, mm)%w(i, j& +& , k, ivy) + dvector(jj, ll, ii+3)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivz) +! update the residual. note the +! multiplication with the density to obtain +! the correct time derivative for the +! momentum variable. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho) + end do + end do + end do + else +! scalar variable. loop over the owned cells to +! add the contribution of wsp to the time +! derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*flowdoms(nn, currentlevel, mm)%w(i, j, k, l) + end do + end do + end do + end if + end do varloopfine + end do timeloopfine + else +! coarse grid level. initialize the owned cells to the +! residual forcing term plus a correction for the +! multigrid treatment of the time derivative term. +! initialization to the residual forcing term. + do l=varstart,varend + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = wr(i, j, k, l) + end do + end do + end do + end do +! loop over the number of terms which contribute +! to the time derivative. +timeloopcoarse:do mm=1,ntimeintervalsspectral +! store the pointer for the variable to be used to +! compute the unsteady source term and the pointers +! for wsp1, the solution when entering this mg level +! and for the volume. +! furthermore store in ii the offset needed for +! vector quantities. + ii = 3*(mm-1) +! loop over the number of variables to be set. +varloopcoarse:do l=varstart,varend +! test for a momentum variable. + if ((l .eq. ivx .or. l .eq. ivy) .or. l .eq. ivz) then +! momentum variable. a special treatment is +! needed because it is a vector and the velocities +! are stored instead of the momentum. set the +! coefficient ll, which defines the row of the +! matrix used later on. + if (l .eq. ivx) ll = 3*sps - 2 + if (l .eq. ivy) ll = 3*sps - 1 + if (l .eq. ivz) ll = 3*sps +! add the contribution of wps to the correction +! of the time derivative. the difference between +! the current time derivative and the one when +! entering this grid level must be added, because +! the residual forcing term only takes the spatial +! part of the coarse grid residual into account. + do k=2,kl + do j=2,jl + do i=2,il +! store the matrix vector product with the +! momentum in tmp. + tmp = dvector(jj, ll, ii+1)*(flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, irho)*flowdoms(nn, & +& currentlevel, mm)%w(i, j, k, ivx)-flowdoms(nn, & +& currentlevel, mm)%w1(i, j, k, irho)*flowdoms(nn& +& , currentlevel, mm)%w1(i, j, k, ivx)) + dvector(& +& jj, ll, ii+2)*(flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w(& +& i, j, k, ivy)-flowdoms(nn, currentlevel, mm)%w1(& +& i, j, k, irho)*flowdoms(nn, currentlevel, mm)%w1& +& (i, j, k, ivy)) + dvector(jj, ll, ii+3)*(& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, irho)*& +& flowdoms(nn, currentlevel, mm)%w(i, j, k, ivz)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, irho)& +& *flowdoms(nn, currentlevel, mm)%w1(i, j, k, ivz)& +& ) +! add tmp to the residual. multiply it by +! the volume to obtain the finite volume +! formulation of the derivative of the +! momentum. + dw(i, j, k, l) = dw(i, j, k, l) + tmp*flowdoms(nn& +& , currentlevel, mm)%vol(i, j, k) + end do + end do + end do + else +! scalar variable. loop over the owned cells +! to add the contribution of wsp to the correction +! of the time derivative. + do k=2,kl + do j=2,jl + do i=2,il + dw(i, j, k, l) = dw(i, j, k, l) + dscalar(jj, sps& +& , mm)*flowdoms(nn, currentlevel, mm)%vol(i, j, k& +& )*(flowdoms(nn, currentlevel, mm)%w(i, j, k, l)-& +& flowdoms(nn, currentlevel, mm)%w1(i, j, k, l)) + end do + end do + end do + end if + end do varloopcoarse + end do timeloopcoarse + end if end select ! set the residual in the halo cells to zero. this is just ! to avoid possible problems. their values do not matter. From f4f35c86879acab00278f274bc49ddc941d5ef75 Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Mon, 9 Mar 2026 12:25:17 -0400 Subject: [PATCH 05/12] add test --- ...djoint_euler_time_spectral_naca64A010.json | 280 ++++++++++++++++++ .../test_adjoint_time_spectral_naca64A010.py | 194 ++++++++++++ 2 files changed, 474 insertions(+) create mode 100644 tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json create mode 100644 tests/reg_tests/test_adjoint_time_spectral_naca64A010.py diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json new file mode 100644 index 000000000..4077b5aeb --- /dev/null +++ b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json @@ -0,0 +1,280 @@ +{ + "Adjoint States": 0.0, + "Eval Functions Sens:": { + "64A010pitchingTS_cd": {}, + "64A010pitchingTS_cl": {}, + "64A010pitchingTS_cmz": {} + }, + "Eval Functions:": { + "64A010pitchingTS_cd": 0.002500060868193104, + "64A010pitchingTS_cl": 0.0037769183216176568, + "64A010pitchingTS_cmz": 0.0011186213900658232 + }, + "metadata": { + "ADPC": false, + "ANKADPC": false, + "ANKAMGLevels": 2, + "ANKAMGNSmooth": 1, + "ANKASMOverlap": 1, + "ANKASMOverlapCoarse": 0, + "ANKCFL0": 5.0, + "ANKCFLCutback": 0.5, + "ANKCFLExponent": 0.5, + "ANKCFLFactor": 10.0, + "ANKCFLLimit": 100000.0, + "ANKCFLMin": 1.0, + "ANKCFLReset": true, + "ANKCharTimeStepType": "None", + "ANKConstCFLStep": 0.4, + "ANKCoupledSwitchTol": 1e-16, + "ANKGlobalPreconditioner": "additive Schwarz", + "ANKInnerPreconIts": 1, + "ANKInnerPreconItsCoarse": 1, + "ANKJacobianLag": 10, + "ANKLinResMax": 0.1, + "ANKLinearSolveBuffer": 0.01, + "ANKLinearSolveTol": 0.05, + "ANKMaxIter": 40, + "ANKNSubiterTurb": 1, + "ANKOuterPreconIts": 1, + "ANKPCILUFill": 2, + "ANKPCILUFillCoarse": 0, + "ANKPCUpdateCutoff": 1e-16, + "ANKPCUpdateTol": 0.5, + "ANKPCUpdateTolAfterCutoff": 0.0001, + "ANKPhysicalLSTol": 0.2, + "ANKPhysicalLSTolTurb": 0.99, + "ANKSecondOrdSwitchTol": 1e-16, + "ANKStepFactor": 1.0, + "ANKStepMin": 0.01, + "ANKSubspaceSize": 50, + "ANKSwitchTol": 0.01, + "ANKTurbCFLScale": 1.0, + "ANKTurbKSPDebug": false, + "ANKUnsteadyLSTol": 1.0, + "ANKUseApproxSA": false, + "ANKUseFullVisc": true, + "ANKUseMatrixFree": true, + "ANKUseTurbDADI": true, + "ASMOverlap": 1, + "ASMOverlapCoarse": 0, + "CFL": 1.7, + "CFLCoarse": 1.0, + "CFLLimit": 1.5, + "GMRESOrthogonalizationType": "modified Gram-Schmidt", + "ILUFill": 2, + "ILUFillCoarse": 0, + "L2Convergence": 1e-13, + "L2ConvergenceCoarse": 0.01, + "L2ConvergenceRel": 1e-16, + "MGCycle": "sg", + "MGStartLevel": -1, + "NKADPC": false, + "NKAMGLevels": 2, + "NKAMGNSmooth": 1, + "NKASMOverlap": 1, + "NKASMOverlapCoarse": 0, + "NKFixedStep": 0.25, + "NKGlobalPreconditioner": "additive Schwarz", + "NKInnerPreconIts": 1, + "NKInnerPreconItsCoarse": 1, + "NKJacobianLag": 20, + "NKLS": "cubic", + "NKLinearSolveTol": 0.3, + "NKOuterPreconIts": 1, + "NKPCILUFill": 2, + "NKPCILUFillCoarse": 0, + "NKSubspaceSize": 400, + "NKSwitchTol": 0.0001, + "NKUseEW": true, + "NKViscPC": false, + "RKReset": false, + "TSStability": false, + "acousticScaleFactor": 1.0, + "adjointAMGLevels": 2, + "adjointAMGNSmooth": 1, + "adjointDivTol": 100000.0, + "adjointL2Convergence": 1e-14, + "adjointL2ConvergenceAbs": 1e-16, + "adjointL2ConvergenceRel": 1e-16, + "adjointMaxIter": 500, + "adjointMaxL2DeviationFactor": 1.0, + "adjointMonitorStep": 10, + "adjointSolver": "GMRES", + "adjointSubspaceSize": 100, + "alphaFollowing": false, + "alphaMode": false, + "altitudeMode": false, + "applyAdjointPCSubspaceSize": 20, + "applyPCSubspaceSize": 400, + "approxPC": true, + "backgroundVolScale": 1.0, + "betaMode": false, + "blockSplitting": true, + "cavExponent": 0, + "cavSensorOffset": 0.0, + "cavSensorSharpness": 10.0, + "cavitationNumber": 1.4, + "closedSurfaceFamilies": null, + "coarseDiscretization": "central plus scalar dissipation", + "computeCavitation": false, + "computeSepSensorKs": false, + "coupledSolution": false, + "cpMinRho": 100.0, + "cutCallback": null, + "debugZipper": false, + "deltaT": 0.01, + "designSurfaceFamily": null, + "discretization": "central plus scalar dissipation", + "dissContMagnitude": 1.0, + "dissContMidpoint": 3.0, + "dissContSharpness": 3.0, + "dissipationLumpingParameter": 6.0, + "dissipationScalingExponent": 0.67, + "eddyVisInfRatio": 0.009, + "equationMode": "time spectral", + "equationType": "Euler", + "eulerWallTreatment": "linear pressure extrapolation", + "explicitSurfaceCallback": null, + "firstRun": true, + "flowType": "external", + "forcesAsTractions": true, + "frozenTurbulence": false, + "globalPreconditioner": "additive Schwarz", + "gridFile": "tests/tests/tests/input_files/naca64A010_euler-L2.cgns", + "gridPrecision": "double", + "gridPrecisionSurface": "single", + "infChangeCorrection": true, + "infChangeCorrectionTol": 1e-12, + "infChangeCorrectionType": "offset", + "innerPreconIts": 1, + "innerPreconItsCoarse": 1, + "isoVariables": [], + "isosurface": {}, + "liftIndex": 2, + "limiter": "van Albada", + "loadBalanceIter": 10, + "loadImbalance": 0.1, + "localPreconditioner": "ILU", + "lowSpeedPreconditioner": false, + "machMode": false, + "matrixOrdering": "RCM", + "maxL2DeviationFactor": 1.0, + "meshMaxSkewness": 1.0, + "meshSurfaceFamily": null, + "monitorVariables": [ + "resrho", + "cl" + ], + "nCycles": 200000, + "nCyclesCoarse": 500, + "nFloodIter": -1, + "nRKReset": 5, + "nRefine": 10, + "nSaveSurface": 1, + "nSaveVolume": 1, + "nSubiter": 1, + "nSubiterTurb": 3, + "nTimeStepsCoarse": 48, + "nTimeStepsFine": 400, + "nearWallDist": 0.1, + "numberSolutions": true, + "outerPreconIts": 3, + "outputDirectory": "tests/tests/tests/tests/output_files", + "outputSurfaceFamily": "allSurfaces", + "overlapFactor": 0.9, + "oversetDebugPrint": false, + "oversetLoadBalance": true, + "oversetPriority": {}, + "oversetProjTol": 1e-12, + "oversetUpdateMode": "frozen", + "pMode": false, + "partitionLikeNProc": -1, + "partitionOnly": false, + "preconditionerSide": "right", + "printAllOptions": true, + "printBadlySkewedCells": false, + "printIntro": true, + "printIterations": true, + "printNegativeVolumes": false, + "printTiming": true, + "printWarnings": true, + "qMode": false, + "rMode": false, + "recomputeOverlapMatrix": true, + "resAveraging": "alternate", + "restartAdjoint": true, + "restartFile": null, + "restrictionRelaxation": 0.8, + "selfZipCutoff": 120.0, + "sepSensorKsOffset": 0.0, + "sepSensorKsPhi": 90.0, + "sepSensorKsRho": 1000.0, + "sepSensorKsSharpness": 25.0, + "sepSensorOffset": 0.0, + "sepSensorSharpness": 10.0, + "setMonitor": true, + "skipAfterFailedAdjoint": true, + "smoothParameter": 1.5, + "smoother": "Runge-Kutta", + "solutionPrecision": "single", + "solutionPrecisionSurface": "single", + "storeConvHist": true, + "storeRindLayer": true, + "surfaceVariables": [ + "cp", + "vx", + "vy", + "vz", + "mach" + ], + "timeAccuracy": 2, + "timeIntegrationScheme": "BDF", + "timeIntervals": 3, + "timeLimit": -1.0, + "turbResScale": 10000.0, + "turbulenceModel": "SA", + "turbulenceOrder": "first order", + "turbulenceProduction": "strain", + "updateWallAssociations": false, + "useALE": true, + "useANKSolver": true, + "useApproxWallDistance": true, + "useBlockettes": false, + "useDiagTSPC": true, + "useDissContinuation": false, + "useExternalDynamicMesh": true, + "useGridMotion": false, + "useLinResMonitor": false, + "useMatrixFreedrdw": true, + "useNKSolver": true, + "useOversetWallScaling": false, + "useQCR": false, + "useRotationSA": false, + "useSkewnessCheck": false, + "useTSInterpolatedGridVelocity": true, + "useWallFunctions": false, + "useZipperMesh": true, + "useft2SA": true, + "verifyExtra": true, + "verifySpatial": true, + "verifyState": true, + "vis2": 0.25, + "vis2Coarse": 0.5, + "vis4": 0.0156, + "viscPC": false, + "viscWallTreatment": "constant pressure extrapolation", + "viscousSurfaceVelocities": true, + "volumeVariables": [ + "resrho" + ], + "wallDistCutoff": 1e+20, + "windAxis": false, + "writeSolutionDigits": 3, + "writeSolutionEachIter": false, + "writeSurfaceSolution": false, + "writeTecplotSurfaceSolution": false, + "writeVolumeSolution": false, + "zipperSurfaceFamily": null + } +} \ No newline at end of file diff --git a/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py b/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py new file mode 100644 index 000000000..8785a1b49 --- /dev/null +++ b/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py @@ -0,0 +1,194 @@ +# built-ins +import unittest +import os +import copy +import numpy + +# MACH testing class +from adflow import ADFLOW +from idwarp import USMesh + +# import the testing utilities +import reg_test_utils as utils + +from reg_default_options import adflowDefOpts + +from reg_aeroproblems import ap_naca64A010_time_spectral +import reg_test_classes + +baseDir = os.path.dirname(os.path.abspath(__file__)) + + +class TestAdjointTimeSpectral(reg_test_classes.RegTest): + """ + Tests for time spectral adjoint on the NACA 64A010 case. + """ + + N_PROCS = 2 + ref_file = "adjoint_euler_time_spectral_naca64A010.json" + + def setUp(self): + # Introduce a transfer class for displacement transfer from struct to aero. + class Transfer: + # simplified transfer class + # converting csd displacement to cfd surface nodes + + def __init__(self, alpha, xRot, aeroSolver): + # takes in displacement history + + self.alpha = alpha + self.ntimeintervalsspectral = len(alpha) + + self.aeroSolver = aeroSolver # a shallow copy of CFD solver + + self.xRot = xRot + + def getUndeformedSurfaceNodes(self): + self.MDGroup = self.aeroSolver.allWallsGroup + + self.cfdPts0 = self.aeroSolver.getSurfaceCoordinates( + self.MDGroup, includeZipper=False + ) + + def setDisplacements(self): + xRot = self.xRot + ntimeintervalsspectral = self.ntimeintervalsspectral + alpha = self.alpha + cfdPoints_init = self.cfdPts0 + + N_pts = cfdPoints_init.shape[0] + + self.cfdPts = [] + + for sps in range(ntimeintervalsspectral): + cfdPoints_deformed = numpy.zeros((N_pts, 3)) + + ptch_loc = alpha[sps] + + cc = numpy.cos(ptch_loc) + ss = numpy.sin(ptch_loc) + + for j in range(N_pts): + cfdPoints_deformed[j, 0] = ( + cc * (cfdPoints_init[j, 0] - xRot) + + ss * cfdPoints_init[j, 1] + + xRot + ) + cfdPoints_deformed[j, 1] = ( + -ss * (cfdPoints_init[j, 0] - xRot) + + cc * cfdPoints_init[j, 1] + ) + cfdPoints_deformed[j, 2] = cfdPoints_init[j, 2] + + self.cfdPts.append(cfdPoints_deformed) + + def setVolumeMesh(self): + ntimeintervalsspectral = self.ntimeintervalsspectral + + for sps in range(ntimeintervalsspectral): + self.aeroSolver.mesh.setSurfaceCoordinates(self.cfdPts[sps]) + self.aeroSolver.mesh.warpMesh() + m = self.aeroSolver.mesh.getSolverGrid() + self.aeroSolver.adflow.warping.setgridforoneinstance(m, sps=sps + 1) + + self.aeroSolver._updateGeomInfo = True + self.aeroSolver.updateGeometryInfo() + + # Set up the test for the parent RegTest class. + super().setUp() + + gridFile = os.path.join(baseDir, "../../input_files/naca64A010_euler-L2.cgns") + + ntimeintervalsspectral = 3 + options = copy.copy(adflowDefOpts) + options.update( + { + "gridfile": gridFile, + "outputDirectory": os.path.join(baseDir, "../output_files"), + "writeVolumeSolution": False, + "writeSurfaceSolution": False, + "blocksplitting": True, + "useblockettes": False, + "equationtype": "Euler", + "equationmode": "time spectral", + "mgcycle": "sg", + "l2convergence": 1e-13, + "ncycles": 200000, + "monitorvariables": ["resrho", "cl"], + "usenksolver": True, + "nkswitchtol": 1e-4, + "NKSubSpaceSize": 400, + "applypcsubspacesize": 400, + "useanksolver": True, + "ankswitchtol": 1e-2, + "anksubspacesize": 50, + "alphafollowing": False, + "timeintervals": ntimeintervalsspectral, + "useexternaldynamicmesh": True, + "usetsinterpolatedgridvelocity": True, + "adjointl2convergence": 1e-14, + } + ) + + # Grid option + meshOptions = { + "gridFile": gridFile, + } + + # Setup aeroproblem + self.ap = copy.copy(ap_naca64A010_time_spectral) + + # Motion history + alpha_0 = 1.01 + deltaAlpha = -alpha_0 * numpy.pi / 180.0 + alpha = numpy.linspace(0.0, 2.0 * numpy.pi, ntimeintervalsspectral + 1)[:-1] + alpha = -numpy.sin(alpha) + alpha *= deltaAlpha + + # Create the solver + self.CFDSolver = ADFLOW(options=options, debug=True) + + # Deform the mesh + mesh = USMesh(options=meshOptions) + self.CFDSolver.setMesh(mesh) + + # deformation + xRot = 0.25 # Hard copied from the reference file. + TSTransfer = Transfer(alpha, xRot, self.CFDSolver) + TSTransfer.getUndeformedSurfaceNodes() + TSTransfer.setDisplacements() + TSTransfer.setVolumeMesh() + + # Solve + self.CFDSolver(self.ap) + + # Check solution + self.assert_solution_failure() + + # Initialize residual for adjoint tests + self.CFDSolver.getResidual(self.ap) + + def test_functions(self): + utils.assert_functions_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-8) + + def test_adjoint(self): + utils.assert_adjoint_sens_allclose( + self.handler, self.CFDSolver, self.ap, tol=1e-10 + ) + self.assert_adjoint_failure() + + def test_adjoint2(self): + utils.assert_adjoint2_sens_allclose( + self.handler, self.CFDSolver, self.ap, tol=1e-10 + ) + self.assert_adjoint_failure() + + def test_adjoint_states(self): + utils.assert_adjoint_states_allclose( + self.handler, self.CFDSolver, self.ap, tol=1e-10 + ) + self.assert_adjoint_failure() + + +if __name__ == "__main__": + unittest.main() From d91f6ca2d82dbf710d09b08847955f35068fc1c9 Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Mon, 9 Mar 2026 12:30:12 -0400 Subject: [PATCH 06/12] reformat --- ...djoint_euler_time_spectral_naca64A010.json | 2 +- .../test_adjoint_time_spectral_naca64A010.py | 27 +++++-------------- 2 files changed, 7 insertions(+), 22 deletions(-) diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json index 4077b5aeb..4a9c6ec3c 100644 --- a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json +++ b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json @@ -277,4 +277,4 @@ "writeVolumeSolution": false, "zipperSurfaceFamily": null } -} \ No newline at end of file +} diff --git a/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py b/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py index 8785a1b49..7311e41c6 100644 --- a/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py +++ b/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py @@ -46,9 +46,7 @@ def __init__(self, alpha, xRot, aeroSolver): def getUndeformedSurfaceNodes(self): self.MDGroup = self.aeroSolver.allWallsGroup - self.cfdPts0 = self.aeroSolver.getSurfaceCoordinates( - self.MDGroup, includeZipper=False - ) + self.cfdPts0 = self.aeroSolver.getSurfaceCoordinates(self.MDGroup, includeZipper=False) def setDisplacements(self): xRot = self.xRot @@ -69,15 +67,8 @@ def setDisplacements(self): ss = numpy.sin(ptch_loc) for j in range(N_pts): - cfdPoints_deformed[j, 0] = ( - cc * (cfdPoints_init[j, 0] - xRot) - + ss * cfdPoints_init[j, 1] - + xRot - ) - cfdPoints_deformed[j, 1] = ( - -ss * (cfdPoints_init[j, 0] - xRot) - + cc * cfdPoints_init[j, 1] - ) + cfdPoints_deformed[j, 0] = cc * (cfdPoints_init[j, 0] - xRot) + ss * cfdPoints_init[j, 1] + xRot + cfdPoints_deformed[j, 1] = -ss * (cfdPoints_init[j, 0] - xRot) + cc * cfdPoints_init[j, 1] cfdPoints_deformed[j, 2] = cfdPoints_init[j, 2] self.cfdPts.append(cfdPoints_deformed) @@ -172,21 +163,15 @@ def test_functions(self): utils.assert_functions_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-8) def test_adjoint(self): - utils.assert_adjoint_sens_allclose( - self.handler, self.CFDSolver, self.ap, tol=1e-10 - ) + utils.assert_adjoint_sens_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-10) self.assert_adjoint_failure() def test_adjoint2(self): - utils.assert_adjoint2_sens_allclose( - self.handler, self.CFDSolver, self.ap, tol=1e-10 - ) + utils.assert_adjoint2_sens_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-10) self.assert_adjoint_failure() def test_adjoint_states(self): - utils.assert_adjoint_states_allclose( - self.handler, self.CFDSolver, self.ap, tol=1e-10 - ) + utils.assert_adjoint_states_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-10) self.assert_adjoint_failure() From 6f355aea5274165128ca89220ba3752ba8c7002c Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Mon, 9 Mar 2026 13:36:53 -0400 Subject: [PATCH 07/12] more --- .../outputReverseFast/residuals_fast_b.f90 | 48 +++++++++---------- ...djoint_euler_time_spectral_naca64A010.json | 2 +- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index fca35d480..221a45f49 100644 --- a/src/adjoint/outputReverseFast/residuals_fast_b.f90 +++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90 @@ -698,8 +698,8 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) do k=kl,2,-1 do j=jl,2,-1 do i=il,2,-1 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & + flowdomsd(nn, 1, mm)%w(i, j, k, l) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, l) + & & dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & & mm)%vol(i, j, k)*dwd(i, j, k, l) end do @@ -714,28 +714,28 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) tempd = dvector(jj, ll, ii+1)*tmpd tempd0 = dvector(jj, ll, ii+2)*tmpd tempd1 = dvector(jj, ll, ii+3)*tmpd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& + flowdomsd(nn, 1, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivz)*tempd1 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivz) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd1 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& + flowdomsd(nn, 1, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivy)*tempd0 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivy) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd0 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& + flowdomsd(nn, 1, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivx)*tempd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivx) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd end do @@ -771,8 +771,8 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) do k=kl,2,-1 do j=jl,2,-1 do i=il,2,-1 - flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & + flowdomsd(nn, 1, mm)%w(i, j, k, l) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, l) + & & dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & & mm)%vol(i, j, k)*dwd(i, j, k, l) end do @@ -786,17 +786,17 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) & )*dwd(i, j, k, l) tmpd = flowdoms(nn, currentlevel, mm)%w(i, j, k, & & irho)*tempd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& + flowdomsd(nn, 1, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& & ) + tmp*tempd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivx) & & + dvector(jj, ll, ii+1)*tmpd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivy) & & + dvector(jj, ll, ii+2)*tmpd - flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & -& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & + flowdomsd(nn, 1, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, 1, mm)%w(i, j, k, ivz) & & + dvector(jj, ll, ii+3)*tmpd end do end do diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json index 4a9c6ec3c..4077b5aeb 100644 --- a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json +++ b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json @@ -277,4 +277,4 @@ "writeVolumeSolution": false, "zipperSurfaceFamily": null } -} +} \ No newline at end of file From d0cb10e2a488baacbb8db5fd28f37b1f211d7a7d Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Tue, 10 Mar 2026 11:56:03 -0400 Subject: [PATCH 08/12] more --- src/adjoint/autoEdit/autoEditReverseFast.py | 5 +++++ .../refs/adjoint_euler_time_spectral_naca64A010.json | 3 ++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/adjoint/autoEdit/autoEditReverseFast.py b/src/adjoint/autoEdit/autoEditReverseFast.py index 5df9ddb9a..87c630da9 100644 --- a/src/adjoint/autoEdit/autoEditReverseFast.py +++ b/src/adjoint/autoEdit/autoEditReverseFast.py @@ -141,6 +141,11 @@ if m: line = line.replace("integer*4", "integer") + # flowdomsd is only allocated for 1 multigrid level, so we must + # replace flowdomsd(nn, currentlevel, mm) with flowdomsd(nn, 1, mm) + if "flowdomsd" in line and "currentlevel" in line: + line = line.replace("flowdomsd(nn, currentlevel, mm)", "flowdomsd(nn, 1, mm)") + # # See if we need to modify the line # m = patt_module_start.match(line) # if m: diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json index 4077b5aeb..18bdb658a 100644 --- a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json +++ b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json @@ -277,4 +277,5 @@ "writeVolumeSolution": false, "zipperSurfaceFamily": null } -} \ No newline at end of file +} + From 468ee51c3770034a0f286e4d3ca25bf373cacf4b Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Tue, 10 Mar 2026 14:17:24 -0400 Subject: [PATCH 09/12] revert the autoedit change --- src/adjoint/autoEdit/autoEditReverseFast.py | 5 -- .../outputReverseFast/residuals_fast_b.f90 | 48 +++++++++---------- 2 files changed, 24 insertions(+), 29 deletions(-) diff --git a/src/adjoint/autoEdit/autoEditReverseFast.py b/src/adjoint/autoEdit/autoEditReverseFast.py index 87c630da9..5df9ddb9a 100644 --- a/src/adjoint/autoEdit/autoEditReverseFast.py +++ b/src/adjoint/autoEdit/autoEditReverseFast.py @@ -141,11 +141,6 @@ if m: line = line.replace("integer*4", "integer") - # flowdomsd is only allocated for 1 multigrid level, so we must - # replace flowdomsd(nn, currentlevel, mm) with flowdomsd(nn, 1, mm) - if "flowdomsd" in line and "currentlevel" in line: - line = line.replace("flowdomsd(nn, currentlevel, mm)", "flowdomsd(nn, 1, mm)") - # # See if we need to modify the line # m = patt_module_start.match(line) # if m: diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index 221a45f49..fca35d480 100644 --- a/src/adjoint/outputReverseFast/residuals_fast_b.f90 +++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90 @@ -698,8 +698,8 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) do k=kl,2,-1 do j=jl,2,-1 do i=il,2,-1 - flowdomsd(nn, 1, mm)%w(i, j, k, l) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, l) + & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & & dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & & mm)%vol(i, j, k)*dwd(i, j, k, l) end do @@ -714,28 +714,28 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) tempd = dvector(jj, ll, ii+1)*tmpd tempd0 = dvector(jj, ll, ii+2)*tmpd tempd1 = dvector(jj, ll, ii+3)*tmpd - flowdomsd(nn, 1, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivz)*tempd1 - flowdomsd(nn, 1, mm)%w(i, j, k, ivz) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivz) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd1 - flowdomsd(nn, 1, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivy)*tempd0 - flowdomsd(nn, 1, mm)%w(i, j, k, ivy) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivy) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd0 - flowdomsd(nn, 1, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& & ) + flowdoms(nn, currentlevel, mm)%w(i, j, k, & & ivx)*tempd - flowdomsd(nn, 1, mm)%w(i, j, k, ivx) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivx) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & & + flowdoms(nn, currentlevel, mm)%w(i, j, k, irho& & )*tempd end do @@ -771,8 +771,8 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) do k=kl,2,-1 do j=jl,2,-1 do i=il,2,-1 - flowdomsd(nn, 1, mm)%w(i, j, k, l) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, l) + & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, l) + & & dscalar(jj, sps, mm)*flowdoms(nn, currentlevel, & & mm)%vol(i, j, k)*dwd(i, j, k, l) end do @@ -786,17 +786,17 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) & )*dwd(i, j, k, l) tmpd = flowdoms(nn, currentlevel, mm)%w(i, j, k, & & irho)*tempd - flowdomsd(nn, 1, mm)%w(i, j, k, irho)& -& = flowdomsd(nn, 1, mm)%w(i, j, k, irho& + flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho)& +& = flowdomsd(nn, currentlevel, mm)%w(i, j, k, irho& & ) + tmp*tempd - flowdomsd(nn, 1, mm)%w(i, j, k, ivx) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivx) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivx) & & + dvector(jj, ll, ii+1)*tmpd - flowdomsd(nn, 1, mm)%w(i, j, k, ivy) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivy) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivy) & & + dvector(jj, ll, ii+2)*tmpd - flowdomsd(nn, 1, mm)%w(i, j, k, ivz) = & -& flowdomsd(nn, 1, mm)%w(i, j, k, ivz) & + flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) = & +& flowdomsd(nn, currentlevel, mm)%w(i, j, k, ivz) & & + dvector(jj, ll, ii+3)*tmpd end do end do From 615d8b0142f9ff53f6a896068b9bf1347e5e7e9f Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Wed, 11 Mar 2026 13:00:58 -0400 Subject: [PATCH 10/12] fix mem overflow issue --- src/modules/constants.F90 | 2 +- .../reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/modules/constants.F90 b/src/modules/constants.F90 index 529cf488c..d4f4aa346 100644 --- a/src/modules/constants.F90 +++ b/src/modules/constants.F90 @@ -296,7 +296,7 @@ module constants integer(kind=intType), parameter :: kMin = 5 integer(kind=intType), parameter :: kMax = 6 - integer(kind=intType) :: myIntStack(32) + integer(kind=intType) :: myIntStack(200) integer(kind=intType) :: myIntPtr = 0 ! BC specific input variable counts diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json index 18bdb658a..4a9c6ec3c 100644 --- a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json +++ b/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json @@ -278,4 +278,3 @@ "zipperSurfaceFamily": null } } - From a758ff5c77ea83cd458af0a58e6e0d89ca8596ea Mon Sep 17 00:00:00 2001 From: Sicheng He Date: Thu, 12 Mar 2026 12:28:37 -0400 Subject: [PATCH 11/12] fix source code issue --- src/adjoint/outputForward/residuals_d.f90 | 16 +--- src/adjoint/outputReverse/residuals_b.f90 | 68 ++--------------- .../outputReverseFast/residuals_fast_b.f90 | 76 +------------------ src/solver/residuals.F90 | 8 +- 4 files changed, 18 insertions(+), 150 deletions(-) diff --git a/src/adjoint/outputForward/residuals_d.f90 b/src/adjoint/outputForward/residuals_d.f90 index 6588e8cac..f330082e0 100644 --- a/src/adjoint/outputForward/residuals_d.f90 +++ b/src/adjoint/outputForward/residuals_d.f90 @@ -625,9 +625,7 @@ subroutine initres_block_d(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -719,9 +717,7 @@ subroutine initres_block_d(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! add the contribution of wps to the correction ! of the time derivative. the difference between ! the current time derivative and the one when @@ -943,9 +939,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -1017,9 +1011,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! add the contribution of wps to the correction ! of the time derivative. the difference between ! the current time derivative and the one when diff --git a/src/adjoint/outputReverse/residuals_b.f90 b/src/adjoint/outputReverse/residuals_b.f90 index 8a831baae..f63db42f7 100644 --- a/src/adjoint/outputReverse/residuals_b.f90 +++ b/src/adjoint/outputReverse/residuals_b.f90 @@ -609,27 +609,8 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) then - call pushinteger4(ll) - ll = 3*sps - 2 - call pushcontrol1b(0) - else - call pushcontrol1b(1) - end if - if (l .eq. ivy) then - call pushinteger4(ll) - ll = 3*sps - 1 - call pushcontrol1b(0) - else - call pushcontrol1b(1) - end if - if (l .eq. ivz) then - call pushinteger4(ll) - ll = 3*sps - call pushcontrol1b(1) - else - call pushcontrol1b(0) - end if + call pushinteger4(ll) + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -678,27 +659,8 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) then - call pushinteger4(ll) - ll = 3*sps - 2 - call pushcontrol1b(0) - else - call pushcontrol1b(1) - end if - if (l .eq. ivy) then - call pushinteger4(ll) - ll = 3*sps - 1 - call pushcontrol1b(0) - else - call pushcontrol1b(1) - end if - if (l .eq. ivz) then - call pushinteger4(ll) - ll = 3*sps - call pushcontrol1b(1) - else - call pushcontrol1b(0) - end if + call pushinteger4(ll) + ll = 3*sps - 2 + (l-ivx) ! add the contribution of wps to the correction ! of the time derivative. the difference between ! the current time derivative and the one when @@ -833,12 +795,7 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do end do - call popcontrol1b(branch) - if (branch .ne. 0) call popinteger4(ll) - call popcontrol1b(branch) - if (branch .eq. 0) call popinteger4(ll) - call popcontrol1b(branch) - if (branch .eq. 0) call popinteger4(ll) + call popinteger4(ll) end if end do call popinteger4(ii) @@ -899,12 +856,7 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do end do - call popcontrol1b(branch) - if (branch .ne. 0) call popinteger4(ll) - call popcontrol1b(branch) - if (branch .eq. 0) call popinteger4(ll) - call popcontrol1b(branch) - if (branch .eq. 0) call popinteger4(ll) + call popinteger4(ll) end if end do call popinteger4(ii) @@ -1041,9 +993,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -1115,9 +1065,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! add the contribution of wps to the correction ! of the time derivative. the difference between ! the current time derivative and the one when diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index fca35d480..1b7da88dd 100644 --- a/src/adjoint/outputReverseFast/residuals_fast_b.f90 +++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90 @@ -552,30 +552,7 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) then - ll = 3*sps - 2 -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - end if - if (l .eq. ivy) then - ll = 3*sps - 1 -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - end if - if (l .eq. ivz) then - ll = 3*sps -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - end if + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -624,30 +601,7 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) then - ll = 3*sps - 2 -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - end if - if (l .eq. ivy) then - ll = 3*sps - 1 -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - end if - if (l .eq. ivz) then - ll = 3*sps -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 1 - else -myIntPtr = myIntPtr + 1 - myIntStack(myIntPtr) = 0 - end if + ll = 3*sps - 2 + (l-ivx) myIntPtr = myIntPtr + 1 myIntStack(myIntPtr) = 1 else @@ -741,15 +695,6 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do end do -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .ne. 0) call popinteger4(ll) -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .eq. 0) call popinteger4(ll) -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .eq. 0) call popinteger4(ll) end if end do end do @@ -801,15 +746,6 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do end do -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .ne. 0) call popinteger4(ll) -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .eq. 0) call popinteger4(ll) -branch = myIntStack(myIntPtr) - myIntPtr = myIntPtr - 1 - if (branch .eq. 0) call popinteger4(ll) end if end do end do @@ -945,9 +881,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! loop over the owned cell centers to add the ! contribution from wsp. do k=2,kl @@ -1019,9 +953,7 @@ subroutine initres_block(varstart, varend, nn, sps) ! are stored instead of the momentum. set the ! coefficient ll, which defines the row of the ! matrix used later on. - if (l .eq. ivx) ll = 3*sps - 2 - if (l .eq. ivy) ll = 3*sps - 1 - if (l .eq. ivz) ll = 3*sps + ll = 3*sps - 2 + (l-ivx) ! add the contribution of wps to the correction ! of the time derivative. the difference between ! the current time derivative and the one when diff --git a/src/solver/residuals.F90 b/src/solver/residuals.F90 index 73c976c84..5b56ccc77 100644 --- a/src/solver/residuals.F90 +++ b/src/solver/residuals.F90 @@ -757,9 +757,7 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! coefficient ll, which defines the row of the ! matrix used later on. - if (l == ivx) ll = 3 * sps - 2 - if (l == ivy) ll = 3 * sps - 1 - if (l == ivz) ll = 3 * sps + ll = 3 * sps - 2 + (l - ivx) ! Loop over the owned cell centers to add the ! contribution from wsp. @@ -879,9 +877,7 @@ subroutine initres_block(varStart, varEnd, nn, sps) ! coefficient ll, which defines the row of the ! matrix used later on. - if (l == ivx) ll = 3 * sps - 2 - if (l == ivy) ll = 3 * sps - 1 - if (l == ivz) ll = 3 * sps + ll = 3 * sps - 2 + (l - ivx) ! Add the contribution of wps to the correction ! of the time derivative. The difference between From e4078103a024e8a2f58b60a8afc8f5d32080de22 Mon Sep 17 00:00:00 2001 From: Alasdair Gray Date: Thu, 16 Apr 2026 17:18:16 -0400 Subject: [PATCH 12/12] Combine time spectral tests --- ...on => euler_time_spectral_naca64A010.json} | 0 .../solve_euler_time_spectral_naca64A010.json | 222 ------------------ ...al_naca64A010.py => test_time_spectral.py} | 10 +- .../test_time_spectral_naca64A010.py | 166 ------------- 4 files changed, 5 insertions(+), 393 deletions(-) rename tests/reg_tests/refs/{adjoint_euler_time_spectral_naca64A010.json => euler_time_spectral_naca64A010.json} (100%) delete mode 100644 tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json rename tests/reg_tests/{test_adjoint_time_spectral_naca64A010.py => test_time_spectral.py} (97%) delete mode 100644 tests/reg_tests/test_time_spectral_naca64A010.py diff --git a/tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/euler_time_spectral_naca64A010.json similarity index 100% rename from tests/reg_tests/refs/adjoint_euler_time_spectral_naca64A010.json rename to tests/reg_tests/refs/euler_time_spectral_naca64A010.json diff --git a/tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json deleted file mode 100644 index e89f229eb..000000000 --- a/tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json +++ /dev/null @@ -1,222 +0,0 @@ -{ - "Eval Functions:": { - "64A010pitchingTS_cd": 0.0025000608682108497, - "64A010pitchingTS_cl": 0.0037769183221013254, - "64A010pitchingTS_cmz": 0.0011186213899540035 - }, - "metadata": { - "ADPC": false, - "AGMGLevels": 1, - "AGMGNSmooth": 3, - "ANKADPC": false, - "ANKASMOverlap": 1, - "ANKCFL0": 5.0, - "ANKCFLCutback": 0.5, - "ANKCFLExponent": 0.5, - "ANKCFLFactor": 10.0, - "ANKCFLLimit": 100000.0, - "ANKCFLMin": 1.0, - "ANKConstCFLStep": 0.4, - "ANKCoupledSwitchTol": 1e-16, - "ANKInnerPreconIts": 1, - "ANKJacobianLag": 10, - "ANKLinResMax": 0.1, - "ANKLinearSolveTol": 0.05, - "ANKMaxIter": 40, - "ANKNSubiterTurb": 1, - "ANKOuterPreconIts": 1, - "ANKPCILUFill": 2, - "ANKPCUpdateTol": 0.5, - "ANKPhysicalLSTol": 0.2, - "ANKPhysicalLSTolTurb": 0.99, - "ANKSecondOrdSwitchTol": 1e-16, - "ANKStepFactor": 1.0, - "ANKStepMin": 0.01, - "ANKSubspaceSize": 50, - "ANKSwitchTol": 0.01, - "ANKTurbCFLScale": 1.0, - "ANKTurbKSPDebug": false, - "ANKUnsteadyLSTol": 1.0, - "ANKUseFullVisc": true, - "ANKUseMatrixFree": true, - "ANKUseTurbDADI": true, - "ASMOverlap": 1, - "CFL": 1.7, - "CFLCoarse": 1.0, - "CFLLimit": 1.5, - "ILUFill": 2, - "L2Convergence": 1e-15, - "L2ConvergenceCoarse": 0.01, - "L2ConvergenceRel": 1e-16, - "MGCycle": "sg", - "MGStartLevel": -1, - "NKADPC": false, - "NKASMOverlap": 1, - "NKFixedStep": 0.25, - "NKInnerPreconIts": 1, - "NKJacobianLag": 20, - "NKLS": "cubic", - "NKLinearSolveTol": 0.3, - "NKOuterPreconIts": 1, - "NKPCILUFill": 2, - "NKSubspaceSize": 400, - "NKSwitchTol": 0.0001, - "NKUseEW": true, - "NKViscPC": false, - "RKReset": false, - "TSStability": false, - "adjointDivTol": 100000.0, - "adjointL2Convergence": 1e-06, - "adjointL2ConvergenceAbs": 1e-16, - "adjointL2ConvergenceRel": 1e-16, - "adjointMaxIter": 500, - "adjointMonitorStep": 10, - "adjointSolver": "GMRES", - "adjointSubspaceSize": 100, - "alphaFollowing": false, - "alphaMode": false, - "altitudeMode": false, - "applyAdjointPCSubspaceSize": 20, - "applyPCSubspaceSize": 400, - "approxPC": true, - "backgroundVolScale": 1.0, - "betaMode": false, - "blockSplitting": true, - "cavitationNumber": 1.4, - "closedSurfaceFamilies": null, - "coarseDiscretization": "central plus scalar dissipation", - "computeCavitation": false, - "coupledSolution": false, - "cutCallback": null, - "debugZipper": false, - "deltaT": 0.01, - "designSurfaceFamily": null, - "discretization": "central plus scalar dissipation", - "dissipationLumpingParameter": 6.0, - "dissipationScalingExponent": 0.67, - "eddyVisInfRatio": 0.009, - "equationMode": "time spectral", - "equationType": "Euler", - "eulerWallTreatment": "linear pressure extrapolation", - "firstRun": true, - "flowType": "external", - "forcesAsTractions": true, - "frozenTurbulence": false, - "globalPreconditioner": "additive Schwarz", - "gridFile": "input_files/naca64A010_euler-L2.cgns", - "gridPrecision": "double", - "gridPrecisionSurface": "single", - "infChangeCorrection": true, - "innerPreconIts": 1, - "isoVariables": [], - "isosurface": {}, - "liftIndex": 2, - "limiter": "van Albada", - "loadBalanceIter": 10, - "loadImbalance": 0.1, - "localPreconditioner": "ILU", - "lowSpeedPreconditioner": false, - "machMode": false, - "matrixOrdering": "RCM", - "maxL2DeviationFactor": 1.0, - "meshSurfaceFamily": null, - "monitorVariables": [ - "resrho", - "cl" - ], - "nCycles": 200000, - "nCyclesCoarse": 500, - "nFloodIter": -1, - "nRKReset": 5, - "nRefine": 10, - "nSaveSurface": 1, - "nSaveVolume": 1, - "nSubiter": 1, - "nSubiterTurb": 3, - "nTimeStepsCoarse": 48, - "nTimeStepsFine": 400, - "nearWallDist": 0.1, - "numberSolutions": true, - "outerPreconIts": 3, - "outputDirectory": "tests/output_files", - "outputSurfaceFamily": "allSurfaces", - "overlapFactor": 0.9, - "oversetLoadBalance": true, - "oversetPriority": {}, - "oversetProjTol": 1e-12, - "oversetUpdateMode": "frozen", - "pMode": false, - "partitionLikeNProc": -1, - "partitionOnly": false, - "preconditionerSide": "right", - "printIterations": true, - "printTiming": true, - "printWarnings": true, - "qMode": false, - "rMode": false, - "resAveraging": "alternate", - "restartAdjoint": true, - "restartFile": null, - "restrictionRelaxation": 0.8, - "selfZipCutoff": 120.0, - "sepSensorOffset": 0.0, - "sepSensorSharpness": 10.0, - "setMonitor": true, - "skipAfterFailedAdjoint": true, - "smoothParameter": 1.5, - "smoother": "Runge-Kutta", - "solutionPrecision": "single", - "solutionPrecisionSurface": "single", - "storeRindLayer": true, - "surfaceVariables": [ - "cp", - "vx", - "vy", - "vz", - "mach" - ], - "timeAccuracy": 2, - "timeIntegrationScheme": "BDF", - "timeIntervals": 3, - "timeLimit": -1.0, - "turbResScale": 10000.0, - "turbulenceModel": "SA", - "turbulenceOrder": "first order", - "turbulenceProduction": "strain", - "useALE": true, - "useANKSolver": true, - "useApproxWallDistance": true, - "useBlockettes": false, - "useDiagTSPC": true, - "useGridMotion": false, - "useLinResMonitor": false, - "useMatrixFreedrdw": true, - "useNKSolver": true, - "useOversetWallScaling": false, - "useQCR": false, - "useRotationSA": false, - "useWallFunctions": false, - "useZipperMesh": true, - "useexternaldynamicmesh": true, - "useft2SA": true, - "usetsinterpolatedgridvelocity": true, - "verifyExtra": true, - "verifySpatial": true, - "verifyState": true, - "vis2": 0.25, - "vis2Coarse": 0.5, - "vis4": 0.0156, - "viscPC": false, - "viscWallTreatment": "constant pressure extrapolation", - "viscousSurfaceVelocities": true, - "volumeVariables": [ - "resrho" - ], - "wallDistCutoff": 1e+20, - "windAxis": false, - "writeSurfaceSolution": false, - "writeTecplotSurfaceSolution": false, - "writeVolumeSolution": false, - "zipperSurfaceFamily": null - } -} diff --git a/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py b/tests/reg_tests/test_time_spectral.py similarity index 97% rename from tests/reg_tests/test_adjoint_time_spectral_naca64A010.py rename to tests/reg_tests/test_time_spectral.py index 7311e41c6..2477e2bbc 100644 --- a/tests/reg_tests/test_adjoint_time_spectral_naca64A010.py +++ b/tests/reg_tests/test_time_spectral.py @@ -25,7 +25,7 @@ class TestAdjointTimeSpectral(reg_test_classes.RegTest): """ N_PROCS = 2 - ref_file = "adjoint_euler_time_spectral_naca64A010.json" + ref_file = "euler_time_spectral_naca64A010.json" def setUp(self): # Introduce a transfer class for displacement transfer from struct to aero. @@ -137,7 +137,7 @@ def setVolumeMesh(self): alpha *= deltaAlpha # Create the solver - self.CFDSolver = ADFLOW(options=options, debug=True) + self.CFDSolver = ADFLOW(options=options, debug=False) # Deform the mesh mesh = USMesh(options=meshOptions) @@ -153,12 +153,12 @@ def setVolumeMesh(self): # Solve self.CFDSolver(self.ap) - # Check solution - self.assert_solution_failure() - # Initialize residual for adjoint tests self.CFDSolver.getResidual(self.ap) + def test_solution_failure(self): + self.assert_solution_failure() + def test_functions(self): utils.assert_functions_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-8) diff --git a/tests/reg_tests/test_time_spectral_naca64A010.py b/tests/reg_tests/test_time_spectral_naca64A010.py deleted file mode 100644 index 7f71acd84..000000000 --- a/tests/reg_tests/test_time_spectral_naca64A010.py +++ /dev/null @@ -1,166 +0,0 @@ -# built-ins -import unittest -import os -import copy -import numpy - -# MACH testing class -from adflow import ADFLOW -from idwarp import USMesh - -# import the testing utilities -import reg_test_utils as utils - -from reg_default_options import adflowDefOpts - -from reg_aeroproblems import ap_naca64A010_time_spectral -import reg_test_classes - - -baseDir = os.path.dirname(os.path.abspath(__file__)) - - -class TestSolve(reg_test_classes.RegTest): - """ - Tests for time spectral case for an airfoil. - """ - - N_PROCS = 2 - ref_file = "solve_euler_time_spectral_naca64A010.json" - - def setUp(self): - # Introduce a transfer class for displacement transfer from struct to aero. - class Transfer: - # simplified transfer class - # converting csd displacement to cfd surface nodes - - def __init__(self, alpha, xRot, aeroSolver): - # takes in displacement history - - self.alpha = alpha - self.ntimeintervalsspectral = len(alpha) - - self.aeroSolver = aeroSolver # a shallow copy of CFD solver - - self.xRot = xRot - - def getUndeformedSurfaceNodes(self): - self.MDGroup = self.aeroSolver.allWallsGroup - - self.cfdPts0 = self.aeroSolver.getSurfaceCoordinates(self.MDGroup, includeZipper=False) - - def setDisplacements(self): - xRot = self.xRot - ntimeintervalsspectral = self.ntimeintervalsspectral - alpha = self.alpha # notice a shallow copy introduced here; dont change the underlying obj! - cfdPoints_init = self.cfdPts0 # notice a shallow copy introduced here; dont change the underlying obj! - - N_pts = cfdPoints_init.shape[0] - - self.cfdPts = [] - - for sps in range(ntimeintervalsspectral): - cfdPoints_deformed = numpy.zeros((N_pts, 3)) - - ptch_loc = alpha[sps] - - cc = numpy.cos(ptch_loc) - ss = numpy.sin(ptch_loc) - - for j in range(N_pts): - cfdPoints_deformed[j, 0] = cc * (cfdPoints_init[j, 0] - xRot) + ss * cfdPoints_init[j, 1] + xRot - cfdPoints_deformed[j, 1] = -ss * (cfdPoints_init[j, 0] - xRot) + cc * cfdPoints_init[j, 1] - cfdPoints_deformed[j, 2] = cfdPoints_init[j, 2] - - self.cfdPts.append(cfdPoints_deformed) - - def setVolumeMesh(self): - ntimeintervalsspectral = self.ntimeintervalsspectral - - for sps in range(ntimeintervalsspectral): - self.aeroSolver.mesh.setSurfaceCoordinates(self.cfdPts[sps]) - self.aeroSolver.mesh.warpMesh() - m = self.aeroSolver.mesh.getSolverGrid() - self.aeroSolver.adflow.warping.setgridforoneinstance(m, sps=sps + 1) - - self.aeroSolver._updateGeomInfo = True - self.aeroSolver.updateGeometryInfo() - - # Set up the test for the parent RegTest class. - super().setUp() - - gridFile = os.path.join(baseDir, "../../input_files/naca64A010_euler-L2.cgns") - - ntimeintervalsspectral = 3 - options = copy.copy(adflowDefOpts) - options.update( - { - "gridfile": gridFile, - "outputDirectory": os.path.join(baseDir, "../output_files"), - "writeVolumeSolution": False, - "writeSurfaceSolution": False, - "blocksplitting": True, - "useblockettes": False, - "equationtype": "Euler", - "equationmode": "time spectral", - "mgcycle": "sg", - "l2convergence": 1e-15, - "ncycles": 200000, - "monitorvariables": ["resrho", "cl"], - "usenksolver": True, - "nkswitchtol": 1e-4, - "NKSubSpaceSize": 400, - "applypcsubspacesize": 400, - "useanksolver": True, - "ankswitchtol": 1e-2, - "anksubspacesize": 50, - "alphafollowing": False, - "timeintervals": ntimeintervalsspectral, - "useexternaldynamicmesh": True, - "usetsinterpolatedgridvelocity": True, - } - ) - - # Grid option - meshOptions = { - "gridFile": gridFile, - } - - # Setup aeroproblem - self.ap = copy.copy(ap_naca64A010_time_spectral) - - # Motion history - alpha_0 = 1.01 - deltaAlpha = -alpha_0 * numpy.pi / 180.0 - alpha = numpy.linspace(0.0, 2.0 * numpy.pi, ntimeintervalsspectral + 1)[:-1] - alpha = -numpy.sin(alpha) - alpha *= deltaAlpha - - # Create the solver - self.CFDSolver = ADFLOW(options=options, debug=False) - # self.CFDSolver.addSlices("z", [0.5]) - - # Deform the mesh - mesh = USMesh(options=meshOptions) - self.CFDSolver.setMesh(mesh) - - # deformation - xRot = 0.25 # Hard copied from the reference file. - TSTransfer = Transfer(alpha, xRot, self.CFDSolver) - TSTransfer.getUndeformedSurfaceNodes() - TSTransfer.setDisplacements() - TSTransfer.setVolumeMesh() - - def test_solve(self): - # do the solve - self.CFDSolver(self.ap) - - # check if the solution failed - self.assert_solution_failure() - - # check its accuracy - utils.assert_functions_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-8) - - -if __name__ == "__main__": - unittest.main()