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 8e12ed2e6..f330082e0 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 -! rw status of diff variables: *dw:in-out -! 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 @@ -535,8 +538,18 @@ subroutine initres_block_d(varstart, varend, nn, sps) ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m 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 @@ -573,6 +586,217 @@ subroutine initres_block_d(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 + dwd(i, j, k, l) = 0.0_8 + 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. + ll = 3*sps - 2 + (l-ivx) +! 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. 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 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 + 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 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. + 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 +! 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. 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 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. @@ -677,6 +901,170 @@ 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. + ll = 3*sps - 2 + (l-ivx) +! 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. + 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 +! 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 b16d723d8..f63db42f7 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 -! rw status of diff variables: *dw:in-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 @@ -558,8 +562,13 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m 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 @@ -569,12 +578,133 @@ subroutine initres_block_b(varstart, varend, nn, sps) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol2b(1) + call pushcontrol3b(3) else - call pushcontrol2b(0) + call pushcontrol3b(2) + 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 +! 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. + 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 + 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 +! 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. + 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 +! 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 pushcontrol2b(2) + call pushcontrol3b(4) end select do l=varend,varstart,-1 do j=jl,2,-1 @@ -602,8 +732,146 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do end do - call popcontrol2b(branch) - if (branch .eq. 0) then + call popcontrol3b(branch) + 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 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 do + end do + 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 + 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 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 do + end do + end if + else if (branch .eq. 2) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -613,7 +881,7 @@ subroutine initres_block_b(varstart, varend, nn, sps) end do end do end do - else if (branch .eq. 1) then + else if (branch .eq. 3) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -687,6 +955,170 @@ 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. + ll = 3*sps - 2 + (l-ivx) +! 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. + 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 +! 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/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90 index edf78f8d1..1b7da88dd 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 -! rw status of diff variables: *dw:in-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 @@ -507,8 +507,12 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! integer(kind=inttype) :: mm, ll, ii, jj, i, j, k, l, m 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 @@ -518,12 +522,98 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) ! steady state computation. ! determine the currently active multigrid level. if (currentlevel .eq. groundlevel) then - call pushcontrol2b(1) + call pushcontrol3b(3) else - call pushcontrol2b(0) + call pushcontrol3b(2) + 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 +! 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. + ll = 3*sps - 2 + (l-ivx) +! 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 +! 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. + ll = 3*sps - 2 + (l-ivx) +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 pushcontrol2b(2) + call pushcontrol3b(4) end select do l=varend,varstart,-1 do j=jl,2,-1 @@ -551,8 +641,125 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do end do - call popcontrol2b(branch) - if (branch .eq. 0) then + call popcontrol3b(branch) + 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 + 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 do + end do + 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 + 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 + 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 do + end do + end if + else if (branch .eq. 2) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -562,7 +769,7 @@ subroutine initres_block_fast_b(varstart, varend, nn, sps) end do end do end do - else if (branch .eq. 1) then + else if (branch .eq. 3) then do l=varend,varstart,-1 do k=kl,2,-1 do j=jl,2,-1 @@ -636,6 +843,170 @@ 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. + ll = 3*sps - 2 + (l-ivx) +! 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. + 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 +! 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/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/src/solver/residuals.F90 b/src/solver/residuals.F90 index 44d8a39a6..5b56ccc77 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. @@ -754,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. @@ -768,17 +769,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 +806,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 +856,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. @@ -855,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 @@ -873,6 +893,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 +919,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 +946,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 +968,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/reg_tests/refs/solve_euler_time_spectral_naca64A010.json b/tests/reg_tests/refs/euler_time_spectral_naca64A010.json similarity index 71% rename from tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json rename to tests/reg_tests/refs/euler_time_spectral_naca64A010.json index e89f229eb..4a9c6ec3c 100644 --- a/tests/reg_tests/refs/solve_euler_time_spectral_naca64A010.json +++ b/tests/reg_tests/refs/euler_time_spectral_naca64A010.json @@ -1,32 +1,47 @@ { + "Adjoint States": 0.0, + "Eval Functions Sens:": { + "64A010pitchingTS_cd": {}, + "64A010pitchingTS_cl": {}, + "64A010pitchingTS_cmz": {} + }, "Eval Functions:": { - "64A010pitchingTS_cd": 0.0025000608682108497, - "64A010pitchingTS_cl": 0.0037769183221013254, - "64A010pitchingTS_cmz": 0.0011186213899540035 + "64A010pitchingTS_cd": 0.002500060868193104, + "64A010pitchingTS_cl": 0.0037769183216176568, + "64A010pitchingTS_cmz": 0.0011186213900658232 }, "metadata": { "ADPC": false, - "AGMGLevels": 1, - "AGMGNSmooth": 3, "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, @@ -37,39 +52,53 @@ "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, - "L2Convergence": 1e-15, + "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-06, + "adjointL2Convergence": 1e-14, "adjointL2ConvergenceAbs": 1e-16, "adjointL2ConvergenceRel": 1e-16, "adjointMaxIter": 500, + "adjointMaxL2DeviationFactor": 1.0, "adjointMonitorStep": 10, "adjointSolver": "GMRES", "adjointSubspaceSize": 100, @@ -82,32 +111,44 @@ "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": "input_files/naca64A010_euler-L2.cgns", + "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, @@ -119,6 +160,7 @@ "machMode": false, "matrixOrdering": "RCM", "maxL2DeviationFactor": 1.0, + "meshMaxSkewness": 1.0, "meshSurfaceFamily": null, "monitorVariables": [ "resrho", @@ -138,9 +180,10 @@ "nearWallDist": 0.1, "numberSolutions": true, "outerPreconIts": 3, - "outputDirectory": "tests/output_files", + "outputDirectory": "tests/tests/tests/tests/output_files", "outputSurfaceFamily": "allSurfaces", "overlapFactor": 0.9, + "oversetDebugPrint": false, "oversetLoadBalance": true, "oversetPriority": {}, "oversetProjTol": 1e-12, @@ -149,16 +192,25 @@ "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, @@ -167,6 +219,7 @@ "smoother": "Runge-Kutta", "solutionPrecision": "single", "solutionPrecisionSurface": "single", + "storeConvHist": true, "storeRindLayer": true, "surfaceVariables": [ "cp", @@ -183,11 +236,14 @@ "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, @@ -195,11 +251,11 @@ "useOversetWallScaling": false, "useQCR": false, "useRotationSA": false, + "useSkewnessCheck": false, + "useTSInterpolatedGridVelocity": true, "useWallFunctions": false, "useZipperMesh": true, - "useexternaldynamicmesh": true, "useft2SA": true, - "usetsinterpolatedgridvelocity": true, "verifyExtra": true, "verifySpatial": true, "verifyState": true, @@ -214,6 +270,8 @@ ], "wallDistCutoff": 1e+20, "windAxis": false, + "writeSolutionDigits": 3, + "writeSolutionEachIter": false, "writeSurfaceSolution": false, "writeTecplotSurfaceSolution": false, "writeVolumeSolution": false, diff --git a/tests/reg_tests/test_time_spectral_naca64A010.py b/tests/reg_tests/test_time_spectral.py similarity index 83% rename from tests/reg_tests/test_time_spectral_naca64A010.py rename to tests/reg_tests/test_time_spectral.py index 7f71acd84..2477e2bbc 100644 --- a/tests/reg_tests/test_time_spectral_naca64A010.py +++ b/tests/reg_tests/test_time_spectral.py @@ -16,17 +16,16 @@ 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): +class TestAdjointTimeSpectral(reg_test_classes.RegTest): """ - Tests for time spectral case for an airfoil. + Tests for time spectral adjoint on the NACA 64A010 case. """ N_PROCS = 2 - ref_file = "solve_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. @@ -52,8 +51,8 @@ def getUndeformedSurfaceNodes(self): 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! + alpha = self.alpha + cfdPoints_init = self.cfdPts0 N_pts = cfdPoints_init.shape[0] @@ -104,7 +103,7 @@ def setVolumeMesh(self): "equationtype": "Euler", "equationmode": "time spectral", "mgcycle": "sg", - "l2convergence": 1e-15, + "l2convergence": 1e-13, "ncycles": 200000, "monitorvariables": ["resrho", "cl"], "usenksolver": True, @@ -118,6 +117,7 @@ def setVolumeMesh(self): "timeintervals": ntimeintervalsspectral, "useexternaldynamicmesh": True, "usetsinterpolatedgridvelocity": True, + "adjointl2convergence": 1e-14, } ) @@ -138,7 +138,6 @@ def setVolumeMesh(self): # Create the solver self.CFDSolver = ADFLOW(options=options, debug=False) - # self.CFDSolver.addSlices("z", [0.5]) # Deform the mesh mesh = USMesh(options=meshOptions) @@ -151,16 +150,30 @@ def setVolumeMesh(self): TSTransfer.setDisplacements() TSTransfer.setVolumeMesh() - def test_solve(self): - # do the solve + # Solve self.CFDSolver(self.ap) - # check if the solution failed + # Initialize residual for adjoint tests + self.CFDSolver.getResidual(self.ap) + + def test_solution_failure(self): self.assert_solution_failure() - # check its accuracy + 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()