GWI_GWJ_SWEEP – FORTVSH

Uses

- KINDS
- VSH
- GLOBALS
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.

Worked example (2e): Geppert-Wiebicke (GWI/GWJ) coefficients swept across a range of mode-number combinations, each checked against an independent formula (not just the closed-form GWI/GWJ code itself evaluated a second way).

GWI check: GWI's own formula is algebraically the standard closed-form Gaunt coefficient -- the textbook expression for the triple-scalar-harmonic integral $I_{j_1m_1j_2m_2}^{\ell m} = \int Y_{j_1}^{m_1}Y_{j_2}^{m_2} (Y_\ell^m)^{*}\,d\Omega.$ This was verified numerically against direct 2D quadrature (SHGLQ in $\theta$ times uniform trapezoidal in $\phi$ , the same technique introduced in vsh_decomposition.f90) before writing this sweep, confirming the identification to $\sim10^{-16}$ across several mode combinations, including a selection-rule-forbidden case that correctly evaluates to exactly zero on both sides.

GWJ check: GWJ's defining relationship is the coefficient in the expansion of a poloidal-toroidal VSH dot product in scalar harmonics (the same pattern src/tests.f90's TEST2_GW/TEST2_VSH already validate for one hardcoded case): $\mathrm{[[PVSH_POL]]}(j_2,m_2)\cdot\mathrm{[[PVSH_TOR]]}(j_1,m_1) = \frac{-i}{\sqrt{j_1(j_1{+}1)j_2(j_2{+}1)}}\sum_\ell J_{j_1m_1j_2m_2}^{\ell,m_1+m_2}\,Y_\ell^{m_1+m_2},$ summed over $\ell=|j_1-j_2|\ldots j_1{+}j_2$ (terms outside the valid triangle/parity range correctly evaluate to zero after the GWJ guard fix). This is a pointwise identity (no integration needed) and was independently verified against several new mode combinations (beyond the one hardcoded in TEST2_GW) before writing this sweep.

Calls

Help

Variables

Type	Attributes		Name		Initial
integer(kind=i4),	parameter	::	JMAX	=	6
integer(kind=i4),	parameter	::	LMAX_QUAD	=	20
integer(kind=i4),	parameter	::	NPHI	=	16
integer(kind=i4)		::	N_FAIL
integer(kind=i4),	parameter	::	N_POINTS	=	2
real(kind=dp),	parameter	::	PHI_PTS(N_POINTS)	=	(/1.3d0, 0.4d0/)
real(kind=dp),	parameter	::	THETA_PTS(N_POINTS)	=	(/0.7d0, 1.9d0/)
real(kind=dp),	parameter	::	TOL	=	1.0E-8_dp

Subroutines

subroutine RUN_GWI_SWEEP(N_FAIL)

GWI vs. the brute-force Gaunt-coefficient integral (via SHGLQ theta-quadrature times a uniform phi sum), swept over $j_1,j_2=0\ldots$ JMAX, all valid $m_1,m_2,\ell,m$ . Writes examples/gwi_gwj_sweep/gwi_check.dat.

Arguments

Type	Intent	Optional	Attributes		Name
integer(kind=i4),	intent(inout)			::	N_FAIL	Running failure counter, incremented if the max error

subroutine RUN_GWJ_SWEEP(N_FAIL)

GWJ vs. the $\mathrm{[[PVSH_POL]]}\cdot\mathrm{[[PVSH_TOR]]}$ pointwise identity, swept over $j_1,j_2=0\ldots$ JMAX, all valid $m_1,m_2$ , at N_POINTS fixed $(\theta,\phi)$ . Writes examples/gwi_gwj_sweep/gwj_check.dat.

Arguments

Type	Intent	Optional	Attributes		Name
integer(kind=i4),	intent(inout)			::	N_FAIL	Running failure counter, incremented if the max error

Source Code

