Update tapenade 3.16

.github/azure-pipelines.yaml

@@ -20,3 +20,4 @@ extends:
     GCC_CONFIG: config/defaults/config.LINUX_GFORTRAN_OPENMPI.mk
     INTEL_CONFIG: config/defaults/config.LINUX_INTEL_OPENMPI.mk
     TAPENADE: true
+    TAPENADE_VERSION: "3.16"

src/adjoint/Makefile_tapenade

@@ -5,10 +5,8 @@ PREPROCESS_DIR = preprocess_temp
 TAPENADE_SRC =	$(TAPENADE_DIRSRC)/warp/getElementProps.F90\
 		$(TAPENADE_DIRSRC)/utils/vectorUtils.f90\
 		$(TAPENADE_DIRSRC)/modules/gridData.F90\
-		$(TAPENADE_DIRSRC)/modules/gridInput.f90\
 		$(TAPENADE_DIRSRC)/modules/constants.F90\
 		$(TAPENADE_DIRSRC)/modules/precision.F90\
-		$(TAPENADE_DIRSRC)/modules/communication.f90\
 		$(TAPENADE_DIRSRC)/modules/kd_tree.F90\
 		$(TAPENADE_DIRSRC)/warp/warpMesh.F90\
 
@@ -22,7 +20,7 @@ preprocess:
 	@echo "Processing input files for AD..."
 	mkdir -p $(PREPROCESS_DIR)
 	@for file in $(TAPENADE_SRC); do \
-		echo Preprocssing $$file; \
+		echo Preprocessing $$file; \
 		cpp -DUSE_TAPENADE -traditional -P $$file $(PREPROCESS_DIR)/"$$(basename $$file)"; \
 	done
 
@@ -32,7 +30,7 @@ ad_reverse:
 	mkdir -p temp_reverse
 	make -f Makefile_tapenade preprocess
 
-# the following ist he actual tapenade command to run
+# the following is the actual tapenade command to run
 	$(TAPENADE_HOME)/bin/tapenade \
 	-head "kd_tree%computeNodalProperties(XsPtr, tp%Mi, tp%Bi)>(tp%Mi, tp%Bi)" \
 	-noisize \
@@ -59,5 +57,3 @@ ad_forward:
 
 # Clean up files
 	make -f Makefile_tapenade cleanpreprocess
-
-

src/adjoint/autoEdit/autoEditForward.py

@@ -33,7 +33,7 @@
 STR_REPLACE_ALL = {"_CD": ""}
 # Also, entirely ignore lines with these strings:
 LINE_IGNORE = []
-FILE_IGNORE = ["griddata_d.F90"]
+FILE_IGNORE = ["gridData_d.f90", "kd_tree_d.f90", "warpMesh_d.f90"]
 # Some feedback
 
 print("Directory of input source files  :", DIR_ORI)

src/adjoint/autoEdit/autoEditReverse.py

@@ -33,7 +33,7 @@
 STR_REPLACE_ALL = {"_CB": ""}
 # Also, entirely ignore lines with these strings:
 LINE_IGNORE = []
-FILE_IGNORE = ["griddata_b.F90"]
+FILE_IGNORE = ["gridData_b.f90", "kd_tree_b.f90", "warpMesh_b.f90"]
 # Some feedback
 
 print("Directory of input source files  :", DIR_ORI)

src/adjoint/outputForward/Makefile

@@ -11,10 +11,8 @@ include ${COMMON_FILE}
 include ${RULES_FILE}
 vpath %.o $(OBJDIR)
 
-FF90_OBJECTS =	cross_product_3d_d.o\
-		getelementprops_d.o\
-		getmag_d.o\
-		getrotationmatrix3d_d.o
+FF90_OBJECTS =	vectorUtils_d.o\
+		getElementProps_d.o
 
 all:
 	$(PMAKE)  src

src/adjoint/outputForward/getelementprops_d.f90 → src/adjoint/outputForward/getElementProps_d.f90

