RUN_PART_B Subroutine

subroutine RUN_PART_B(N_FAIL)

Part B: DECOMPOSE+RECONSTRUCT FIELD_GAUSSIAN at increasing (quadrature order scaled with each: LMAX_QUAD = 3*LMAX_TEST+10, NPHI = 4*LMAX_TEST+16), tracking max reconstruction error against the true closed-form field on an independent evaluation grid. Self-check: error at the largest must be much smaller than at the smallest (genuine convergence, not a flat/noisy curve -- which would indicate a quadrature-order-vs- scaling bug). Writes examples/vsh_decomposition/convergence.dat.

largest is not far smaller than at the smallest.

Arguments

Type IntentOptional Attributes Name
integer(kind=i4), intent(inout) :: N_FAIL

Running failure counter, incremented if the error at the

Calls

proc~~run_part_b~~CallsGraph proc~run_part_b RUN_PART_B proc~decompose DECOMPOSE proc~run_part_b->proc~decompose proc~field_gaussian FIELD_GAUSSIAN proc~run_part_b->proc~field_gaussian proc~reconstruct RECONSTRUCT proc~run_part_b->proc~reconstruct dacos dacos proc~decompose->dacos interface~dot DOT proc~decompose->interface~dot interface~pvsh_pol_all PVSH_POL_ALL proc~decompose->interface~pvsh_pol_all interface~pvsh_rad_all PVSH_RAD_ALL proc~decompose->interface~pvsh_rad_all interface~pvsh_tor_all PVSH_TOR_ALL proc~decompose->interface~pvsh_tor_all interface~shglq SHGLQ proc~decompose->interface~shglq dcos dcos proc~field_gaussian->dcos dexp dexp proc~field_gaussian->dexp dsin dsin proc~field_gaussian->dsin proc~reconstruct->interface~pvsh_pol_all proc~reconstruct->interface~pvsh_rad_all proc~reconstruct->interface~pvsh_tor_all proc~dot_dp DOT_DP interface~dot->proc~dot_dp proc~pvsh_pol_all_dp PVSH_POL_ALL_DP interface~pvsh_pol_all->proc~pvsh_pol_all_dp proc~pvsh_rad_all_dp PVSH_RAD_ALL_DP interface~pvsh_rad_all->proc~pvsh_rad_all_dp proc~pvsh_tor_all_dp PVSH_TOR_ALL_DP interface~pvsh_tor_all->proc~pvsh_tor_all_dp proc~shglq_dp SHGLQ_DP interface~shglq->proc~shglq_dp interface~vsh_core VSH_CORE proc~pvsh_pol_all_dp->interface~vsh_core proc~ylm_index YLM_INDEX proc~pvsh_pol_all_dp->proc~ylm_index proc~pvsh_rad_all_dp->interface~vsh_core proc~pvsh_tor_all_dp->interface~vsh_core proc~pvsh_tor_all_dp->proc~ylm_index proc~shglq_dp->dcos proc~vsh_core_dp VSH_CORE_DP interface~vsh_core->proc~vsh_core_dp proc~vsh_core_dp->dcos proc~vsh_core_dp->dsin proc~vsh_core_dp->proc~ylm_index interface~assoc_legendre_norm_all ASSOC_LEGENDRE_NORM_ALL proc~vsh_core_dp->interface~assoc_legendre_norm_all interface~ddx_assoc_legendre_norm_all DDX_ASSOC_LEGENDRE_NORM_ALL proc~vsh_core_dp->interface~ddx_assoc_legendre_norm_all proc~plm_index PLM_INDEX proc~vsh_core_dp->proc~plm_index proc~assoc_legendre_norm_all_dp ASSOC_LEGENDRE_NORM_ALL_DP interface~assoc_legendre_norm_all->proc~assoc_legendre_norm_all_dp proc~ddx_assoc_legendre_norm_all_dp DDX_ASSOC_LEGENDRE_NORM_ALL_DP interface~ddx_assoc_legendre_norm_all->proc~ddx_assoc_legendre_norm_all_dp proc~assoc_legendre_norm_all_dp->proc~plm_index proc~ddx_assoc_legendre_norm_all_dp->proc~plm_index

Called by

proc~~run_part_b~~CalledByGraph proc~run_part_b RUN_PART_B program~vsh_decomposition VSH_DECOMPOSITION program~vsh_decomposition->proc~run_part_b

Source Code

  SUBROUTINE RUN_PART_B(N_FAIL)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(INOUT) :: N_FAIL
  INTEGER(KIND=i4), PARAMETER :: NEVAL = 40
  REAL(KIND=dp),    PARAMETER :: CONVERGENCE_RATIO_TOL = 0.01_dp
  INTEGER(KIND=i4) :: IB, LMAX_TEST, LMAX_QUAD, NPHI, ITH, IPH
  COMPLEX(KIND=dp), ALLOCATABLE :: A_RAD(:), A_POL(:), A_TOR(:)
  COMPLEX(KIND=dp), DIMENSION(3) :: F_TRUE, F_RECON
  REAL(KIND=dp) :: THETA, PHI, MAX_ERR, FIRST_ERR, LAST_ERR

  OPEN(UNIT=72, FILE="./examples/vsh_decomposition/convergence.dat")
  WRITE(72,'(A)') '# Lmax  max_reconstruction_error'

  DO IB = 1, N_LMAX_B
    LMAX_TEST = LMAX_B_LIST(IB)
    LMAX_QUAD = 3*LMAX_TEST + 10
    NPHI = 4*LMAX_TEST + 16

    ALLOCATE(A_RAD((LMAX_TEST+1)**2), A_POL((LMAX_TEST+1)**2), &
             A_TOR((LMAX_TEST+1)**2))
    CALL DECOMPOSE(FIELD_GAUSSIAN, LMAX_TEST, LMAX_QUAD, NPHI, &
                   A_RAD, A_POL, A_TOR)

    MAX_ERR = 0.0_dp
    DO ITH = 0, NEVAL-1
      THETA = ITH*pi/(NEVAL-1)
      DO IPH = 0, NEVAL-1
        PHI = IPH*2.d0*pi/NEVAL
        F_TRUE  = FIELD_GAUSSIAN(THETA, PHI)
        F_RECON = RECONSTRUCT(LMAX_TEST, A_RAD, A_POL, A_TOR, THETA, PHI)
        MAX_ERR = MAX(MAX_ERR, ABS(F_RECON(1)-F_TRUE(1)), &
                      ABS(F_RECON(2)-F_TRUE(2)), ABS(F_RECON(3)-F_TRUE(3)))
      ENDDO
    ENDDO

    IF (IB == 1) FIRST_ERR = MAX_ERR
    IF (IB == N_LMAX_B) LAST_ERR = MAX_ERR

    WRITE(72,*) LMAX_TEST, MAX_ERR
    WRITE(*,'(A,I0,A,ES10.3)') "  Lmax=", LMAX_TEST, &
      "  max_reconstruction_error=", MAX_ERR

    DEALLOCATE(A_RAD, A_POL, A_TOR)
  ENDDO
  CLOSE(72)

  IF (LAST_ERR > CONVERGENCE_RATIO_TOL*FIRST_ERR) THEN
    N_FAIL = N_FAIL + 1
    WRITE(*,'(A,ES10.3,A,ES10.3)') "FAIL  spectral_convergence  first=", &
      FIRST_ERR, "  last=", LAST_ERR
  ELSE
    WRITE(*,'(A,ES10.3,A,ES10.3)') "PASS  spectral_convergence  first=", &
      FIRST_ERR, "  last=", LAST_ERR
  END IF
  END SUBROUTINE RUN_PART_B