PROGRAM GWI_GWJ_SWEEP
USE KINDS,   ONLY: dp, i4
USE GLOBALS, ONLY: pi, j
USE VSH,     ONLY: GWI, GWJ, SSH, PVSH_POL, PVSH_TOR, DOT, SHGLQ
IMPLICIT NONE

INTEGER(KIND=i4), PARAMETER :: JMAX = 6
INTEGER(KIND=i4), PARAMETER :: LMAX_QUAD = 20
INTEGER(KIND=i4), PARAMETER :: NPHI = 16
REAL(KIND=dp),    PARAMETER :: TOL = 1.0E-8_dp

INTEGER(KIND=i4), PARAMETER :: N_POINTS = 2
REAL(KIND=dp), PARAMETER :: THETA_PTS(N_POINTS) = (/0.7d0, 1.9d0/)
REAL(KIND=dp), PARAMETER :: PHI_PTS(N_POINTS)   = (/1.3d0, 0.4d0/)

INTEGER(KIND=i4) :: N_FAIL

N_FAIL = 0
CALL RUN_GWI_SWEEP(N_FAIL)
CALL RUN_GWJ_SWEEP(N_FAIL)

IF (N_FAIL > 0) THEN
  WRITE(*,'(A,I0,A)') "RESULT: FAILED - ", N_FAIL, " check(s) failed"
  STOP 1
ELSE
  WRITE(*,'(A)') "RESULT: PASSED - gwi_gwj_sweep"
END IF

CONTAINS

!> [[GWI]] vs. the brute-force Gaunt-coefficient integral (via [[SHGLQ]]
!> theta-quadrature times a uniform phi sum), swept over
!> \( j_1,j_2=0\ldots\)`JMAX`, all valid \( m_1,m_2,\ell,m \). Writes
!> `examples/gwi_gwj_sweep/gwi_check.dat`.
!>
!> @param N_FAIL Running failure counter, incremented if the max error
!>   across the sweep exceeds `TOL`.
  SUBROUTINE RUN_GWI_SWEEP(N_FAIL)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(INOUT) :: N_FAIL
  REAL(KIND=dp) :: ZERO(LMAX_QUAD+1), WTHETA(LMAX_QUAD+1)
  REAL(KIND=dp) :: DPHI, THETA, PHI, MAX_ERR, ABSDIFF
  COMPLEX(KIND=dp) :: CLOSED, BRUTE
  INTEGER(KIND=i4) :: J1, J2, M1, M2, L, M, IQ, IP

  CALL SHGLQ(ZERO, WTHETA, LMAX_QUAD)
  DPHI = 2.d0*pi/NPHI
  MAX_ERR = 0.0_dp

  OPEN(UNIT=81, FILE="./examples/gwi_gwj_sweep/gwi_check.dat")
  WRITE(81,'(A)') '# J1  M1  J2  M2  L  M  '// &
                  'Re(closed)  Im(closed)  Re(brute)  Im(brute)  absdiff'

  DO J1 = 0, JMAX
    DO J2 = 0, JMAX
      DO M1 = -J1, J1
        DO M2 = -J2, J2
          DO L = ABS(J1-J2), J1+J2
            M = M1+M2
            IF (ABS(M) > L) CYCLE

            CLOSED = GWI(J1,M1,J2,M2,L,M)

            BRUTE = DCMPLX(0.d0, 0.d0)
            DO IQ = 1, LMAX_QUAD+1
              THETA = DACOS(ZERO(IQ))
              DO IP = 0, NPHI-1
                PHI = IP*DPHI
                BRUTE = BRUTE + WTHETA(IQ)*DPHI * SSH(J1,M1,THETA,PHI) * &
                        SSH(J2,M2,THETA,PHI) * CONJG(SSH(L,M,THETA,PHI))
              ENDDO
            ENDDO

            ABSDIFF = ABS(CLOSED-BRUTE)
            MAX_ERR = MAX(MAX_ERR, ABSDIFF)
            WRITE(81,*) J1, M1, J2, M2, L, M, DBLE(CLOSED), AIMAG(CLOSED), &
                        DBLE(BRUTE), AIMAG(BRUTE), ABSDIFF
          ENDDO
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  CLOSE(81)

  IF (MAX_ERR > TOL) THEN
    N_FAIL = N_FAIL + 1
    WRITE(*,'(A,ES10.3)') "FAIL  gwi_gaunt_sweep  max_err=", MAX_ERR
  ELSE
    WRITE(*,'(A,ES10.3)') "PASS  gwi_gaunt_sweep  max_err=", MAX_ERR
  END IF
  END SUBROUTINE RUN_GWI_SWEEP