@@ -1,5 +1,5 @@
-   !        Generated by TAPENADE     (INRIA, Tropics team)
-   !  Tapenade 3.10 (r5363) -  9 Sep 2014 09:53
+   !        Generated by TAPENADE     (INRIA, Ecuador team)
+   !  Tapenade 3.16 (master) -  9 Oct 2020 17:47
    !
    !  Differentiation of getelementprops in forward (tangent) mode (with options noISIZE i4 dr8 r8):
    !   variations   of useful results: area normal
@@ -29,6 +29,7 @@ SUBROUTINE GETELEMENTPROPS_D(points, pointsb, npts, area, areab, normal&
    REAL(kind=realtype) :: arg1b
    REAL(kind=realtype) :: result1
    REAL(kind=realtype) :: result1b
+   REAL(kind=realtype) :: temp
    ! Compute the element center:
    center(:) = zero
    centerb = 0.0_8
@@ -64,28 +65,31 @@ SUBROUTINE GETELEMENTPROPS_D(points, pointsb, npts, area, areab, normal&
    arg1b = 2*cross(1)*crossb(1) + 2*cross(2)*crossb(2) + 2*cross(3)*&
    &     crossb(3)
    arg1 = cross(1)**2 + cross(2)**2 + cross(3)**2 + 1e-15
+   temp = SQRT(arg1)
    IF (arg1 .EQ. 0.0_8) THEN
    crossnormb(1) = 0.0_8
    ELSE
-   crossnormb(1) = arg1b/(2.0*SQRT(arg1))
+   crossnormb(1) = arg1b/(2.0*temp)
    END IF
-   crossnorm(1) = SQRT(arg1)
+   crossnorm(1) = temp
    areab = areab + half*crossnormb(1)
    area = area + half*crossnorm(1)
-   normalb = normalb + (crossb*crossnorm(1)-cross*crossnormb(1))/&
-   &     crossnorm(1)**2
+   normalb = normalb + (crossb-cross*crossnormb(1)/crossnorm(1))/&
+   &     crossnorm(1)
    normal = normal + cross/crossnorm(1)
    END DO
    ! Also make the normal a unit lenth
    arg1b = 2*normal(1)*normalb(1) + 2*normal(2)*normalb(2) + 2*normal(3)*&
    &   normalb(3)
    arg1 = normal(1)**2 + normal(2)**2 + normal(3)**2
+   temp = SQRT(arg1)
    IF (arg1 .EQ. 0.0_8) THEN
    result1b = 0.0_8
    ELSE
-   result1b = arg1b/(2.0*SQRT(arg1))
+   result1b = arg1b/(2.0*temp)
    END IF
-   result1 = SQRT(arg1)
-   normalb = (normalb*result1-normal*result1b)/result1**2
+   result1 = temp
+   normalb = (normalb-normal*result1b/result1)/result1
    normal = normal/result1
    END SUBROUTINE GETELEMENTPROPS_D
+

src/adjoint/outputForward/getrotationmatrix3d_d.f90 → src/adjoint/outputForward/vectorUtils_d.f90

@@ -1,5 +1,5 @@
-   !        Generated by TAPENADE     (INRIA, Tropics team)
-   !  Tapenade 3.10 (r5363) -  9 Sep 2014 09:53
+   !        Generated by TAPENADE     (INRIA, Ecuador team)
+   !  Tapenade 3.16 (master) -  9 Oct 2020 17:47
    !
    !  Differentiation of getrotationmatrix3d in forward (tangent) mode (with options noISIZE i4 dr8 r8):
    !   variations   of useful results: mi
@@ -25,6 +25,8 @@ SUBROUTINE GETROTATIONMATRIX3D_D(v1, v2, v2b, mi, mib)
    INTRINSIC ACOS
    INTRINSIC SIN
    INTRINSIC COS
+   REAL(kind=realtype) :: temp
+   REAL(kind=realtype) :: temp0
    REAL(kind=realtype) :: v1b(3)
    CALL GETMAG(v1, magv1)
    CALL GETMAG_D(v2, v2b, magv2, magv2b)
@@ -46,11 +48,11 @@ SUBROUTINE GETROTATIONMATRIX3D_D(v1, v2, v2b, mi, mib)
    axisb = 0.0_8
    angleb = 0.0_8
    ELSE
-   axisb = (axisb*axismag-axis*axismagb)/axismag**2
+   axisb = (axisb-axis*axismagb/axismag)/axismag
    axis = axis/axismag
    ! Now compute the rotation angle about that axis
    vv1 = v1/magv1
-   vv2b = (v2b*magv2-v2*magv2b)/magv2**2
+   vv2b = (v2b-v2*magv2b/magv2)/magv2
    vv2 = v2/magv2
    IF (one .GT. vv1(1)*vv2(1) + vv1(2)*vv2(2) + vv1(3)*vv2(3)) THEN
    argb = vv1(1)*vv2b(1) + vv1(2)*vv2b(2) + vv1(3)*vv2b(3)
@@ -88,39 +90,88 @@ SUBROUTINE GETROTATIONMATRIX3D_D(v1, v2, v2b, mi, mib)
    a(3, 3) = zero
    !C = A*A
    cb = 0.0_8
-   cb(1, 1) = ab(1, 1)*a(1, 1) + a(1, 1)*ab(1, 1) + ab(1, 2)*a(2, 1) + a(&
-   &   1, 2)*ab(2, 1) + ab(1, 3)*a(3, 1) + a(1, 3)*ab(3, 1)
+   cb(1, 1) = 2*a(1, 1)*ab(1, 1) + a(2, 1)*ab(1, 2) + a(1, 2)*ab(2, 1) + &
+   &   a(3, 1)*ab(1, 3) + a(1, 3)*ab(3, 1)
    c(1, 1) = a(1, 1)*a(1, 1) + a(1, 2)*a(2, 1) + a(1, 3)*a(3, 1)
-   cb(1, 2) = ab(1, 1)*a(1, 2) + a(1, 1)*ab(1, 2) + ab(1, 2)*a(2, 2) + a(&
-   &   1, 2)*ab(2, 2) + ab(1, 3)*a(3, 2) + a(1, 3)*ab(3, 2)
+   cb(1, 2) = a(1, 2)*ab(1, 1) + (a(1, 1)+a(2, 2))*ab(1, 2) + a(1, 2)*ab(&
+   &   2, 2) + a(3, 2)*ab(1, 3) + a(1, 3)*ab(3, 2)
    c(1, 2) = a(1, 1)*a(1, 2) + a(1, 2)*a(2, 2) + a(1, 3)*a(3, 2)
-   cb(1, 3) = ab(1, 1)*a(1, 3) + a(1, 1)*ab(1, 3) + ab(1, 2)*a(2, 3) + a(&
-   &   1, 2)*ab(2, 3) + ab(1, 3)*a(3, 3) + a(1, 3)*ab(3, 3)
+   cb(1, 3) = a(1, 3)*ab(1, 1) + (a(1, 1)+a(3, 3))*ab(1, 3) + a(2, 3)*ab(&
+   &   1, 2) + a(1, 2)*ab(2, 3) + a(1, 3)*ab(3, 3)
    c(1, 3) = a(1, 1)*a(1, 3) + a(1, 2)*a(2, 3) + a(1, 3)*a(3, 3)
-   cb(2, 1) = ab(2, 1)*a(1, 1) + a(2, 1)*ab(1, 1) + ab(2, 2)*a(2, 1) + a(&
-   &   2, 2)*ab(2, 1) + ab(2, 3)*a(3, 1) + a(2, 3)*ab(3, 1)
+   cb(2, 1) = (a(1, 1)+a(2, 2))*ab(2, 1) + a(2, 1)*ab(1, 1) + a(2, 1)*ab(&
+   &   2, 2) + a(3, 1)*ab(2, 3) + a(2, 3)*ab(3, 1)
    c(2, 1) = a(2, 1)*a(1, 1) + a(2, 2)*a(2, 1) + a(2, 3)*a(3, 1)
-   cb(2, 2) = ab(2, 1)*a(1, 2) + a(2, 1)*ab(1, 2) + ab(2, 2)*a(2, 2) + a(&
-   &   2, 2)*ab(2, 2) + ab(2, 3)*a(3, 2) + a(2, 3)*ab(3, 2)
+   cb(2, 2) = a(1, 2)*ab(2, 1) + a(2, 1)*ab(1, 2) + 2*a(2, 2)*ab(2, 2) + &
+   &   a(3, 2)*ab(2, 3) + a(2, 3)*ab(3, 2)
    c(2, 2) = a(2, 1)*a(1, 2) + a(2, 2)*a(2, 2) + a(2, 3)*a(3, 2)
-   cb(2, 3) = ab(2, 1)*a(1, 3) + a(2, 1)*ab(1, 3) + ab(2, 2)*a(2, 3) + a(&
-   &   2, 2)*ab(2, 3) + ab(2, 3)*a(3, 3) + a(2, 3)*ab(3, 3)
+   cb(2, 3) = a(1, 3)*ab(2, 1) + a(2, 1)*ab(1, 3) + a(2, 3)*ab(2, 2) + (a&
+   &   (2, 2)+a(3, 3))*ab(2, 3) + a(2, 3)*ab(3, 3)
    c(2, 3) = a(2, 1)*a(1, 3) + a(2, 2)*a(2, 3) + a(2, 3)*a(3, 3)
-   cb(3, 1) = ab(3, 1)*a(1, 1) + a(3, 1)*ab(1, 1) + ab(3, 2)*a(2, 1) + a(&
-   &   3, 2)*ab(2, 1) + ab(3, 3)*a(3, 1) + a(3, 3)*ab(3, 1)
+   cb(3, 1) = (a(1, 1)+a(3, 3))*ab(3, 1) + a(3, 1)*ab(1, 1) + a(2, 1)*ab(&
+   &   3, 2) + a(3, 2)*ab(2, 1) + a(3, 1)*ab(3, 3)
    c(3, 1) = a(3, 1)*a(1, 1) + a(3, 2)*a(2, 1) + a(3, 3)*a(3, 1)
-   cb(3, 2) = ab(3, 1)*a(1, 2) + a(3, 1)*ab(1, 2) + ab(3, 2)*a(2, 2) + a(&
-   &   3, 2)*ab(2, 2) + ab(3, 3)*a(3, 2) + a(3, 3)*ab(3, 2)
+   cb(3, 2) = a(1, 2)*ab(3, 1) + a(3, 1)*ab(1, 2) + (a(2, 2)+a(3, 3))*ab(&
+   &   3, 2) + a(3, 2)*ab(2, 2) + a(3, 2)*ab(3, 3)
    c(3, 2) = a(3, 1)*a(1, 2) + a(3, 2)*a(2, 2) + a(3, 3)*a(3, 2)
-   cb(3, 3) = ab(3, 1)*a(1, 3) + a(3, 1)*ab(1, 3) + ab(3, 2)*a(2, 3) + a(&
-   &   3, 2)*ab(2, 3) + ab(3, 3)*a(3, 3) + a(3, 3)*ab(3, 3)
+   cb(3, 3) = a(1, 3)*ab(3, 1) + a(3, 1)*ab(1, 3) + a(2, 3)*ab(3, 2) + a(&
+   &   3, 2)*ab(2, 3) + 2*a(3, 3)*ab(3, 3)
    c(3, 3) = a(3, 1)*a(1, 3) + a(3, 2)*a(2, 3) + a(3, 3)*a(3, 3)
    ! Rodrigues formula for the rotation matrix 
    mi = zero
    mi(1, 1) = one
    mi(2, 2) = one
    mi(3, 3) = one
-   mib = angleb*COS(angle)*a + SIN(angle)*ab + angleb*SIN(angle)*c + (one&
-   &   -COS(angle))*cb
-   mi = mi + SIN(angle)*a + (one-COS(angle))*c
+   temp = SIN(angle)
+   temp0 = one - COS(angle)
+   mib = (a*COS(angle)+c*SIN(angle))*angleb + temp*ab + temp0*cb
+   mi = mi + temp*a + temp0*c
    END SUBROUTINE GETROTATIONMATRIX3D_D
+      !  Differentiation of cross_product_3d in forward (tangent) mode (with options noISIZE i4 dr8 r8):
+   !   variations   of useful results: cross
+   !   with respect to varying inputs: cross v1 v2
+   ! ====================================================================
+   ! File: vectorUtils.f90
+   ! Author: C.A.(Sandy) Mader
+   ! Date Started: July 14, 2014
+   ! Date Modified:
+   SUBROUTINE CROSS_PRODUCT_3D_D(v1, v1b, v2, v2b, cross, crossb)
+   USE CONSTANTS
+   IMPLICIT NONE
+   REAL(kind=realtype), DIMENSION(3), INTENT(IN) :: v1, v2
+   REAL(kind=realtype), DIMENSION(3), INTENT(IN) :: v1b, v2b
+   REAL(kind=realtype), DIMENSION(3), INTENT(OUT) :: cross
+   REAL(kind=realtype), DIMENSION(3), INTENT(OUT) :: crossb
+   crossb(1) = v2(3)*v1b(2) + v1(2)*v2b(3) - v2(2)*v1b(3) - v1(3)*v2b(2)
+   cross(1) = v1(2)*v2(3) - v1(3)*v2(2)
+   crossb(2) = v2(1)*v1b(3) + v1(3)*v2b(1) - v2(3)*v1b(1) - v1(1)*v2b(3)
+   cross(2) = v1(3)*v2(1) - v1(1)*v2(3)
+   crossb(3) = v2(2)*v1b(1) + v1(1)*v2b(2) - v2(1)*v1b(2) - v1(2)*v2b(1)
+   cross(3) = v1(1)*v2(2) - v1(2)*v2(1)
+   END SUBROUTINE CROSS_PRODUCT_3D_D
+      !  Differentiation of getmag in forward (tangent) mode (with options noISIZE i4 dr8 r8):
+   !   variations   of useful results: mag
+   !   with respect to varying inputs: v
+   SUBROUTINE GETMAG_D(v, vb, mag, magb)
+   USE CONSTANTS
+   IMPLICIT NONE
+   ! Subroutine Variables
+   REAL(kind=realtype), DIMENSION(3), INTENT(IN) :: v
+   REAL(kind=realtype), DIMENSION(3), INTENT(IN) :: vb
+   REAL(kind=realtype), INTENT(OUT) :: mag
+   REAL(kind=realtype), INTENT(OUT) :: magb
+   INTRINSIC SQRT
+   REAL(kind=realtype) :: arg1
+   REAL(kind=realtype) :: arg1b
+   REAL(kind=realtype) :: temp
+   arg1b = 2*v(1)*vb(1) + 2*v(2)*vb(2) + 2*v(3)*vb(3)
+   arg1 = v(1)**2 + v(2)**2 + v(3)**2 + 1e-30
+   temp = SQRT(arg1)
+   IF (arg1 .EQ. 0.0_8) THEN
+   magb = 0.0_8
+   ELSE
+   magb = arg1b/(2.0*temp)
+   END IF
+   mag = temp
+   END SUBROUTINE GETMAG_D
+

src/adjoint/outputReverse/Makefile

@@ -11,10 +11,8 @@ include ${COMMON_FILE}
 include ${RULES_FILE}
 vpath %.o $(OBJDIR)
 
-FF90_OBJECTS =	cross_product_3d_b.o\
-		getelementprops_b.o\
-		getmag_b.o\
-		getrotationmatrix3d_b.o
+FF90_OBJECTS =	vectorUtils_b.o\
+		getElementProps_b.o
 
 all:
 	$(PMAKE)  src