diff --git a/src/adjoint/Makefile_tapenade b/src/adjoint/Makefile_tapenade index d0077a2ea..01ba4405d 100644 --- a/src/adjoint/Makefile_tapenade +++ b/src/adjoint/Makefile_tapenade @@ -138,6 +138,12 @@ BCRoutines%bcSymm1stHalo(ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, bcData%norm BCRoutines%bcSymm2ndHalo(ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bcData%norm) > \ (ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bcData%norm) \ \ +BCRoutines%bcAntiSymm1stHalo(ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, pInf, wInf, bcData%norm) > \ + (ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, pInf, wInf, bcData%norm) \ +\ +BCRoutines%bcAntiSymm2ndHalo(ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, pInf, wInf, bcData%norm) > \ + (ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, pInf, wInf, bcData%norm) \ +\ BCRoutines%bcSymmPolar1stHalo(xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2) > \ (xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2) \ \ diff --git a/src/adjoint/adjointExtra.F90 b/src/adjoint/adjointExtra.F90 index 71ef9e328..b6dc9e310 100644 --- a/src/adjoint/adjointExtra.F90 +++ b/src/adjoint/adjointExtra.F90 @@ -436,7 +436,7 @@ subroutine xhalo_block ! The actual correction of the coordinates only takes ! place for symmetry planes. - testSymmetry: if (BCType(mm) == Symm) then + testSymmetry: if (BCType(mm) == Symm .or. BCType(mm) == AntiSymm) then ! Set some variables, depending on the block face on ! which the subface is located. diff --git a/src/adjoint/outputForward/BCExtra_d.F90 b/src/adjoint/outputForward/BCExtra_d.F90 index faf616df6..dfec30305 100644 --- a/src/adjoint/outputForward/BCExtra_d.F90 +++ b/src/adjoint/outputForward/BCExtra_d.F90 @@ -54,6 +54,25 @@ subroutine applyAllBC_block_d(secondHalo) end do end if + ! ------------------------------------ + ! AntiSymmetry Boundary Condition + ! ------------------------------------ + do nn=1, nBocos + if (bcType(nn) == symm) then + call setBCPointers_d(nn, .False.) + call bcAntiSymm1stHalo_d(nn) + end if + end do + + if (secondHalo) then + do nn=1, nBocos + if (bcType(nn) == symm) then + call setBCPointers_d(nn, .False.) + call bcAntiSymm2ndHalo_d(nn) + end if + end do + end if + ! ------------------------------------ ! Symmetry Polar Boundary Condition ! ------------------------------------ diff --git a/src/adjoint/outputForward/BCRoutines_d.f90 b/src/adjoint/outputForward/BCRoutines_d.f90 index 3d28d0bbf..56cc0d292 100644 --- a/src/adjoint/outputForward/BCRoutines_d.f90 +++ b/src/adjoint/outputForward/BCRoutines_d.f90 @@ -52,6 +52,25 @@ subroutine applyallbc_block(secondhalo) end if !$ad ii-loop ! ------------------------------------ +! antisymmetry boundary condition +! ------------------------------------ + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm1sthalo(nn) + end if + end do + if (secondhalo) then +!$ad ii-loop + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm2ndhalo(nn) + end if + end do + end if +!$ad ii-loop +! ------------------------------------ ! symmetry polar boundary condition ! ------------------------------------ do nn=1,nbocos @@ -411,6 +430,312 @@ subroutine bcsymm2ndhalo(nn) end do end subroutine bcsymm2ndhalo +! differentiation of bcantisymm1sthalo in forward (tangent) mode (with options i4 dr8 r8): +! variations of useful results: *rev1 *pp1 *rlv1 *ww1 +! with respect to varying inputs: pinf winf *rev1 *rev2 *pp1 +! *pp2 *rlv1 *rlv2 *ww1 *ww2 *(*bcdata.norm) +! rw status of diff variables: pinf:in winf:in *rev1:in-out *rev2:in +! *pp1:in-out *pp2:in *rlv1:in-out *rlv2:in *ww1:in-out +! *ww2:in *(*bcdata.norm):in +! plus diff mem management of: rev1:in rev2:in pp1:in pp2:in +! rlv1:in rlv2:in ww1:in ww2:in bcdata:in *bcdata.norm:in + subroutine bcantisymm1sthalo_d(nn) +! bcantisymm1sthalo applies the antisymmetry boundary conditions to a +! block. * it is assumed that the pointers in blockpointers are +! already set to the correct block on the correct grid level. +! it should only be applied for high-speed near free-face operation, +! unless you are clear what you use it for +! +! in case also the second halo must be set, a second loop is +! execulted calling bcsymm2ndhalo. this is the only correct way +! in case the block contains only 1 cell between two symmetry +! planes, i.e. a 2d problem. + use constants + use blockpointers, only : bcdata, bcdatad + use flowvarrefstate, only : viscous, eddymodel, pinf, pinfd, winf,& +& winfd + use bcpointers_d, only : gamma1, gamma2, ww1, ww1d, ww2, ww2d, pp1, & +& pp1d, pp2, pp2d, rlv1, rlv1d, rlv2, rlv2d, istart, jstart, isize, & +& jsize, rev1, rev1d, rev2, rev2d + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype) :: vnd + real(kind=realtype), dimension(5) :: dww2 + real(kind=realtype), dimension(5) :: dww2d + intrinsic mod + real(kind=realtype) :: temp + real(kind=realtype) :: temp0 + real(kind=realtype) :: temp1 + dww2d = 0.0_8 +! loop over the generic subface to set the state in the +! 1-st level halos + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart +! determine twice the normal velocity component, +! which must be substracted from the donor velocity +! to obtain the halo velocity. + dww2d(irho) = ww2d(i, j, irho) - winfd(irho) + dww2(irho) = ww2(i, j, irho) - winf(irho) + dww2d(ivx) = ww2d(i, j, ivx) - winfd(ivx) + dww2(ivx) = ww2(i, j, ivx) - winf(ivx) + dww2d(ivy) = ww2d(i, j, ivy) - winfd(ivy) + dww2(ivy) = ww2(i, j, ivy) - winf(ivy) + dww2d(ivz) = ww2d(i, j, ivz) - winfd(ivz) + dww2(ivz) = ww2(i, j, ivz) - winf(ivz) + dww2d(irhoe) = ww2d(i, j, irhoe) - winfd(irhoe) + dww2(irhoe) = ww2(i, j, irhoe) - winf(irhoe) + temp = bcdata(nn)%norm(i, j, 1) + temp0 = bcdata(nn)%norm(i, j, 2) + temp1 = bcdata(nn)%norm(i, j, 3) + vnd = two*(temp*dww2d(ivx)+dww2(ivx)*bcdatad(nn)%norm(i, j, 1)+& +& temp0*dww2d(ivy)+dww2(ivy)*bcdatad(nn)%norm(i, j, 2)+temp1*dww2d& +& (ivz)+dww2(ivz)*bcdatad(nn)%norm(i, j, 3)) + vn = two*(dww2(ivx)*temp+dww2(ivy)*temp0+dww2(ivz)*temp1) +! determine the flow variables in the halo cell. + temp1 = bcdata(nn)%norm(i, j, 1) + ww1d(i, j, ivx) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 1) - dww2d& +& (ivx) + ww1(i, j, ivx) = vn*temp1 - dww2(ivx) + temp1 = bcdata(nn)%norm(i, j, 2) + ww1d(i, j, ivy) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 2) - dww2d& +& (ivy) + ww1(i, j, ivy) = vn*temp1 - dww2(ivy) + temp1 = bcdata(nn)%norm(i, j, 3) + ww1d(i, j, ivz) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 3) - dww2d& +& (ivz) + ww1(i, j, ivz) = vn*temp1 - dww2(ivz) + ww1d(i, j, irho) = winfd(irho) - dww2d(irho) + ww1(i, j, irho) = -dww2(irho) + winf(irho) + ww1d(i, j, ivx) = winfd(ivx) + ww1d(i, j, ivx) + ww1(i, j, ivx) = winf(ivx) + ww1(i, j, ivx) + ww1d(i, j, ivy) = winfd(ivy) + ww1d(i, j, ivy) + ww1(i, j, ivy) = winf(ivy) + ww1(i, j, ivy) + ww1d(i, j, ivz) = winfd(ivz) + ww1d(i, j, ivz) + ww1(i, j, ivz) = winf(ivz) + ww1(i, j, ivz) + ww1d(i, j, irhoe) = winfd(irhoe) - dww2d(irhoe) + ww1(i, j, irhoe) = -dww2(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma1(i, j) = gamma2(i, j) + pp1d(i, j) = 2*pinfd - pp2d(i, j) + pp1(i, j) = pinf - (pp2(i, j)-pinf) + if (viscous) then + rlv1d(i, j) = -rlv2d(i, j) + rlv1(i, j) = -rlv2(i, j) + end if + if (eddymodel) then + rev1d(i, j) = -rev2d(i, j) + rev1(i, j) = -rev2(i, j) + end if + end do + end subroutine bcantisymm1sthalo_d + + subroutine bcantisymm1sthalo(nn) +! bcantisymm1sthalo applies the antisymmetry boundary conditions to a +! block. * it is assumed that the pointers in blockpointers are +! already set to the correct block on the correct grid level. +! it should only be applied for high-speed near free-face operation, +! unless you are clear what you use it for +! +! in case also the second halo must be set, a second loop is +! execulted calling bcsymm2ndhalo. this is the only correct way +! in case the block contains only 1 cell between two symmetry +! planes, i.e. a 2d problem. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers_d, only : gamma1, gamma2, ww1, ww2, pp1, pp2, rlv1, & +& rlv2, istart, jstart, isize, jsize, rev1, rev2 + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww2 + intrinsic mod +!$ad ii-loop +! loop over the generic subface to set the state in the +! 1-st level halos + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart +! determine twice the normal velocity component, +! which must be substracted from the donor velocity +! to obtain the halo velocity. + dww2(irho) = ww2(i, j, irho) - winf(irho) + dww2(ivx) = ww2(i, j, ivx) - winf(ivx) + dww2(ivy) = ww2(i, j, ivy) - winf(ivy) + dww2(ivz) = ww2(i, j, ivz) - winf(ivz) + dww2(irhoe) = ww2(i, j, irhoe) - winf(irhoe) + vn = two*(dww2(ivx)*bcdata(nn)%norm(i, j, 1)+dww2(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww2(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww1(i, j, ivx) = -dww2(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww1(i, j, ivy) = -dww2(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww1(i, j, ivz) = -dww2(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww1(i, j, irho) = -dww2(irho) + winf(irho) + ww1(i, j, ivx) = winf(ivx) + ww1(i, j, ivx) + ww1(i, j, ivy) = winf(ivy) + ww1(i, j, ivy) + ww1(i, j, ivz) = winf(ivz) + ww1(i, j, ivz) + ww1(i, j, irhoe) = -dww2(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma1(i, j) = gamma2(i, j) + pp1(i, j) = pinf - (pp2(i, j)-pinf) + if (viscous) rlv1(i, j) = -rlv2(i, j) + if (eddymodel) rev1(i, j) = -rev2(i, j) + end do + end subroutine bcantisymm1sthalo + +! differentiation of bcantisymm2ndhalo in forward (tangent) mode (with options i4 dr8 r8): +! variations of useful results: *rev0 *pp0 *rlv0 *ww0 +! with respect to varying inputs: pinf winf *rev0 *rev3 *pp0 +! *pp3 *rlv0 *rlv3 *ww0 *ww3 *(*bcdata.norm) +! rw status of diff variables: pinf:in winf:in *rev0:in-out *rev3:in +! *pp0:in-out *pp3:in *rlv0:in-out *rlv3:in *ww0:in-out +! *ww3:in *(*bcdata.norm):in +! plus diff mem management of: rev0:in rev3:in pp0:in pp3:in +! rlv0:in rlv3:in ww0:in ww3:in bcdata:in *bcdata.norm:in + subroutine bcantisymm2ndhalo_d(nn) +! bcsymm2ndhalo applies the symmetry boundary conditions to a +! block for the 2nd halo. this routine is separate as it makes +! ad slightly easier. + use constants + use blockpointers, only : bcdata, bcdatad + use flowvarrefstate, only : viscous, eddymodel, pinf, pinfd, winf,& +& winfd + use bcpointers_d, only : gamma0, gamma3, ww0, ww0d, ww3, ww3d, pp0, & +& pp0d, pp3, pp3d, rlv0, rlv0d, rlv3, rlv3d, rev0, rev0d, rev3, rev3d,& +& istart, jstart, isize, jsize + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype) :: vnd + real(kind=realtype), dimension(5) :: dww3 + real(kind=realtype), dimension(5) :: dww3d + intrinsic mod + real(kind=realtype) :: temp + real(kind=realtype) :: temp0 + real(kind=realtype) :: temp1 + dww3d = 0.0_8 +! if we need the second halo, do everything again, but using ww0, +! ww3 etc instead of ww2 and ww1. + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart + dww3d(irho) = ww3d(i, j, irho) - winfd(irho) + dww3(irho) = ww3(i, j, irho) - winf(irho) + dww3d(ivx) = ww3d(i, j, ivx) - winfd(ivx) + dww3(ivx) = ww3(i, j, ivx) - winf(ivx) + dww3d(ivy) = ww3d(i, j, ivy) - winfd(ivy) + dww3(ivy) = ww3(i, j, ivy) - winf(ivy) + dww3d(ivz) = ww3d(i, j, ivz) - winfd(ivz) + dww3(ivz) = ww3(i, j, ivz) - winf(ivz) + dww3d(irhoe) = ww3d(i, j, irhoe) - winfd(irhoe) + dww3(irhoe) = ww3(i, j, irhoe) - winf(irhoe) + temp = bcdata(nn)%norm(i, j, 1) + temp0 = bcdata(nn)%norm(i, j, 2) + temp1 = bcdata(nn)%norm(i, j, 3) + vnd = two*(temp*dww3d(ivx)+dww3(ivx)*bcdatad(nn)%norm(i, j, 1)+& +& temp0*dww3d(ivy)+dww3(ivy)*bcdatad(nn)%norm(i, j, 2)+temp1*dww3d& +& (ivz)+dww3(ivz)*bcdatad(nn)%norm(i, j, 3)) + vn = two*(dww3(ivx)*temp+dww3(ivy)*temp0+dww3(ivz)*temp1) +! determine the flow variables in the halo cell. + temp1 = bcdata(nn)%norm(i, j, 1) + ww0d(i, j, ivx) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 1) - dww3d& +& (ivx) + ww0(i, j, ivx) = vn*temp1 - dww3(ivx) + temp1 = bcdata(nn)%norm(i, j, 2) + ww0d(i, j, ivy) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 2) - dww3d& +& (ivy) + ww0(i, j, ivy) = vn*temp1 - dww3(ivy) + temp1 = bcdata(nn)%norm(i, j, 3) + ww0d(i, j, ivz) = temp1*vnd + vn*bcdatad(nn)%norm(i, j, 3) - dww3d& +& (ivz) + ww0(i, j, ivz) = vn*temp1 - dww3(ivz) + ww0d(i, j, irho) = winfd(irho) - dww3d(irho) + ww0(i, j, irho) = -dww3(irho) + winf(irho) + ww0d(i, j, ivx) = winfd(ivx) + ww0d(i, j, ivx) + ww0(i, j, ivx) = winf(ivx) + ww0(i, j, ivx) + ww0d(i, j, ivy) = winfd(ivy) + ww0d(i, j, ivy) + ww0(i, j, ivy) = winf(ivy) + ww0(i, j, ivy) + ww0d(i, j, ivz) = winfd(ivz) + ww0d(i, j, ivz) + ww0(i, j, ivz) = winf(ivz) + ww0(i, j, ivz) + ww0d(i, j, irhoe) = winfd(irhoe) - dww3d(irhoe) + ww0(i, j, irhoe) = -dww3(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma0(i, j) = gamma3(i, j) + pp0d(i, j) = 2*pinfd - pp3d(i, j) + pp0(i, j) = pinf - (pp3(i, j)-pinf) + if (viscous) then + rlv0d(i, j) = -rlv3d(i, j) + rlv0(i, j) = -rlv3(i, j) + end if + if (eddymodel) then + rev0d(i, j) = -rev3d(i, j) + rev0(i, j) = -rev3(i, j) + end if + end do + end subroutine bcantisymm2ndhalo_d + + subroutine bcantisymm2ndhalo(nn) +! bcsymm2ndhalo applies the symmetry boundary conditions to a +! block for the 2nd halo. this routine is separate as it makes +! ad slightly easier. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers_d, only : gamma0, gamma3, ww0, ww3, pp0, pp3, rlv0, & +& rlv3, rev0, rev3, istart, jstart, isize, jsize + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww3 + intrinsic mod +!$ad ii-loop +! if we need the second halo, do everything again, but using ww0, +! ww3 etc instead of ww2 and ww1. + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart + dww3(irho) = ww3(i, j, irho) - winf(irho) + dww3(ivx) = ww3(i, j, ivx) - winf(ivx) + dww3(ivy) = ww3(i, j, ivy) - winf(ivy) + dww3(ivz) = ww3(i, j, ivz) - winf(ivz) + dww3(irhoe) = ww3(i, j, irhoe) - winf(irhoe) + vn = two*(dww3(ivx)*bcdata(nn)%norm(i, j, 1)+dww3(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww3(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww0(i, j, ivx) = -dww3(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww0(i, j, ivy) = -dww3(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww0(i, j, ivz) = -dww3(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww0(i, j, irho) = -dww3(irho) + winf(irho) + ww0(i, j, ivx) = winf(ivx) + ww0(i, j, ivx) + ww0(i, j, ivy) = winf(ivy) + ww0(i, j, ivy) + ww0(i, j, ivz) = winf(ivz) + ww0(i, j, ivz) + ww0(i, j, irhoe) = -dww3(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma0(i, j) = gamma3(i, j) + pp0(i, j) = pinf - (pp3(i, j)-pinf) + if (viscous) rlv0(i, j) = -rlv3(i, j) + if (eddymodel) rev0(i, j) = -rev3(i, j) + end do + end subroutine bcantisymm2ndhalo + ! differentiation of bcsymmpolar1sthalo in forward (tangent) mode (with options i4 dr8 r8): ! variations of useful results: *rev1 *pp1 *rlv1 *ww1 ! with respect to varying inputs: *xx *rev1 *rev2 *pp1 *pp2 *rlv1 diff --git a/src/adjoint/outputForward/adjointExtra_d.f90 b/src/adjoint/outputForward/adjointExtra_d.f90 index d723612dc..a412f404a 100644 --- a/src/adjoint/outputForward/adjointExtra_d.f90 +++ b/src/adjoint/outputForward/adjointExtra_d.f90 @@ -951,7 +951,7 @@ subroutine xhalo_block_d() loopbocos:do mm=1,nbocos ! the actual correction of the coordinates only takes ! place for symmetry planes. - if (bctype(mm) .eq. symm) then + if (bctype(mm) .eq. symm .or. bctype(mm) .eq. antisymm) then ! set some variables, depending on the block face on ! which the subface is located. norm(1) = bcdata(mm)%symnorm(1) @@ -1216,7 +1216,7 @@ subroutine xhalo_block() loopbocos:do mm=1,nbocos ! the actual correction of the coordinates only takes ! place for symmetry planes. - if (bctype(mm) .eq. symm) then + if (bctype(mm) .eq. symm .or. bctype(mm) .eq. antisymm) then ! set some variables, depending on the block face on ! which the subface is located. norm(1) = bcdata(mm)%symnorm(1) diff --git a/src/adjoint/outputForward/turbBCRoutines_d.f90 b/src/adjoint/outputForward/turbBCRoutines_d.f90 index 6f41ad314..e8264436e 100644 --- a/src/adjoint/outputForward/turbBCRoutines_d.f90 +++ b/src/adjoint/outputForward/turbBCRoutines_d.f90 @@ -1046,7 +1046,7 @@ subroutine bcturbtreatment_d() call bcturbwall_d(nn) !============================================================= !============================================================= - case (symm, symmpolar, eulerwall) + case (symm, antisymm, symmpolar, eulerwall) ! symmetry, polar symmetry or inviscid wall. treatment of ! the turbulent equations is identical. call bcturbsymm(nn) @@ -1138,7 +1138,7 @@ subroutine bcturbtreatment() call bcturbwall(nn) !============================================================= !============================================================= - case (symm, symmpolar, eulerwall) + case (symm, antisymm, symmpolar, eulerwall) ! symmetry, polar symmetry or inviscid wall. treatment of ! the turbulent equations is identical. call bcturbsymm(nn) diff --git a/src/adjoint/outputReverse/BCExtra_b.F90 b/src/adjoint/outputReverse/BCExtra_b.F90 index e7bb4df83..8af390295 100644 --- a/src/adjoint/outputReverse/BCExtra_b.F90 +++ b/src/adjoint/outputReverse/BCExtra_b.F90 @@ -126,5 +126,24 @@ subroutine applyAllBC_block_b(secondHalo) END IF END DO END if + + ! ------------------------------------ + ! AntiSymmetry Boundary Condition + ! ------------------------------------ + if (secondHalo) then + DO mm=1,nbocos + IF (bctype(mm) .EQ. antiSymm) THEN + CALL setBCPointers_d(mm, .false.) + CALL BCAntiSymm2ndhalo_b(mm) + END IF + END DO + END if + + DO mm=1,nbocos + IF (bctype(mm) .EQ. antiSymm) THEN + CALL setBCPointers_d(mm, .false.) + CALL BCAntiSymm1sthalo_b(mm) + END IF + END DO end subroutine applyAllBC_block_b end module BCExtra_b diff --git a/src/adjoint/outputReverse/BCRoutines_b.f90 b/src/adjoint/outputReverse/BCRoutines_b.f90 index cdbf2ebf7..d31206b77 100644 --- a/src/adjoint/outputReverse/BCRoutines_b.f90 +++ b/src/adjoint/outputReverse/BCRoutines_b.f90 @@ -52,6 +52,25 @@ subroutine applyallbc_block(secondhalo) end if !$ad ii-loop ! ------------------------------------ +! antisymmetry boundary condition +! ------------------------------------ + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm1sthalo(nn) + end if + end do + if (secondhalo) then +!$ad ii-loop + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm2ndhalo(nn) + end if + end do + end if +!$ad ii-loop +! ------------------------------------ ! symmetry polar boundary condition ! ------------------------------------ do nn=1,nbocos @@ -429,6 +448,352 @@ subroutine bcsymm2ndhalo(nn) end do end subroutine bcsymm2ndhalo +! differentiation of bcantisymm1sthalo in reverse (adjoint) mode (with options noisize i4 dr8 r8): +! gradient of useful results: pinf winf *rev1 *rev2 *pp1 +! *pp2 *rlv1 *rlv2 *ww1 *ww2 *(*bcdata.norm) +! with respect to varying inputs: pinf winf *rev1 *rev2 *pp1 +! *pp2 *rlv1 *rlv2 *ww1 *ww2 *(*bcdata.norm) +! rw status of diff variables: pinf:incr winf:incr *rev1:in-out +! *rev2:incr *pp1:in-out *pp2:incr *rlv1:in-out +! *rlv2:incr *ww1:in-out *ww2:incr *(*bcdata.norm):incr +! plus diff mem management of: rev1:in rev2:in pp1:in pp2:in +! rlv1:in rlv2:in ww1:in ww2:in bcdata:in *bcdata.norm:in + subroutine bcantisymm1sthalo_b(nn) +! bcantisymm1sthalo applies the antisymmetry boundary conditions to a +! block. * it is assumed that the pointers in blockpointers are +! already set to the correct block on the correct grid level. +! it should only be applied for high-speed near free-face operation, +! unless you are clear what you use it for +! +! in case also the second halo must be set, a second loop is +! execulted calling bcsymm2ndhalo. this is the only correct way +! in case the block contains only 1 cell between two symmetry +! planes, i.e. a 2d problem. + use constants + use blockpointers, only : bcdata, bcdatad + use flowvarrefstate, only : viscous, eddymodel, pinf, pinfd, winf,& +& winfd + use bcpointers_b, only : gamma1, gamma2, ww1, ww1d, ww2, ww2d, pp1, & +& pp1d, pp2, pp2d, rlv1, rlv1d, rlv2, rlv2d, istart, jstart, isize, & +& jsize, rev1, rev1d, rev2, rev2d + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype) :: vnd + real(kind=realtype), dimension(5) :: dww2 + real(kind=realtype), dimension(5) :: dww2d + intrinsic mod + real(kind=realtype) :: tempd + integer :: branch + dww2d = 0.0_8 +!$bwd-of ii-loop + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart +! determine twice the normal velocity component, +! which must be substracted from the donor velocity +! to obtain the halo velocity. + dww2(irho) = ww2(i, j, irho) - winf(irho) + dww2(ivx) = ww2(i, j, ivx) - winf(ivx) + dww2(ivy) = ww2(i, j, ivy) - winf(ivy) + dww2(ivz) = ww2(i, j, ivz) - winf(ivz) + dww2(irhoe) = ww2(i, j, irhoe) - winf(irhoe) + vn = two*(dww2(ivx)*bcdata(nn)%norm(i, j, 1)+dww2(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww2(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + if (viscous) then + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (eddymodel) then + rev2d(i, j) = rev2d(i, j) - rev1d(i, j) + rev1d(i, j) = 0.0_8 + end if + call popcontrol1b(branch) + if (branch .eq. 0) then + rlv2d(i, j) = rlv2d(i, j) - rlv1d(i, j) + rlv1d(i, j) = 0.0_8 + end if + pinfd = pinfd + 2*pp1d(i, j) + pp2d(i, j) = pp2d(i, j) - pp1d(i, j) + pp1d(i, j) = 0.0_8 + winfd(irhoe) = winfd(irhoe) + ww1d(i, j, irhoe) + dww2d(irhoe) = dww2d(irhoe) - ww1d(i, j, irhoe) + ww1d(i, j, irhoe) = 0.0_8 + winfd(ivz) = winfd(ivz) + ww1d(i, j, ivz) + winfd(ivy) = winfd(ivy) + ww1d(i, j, ivy) + winfd(ivx) = winfd(ivx) + ww1d(i, j, ivx) + winfd(irho) = winfd(irho) + ww1d(i, j, irho) + dww2d(irho) = dww2d(irho) - ww1d(i, j, irho) + ww1d(i, j, irho) = 0.0_8 + vnd = bcdata(nn)%norm(i, j, 3)*ww1d(i, j, ivz) + bcdatad(nn)%norm(i, j, 3) = bcdatad(nn)%norm(i, j, 3) + vn*ww1d(i& +& , j, ivz) + dww2d(ivz) = dww2d(ivz) - ww1d(i, j, ivz) + ww1d(i, j, ivz) = 0.0_8 + vnd = vnd + bcdata(nn)%norm(i, j, 2)*ww1d(i, j, ivy) + bcdatad(nn)%norm(i, j, 2) = bcdatad(nn)%norm(i, j, 2) + vn*ww1d(i& +& , j, ivy) + dww2d(ivy) = dww2d(ivy) - ww1d(i, j, ivy) + ww1d(i, j, ivy) = 0.0_8 + vnd = vnd + bcdata(nn)%norm(i, j, 1)*ww1d(i, j, ivx) + tempd = two*vnd + bcdatad(nn)%norm(i, j, 1) = bcdatad(nn)%norm(i, j, 1) + vn*ww1d(i& +& , j, ivx) + dww2(ivx)*tempd + dww2d(ivx) = dww2d(ivx) + bcdata(nn)%norm(i, j, 1)*tempd - ww1d(i& +& , j, ivx) + ww1d(i, j, ivx) = 0.0_8 + dww2d(ivy) = dww2d(ivy) + bcdata(nn)%norm(i, j, 2)*tempd + bcdatad(nn)%norm(i, j, 2) = bcdatad(nn)%norm(i, j, 2) + dww2(ivy)*& +& tempd + dww2d(ivz) = dww2d(ivz) + bcdata(nn)%norm(i, j, 3)*tempd + bcdatad(nn)%norm(i, j, 3) = bcdatad(nn)%norm(i, j, 3) + dww2(ivz)*& +& tempd + ww2d(i, j, irhoe) = ww2d(i, j, irhoe) + dww2d(irhoe) + winfd(irhoe) = winfd(irhoe) - dww2d(irhoe) + dww2d(irhoe) = 0.0_8 + ww2d(i, j, ivz) = ww2d(i, j, ivz) + dww2d(ivz) + winfd(ivz) = winfd(ivz) - dww2d(ivz) + dww2d(ivz) = 0.0_8 + ww2d(i, j, ivy) = ww2d(i, j, ivy) + dww2d(ivy) + winfd(ivy) = winfd(ivy) - dww2d(ivy) + dww2d(ivy) = 0.0_8 + ww2d(i, j, ivx) = ww2d(i, j, ivx) + dww2d(ivx) + winfd(ivx) = winfd(ivx) - dww2d(ivx) + dww2d(ivx) = 0.0_8 + ww2d(i, j, irho) = ww2d(i, j, irho) + dww2d(irho) + winfd(irho) = winfd(irho) - dww2d(irho) + dww2d(irho) = 0.0_8 + end do + end subroutine bcantisymm1sthalo_b + + subroutine bcantisymm1sthalo(nn) +! bcantisymm1sthalo applies the antisymmetry boundary conditions to a +! block. * it is assumed that the pointers in blockpointers are +! already set to the correct block on the correct grid level. +! it should only be applied for high-speed near free-face operation, +! unless you are clear what you use it for +! +! in case also the second halo must be set, a second loop is +! execulted calling bcsymm2ndhalo. this is the only correct way +! in case the block contains only 1 cell between two symmetry +! planes, i.e. a 2d problem. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers_b, only : gamma1, gamma2, ww1, ww2, pp1, pp2, rlv1, & +& rlv2, istart, jstart, isize, jsize, rev1, rev2 + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww2 + intrinsic mod +!$ad ii-loop +! loop over the generic subface to set the state in the +! 1-st level halos + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart +! determine twice the normal velocity component, +! which must be substracted from the donor velocity +! to obtain the halo velocity. + dww2(irho) = ww2(i, j, irho) - winf(irho) + dww2(ivx) = ww2(i, j, ivx) - winf(ivx) + dww2(ivy) = ww2(i, j, ivy) - winf(ivy) + dww2(ivz) = ww2(i, j, ivz) - winf(ivz) + dww2(irhoe) = ww2(i, j, irhoe) - winf(irhoe) + vn = two*(dww2(ivx)*bcdata(nn)%norm(i, j, 1)+dww2(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww2(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww1(i, j, ivx) = -dww2(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww1(i, j, ivy) = -dww2(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww1(i, j, ivz) = -dww2(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww1(i, j, irho) = -dww2(irho) + winf(irho) + ww1(i, j, ivx) = winf(ivx) + ww1(i, j, ivx) + ww1(i, j, ivy) = winf(ivy) + ww1(i, j, ivy) + ww1(i, j, ivz) = winf(ivz) + ww1(i, j, ivz) + ww1(i, j, irhoe) = -dww2(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma1(i, j) = gamma2(i, j) + pp1(i, j) = pinf - (pp2(i, j)-pinf) + if (viscous) rlv1(i, j) = -rlv2(i, j) + if (eddymodel) rev1(i, j) = -rev2(i, j) + end do + end subroutine bcantisymm1sthalo + +! differentiation of bcantisymm2ndhalo in reverse (adjoint) mode (with options noisize i4 dr8 r8): +! gradient of useful results: pinf winf *rev0 *rev3 *pp0 +! *pp3 *rlv0 *rlv3 *ww0 *ww3 *(*bcdata.norm) +! with respect to varying inputs: pinf winf *rev0 *rev3 *pp0 +! *pp3 *rlv0 *rlv3 *ww0 *ww3 *(*bcdata.norm) +! rw status of diff variables: pinf:incr winf:incr *rev0:in-out +! *rev3:incr *pp0:in-out *pp3:incr *rlv0:in-out +! *rlv3:incr *ww0:in-out *ww3:incr *(*bcdata.norm):incr +! plus diff mem management of: rev0:in rev3:in pp0:in pp3:in +! rlv0:in rlv3:in ww0:in ww3:in bcdata:in *bcdata.norm:in + subroutine bcantisymm2ndhalo_b(nn) +! bcsymm2ndhalo applies the symmetry boundary conditions to a +! block for the 2nd halo. this routine is separate as it makes +! ad slightly easier. + use constants + use blockpointers, only : bcdata, bcdatad + use flowvarrefstate, only : viscous, eddymodel, pinf, pinfd, winf,& +& winfd + use bcpointers_b, only : gamma0, gamma3, ww0, ww0d, ww3, ww3d, pp0, & +& pp0d, pp3, pp3d, rlv0, rlv0d, rlv3, rlv3d, rev0, rev0d, rev3, rev3d,& +& istart, jstart, isize, jsize + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype) :: vnd + real(kind=realtype), dimension(5) :: dww3 + real(kind=realtype), dimension(5) :: dww3d + intrinsic mod + real(kind=realtype) :: tempd + integer :: branch + dww3d = 0.0_8 +!$bwd-of ii-loop + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart + dww3(irho) = ww3(i, j, irho) - winf(irho) + dww3(ivx) = ww3(i, j, ivx) - winf(ivx) + dww3(ivy) = ww3(i, j, ivy) - winf(ivy) + dww3(ivz) = ww3(i, j, ivz) - winf(ivz) + dww3(irhoe) = ww3(i, j, irhoe) - winf(irhoe) + vn = two*(dww3(ivx)*bcdata(nn)%norm(i, j, 1)+dww3(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww3(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + if (viscous) then + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if + if (eddymodel) then + rev3d(i, j) = rev3d(i, j) - rev0d(i, j) + rev0d(i, j) = 0.0_8 + end if + call popcontrol1b(branch) + if (branch .eq. 0) then + rlv3d(i, j) = rlv3d(i, j) - rlv0d(i, j) + rlv0d(i, j) = 0.0_8 + end if + pinfd = pinfd + 2*pp0d(i, j) + pp3d(i, j) = pp3d(i, j) - pp0d(i, j) + pp0d(i, j) = 0.0_8 + winfd(irhoe) = winfd(irhoe) + ww0d(i, j, irhoe) + dww3d(irhoe) = dww3d(irhoe) - ww0d(i, j, irhoe) + ww0d(i, j, irhoe) = 0.0_8 + winfd(ivz) = winfd(ivz) + ww0d(i, j, ivz) + winfd(ivy) = winfd(ivy) + ww0d(i, j, ivy) + winfd(ivx) = winfd(ivx) + ww0d(i, j, ivx) + winfd(irho) = winfd(irho) + ww0d(i, j, irho) + dww3d(irho) = dww3d(irho) - ww0d(i, j, irho) + ww0d(i, j, irho) = 0.0_8 + vnd = bcdata(nn)%norm(i, j, 3)*ww0d(i, j, ivz) + bcdatad(nn)%norm(i, j, 3) = bcdatad(nn)%norm(i, j, 3) + vn*ww0d(i& +& , j, ivz) + dww3d(ivz) = dww3d(ivz) - ww0d(i, j, ivz) + ww0d(i, j, ivz) = 0.0_8 + vnd = vnd + bcdata(nn)%norm(i, j, 2)*ww0d(i, j, ivy) + bcdatad(nn)%norm(i, j, 2) = bcdatad(nn)%norm(i, j, 2) + vn*ww0d(i& +& , j, ivy) + dww3d(ivy) = dww3d(ivy) - ww0d(i, j, ivy) + ww0d(i, j, ivy) = 0.0_8 + vnd = vnd + bcdata(nn)%norm(i, j, 1)*ww0d(i, j, ivx) + tempd = two*vnd + bcdatad(nn)%norm(i, j, 1) = bcdatad(nn)%norm(i, j, 1) + vn*ww0d(i& +& , j, ivx) + dww3(ivx)*tempd + dww3d(ivx) = dww3d(ivx) + bcdata(nn)%norm(i, j, 1)*tempd - ww0d(i& +& , j, ivx) + ww0d(i, j, ivx) = 0.0_8 + dww3d(ivy) = dww3d(ivy) + bcdata(nn)%norm(i, j, 2)*tempd + bcdatad(nn)%norm(i, j, 2) = bcdatad(nn)%norm(i, j, 2) + dww3(ivy)*& +& tempd + dww3d(ivz) = dww3d(ivz) + bcdata(nn)%norm(i, j, 3)*tempd + bcdatad(nn)%norm(i, j, 3) = bcdatad(nn)%norm(i, j, 3) + dww3(ivz)*& +& tempd + ww3d(i, j, irhoe) = ww3d(i, j, irhoe) + dww3d(irhoe) + winfd(irhoe) = winfd(irhoe) - dww3d(irhoe) + dww3d(irhoe) = 0.0_8 + ww3d(i, j, ivz) = ww3d(i, j, ivz) + dww3d(ivz) + winfd(ivz) = winfd(ivz) - dww3d(ivz) + dww3d(ivz) = 0.0_8 + ww3d(i, j, ivy) = ww3d(i, j, ivy) + dww3d(ivy) + winfd(ivy) = winfd(ivy) - dww3d(ivy) + dww3d(ivy) = 0.0_8 + ww3d(i, j, ivx) = ww3d(i, j, ivx) + dww3d(ivx) + winfd(ivx) = winfd(ivx) - dww3d(ivx) + dww3d(ivx) = 0.0_8 + ww3d(i, j, irho) = ww3d(i, j, irho) + dww3d(irho) + winfd(irho) = winfd(irho) - dww3d(irho) + dww3d(irho) = 0.0_8 + end do + end subroutine bcantisymm2ndhalo_b + + subroutine bcantisymm2ndhalo(nn) +! bcsymm2ndhalo applies the symmetry boundary conditions to a +! block for the 2nd halo. this routine is separate as it makes +! ad slightly easier. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers_b, only : gamma0, gamma3, ww0, ww3, pp0, pp3, rlv0, & +& rlv3, rev0, rev3, istart, jstart, isize, jsize + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww3 + intrinsic mod +!$ad ii-loop +! if we need the second halo, do everything again, but using ww0, +! ww3 etc instead of ww2 and ww1. + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart + dww3(irho) = ww3(i, j, irho) - winf(irho) + dww3(ivx) = ww3(i, j, ivx) - winf(ivx) + dww3(ivy) = ww3(i, j, ivy) - winf(ivy) + dww3(ivz) = ww3(i, j, ivz) - winf(ivz) + dww3(irhoe) = ww3(i, j, irhoe) - winf(irhoe) + vn = two*(dww3(ivx)*bcdata(nn)%norm(i, j, 1)+dww3(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww3(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww0(i, j, ivx) = -dww3(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww0(i, j, ivy) = -dww3(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww0(i, j, ivz) = -dww3(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww0(i, j, irho) = -dww3(irho) + winf(irho) + ww0(i, j, ivx) = winf(ivx) + ww0(i, j, ivx) + ww0(i, j, ivy) = winf(ivy) + ww0(i, j, ivy) + ww0(i, j, ivz) = winf(ivz) + ww0(i, j, ivz) + ww0(i, j, irhoe) = -dww3(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma0(i, j) = gamma3(i, j) + pp0(i, j) = pinf - (pp3(i, j)-pinf) + if (viscous) rlv0(i, j) = -rlv3(i, j) + if (eddymodel) rev0(i, j) = -rev3(i, j) + end do + end subroutine bcantisymm2ndhalo + ! differentiation of bcsymmpolar1sthalo in reverse (adjoint) mode (with options noisize i4 dr8 r8): ! gradient of useful results: *xx *rev1 *rev2 *pp1 *pp2 *rlv1 ! *rlv2 *ww1 *ww2 diff --git a/src/adjoint/outputReverse/adjointExtra_b.f90 b/src/adjoint/outputReverse/adjointExtra_b.f90 index 84ee6e537..b93acc5ee 100644 --- a/src/adjoint/outputReverse/adjointExtra_b.f90 +++ b/src/adjoint/outputReverse/adjointExtra_b.f90 @@ -1279,7 +1279,7 @@ subroutine xhalo_block_b() loopbocos:do mm=1,nbocos ! the actual correction of the coordinates only takes ! place for symmetry planes. - if (bctype(mm) .eq. symm) then + if (bctype(mm) .eq. symm .or. bctype(mm) .eq. antisymm) then ! set some variables, depending on the block face on ! which the subface is located. call pushreal8(norm(1)) @@ -1781,7 +1781,7 @@ subroutine xhalo_block() loopbocos:do mm=1,nbocos ! the actual correction of the coordinates only takes ! place for symmetry planes. - if (bctype(mm) .eq. symm) then + if (bctype(mm) .eq. symm .or. bctype(mm) .eq. antisymm) then ! set some variables, depending on the block face on ! which the subface is located. norm(1) = bcdata(mm)%symnorm(1) diff --git a/src/adjoint/outputReverse/turbBCRoutines_b.f90 b/src/adjoint/outputReverse/turbBCRoutines_b.f90 index 589ef941c..8d9ec6d47 100644 --- a/src/adjoint/outputReverse/turbBCRoutines_b.f90 +++ b/src/adjoint/outputReverse/turbBCRoutines_b.f90 @@ -1195,7 +1195,7 @@ subroutine bcturbtreatment_b() select case (bctype(nn)) case (nswalladiabatic, nswallisothermal) call pushcontrol2b(2) - case (symm, symmpolar, eulerwall) + case (symm, antisymm, symmpolar, eulerwall) call pushcontrol2b(3) case (farfield) call pushcontrol2b(1) @@ -1316,7 +1316,7 @@ subroutine bcturbtreatment() call bcturbwall(nn) !============================================================= !============================================================= - case (symm, symmpolar, eulerwall) + case (symm, antisymm, symmpolar, eulerwall) ! symmetry, polar symmetry or inviscid wall. treatment of ! the turbulent equations is identical. call bcturbsymm(nn) diff --git a/src/adjoint/outputReverseFast/BCRoutines_fast_b.f90 b/src/adjoint/outputReverseFast/BCRoutines_fast_b.f90 index 37f0e540f..74f9fe1a4 100644 --- a/src/adjoint/outputReverseFast/BCRoutines_fast_b.f90 +++ b/src/adjoint/outputReverseFast/BCRoutines_fast_b.f90 @@ -52,6 +52,25 @@ subroutine applyallbc_block(secondhalo) end if !$ad ii-loop ! ------------------------------------ +! antisymmetry boundary condition +! ------------------------------------ + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm1sthalo(nn) + end if + end do + if (secondhalo) then +!$ad ii-loop + do nn=1,nbocos + if (bctype(nn) .eq. antisymm) then + call setbcpointers(nn, .false.) + call bcantisymm2ndhalo(nn) + end if + end do + end if +!$ad ii-loop +! ------------------------------------ ! symmetry polar boundary condition ! ------------------------------------ do nn=1,nbocos @@ -246,6 +265,112 @@ subroutine bcsymm2ndhalo(nn) end do end subroutine bcsymm2ndhalo + subroutine bcantisymm1sthalo(nn) +! bcantisymm1sthalo applies the antisymmetry boundary conditions to a +! block. * it is assumed that the pointers in blockpointers are +! already set to the correct block on the correct grid level. +! it should only be applied for high-speed near free-face operation, +! unless you are clear what you use it for +! +! in case also the second halo must be set, a second loop is +! execulted calling bcsymm2ndhalo. this is the only correct way +! in case the block contains only 1 cell between two symmetry +! planes, i.e. a 2d problem. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers, only : gamma1, gamma2, ww1, ww2, pp1, pp2, rlv1, & +& rlv2, istart, jstart, isize, jsize, rev1, rev2 + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww2 + intrinsic mod +!$ad ii-loop +! loop over the generic subface to set the state in the +! 1-st level halos + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart +! determine twice the normal velocity component, +! which must be substracted from the donor velocity +! to obtain the halo velocity. + dww2(irho) = ww2(i, j, irho) - winf(irho) + dww2(ivx) = ww2(i, j, ivx) - winf(ivx) + dww2(ivy) = ww2(i, j, ivy) - winf(ivy) + dww2(ivz) = ww2(i, j, ivz) - winf(ivz) + dww2(irhoe) = ww2(i, j, irhoe) - winf(irhoe) + vn = two*(dww2(ivx)*bcdata(nn)%norm(i, j, 1)+dww2(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww2(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww1(i, j, ivx) = -dww2(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww1(i, j, ivy) = -dww2(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww1(i, j, ivz) = -dww2(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww1(i, j, irho) = -dww2(irho) + winf(irho) + ww1(i, j, ivx) = winf(ivx) + ww1(i, j, ivx) + ww1(i, j, ivy) = winf(ivy) + ww1(i, j, ivy) + ww1(i, j, ivz) = winf(ivz) + ww1(i, j, ivz) + ww1(i, j, irhoe) = -dww2(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma1(i, j) = gamma2(i, j) + pp1(i, j) = pinf - (pp2(i, j)-pinf) + if (viscous) rlv1(i, j) = -rlv2(i, j) + if (eddymodel) rev1(i, j) = -rev2(i, j) + end do + end subroutine bcantisymm1sthalo + + subroutine bcantisymm2ndhalo(nn) +! bcsymm2ndhalo applies the symmetry boundary conditions to a +! block for the 2nd halo. this routine is separate as it makes +! ad slightly easier. + use constants + use blockpointers, only : bcdata + use flowvarrefstate, only : viscous, eddymodel, pinf, winf + use bcpointers, only : gamma0, gamma3, ww0, ww3, pp0, pp3, rlv0, & +& rlv3, rev0, rev3, istart, jstart, isize, jsize + implicit none +! subroutine arguments. + integer(kind=inttype), intent(in) :: nn +! local variables. + integer(kind=inttype) :: i, j, l, ii + real(kind=realtype) :: vn, nnx, nny, nnz + real(kind=realtype), dimension(5) :: dww3 + intrinsic mod +!$ad ii-loop +! if we need the second halo, do everything again, but using ww0, +! ww3 etc instead of ww2 and ww1. + do ii=0,isize*jsize-1 + i = mod(ii, isize) + istart + j = ii/isize + jstart + dww3(irho) = ww3(i, j, irho) - winf(irho) + dww3(ivx) = ww3(i, j, ivx) - winf(ivx) + dww3(ivy) = ww3(i, j, ivy) - winf(ivy) + dww3(ivz) = ww3(i, j, ivz) - winf(ivz) + dww3(irhoe) = ww3(i, j, irhoe) - winf(irhoe) + vn = two*(dww3(ivx)*bcdata(nn)%norm(i, j, 1)+dww3(ivy)*bcdata(nn)%& +& norm(i, j, 2)+dww3(ivz)*bcdata(nn)%norm(i, j, 3)) +! determine the flow variables in the halo cell. + ww0(i, j, ivx) = -dww3(ivx) + vn*bcdata(nn)%norm(i, j, 1) + ww0(i, j, ivy) = -dww3(ivy) + vn*bcdata(nn)%norm(i, j, 2) + ww0(i, j, ivz) = -dww3(ivz) + vn*bcdata(nn)%norm(i, j, 3) + ww0(i, j, irho) = -dww3(irho) + winf(irho) + ww0(i, j, ivx) = winf(ivx) + ww0(i, j, ivx) + ww0(i, j, ivy) = winf(ivy) + ww0(i, j, ivy) + ww0(i, j, ivz) = winf(ivz) + ww0(i, j, ivz) + ww0(i, j, irhoe) = -dww3(irhoe) + winf(irhoe) +! set the pressure and gamma and possibly the +! laminar and eddy viscosity in the halo. + gamma0(i, j) = gamma3(i, j) + pp0(i, j) = pinf - (pp3(i, j)-pinf) + if (viscous) rlv0(i, j) = -rlv3(i, j) + if (eddymodel) rev0(i, j) = -rev3(i, j) + end do + end subroutine bcantisymm2ndhalo + subroutine bcsymmpolar1sthalo(nn) ! bcsymmpolar applies the polar symmetry boundary conditions to a ! singular line of a block. it is assumed that the pointers in diff --git a/src/adjoint/outputReverseFast/turbBCRoutines_fast_b.f90 b/src/adjoint/outputReverseFast/turbBCRoutines_fast_b.f90 index aacda3694..2daba532d 100644 --- a/src/adjoint/outputReverseFast/turbBCRoutines_fast_b.f90 +++ b/src/adjoint/outputReverseFast/turbBCRoutines_fast_b.f90 @@ -589,7 +589,7 @@ subroutine bcturbtreatment() call bcturbwall(nn) !============================================================= !============================================================= - case (symm, symmpolar, eulerwall) + case (symm, antisymm, symmpolar, eulerwall) ! symmetry, polar symmetry or inviscid wall. treatment of ! the turbulent equations is identical. call bcturbsymm(nn) diff --git a/src/bcdata/BCData.F90 b/src/bcdata/BCData.F90 index f8ebb836d..412df983c 100644 --- a/src/bcdata/BCData.F90 +++ b/src/bcdata/BCData.F90 @@ -2315,7 +2315,7 @@ subroutine allocMemBCData &a farfield") !======================================================= - case (symm, symmPolar) + case (symm, antiSymm, symmPolar) ! Allocate for symm as well. This is not necessary ! but we need it for the reverse AD. diff --git a/src/modules/constants.F90 b/src/modules/constants.F90 index 086240822..d70ed0255 100644 --- a/src/modules/constants.F90 +++ b/src/modules/constants.F90 @@ -276,7 +276,8 @@ module constants integer(kind=intType), parameter :: DomainInterfaceP = -22 integer(kind=intType), parameter :: DomainInterfaceRho = -23 integer(kind=intType), parameter :: DomainInterfaceTotal = -24 - integer(kind=intType), parameter :: BCNotValid = -25 + integer(kind=intType), parameter :: AntiSymm = -25 + integer(kind=intType), parameter :: BCNotValid = -26 ! ! Number of actual boundary conditions supported by the code ! This number refers to bocos, not flow-through BCs diff --git a/src/output/writeCGNSSurface.F90 b/src/output/writeCGNSSurface.F90 index 8a6bcf3b6..5fa8bc312 100644 --- a/src/output/writeCGNSSurface.F90 +++ b/src/output/writeCGNSSurface.F90 @@ -916,6 +916,8 @@ subroutine createSurfaceZone select case (cgnsDoms(zone)%bocoInfo(subface)%BCType) case (Symm) zonename = "Symmetry" + case (AntiSymm) + zonename = "AntiSymmetry" case (SymmPolar) zonename = "SymmetryPolar" case (NSWallAdiabatic) diff --git a/src/partitioning/loadBalance.F90 b/src/partitioning/loadBalance.F90 index 6b1fe989c..46394e8cf 100644 --- a/src/partitioning/loadBalance.F90 +++ b/src/partitioning/loadBalance.F90 @@ -1262,59 +1262,63 @@ subroutine sortSubfaces(oldSubfaceID, blockID) case (Symm) bcPrior(i) = 4 - case (SymmPolar) + case (AntiSymm) bcPrior(i) = 5 - case (FarField) + case (SymmPolar) bcPrior(i) = 6 - case (SupersonicInflow) + case (FarField) bcPrior(i) = 7 - case (SubsonicInflow) + case (SupersonicInflow) bcPrior(i) = 8 - case (SupersonicOutflow) + case (SubsonicInflow) bcPrior(i) = 9 - case (SubsonicOutflow) + case (SupersonicOutflow) bcPrior(i) = 10 - case (MassBleedInflow) + case (SubsonicOutflow) bcPrior(i) = 11 - case (MassBleedOutflow) + case (MassBleedInflow) bcPrior(i) = 12 - case (mDot) + case (MassBleedOutflow) bcPrior(i) = 13 - case (bcThrust) + case (mDot) bcPrior(i) = 14 - case (Extrap) + case (bcThrust) bcPrior(i) = 15 + case (Extrap) + bcPrior(i) = 16 + + ! TODO: original code was also skipping 4 values. check why case (SlidingInterface) - bcPrior(i) = 19 + bcPrior(i) = 20 case (OversetOuterBound) - bcPrior(i) = 20 + bcPrior(i) = 21 case (DomainInterfaceAll) - bcPrior(i) = 21 + bcPrior(i) = 22 case (DomainInterfaceRhoUVW) - bcPrior(i) = 22 + bcPrior(i) = 23 case (DomainInterfaceP) - bcPrior(i) = 23 + bcPrior(i) = 24 case (DomainInterfaceRho) - bcPrior(i) = 24 + bcPrior(i) = 25 case (DomainInterfaceTotal) - bcPrior(i) = 25 + bcPrior(i) = 26 end select diff --git a/src/partitioning/readCGNSGrid.F90 b/src/partitioning/readCGNSGrid.F90 index be149fc2d..d11657ae7 100644 --- a/src/partitioning/readCGNSGrid.F90 +++ b/src/partitioning/readCGNSGrid.F90 @@ -2721,6 +2721,9 @@ function internalBC(cgnsBocoType, userDefinedName) case ("BCDomainInterfaceTotal") internalBC = DomainInterfaceTotal + case ("BCAntiSymmetry") + internalBC = AntiSymm + case default internalBC = BCNull diff --git a/src/preprocessing/preprocessingAPI.F90 b/src/preprocessing/preprocessingAPI.F90 index 6516357ec..591c73c9a 100644 --- a/src/preprocessing/preprocessingAPI.F90 +++ b/src/preprocessing/preprocessingAPI.F90 @@ -904,7 +904,7 @@ subroutine checkSymmetry(level) ! Check for symmetry boundary condition. - symmetry: if (BCType(mm) == symm) then + symmetry: if (BCType(mm) == symm .or. BCType(mm) == antisymm) then ! Determine the block face on which this subface is ! located and set some variables accordingly. @@ -1139,7 +1139,7 @@ subroutine xhalo(level) ! The actual correction of the coordinates only takes ! place for symmetry planes. - testSymmetry: if (BCType(mm) == Symm) then + testSymmetry: if (BCType(mm) == Symm .or. BCType(mm) == antisymm) then ! Set some variables, depending on the block face on ! which the subface is located. @@ -1453,7 +1453,7 @@ subroutine setSurfaceFamilyInfo ! scatterd based on groups of BC types. Specifically the following groups: ! 1. Walls : EulerWall, NSWallAdiabatic, NSWallIsothermal - ! 2. Symm : Symm, SymmPolar + ! 2. Symm : Symm, AntiSymm, SymmPolar ! 3. Inflow/Outflow : subSonicInflow, subSonicOutflow, supersonicInflow, superSonicOutflow ! 4. Farfield : Farfield ! 5. Overset : OversetouterBound @@ -1494,7 +1494,7 @@ subroutine setSurfaceFamilyInfo end if case (iBCGroupSymm) - if (BCType(mm) == Symm .or. BCType(mm) == SymmPolar) then + if (BCType(mm) == Symm .or. BCType(mm) == SymmPolar .or. BCType(mm) == AntiSymm) then localFlag(iFam) = 1 end if diff --git a/src/solver/BCRoutines.F90 b/src/solver/BCRoutines.F90 index 9ce054802..25d347011 100644 --- a/src/solver/BCRoutines.F90 +++ b/src/solver/BCRoutines.F90 @@ -102,6 +102,27 @@ subroutine applyAllBC_block(secondHalo) end do end if + ! ------------------------------------ + ! AntiSymmetry Boundary Condition + ! ------------------------------------ + !$AD II-LOOP + do nn = 1, nBocos + if (BCType(nn) == antiSymm) then + call setBCPointers(nn, .False.) + call bcAntiSymm1stHalo(nn) + end if + end do + + if (secondHalo) then + !$AD II-LOOP + do nn = 1, nBocos + if (BCType(nn) == antiSymm) then + call setBCPointers(nn, .False.) + call bcAntiSymm2ndHalo(nn) + end if + end do + end if + ! ------------------------------------ ! Symmetry Polar Boundary Condition ! ------------------------------------ @@ -329,6 +350,138 @@ subroutine bcSymm2ndHalo(nn) end subroutine bcSymm2ndHalo + subroutine bcAntiSymm1stHalo(nn) + + ! bcAntiSymm1stHalo applies the antisymmetry boundary conditions to a + ! block. * It is assumed that the pointers in blockPointers are + ! already set to the correct block on the correct grid level. + ! It should only be applied for high-speed near free-face operation, + ! unless you are clear what you use it for + ! + ! In case also the second halo must be set, a second loop is + ! execulted calling bcSymm2ndhalo. This is the only correct way + ! in case the block contains only 1 cell between two symmetry + ! planes, i.e. a 2D problem. + + use constants + use blockPointers, only: BCdata + use flowVarRefState, only: viscous, eddyModel, pinf, wInf + use BCPointers, only: gamma1, gamma2, ww1, ww2, pp1, pp2, rlv1, rlv2, & + iStart, jStart, iSize, jSize, rev1, rev2 + implicit none + + ! Subroutine arguments. + integer(kind=intType), intent(in) :: nn + + ! Local variables. + integer(kind=intType) :: i, j, l, ii + real(kind=realType) :: vn, nnx, nny, nnz + real(kind=realType), dimension(5) :: dww2 + + ! Loop over the generic subface to set the state in the + ! 1-st level halos + + !$AD II-LOOP + do ii = 0, isize * jsize - 1 + i = mod(ii, isize) + iStart + j = ii / isize + jStart + + ! Determine twice the normal velocity component, + ! which must be substracted from the donor velocity + ! to obtain the halo velocity. + + dww2(irho) = ww2(i, j, irho) - wInf(irho) + dww2(ivx) = ww2(i, j, ivx) - wInf(ivx) + dww2(ivy) = ww2(i, j, ivy) - wInf(ivy) + dww2(ivz) = ww2(i, j, ivz) - wInf(ivz) + dww2(irhoE) = ww2(i, j, irhoE) - wInf(irhoE) + + vn = two * (dww2(ivx) * BCData(nn)%norm(i, j, 1) + & + dww2(ivy) * BCData(nn)%norm(i, j, 2) + & + dww2(ivz) * BCData(nn)%norm(i, j, 3)) + + ! Determine the flow variables in the halo cell. + + ww1(i, j, ivx) = -dww2(ivx) + vn * BCData(nn)%norm(i, j, 1) + ww1(i, j, ivy) = -dww2(ivy) + vn * BCData(nn)%norm(i, j, 2) + ww1(i, j, ivz) = -dww2(ivz) + vn * BCData(nn)%norm(i, j, 3) + + ww1(i, j, irho) = -dww2(irho) + wInf(irho) + ww1(i, j, ivx) = wInf(ivx) + ww1(i, j, ivx) + ww1(i, j, ivy) = wInf(ivy) + ww1(i, j, ivy) + ww1(i, j, ivz) = wInf(ivz) + ww1(i, j, ivz) + ww1(i, j, irhoE) = -dww2(irhoE) + wInf(irhoE) + + ! Set the pressure and gamma and possibly the + ! laminar and eddy viscosity in the halo. + + gamma1(i, j) = gamma2(i, j) + pp1(i, j) = pInf - (pp2(i, j) - pInf) + if (viscous) rlv1(i, j) = -rlv2(i, j) + if (eddyModel) rev1(i, j) = -rev2(i, j) + end do + end subroutine bcAntiSymm1stHalo + + subroutine bcAntiSymm2ndHalo(nn) + + ! bcSymm2ndHalo applies the symmetry boundary conditions to a + ! block for the 2nd halo. This routine is separate as it makes + ! AD slightly easier. + use constants + use blockPointers, only: BCdata + use flowVarRefState, only: viscous, eddyModel, pinf, wInf + use BCPointers, only: gamma0, gamma3, ww0, ww3, pp0, pp3, rlv0, rlv3, & + rev0, rev3, iStart, jStart, iSize, jSize + implicit none + + ! Subroutine arguments. + integer(kind=intType), intent(in) :: nn + + ! Local variables. + integer(kind=intType) :: i, j, l, ii + real(kind=realType) :: vn, nnx, nny, nnz + real(kind=realType), dimension(5) :: dww3 + + ! If we need the second halo, do everything again, but using ww0, + ! ww3 etc instead of ww2 and ww1. + + !$AD II-LOOP + do ii = 0, isize * jsize - 1 + i = mod(ii, isize) + iStart + j = ii / isize + jStart + + dww3(irho) = ww3(i, j, irho) - wInf(irho) + dww3(ivx) = ww3(i, j, ivx) - wInf(ivx) + dww3(ivy) = ww3(i, j, ivy) - wInf(ivy) + dww3(ivz) = ww3(i, j, ivz) - wInf(ivz) + dww3(irhoE) = ww3(i, j, irhoE) - wInf(irhoE) + + vn = two * (dww3(ivx) * BCData(nn)%norm(i, j, 1) + & + dww3(ivy) * BCData(nn)%norm(i, j, 2) + & + dww3(ivz) * BCData(nn)%norm(i, j, 3)) + + ! Determine the flow variables in the halo cell. + ww0(i, j, ivx) = -dww3(ivx) + vn * BCData(nn)%norm(i, j, 1) + ww0(i, j, ivy) = -dww3(ivy) + vn * BCData(nn)%norm(i, j, 2) + ww0(i, j, ivz) = -dww3(ivz) + vn * BCData(nn)%norm(i, j, 3) + + ww0(i, j, irho) = -dww3(irho) + wInf(irho) + ww0(i, j, ivx) = wInf(ivx) + ww0(i, j, ivx) + ww0(i, j, ivy) = wInf(ivy) + ww0(i, j, ivy) + ww0(i, j, ivz) = wInf(ivz) + ww0(i, j, ivz) + ww0(i, j, irhoE) = -dww3(irhoE) + wInf(irhoE) + + ! Set the pressure and gamma and possibly the + ! laminar and eddy viscosity in the halo. + + gamma0(i, j) = gamma3(i, j) + pp0(i, j) = pInf - (pp3(i, j) - pInf) + if (viscous) rlv0(i, j) = -rlv3(i, j) + if (eddyModel) rev0(i, j) = -rev3(i, j) + end do + + end subroutine bcAntiSymm2ndHalo + subroutine bcSymmPolar1stHalo(nn) ! bcSymmPolar applies the polar symmetry boundary conditions to a diff --git a/src/turbulence/turbBCRoutines.F90 b/src/turbulence/turbBCRoutines.F90 index 24c7cadc9..4f1328fe0 100644 --- a/src/turbulence/turbBCRoutines.F90 +++ b/src/turbulence/turbBCRoutines.F90 @@ -762,7 +762,7 @@ subroutine bcTurbTreatment #endif !============================================================= - case (Symm, SymmPolar, EulerWall) + case (Symm, AntiSymm, SymmPolar, EulerWall) ! Symmetry, polar symmetry or inviscid wall. Treatment of ! the turbulent equations is identical. diff --git a/src/utils/utils.F90 b/src/utils/utils.F90 index b50cc8c17..9e10395a1 100644 --- a/src/utils/utils.F90 +++ b/src/utils/utils.F90 @@ -4029,7 +4029,7 @@ subroutine getLiftDirFromSymmetry(liftDir) do nn = 1, nDom call setPointers(nn, 1, 1) do mm = 1, nBocos - if (bcType(mm) == symm) then + if (bcType(mm) == symm .or. bcType(mm) == antisymm) then select case (BCFaceID(mm)) case (iMin)