!> [[GWJ]] vs. the \( \mathrm{[[PVSH_POL]]}\cdot\mathrm{[[PVSH_TOR]]} \)
!> pointwise identity, swept over \( j_1,j_2=0\ldots\)`JMAX`, all valid
!> \( m_1,m_2 \), at `N_POINTS` fixed \( (\theta,\phi) \). Writes
!> `examples/gwi_gwj_sweep/gwj_check.dat`.
!>
!> @param N_FAIL Running failure counter, incremented if the max error
!>   across the sweep exceeds `TOL`.
  SUBROUTINE RUN_GWJ_SWEEP(N_FAIL)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(INOUT) :: N_FAIL
  COMPLEX(KIND=dp) :: LHS, RHS, PREFACTOR
  REAL(KIND=dp) :: MAX_ERR, ABSDIFF, THETA, PHI
  INTEGER(KIND=i4) :: J1, J2, M1, M2, L, M, IP

  MAX_ERR = 0.0_dp

  OPEN(UNIT=82, FILE="./examples/gwi_gwj_sweep/gwj_check.dat")
  WRITE(82,'(A)') '# J1  M1  J2  M2  theta  phi  '// &
                  'Re(lhs)  Im(lhs)  Re(rhs)  Im(rhs)  absdiff'

  DO J1 = 1, JMAX
    DO J2 = 1, JMAX
      PREFACTOR = -j/DSQRT(DBLE(J1)*(J1+1.d0)*DBLE(J2)*(J2+1.d0))
      DO M1 = -J1, J1
        DO M2 = -J2, J2
          M = M1+M2
          DO IP = 1, N_POINTS
            THETA = THETA_PTS(IP)
            PHI   = PHI_PTS(IP)

            LHS = DOT(PVSH_POL(J2,M2,THETA,PHI), PVSH_TOR(J1,M1,THETA,PHI))

            RHS = DCMPLX(0.d0, 0.d0)
            DO L = ABS(J1-J2), J1+J2
              IF (ABS(M) <= L) THEN
                RHS = RHS + GWJ(J1,M1,J2,M2,L,M)*SSH(L,M,THETA,PHI)
              END IF
            ENDDO
            RHS = PREFACTOR * RHS

            ABSDIFF = ABS(LHS-RHS)
            MAX_ERR = MAX(MAX_ERR, ABSDIFF)
            WRITE(82,*) J1, M1, J2, M2, THETA, PHI, DBLE(LHS), AIMAG(LHS), &
                        DBLE(RHS), AIMAG(RHS), ABSDIFF
          ENDDO
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  CLOSE(82)

  IF (MAX_ERR > TOL) THEN
    N_FAIL = N_FAIL + 1
    WRITE(*,'(A,ES10.3)') "FAIL  gwj_identity_sweep  max_err=", MAX_ERR
  ELSE
    WRITE(*,'(A,ES10.3)') "PASS  gwj_identity_sweep  max_err=", MAX_ERR
  END IF
  END SUBROUTINE RUN_GWJ_SWEEP

END PROGRAM GWI_GWJ_SWEEP

GWI_GWJ_SWEEP Program

Contents

Variables

Subroutines

Source Code

Uses

Calls

Variables

Subroutines

subroutine RUN_GWI_SWEEP(N_FAIL)

Arguments

subroutine RUN_GWJ_SWEEP(N_FAIL)

Arguments

Source Code