Antisymmetry BC with surface-piercing topology

src/adjoint/Makefile_tapenade

@@ -134,6 +134,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) \
 \

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.

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
     ! ------------------------------------

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