RUN_PART_B Subroutine

subroutine RUN_PART_B(N_FAIL)

Part B: DECOMPOSE+RECONSTRUCT FIELD_QUAD_BUMP at increasing (quadrature order scaled with each: LMAX_QUAD = 3*LMAX_TEST+30, NPHI = 4*LMAX_TEST+32 -- a higher floor than example 2d needs, since the equator-centered bump's sharper angular structure otherwise leaves a small quadrature-truncation leak in a_rad/a_pol at the smallest LMAX_TEST, shrinking with quadrature order rather than indicating an impure construction), tracking max reconstruction error against the true closed-form field on an independent evaluation grid, plus max a_rad/a_pol leakage across the sweep. Self-check: error at the largest must be much smaller than at the smallest (genuine convergence). Writes examples/vsh_decomposition_tor/convergence.dat.

Note

The evaluation grid still avoids landing exactly on the poles (THETA = (ITH+0.5)*pi/NEVAL, not ITH*pi/(NEVAL-1) as example 2d uses), a holdover from an earlier version of this example where it was load-bearing: PVSH_TOR_ALL/PVSH_POL_ALL guard against dividing by by returning exactly zero for every basis function at , and that earlier bump had a nonzero, perfectly well-defined true limit right at the pole (its cancelled analytically before being coded up), so comparing exactly at the pole produced a hard plateau at reconstruction error that never shrank with -- not an unphysical field, but a coordinate-singular comparison the library's own safety guard could never match. The current FIELD_QUAD_BUMP no longer has that problem -- its true value is exactly zero at both poles (see module docstring), matching what PVSH_TOR_ALL/PVSH_POL_ALL's guard already returns there -- but the off-pole grid is kept anyway as ordinary numerical hygiene near a coordinate singularity.

largest is not far smaller than at the smallest, or if a_rad/a_pol leakage exceeds TOL_A anywhere in the sweep.

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~2~~CallsGraph proc~run_part_b~2 RUN_PART_B proc~decompose~2 DECOMPOSE proc~run_part_b~2->proc~decompose~2 proc~field_quad_bump FIELD_QUAD_BUMP proc~run_part_b~2->proc~field_quad_bump proc~reconstruct~2 RECONSTRUCT proc~run_part_b~2->proc~reconstruct~2 dacos dacos proc~decompose~2->dacos interface~dot DOT proc~decompose~2->interface~dot interface~pvsh_pol_all PVSH_POL_ALL proc~decompose~2->interface~pvsh_pol_all interface~pvsh_rad_all PVSH_RAD_ALL proc~decompose~2->interface~pvsh_rad_all interface~pvsh_tor_all PVSH_TOR_ALL proc~decompose~2->interface~pvsh_tor_all interface~shglq SHGLQ proc~decompose~2->interface~shglq dcos dcos proc~field_quad_bump->dcos dexp dexp proc~field_quad_bump->dexp dsin dsin proc~field_quad_bump->dsin interface~pvsh_tor PVSH_TOR proc~field_quad_bump->interface~pvsh_tor proc~reconstruct~2->interface~pvsh_pol_all proc~reconstruct~2->interface~pvsh_rad_all proc~reconstruct~2->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_dp PVSH_TOR_DP interface~pvsh_tor->proc~pvsh_tor_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~pvsh_tor_dp->dcos proc~pvsh_tor_dp->dsin interface~assoc_legendre_and_deriv ASSOC_LEGENDRE_AND_DERIV proc~pvsh_tor_dp->interface~assoc_legendre_and_deriv interface~log_fact LOG_FACT proc~pvsh_tor_dp->interface~log_fact proc~shglq_dp->dcos proc~assoc_legendre_and_deriv_dp ASSOC_LEGENDRE_AND_DERIV_DP interface~assoc_legendre_and_deriv->proc~assoc_legendre_and_deriv_dp proc~log_fact_dp LOG_FACT_DP interface~log_fact->proc~log_fact_dp proc~vsh_core_dp VSH_CORE_DP interface~vsh_core->proc~vsh_core_dp proc~assoc_legendre_and_deriv_dp->interface~log_fact 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~2~~CalledByGraph proc~run_part_b~2 RUN_PART_B program~vsh_decomposition_tor VSH_DECOMPOSITION_TOR program~vsh_decomposition_tor->proc~run_part_b~2

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, MAX_PURITY

  MAX_PURITY = 0.0_dp

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

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

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

    MAX_PURITY = MAX(MAX_PURITY, MAXVAL(ABS(A_RAD)), MAXVAL(ABS(A_POL)))

    MAX_ERR = 0.0_dp
    DO ITH = 0, NEVAL-1
      THETA = (ITH+0.5d0)*pi/NEVAL
      DO IPH = 0, NEVAL-1
        PHI = IPH*2.d0*pi/NEVAL
        F_TRUE  = FIELD_QUAD_BUMP(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(74,*) 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(74)

  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

  IF (MAX_PURITY > TOL_A) THEN
    N_FAIL = N_FAIL + 1
    WRITE(*,'(A,ES10.3)') "FAIL  toroidal_purity_sweep  max_leakage=", MAX_PURITY
  ELSE
    WRITE(*,'(A,ES10.3)') "PASS  toroidal_purity_sweep  max_leakage=", MAX_PURITY
  END IF
  END SUBROUTINE RUN_PART_B