From afccdba6d5a61d813337959761296b7df1dc106c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 25 Mar 2016 16:04:14 -0400 Subject: [PATCH] Improve multipole readability --- src/cross_section.F90 | 145 +++++++++++++++++++++--------------------- src/initialize.F90 | 2 +- src/math.F90 | 29 +++++---- 3 files changed, 88 insertions(+), 88 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1717b733d7..930d1a0529 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -6,7 +6,7 @@ module cross_section use global use list_header, only: ListElemInt use material_header, only: Material - use math, only: w, broaden_n_polynomials + use math, only: faddeeva, broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, & MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,& FIT_F, MultipoleArray, max_poly, max_L, max_poles @@ -588,123 +588,119 @@ contains real(8), intent(out) :: sigT ! Total cross section real(8), intent(out) :: sigA ! Absorption cross section real(8), intent(out) :: sigF ! Fission cross section - complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic - ! form + complex(8) :: psi_chi ! The value of the psi-chi function for the + ! asymptotic form complex(8) :: c_temp ! complex temporary variable complex(8) :: w_val ! The faddeeva function evaluated at Z complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole) real(8) :: sqrtE ! sqrt(E), eV real(8) :: invE ! 1/E, eV - real(8) :: dopp ! sqrt(atomic weight ratio / kT) - real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E + real(8) :: dopp ! sqrt(atomic weight ratio / kT) = 1 / (2 sqrt(xi)) real(8) :: temp ! real temporary value real(8) :: E ! energy, eV real(8) :: sqrtkT ! sqrt(kT (in eV)) - integer :: iP ! index of pole - integer :: iC ! index of curvefit - integer :: iW ! index of window + integer :: i_pole ! index of pole + integer :: i_poly ! index of curvefit + integer :: i_window ! index of window integer :: startw ! window start pointer (for poles) - integer :: startw_1 ! window start pointer - 1 - integer :: startw_endw ! window start pointer - window end pointer integer :: endw ! window end pointer - ! Convert to eV + ! ========================================================================== + ! Bookkeeping + + ! Convert to eV. E = Emev * 1.0e6_8 sqrtkT = sqrtkT_ * 1.0e3_8 + ! Define some frequently used variables. sqrtE = sqrt(E) - invE = ONE/E + invE = ONE / E + dopp = multipole % sqrtAWR / sqrtkT - if(.not. mp_already_alloc) then + if (.not. mp_already_alloc) then call multipole_eval_allocate() end if - ! Locate us - iW = floor((sqrtE - sqrt(multipole % start_E))/multipole % spacing + ONE) + ! Locate us. + i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing & + + ONE) + startw = multipole % w_start(i_window) + endw = multipole % w_end(i_window) - startw = multipole % w_start(iW) - startw_1 = startw - 1 ! This is an index shift parameter. - endw = multipole % w_end(iW) - startw_endw = endw - startw + 1 - - ! Fill in factors + ! Fill in factors. if (startw <= endw) then - call fill_factors(multipole, sqrtE, sigT_factor, twophi, multipole % num_l) + call fill_factors(multipole, sqrtE, sigT_factor, twophi, & + multipole % num_l) end if - ! Generate some doppler broadening parameters - - ! dopp_ecoef is inverse of dopp, divided by E, multiplied by sqrt(pi). - dopp = multipole % sqrtAWR / sqrtKT - dopp_ecoef = dopp * invE * SQRT_PI - + ! Initialize the ouptut cross sections. sigT = ZERO sigA = ZERO sigF = ZERO - ! Evaluate linefit first - if(sqrtkT /= 0 .and. multipole % broaden_poly(iW) == 1) then ! Broaden the curvefit. - call broaden_n_polynomials(E, dopp, multipole % fit_order + 1, broadened_polynomials) + ! ========================================================================== + ! Add the contribution from the curvefit polynomial. - do iC = 1, multipole % fit_order+1 - sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*broadened_polynomials(iC) - sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*broadened_polynomials(iC) - if (multipole % fissionable) then - sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*broadened_polynomials(iC) - end if + if (sqrtkT /= ZERO .and. multipole % broaden_poly(i_window) == 1) then + ! Broaden the curvefit. + call broaden_wmp_polynomials(E, dopp, multipole % fit_order + 1, & + broadened_polynomials) + do i_poly = 1, multipole % fit_order+1 + sigT = sigT + multipole % curvefit(FIT_T, i_poly, i_window) & + * broadened_polynomials(i_poly) + sigA = sigA + multipole % curvefit(FIT_A, i_poly, i_window) & + * broadened_polynomials(i_poly) + sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) & + * broadened_polynomials(i_poly) end do else ! Evaluate as if it were a polynomial temp = invE - do iC = 1, multipole % fit_order+1 - - sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*temp - sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*temp - if (multipole % fissionable) then - sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*temp - end if - + do i_poly = 1, multipole % fit_order+1 + sigT = sigT + multipole % curvefit(FIT_T, i_poly, i_window) * temp + sigA = sigA + multipole % curvefit(FIT_A, i_poly, i_window) * temp + sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) * temp temp = temp * sqrtE end do end if - ! Then get the poles we want and broaden them. + ! ========================================================================== + ! Add the contribution from the poles in this window. if (sqrtkT == ZERO) then ! If at 0K, use asymptotic form. - do iP = startw, endw - psi_ki = -ONEI/(multipole % data(MP_EA, iP) - sqrtE) - c_temp = psi_ki/E + do i_pole = startw, endw + psi_chi = -ONEI / (multipole % data(MP_EA, i_pole) - sqrtE) + c_temp = psi_chi / E if (multipole % formalism == FORM_MLBW) then - sigT = sigT + real(multipole % data(MLBW_RT, iP) * c_temp * & - sigT_factor(multipole % l_value(iP))) & - + real(multipole % data(MLBW_RX, iP) * c_temp) - sigA = sigA + real(multipole % data(MLBW_RA, iP) * c_temp) - sigF = sigF + real(multipole % data(MLBW_RF, iP) * c_temp) + sigT = sigT + real(multipole % data(MLBW_RT, i_pole) * c_temp * & + sigT_factor(multipole % l_value(i_pole))) & + + real(multipole % data(MLBW_RX, i_pole) * c_temp) + sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * c_temp) + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * c_temp) else if (multipole % formalism == FORM_RM) then - sigT = sigT + real(multipole % data(RM_RT, iP) * c_temp* & - sigT_factor(multipole % l_value(iP))) - sigA = sigA + real(multipole % data(RM_RA, iP) * c_temp) - sigF = sigF + real(multipole % data(RM_RF, iP) * c_temp) + sigT = sigT + real(multipole % data(RM_RT, i_pole) * c_temp * & + sigT_factor(multipole % l_value(i_pole))) + sigA = sigA + real(multipole % data(RM_RA, i_pole) * c_temp) + sigF = sigF + real(multipole % data(RM_RF, i_pole) * c_temp) end if end do else ! At temperature, use Faddeeva function-based form. - if(endw >= startw) then - do iP = startw, endw - Z = (sqrtE - multipole % data(MP_EA, iP)) * dopp - w_val = w(Z) * dopp_ecoef - + if (endw >= startw) then + do i_pole = startw, endw + Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp + w_val = faddeeva(Z) * dopp * invE * SQRT_PI if (multipole % formalism == FORM_MLBW) then - sigT = sigT + real((multipole % data(MLBW_RT, iP) * & - sigT_factor(multipole%l_value(iP)) + & - multipole % data(MLBW_RX, iP)) * w_val) - sigA = sigA + real(multipole % data(MLBW_RA, iP) * w_val) - sigF = sigF + real(multipole % data(MLBW_RF, iP) * w_val) + sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * & + sigT_factor(multipole % l_value(i_pole)) + & + multipole % data(MLBW_RX, i_pole)) * w_val) + sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) else if (multipole % formalism == FORM_RM) then - sigT = sigT + real(multipole % data(RM_RT, iP) * w_val * & - sigT_factor(multipole % l_value(iP))) - sigA = sigA + real(multipole % data(RM_RA, iP) * w_val) - sigF = sigF + real(multipole % data(RM_RF, iP) * w_val) + sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * & + sigT_factor(multipole % l_value(i_pole))) + sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) end if end do end if @@ -735,13 +731,14 @@ contains arg = 3.0_8 * twophi(iL) / (3.0_8 - twophi(iL)**2) twophi(iL) = twophi(iL) - atan(arg) else if (iL == 4) then - arg = twophi(iL) * (15.0_8 - twophi(iL)**2) / (15.0_8 - 6.0_8 * twophi(iL)**2) + arg = twophi(iL) * (15.0_8 - twophi(iL)**2) & + / (15.0_8 - 6.0_8 * twophi(iL)**2) twophi(iL) = twophi(iL) - atan(arg) end if end do twophi = 2.0_8 * twophi - sigT_factor = cmplx(cos(twophi),-sin(twophi), KIND=8) + sigT_factor = cmplx(cos(twophi), -sin(twophi), KIND=8) end subroutine !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index 8502a361b3..f34215fbab 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1178,7 +1178,7 @@ contains do k = 2, mat % n_nuclides ! Warn the user if the nuclides don't have identical temperatues. if (nuclides(mat % nuclide(k)) % kT /= min_temp & - .and. .not. warning_given) then + .and. .not. warning_given .and. multipole_active) then call warning("OpenMC cannot & &identify the temperature of at least one cell. For the & &purposes of multipole cross section evaluations, all cells & diff --git a/src/math.F90 b/src/math.F90 index 57e238afa0..fce57e6fc2 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -719,13 +719,14 @@ contains end function watt_spectrum !=============================================================================== -! W acts as a front end to the MIT Faddeeva function, Faddeeva_w. +! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation !=============================================================================== - function w(z) result(wv) + function faddeeva(z) result(wv) complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at complex(8) :: wv ! The resulting w(z) value - real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT Faddeeva + real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT + ! Faddeeva ! Technically, the value we want is given by the equation: ! w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}] @@ -747,23 +748,24 @@ contains wv = -conjg(faddeeva_w(conjg(z), relerr)) end if - end function w + end function faddeeva !=============================================================================== -! BROADEN_N_POLYNOMIALS doppler broadens polynomials of the form +! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. The +! curvefit is a polynomial of the form ! a/En + b/sqrt(En) + c + d sqrt(En) ... -! exactly and quickly. !=============================================================================== - subroutine broaden_n_polynomials(En, dopp, n, factors) + subroutine broaden_wmp_polynomials(En, dopp, n, factors) real(8), intent(in) :: En ! Energy to evaluate at - real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), kT given in eV. + real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), + ! kT given in eV. integer, intent(in) :: n ! number of components to polynomial real(8), intent(out):: factors(n) ! output leading coefficient integer :: i - real(8) :: sqrtE ! Sqrt(energy) + real(8) :: sqrtE ! sqrt(energy) real(8) :: beta ! sqrt(atomic weight ratio * E / kT) real(8) :: half_inv_dopp2 ! 0.5 / dopp**2 real(8) :: quarter_inv_dopp4 ! 0.25 / dopp**4 @@ -790,20 +792,21 @@ contains factors(1) = erfbeta / En factors(2) = ONE / sqrtE - factors(3) = factors(1) * (half_inv_dopp2 + En) + exp_m_beta2 / (beta * SQRT_PI) + factors(3) = factors(1) * (half_inv_dopp2 + En) & + + exp_m_beta2 / (beta * SQRT_PI) ! Perform recursive broadening of high order components do i = 1, n-3 if (i /= 1) then - factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 + & - factors(i+1) * (En + (ONE + TWO * i) * half_inv_dopp2) + factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 & + + factors(i+1) * (En + (ONE + TWO * i) * half_inv_dopp2) else ! Although it's mathematically identical, factors(0) will contain ! nothing, and we don't want to have to worry about memory. factors(i+3) = factors(i+1)*(En + (ONE + TWO * i) * half_inv_dopp2) end if end do - end subroutine broaden_n_polynomials + end subroutine broaden_wmp_polynomials !=============================================================================== ! find_angle finds the closest angle on the data grid and returns that index