From 39718fe2cdc984c7b3927341382aa10b8370d94b Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Thu, 16 May 2024 16:11:29 +0100 Subject: [PATCH] Convert neoclassics_module.f90 to Python * Remove unused neoclassic routine * Add tests for neoclassics module * Convert init_neoclassics to Python * Convert init_profile_values_from_PROCESS to Python * Convert neoclassics roots and weights functions to Python * Convert neoclassics_calc_KT to Python * Convert neoclassics_calc_nu to Python * Convert neoclassics_calc_nu_star to Python * Convert neoclassics_calc_nu_star_fromT to Python * Convert neoclassics_calc_vd to Python * Convert neoclassics_calc_D11_plateau to Python * Convert neoclassics_calc_d11_mono to Python * Convert neoclassics_calc_D11(1/2/3) to Python * Convert neoclassics_calc_Gamma_flux to Python * Convert neoclassics_calc_q_flux to Python --- process/main.py | 4 +- process/stellarator.py | 571 +++++++++++++- source/fortran/constants.f90 | 3 + source/fortran/neoclassics_module.f90 | 538 ------------- tests/unit/test_neoclassics.py | 1053 +++++++++++++++++++++++++ tests/unit/test_stellarator.py | 3 +- 6 files changed, 1625 insertions(+), 547 deletions(-) create mode 100644 tests/unit/test_neoclassics.py diff --git a/process/main.py b/process/main.py index 284bea95..ddf5acac 100644 --- a/process/main.py +++ b/process/main.py @@ -50,7 +50,7 @@ from process.plasma_geometry import PlasmaGeom from process.pulse import Pulse from process.scan import Scan -from process.stellarator import Stellarator +from process.stellarator import Stellarator, Neoclassics from process.structure import Structure from process.build import Build from process.utilities.f2py_string_patch import string_to_f2py_compatible @@ -629,6 +629,7 @@ def __init__(self): self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive ) + self.neoclassics = Neoclassics() self.stellarator = Stellarator( availability=self.availability, buildings=self.buildings, @@ -639,6 +640,7 @@ def __init__(self): hcpb=self.ccfe_hcpb, current_drive=self.current_drive, physics=self.physics, + neoclassics=self.neoclassics, ) self.dcll = DCLL(blanket_library=self.blanket_library) diff --git a/process/stellarator.py b/process/stellarator.py index c11f7863..fd79eaed 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -39,6 +39,11 @@ s_handler.setLevel(logging.ERROR) logger.addHandler(s_handler) +# NOTE: a different value of echarge was used in the original implementation +# making the post-Python results slightly different. As a result, there is a +# relative tolerance on the neoclassics tests of 1e-3 +KEV = 1e3 * constants.echarge # Kiloelectron-volt (keV) + class Stellarator: """Module containing stellarator routines @@ -61,6 +66,7 @@ def __init__( hcpb, current_drive, physics, + neoclassics, ) -> None: """Initialises the Stellarator model's variables @@ -80,6 +86,8 @@ def __init__( :type current_drive: process.current_drive.CurrentDrive :param physics: a pointer to the Physics model, allowing use of Physics's variables/methods :type physics: process.physics.Physics + :param neoclassics: a pointer to the Neoclassics model, allowing use of neoclassics's variables/methods + :type neoclassics: process.stellarator.Neoclassics """ self.outfile: int = constants.nout @@ -94,6 +102,7 @@ def __init__( self.hcpb = hcpb self.current_drive = current_drive self.physics = physics + self.neoclassics = neoclassics def run(self, output: bool): """Routine to call the physics and engineering modules @@ -3238,7 +3247,7 @@ def stopt_output( ) po.ovarre( self.outfile, - "Maximum reachable ECRH temperature (pseudo) (keV)", + "Maximum reachable ECRH temperature (pseudo) (KEV)", "(te0_ecrh_achievable)", te0_ecrh_achievable, ) @@ -3269,7 +3278,7 @@ def power_at_ignition_point(self, gyro_frequency_max, te0_available): current ECRH achievable peak temperature (which is inaccurate as the cordey pass should be calculated) author: J Lion, IPP Greifswald gyro_frequency_max : input real : Maximal available Gyrotron frequency (1/s) NOT (rad/s) - te0_available : input real : Reachable peak electron temperature, reached by ECRH (keV) + te0_available : input real : Reachable peak electron temperature, reached by ECRH (KEV) powerht_out : output real: Heating Power at ignition point (MW) pscalingmw_out : output real: Heating Power loss at ignition point (MW) This routine calculates the density limit due to an ECRH heating scheme on axis @@ -4472,7 +4481,7 @@ def stphys_output( ) def calc_neoclassics(self): - neoclassics_module.init_neoclassics( + self.neoclassics.init_neoclassics( 0.6, stellarator_configuration.stella_config_epseff, stellarator_variables.iotabar, @@ -4632,7 +4641,7 @@ def st_calc_eff_chi(self): - physics_variables.pcoreradpv ) * volscaling - # in fortran there was a 0d0*alphan term which I have removed for obvious reasons + # in fortran there was a 0*alphan term which I have removed for obvious reasons # the following comment seems to describe this? # "include alphan if chi should be incorporate density gradients too" # but the history can be consulted if required (23/11/22 TN) @@ -4773,9 +4782,9 @@ def stheat(self, output: bool): if abs(current_drive_variables.pnbeam) > 1e-8: po.ovarre( self.outfile, - "Neutral beam energy (keV)", - "(beam_energy)", - current_drive_variables.beam_energy, + "Neutral beam energy (KEV)", + "(enbeam)", + current_drive_variables.enbeam, ) po.ovarre( self.outfile, @@ -4828,3 +4837,551 @@ def stheat(self, output: bool): "(taubeam)", current_drive_variables.taubeam, ) + + +class Neoclassics: + @property + def no_roots(self): + return neoclassics_module.roots.shape[0] + + def init_neoclassics(self, r_effin, eps_effin, iotain): + """Constructor of the neoclassics object from the effective radius, + epsilon effective and iota only. + """ + ( + neoclassics_module.densities, + neoclassics_module.temperatures, + neoclassics_module.dr_densities, + neoclassics_module.dr_temperatures, + ) = self.init_profile_values_from_PROCESS(r_effin) + neoclassics_module.roots = np.array( + [ + 4.740718054080526184e-2, + 2.499239167531593919e-1, + 6.148334543927683749e-1, + 1.143195825666101451, + 1.836454554622572344, + 2.696521874557216147, + 3.725814507779509288, + 4.927293765849881879, + 6.304515590965073635, + 7.861693293370260349, + 9.603775985479263255, + 1.153654659795613924e1, + 1.366674469306423489e1, + 1.600222118898106771e1, + 1.855213484014315029e1, + 2.132720432178312819e1, + 2.434003576453269346e1, + 2.760555479678096091e1, + 3.114158670111123683e1, + 3.496965200824907072e1, + 3.911608494906788991e1, + 4.361365290848483056e1, + 4.850398616380419980e1, + 5.384138540650750571e1, + 5.969912185923549686e1, + 6.618061779443848991e1, + 7.344123859555988076e1, + 8.173681050672767867e1, + 9.155646652253683726e1, + 1.041575244310588886e2, + ] + ) + neoclassics_module.weights = np.array( + [ + 1.160440860204388913e-1, + 2.208511247506771413e-1, + 2.413998275878537214e-1, + 1.946367684464170855e-1, + 1.237284159668764899e-1, + 6.367878036898660943e-2, + 2.686047527337972682e-2, + 9.338070881603925677e-3, + 2.680696891336819664e-3, + 6.351291219408556439e-4, + 1.239074599068830081e-4, + 1.982878843895233056e-5, + 2.589350929131392509e-6, + 2.740942840536013206e-7, + 2.332831165025738197e-8, + 1.580745574778327984e-9, + 8.427479123056716393e-11, + 3.485161234907855443e-12, + 1.099018059753451500e-13, + 2.588312664959080167e-15, + 4.437838059840028968e-17, + 5.365918308212045344e-19, + 4.393946892291604451e-21, + 2.311409794388543236e-23, + 7.274588498292248063e-26, + 1.239149701448267877e-28, + 9.832375083105887477e-32, + 2.842323553402700938e-35, + 1.878608031749515392e-39, + 8.745980440465011553e-45, + ] + ) + + neoclassics_module.kt = self.neoclassics_calc_KT() + neoclassics_module.nu = self.neoclassics_calc_nu() + neoclassics_module.nu_star = self.neoclassics_calc_nu_star() + neoclassics_module.nu_star_averaged = self.neoclassics_calc_nu_star_fromT( + iotain + ) + neoclassics_module.vd = self.neoclassics_calc_vd() + + neoclassics_module.d11_plateau = self.neoclassics_calc_D11_plateau() + + neoclassics_module.d11_mono = self.neoclassics_calc_d11_mono( + eps_effin + ) # for using epseff + + neoclassics_module.d111 = self.calc_integrated_radial_transport_coeffs(index=1) + neoclassics_module.d112 = self.calc_integrated_radial_transport_coeffs(index=2) + neoclassics_module.d113 = self.calc_integrated_radial_transport_coeffs(index=3) + + neoclassics_module.gamma_flux = self.neoclassics_calc_Gamma_flux( + neoclassics_module.densities, + neoclassics_module.temperatures, + neoclassics_module.dr_densities, + neoclassics_module.dr_temperatures, + ) + neoclassics_module.q_flux = self.neoclassics_calc_q_flux() + + def init_profile_values_from_PROCESS(self, rho): + """Initializes the profile_values object from PROCESS' parabolic profiles""" + tempe = physics_variables.te0 * (1 - rho**2) ** physics_variables.alphat * KEV + tempT = physics_variables.ti0 * (1 - rho**2) ** physics_variables.alphat * KEV + tempD = physics_variables.ti0 * (1 - rho**2) ** physics_variables.alphat * KEV + tempa = physics_variables.ti0 * (1 - rho**2) ** physics_variables.alphat * KEV + + dense = physics_variables.ne0 * (1 - rho**2) ** physics_variables.alphan + densT = ( + (1 - physics_variables.f_deuterium) + * physics_variables.ni0 + * (1 - rho**2) ** physics_variables.alphan + ) + densD = ( + physics_variables.f_deuterium + * physics_variables.ni0 + * (1 - rho**2) ** physics_variables.alphan + ) + densa = ( + physics_variables.dnalp + * (1 + physics_variables.alphan) + * (1 - rho**2) ** physics_variables.alphan + ) + + # Derivatives in real space + dr_tempe = ( + -2.0 + * 1.0 + / physics_variables.rminor + * physics_variables.te0 + * rho + * (1.0 - rho**2) ** (physics_variables.alphat - 1.0) + * physics_variables.alphat + * KEV + ) + dr_tempT = ( + -2.0 + * 1.0 + / physics_variables.rminor + * physics_variables.ti0 + * rho + * (1.0 - rho**2) ** (physics_variables.alphat - 1.0) + * physics_variables.alphat + * KEV + ) + dr_tempD = ( + -2.0 + * 1.0 + / physics_variables.rminor + * physics_variables.ti0 + * rho + * (1.0 - rho**2) ** (physics_variables.alphat - 1.0) + * physics_variables.alphat + * KEV + ) + dr_tempa = ( + -2.0 + * 1.0 + / physics_variables.rminor + * physics_variables.ti0 + * rho + * (1.0 - rho**2) ** (physics_variables.alphat - 1.0) + * physics_variables.alphat + * KEV + ) + + dr_dense = ( + -2.0 + * 1.0 + / physics_variables.rminor + * rho + * physics_variables.ne0 + * (1.0 - rho**2) ** (physics_variables.alphan - 1.0) + * physics_variables.alphan + ) + dr_densT = ( + -2.0 + * 1.0 + / physics_variables.rminor + * rho + * (1 - physics_variables.f_deuterium) + * physics_variables.ni0 + * (1.0 - rho**2) ** (physics_variables.alphan - 1.0) + * physics_variables.alphan + ) + dr_densD = ( + -2.0 + * 1.0 + / physics_variables.rminor + * rho + * physics_variables.f_deuterium + * physics_variables.ni0 + * (1.0 - rho**2) ** (physics_variables.alphan - 1.0) + * physics_variables.alphan + ) + dr_densa = ( + -2.0 + * 1.0 + / physics_variables.rminor + * rho + * physics_variables.dnalp + * (1 + physics_variables.alphan) + * (1.0 - rho**2) ** (physics_variables.alphan - 1.0) + * physics_variables.alphan + ) + + dens = np.array([dense, densD, densT, densa]) + temp = np.array([tempe, tempD, tempT, tempa]) + dr_dens = np.array([dr_dense, dr_densD, dr_densT, dr_densa]) + dr_temp = np.array([dr_tempe, dr_tempD, dr_tempT, dr_tempa]) + + return dens, temp, dr_dens, dr_temp + + def neoclassics_calc_KT(self): + """Calculates the energy on the given grid + which is given by the gauss laguerre roots. + """ + K = np.repeat((neoclassics_module.roots / KEV)[:, np.newaxis], 4, axis=1) + + return (K * neoclassics_module.temperatures).T + + def neoclassics_calc_nu(self): + """Calculates the collision frequency""" + mass = np.array( + [ + constants.emass, + constants.mproton * 2.0, + constants.mproton * 3.0, + constants.mproton * 4.0, + ] + ) + z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.echarge + + # transform the temperature back in eV + # Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962 + lnlambda = ( + 32.2 + - 1.15 * np.log10(neoclassics_module.densities[0]) + + 2.3 * np.log10(neoclassics_module.temperatures[0] / constants.echarge) + ) + + neoclassics_calc_nu = np.zeros((4, self.no_roots), order="F") + + for j in range(4): + for i in range(self.no_roots): + x = neoclassics_module.roots[i] + for k in range(4): + xk = ( + (mass[k] / mass[j]) + * ( + neoclassics_module.temperatures[j] + / neoclassics_module.temperatures[k] + ) + * x + ) + expxk = np.exp(-xk) + t = 1.0 / (1.0 + 0.3275911 * np.sqrt(xk)) + erfn = ( + 1.0 + - t + * ( + 0.254829592 + + t + * ( + -0.284496736 + + t + * (1.421413741 + t * (-1.453152027 + t * 1.061405429)) + ) + ) + * expxk + ) + phixmgx = (1.0 - 0.5 / xk) * erfn + expxk / np.sqrt(np.pi * xk) + v = np.sqrt(2.0 * x * neoclassics_module.temperatures[j] / mass[j]) + neoclassics_calc_nu[j, i] = neoclassics_calc_nu[ + j, i + ] + neoclassics_module.densities[k] * ( + z[j] * z[k] + ) ** 2 * lnlambda * phixmgx / ( + 4.0 * np.pi * constants.epsilon0**2 * mass[j] ** 2 * v**3 + ) + + return neoclassics_calc_nu + + def neoclassics_calc_nu_star(self): + """Calculates the normalized collision frequency""" + k = np.repeat(neoclassics_module.roots[:, np.newaxis], 4, axis=1) + kk = (k * neoclassics_module.temperatures).T + + mass = np.array( + [ + constants.emass, + constants.mproton * 2.0, + constants.mproton * 3.0, + constants.mproton * 4.0, + ] + ) + + v = np.empty((4, self.no_roots)) + v[0, :] = constants.speed_light * np.sqrt( + 1.0 - (kk[0, :] / (mass[0] * constants.speed_light**2) + 1) ** (-1) + ) + v[1, :] = constants.speed_light * np.sqrt( + 1.0 - (kk[1, :] / (mass[1] * constants.speed_light**2) + 1) ** (-1) + ) + v[2, :] = constants.speed_light * np.sqrt( + 1.0 - (kk[2, :] / (mass[2] * constants.speed_light**2) + 1) ** (-1) + ) + v[3, :] = constants.speed_light * np.sqrt( + 1.0 - (kk[3, :] / (mass[3] * constants.speed_light**2) + 1) ** (-1) + ) + + return ( + physics_variables.rmajor + * neoclassics_module.nu + / (neoclassics_module.iota * v) + ) + + def neoclassics_calc_nu_star_fromT(self, iota): + """Calculates the collision frequency""" + temp = ( + np.array( + [ + physics_variables.te, + physics_variables.ti, + physics_variables.ti, + physics_variables.ti, + ] + ) + * KEV + ) + density = np.array( + [ + physics_variables.dene, + physics_variables.deni * physics_variables.f_deuterium, + physics_variables.deni * (1 - physics_variables.f_deuterium), + physics_variables.dnalp, + ] + ) + + mass = np.array( + [ + constants.emass, + constants.mproton * 2.0, + constants.mproton * 3.0, + constants.mproton * 4.0, + ] + ) + z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.echarge + + # transform the temperature back in eV + # Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962 + lnlambda = ( + 32.2 + - 1.15 * np.log10(density[0]) + + 2.3 * np.log10(temp[0] / constants.echarge) + ) + + neoclassics_calc_nu_star_fromT = np.zeros((4,)) + + for j in range(4): + v = np.sqrt(2.0 * temp[j] / mass[j]) + for k in range(4): + xk = (mass[k] / mass[j]) * (temp[j] / temp[k]) + + expxk = 0.0 + if xk < 200.0: + expxk = np.exp(-xk) + + t = 1.0 / (1.0 + 0.3275911 * np.sqrt(xk)) + erfn = ( + 1.0 + - t + * ( + 0.254829592 + + t + * ( + -0.284496736 + + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429)) + ) + ) + * expxk + ) + phixmgx = (1.0 - 0.5 / xk) * erfn + expxk / np.sqrt(np.pi * xk) + neoclassics_calc_nu_star_fromT[j] = ( + neoclassics_calc_nu_star_fromT[j] + + density[k] + * (z[j] * z[k]) ** 2 + * lnlambda + * phixmgx + / (4.0 * np.pi * constants.epsilon0**2 * mass[j] ** 2 * v**4) + * physics_variables.rmajor + / iota + ) + return neoclassics_calc_nu_star_fromT + + def neoclassics_calc_vd(self): + vde = ( + neoclassics_module.roots + * neoclassics_module.temperatures[0] + / (constants.echarge * physics_variables.rmajor * physics_variables.bt) + ) + vdD = ( + neoclassics_module.roots + * neoclassics_module.temperatures[1] + / (constants.echarge * physics_variables.rmajor * physics_variables.bt) + ) + vdT = ( + neoclassics_module.roots + * neoclassics_module.temperatures[2] + / (constants.echarge * physics_variables.rmajor * physics_variables.bt) + ) + vda = ( + neoclassics_module.roots + * neoclassics_module.temperatures[3] + / ( + 2.0 + * constants.echarge + * physics_variables.rmajor + * physics_variables.bt + ) + ) + + vd = np.empty((4, self.no_roots)) + + vd[0, :] = vde + vd[1, :] = vdD + vd[2, :] = vdT + vd[3, :] = vda + + return vd + + def neoclassics_calc_D11_plateau(self): + """Calculates the plateau transport coefficients (D11_star sometimes)""" + mass = np.array( + [ + constants.emass, + constants.mproton * 2.0, + constants.mproton * 3.0, + constants.mproton * 4.0, + ] + ) + + v = np.empty((4, self.no_roots)) + v[0, :] = constants.speed_light * np.sqrt( + 1.0 + - (neoclassics_module.kt[0, :] / (mass[0] * constants.speed_light**2) + 1) + ** (-1) + ) + v[1, :] = constants.speed_light * np.sqrt( + 1.0 + - (neoclassics_module.kt[1, :] / (mass[1] * constants.speed_light**2) + 1) + ** (-1) + ) + v[2, :] = constants.speed_light * np.sqrt( + 1.0 + - (neoclassics_module.kt[2, :] / (mass[2] * constants.speed_light**2) + 1) + ** (-1) + ) + v[3, :] = constants.speed_light * np.sqrt( + 1.0 + - (neoclassics_module.kt[3, :] / (mass[3] * constants.speed_light**2) + 1) + ** (-1) + ) + + return ( + np.pi + / 4.0 + * neoclassics_module.vd**2 + * physics_variables.rmajor + / neoclassics_module.iota + / v + ) + + def neoclassics_calc_d11_mono(self, eps_eff): + """Calculates the monoenergetic radial transport coefficients + using epsilon effective + """ + return ( + 4.0 + / (9.0 * np.pi) + * (2.0 * eps_eff) ** (3.0 / 2.0) + * neoclassics_module.vd**2 + / neoclassics_module.nu + ) + + def calc_integrated_radial_transport_coeffs(self, index: int): + """Calculates the integrated radial transport coefficients (index `index`) + It uses Gauss laguerre integration + https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature + """ + return np.sum( + 2.0 + / np.sqrt(np.pi) + * neoclassics_module.d11_mono + * neoclassics_module.roots ** (index - 0.5) + * neoclassics_module.weights, + axis=1, + ) + + def neoclassics_calc_Gamma_flux( + self, densities, temperatures, dr_densities, dr_temperatures + ): + """Calculates the Energy flux by neoclassical particle transport""" + + z = np.array([-1.0, 1.0, 1.0, 2.0]) + + return ( + -densities + * neoclassics_module.d111 + * ( + (dr_densities / densities - z * neoclassics_module.er / temperatures) + + (neoclassics_module.d112 / neoclassics_module.d111 - 3.0 / 2.0) + * dr_temperatures + / temperatures + ) + ) + + def neoclassics_calc_q_flux(self): + """Calculates the Energy flux by neoclassicsal energy transport""" + + z = np.array([-1.0, 1.0, 1.0, 2.0]) + + return ( + -neoclassics_module.densities + * neoclassics_module.temperatures + * neoclassics_module.d112 + * ( + ( + neoclassics_module.dr_densities / neoclassics_module.densities + - z * neoclassics_module.er / neoclassics_module.temperatures + ) + + (neoclassics_module.d113 / neoclassics_module.d112 - 3.0 / 2.0) + * neoclassics_module.dr_temperatures + / neoclassics_module.temperatures + ) + ) diff --git a/source/fortran/constants.f90 b/source/fortran/constants.f90 index da629e98..64b7b4de 100644 --- a/source/fortran/constants.f90 +++ b/source/fortran/constants.f90 @@ -191,6 +191,9 @@ module constants ! Average number of days in a year real(dp), parameter :: n_day_year = 365.2425D0 + real(dp), parameter :: speed_light = 2.997925d+08 + !! The speed of light (m/s) + contains subroutine init_constants diff --git a/source/fortran/neoclassics_module.f90 b/source/fortran/neoclassics_module.f90 index fefc5946..afb9576e 100644 --- a/source/fortran/neoclassics_module.f90 +++ b/source/fortran/neoclassics_module.f90 @@ -60,14 +60,6 @@ module neoclassics_module real(dp), dimension(4,no_roots) :: D11_mono = 0 ! Radial monoenergetic transport coefficient on GL roots (species dependent) real(dp), dimension(4,no_roots) :: D11_plateau = 0 - ! Radial monoenergetic transport coefficient on GL roots (species dependent) - real(dp), dimension(:), allocatable :: nu_star_mono_input - ! Radial monoenergetic transport coefficient as given by the stellarator input json - ! on GL roots (species dependent) - real(dp), dimension(:), allocatable :: D11_star_mono_input - ! Radial monoenergetic transport coefficient as given by the stellarator input json - ! as function of nu_star, normalized by the plateau value. - real(dp), dimension(:), allocatable :: D13_star_mono_input ! Toroidal monoenergetic transport coefficient as given by the stellarator ! input json file as function of nu_star, normalized by the banana value. real(dp), dimension(4) :: D111 = 0 @@ -109,9 +101,6 @@ subroutine init_neoclassics_module iota = 1.0d0 D11_mono = 0 D11_plateau = 0 - nu_star_mono_input = 0 - D11_star_mono_input = 0 - D13_star_mono_input = 0 D111 = 0 D112 = 0 D113 = 0 @@ -121,531 +110,4 @@ subroutine init_neoclassics_module eps_eff = 1d-5 end subroutine init_neoclassics_module - subroutine init_neoclassics(r_eff, eps_eff, iota) - !! Constructor of the neoclassics object from the effective radius, - !! epsilon effective and iota only. - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp), intent(in) :: r_eff - real(dp), dimension(4,no_roots) :: mynu - real(dp), intent(in)::eps_eff, iota - - !! This should be called as the standard constructor - call init_profile_values_from_PROCESS(r_eff, densities, temperatures, dr_densities, dr_temperatures) - call gauss_laguerre_30_roots(roots) - call gauss_laguerre_30_weights(weights) - - - KT = neoclassics_calc_KT() - nu = neoclassics_calc_nu() - nu_star = neoclassics_calc_nu_star() - nu_star_averaged = neoclassics_calc_nu_star_fromT(iota) - vd = neoclassics_calc_vd() - - D11_plateau = neoclassics_calc_D11_plateau() - - D11_mono = neoclassics_calc_D11_mono(eps_eff) !for using epseff - - !alternatively use: = myneo%interpolate_D11_mono() ! - - D111 = neoclassics_calc_D111() - - D112 = neoclassics_calc_D112() - D113 = neoclassics_calc_D113() - - Gamma_flux = neoclassics_calc_Gamma_flux(densities, temperatures, dr_densities, dr_temperatures) - q_flux = neoclassics_calc_q_flux() - - ! Return: - - end subroutine init_neoclassics - - - function neoclassics_calc_KT() result(KK) - !! Calculates the energy on the given grid - !! which is given by the gauss laguerre roots. - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: e_,c_,me_,mp_,keV_ - - real(dp), dimension(no_roots) :: K - real(dp), dimension(4,no_roots) :: KK - - K = roots/keV_ - - KK(1,:) = K * temperatures(1) ! electrons - KK(2,:) = K * temperatures(2) ! deuterium - KK(3,:) = K * temperatures(3) ! tritium - KK(4,:) = K * temperatures(4) ! helium - - end function neoclassics_calc_KT - - function neoclassics_calc_Gamma_flux(densities, temperatures, dr_densities, dr_temperatures) - !! Calculates the Energy flux by neoclassical particle transport - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp),dimension(4) :: neoclassics_calc_Gamma_flux, densities, dr_densities, z, temperatures, dr_temperatures - - - z = (/-1.0,1.0,1.0,2.0/) - - neoclassics_calc_Gamma_flux = - densities * D111 * ((dr_densities/densities - z * Er/temperatures)+ & - (D112/D111-3.0/2.0) * dr_temperatures/temperatures ) - - end function neoclassics_calc_Gamma_flux - - function neoclassics_calc_q_flux() - !! Calculates the Energy flux by neoclassicsal energy transport - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp),dimension(4) :: z, neoclassics_calc_q_flux - - - ! densities = densities - ! temps = temperatures - ! dr_densities = dr_densities - ! dr_temps = dr_temperatures - - z = (/-1.0,1.0,1.0,2.0/) - - q_flux = - densities * temperatures * D112 * ((dr_densities/densities - z * Er/temperatures) + & - (D113/D112-3.0/2.0) * dr_temperatures/temperatures ) - - neoclassics_calc_q_flux = q_flux - end function neoclassics_calc_q_flux - - function neoclassics_calc_D11_mono(eps_eff) result(D11_mono) - !! Calculates the monoenergetic radial transport coefficients - !! using epsilon effective. - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi - - real(dp),dimension(4,no_roots) :: D11_mono - real(dp), intent(in):: eps_eff - - D11_mono = 4.0d0/(9.0d0*pi) * (2.0d0 * eps_eff)**(3.0d0/2.0d0) & - * vd**2/nu - - end function neoclassics_calc_D11_mono - - function neoclassics_calc_D11_plateau() result(D11_plateau) - !! Calculates the plateau transport coefficients (D11_star sometimes) - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi, me_, mp_, c_ - use physics_variables, only: rmajor - - real(dp),dimension(4,no_roots) :: D11_plateau, v - real(dp),dimension(4) :: mass - - mass = (/me_,mp_*2.0d0,mp_*3.0d0,mp_*4.0d0/) - - v(1,:) = c_ * sqrt(1.0d0-(KT(1,:)/(mass(1) * c_**2)+1)**(-1)) - v(2,:) = c_ * sqrt(1.0d0-(KT(2,:)/(mass(2) * c_**2)+1)**(-1)) - v(3,:) = c_ * sqrt(1.0d0-(KT(3,:)/(mass(3) * c_**2)+1)**(-1)) - v(4,:) = c_ * sqrt(1.0d0-(KT(4,:)/(mass(4) * c_**2)+1)**(-1)) - - D11_plateau = pi/4.0 * vd**2 * rmajor/ iota / v - - end function neoclassics_calc_D11_plateau - - function neoclassics_interpolate_D11_mono() result(D11_mono) - !! Interpolates the D11 coefficients on the Gauss laguerre grid - !! (This method is unused as of now, but is needed when taking D11 explicitely as input) - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use const_and_precisions, only: pi - use maths_library, only: find_y_nonuniform_x - ! use grad_func, only: interp1_ef - - real(dp),dimension(4,no_roots) :: D11_mono - integer :: ii,jj - - do ii = 1,4 - do jj = 1,no_roots - D11_mono(ii,jj) = find_y_nonuniform_x(nu_star(ii,jj),nu_star_mono_input, & - D11_star_mono_input,size(nu_star_mono_input)) * & - D11_plateau(ii,jj) - end do - end do - - end function neoclassics_interpolate_D11_mono - - function neoclassics_calc_vd() - !! Calculates the drift velocities - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use const_and_precisions, only: e_ - use physics_variables, only: rmajor, bt - - real(dp), dimension(no_roots) :: vde,vdT,vdD,vda, K - real(dp), dimension(4,no_roots) :: vd,neoclassics_calc_vd - - K = roots - - vde = K * temperatures(1)/(e_ * rmajor * bt) - vdD = K * temperatures(2)/(e_ * rmajor * bt) - vdT = K * temperatures(3)/(e_ * rmajor * bt) - vda = K * temperatures(4)/(2.0*e_ * rmajor * bt) - - vd(1,:) = vde - vd(2,:) = vdD - vd(3,:) = vdT - vd(4,:) = vda - - neoclassics_calc_vd = vd - end function neoclassics_calc_vd - - function neoclassics_calc_nu_star() result(nu_star) - !! Calculates the normalized collision frequency - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use const_and_precisions, only: e_,c_,me_,mp_ - use physics_variables, only: rmajor - - real(dp), dimension(no_roots) :: K - real(dp), dimension(4,no_roots) :: v,nu_star,KK - real(dp), dimension(4) :: mass - - K = roots - - KK(1,:) = K * temperatures(1) - KK(2,:) = K * temperatures(2) - KK(3,:) = K * temperatures(3) - KK(4,:) = K * temperatures(4) - - mass = (/me_,mp_*2.0d0,mp_*3.0d0,mp_*4.0d0/) - - v(1,:) = c_ * sqrt(1.0d0-(KK(1,:)/(mass(1) * c_**2)+1)**(-1)) - v(2,:) = c_ * sqrt(1.0d0-(KK(2,:)/(mass(2) * c_**2)+1)**(-1)) - v(3,:) = c_ * sqrt(1.0d0-(KK(3,:)/(mass(3) * c_**2)+1)**(-1)) - v(4,:) = c_ * sqrt(1.0d0-(KK(4,:)/(mass(4) * c_**2)+1)**(-1)) - - nu_star = rmajor * nu/(iota*v) - - end function neoclassics_calc_nu_star - - - function neoclassics_calc_nu_star_fromT(iota) - !! Calculates the collision frequency - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi, me_, mp_, eps0_,e_, keV_ - use physics_variables, only: rmajor, te,ti, dene,deni, dnalp, f_deuterium - - real(dp),dimension(4) :: neoclassics_calc_nu_star_fromT - real(dp) :: t,erfn,phixmgx,expxk,xk, lnlambda,x,v - real(dp),dimension(4) :: temp, mass,density,z - real(dp) :: iota - - integer :: jj,kk - - temp = (/te,ti,ti,ti /) * keV_ - density = (/dene,deni * f_deuterium,deni*(1-f_deuterium),dnalp /) - - ! e D T a (He) - mass = (/me_,mp_*2.0d0,mp_*3.0d0,mp_*4.0d0/) - z = (/-1.0d0,1.0d0,1.0d0,2.0d0/) * e_ - - ! transform the temperature back in eV - ! Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962 - lnlambda = 32.2d0 - 1.15d0*log10(density(1)) + 2.3d0*log10(temp(1)/e_) - - neoclassics_calc_nu_star_fromT(:) = 0.0d0 - - do jj = 1, 4 - v = sqrt(2d0 * temp(jj)/mass(jj)) - do kk = 1,4 - xk = (mass(kk)/mass(jj))*(temp(jj)/temp(kk)) - - if (xk < 200.d0) then - expxk = exp(-xk) - else - expxk = 0.0d0 - endif - - t = 1.0d0/(1.0d0+0.3275911d0*sqrt(xk)) - erfn = 1.0d0-t*(.254829592d0 + t*(-.284496736d0 + t*(1.421413741d0 & - + t*(-1.453152027d0 +t*1.061405429d0))))*expxk - phixmgx = (1.0d0-0.5d0/xk)*erfn + expxk/sqrt(pi*xk) - neoclassics_calc_nu_star_fromT(jj) = neoclassics_calc_nu_star_fromT(jj) + density(kk)*(z(jj)*z(kk))**2 & - *lnlambda*phixmgx/(4.0d0*pi*eps0_**2*mass(jj)**2*v**4) * rmajor/iota - enddo - enddo - - end function neoclassics_calc_nu_star_fromT - - function neoclassics_calc_D111() - !! Calculates the integrated radial transport coefficients (index 1) - !! It uses Gauss laguerre integration - !! https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi - - real(dp),dimension(4) :: D111, neoclassics_calc_D111 - - real(dp),dimension(no_roots) :: xi,wi - - xi = roots - wi = weights - - D111(1) = sum(2.0d0/sqrt(pi) * D11_mono(1,:) * xi**(1.0d0-0.5d0) * wi) - D111(2) = sum(2.0d0/sqrt(pi) * D11_mono(2,:) * xi**(1.0d0-0.5d0) * wi) - D111(3) = sum(2.0d0/sqrt(pi) * D11_mono(3,:) * xi**(1.0d0-0.5d0) * wi) - D111(4) = sum(2.0d0/sqrt(pi) * D11_mono(4,:) * xi**(1.0d0-0.5d0) * wi) - - neoclassics_calc_D111 = D111 - - end function neoclassics_calc_D111 - - function neoclassics_calc_D112() - !! Calculates the integrated radial transport coefficients (index 2) - !! It uses Gauss laguerre integration - !! https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi - - real(dp),dimension(4) :: D112, neoclassics_calc_D112 - - real(dp),dimension(no_roots) :: xi,wi - - xi = roots - wi = weights - - D112(1) = sum(2.0d0/sqrt(pi) * D11_mono(1,:) * xi**(2.0d0-0.5d0) * wi) - D112(2) = sum(2.0d0/sqrt(pi) * D11_mono(2,:) * xi**(2.0d0-0.5d0) * wi) - D112(3) = sum(2.0d0/sqrt(pi) * D11_mono(3,:) * xi**(2.0d0-0.5d0) * wi) - D112(4) = sum(2.0d0/sqrt(pi) * D11_mono(4,:) * xi**(2.0d0-0.5d0) * wi) - - neoclassics_calc_D112 = D112 - - end function neoclassics_calc_D112 - - function neoclassics_calc_D113() - !! Calculates the integrated radial transport coefficients (index 3) - !! It uses Gauss laguerre integration - !! https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi - - real(dp),dimension(4) :: D113, neoclassics_calc_D113 - - real(dp),dimension(no_roots) :: xi,wi - - xi = roots - wi = weights - - D113(1) = sum(2.0d0/sqrt(pi) * D11_mono(1,:) * xi**(3.0d0-0.5d0) * wi) - D113(2) = sum(2.0d0/sqrt(pi) * D11_mono(2,:) * xi**(3.0d0-0.5d0) * wi) - D113(3) = sum(2.0d0/sqrt(pi) * D11_mono(3,:) * xi**(3.0d0-0.5d0) * wi) - D113(4) = sum(2.0d0/sqrt(pi) * D11_mono(4,:) * xi**(3.0d0-0.5d0) * wi) - - neoclassics_calc_D113 = D113 - end function neoclassics_calc_D113 - - function neoclassics_calc_nu() - !! Calculates the collision frequency - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use const_and_precisions, only: pi, me_, mp_, eps0_,e_ - - real(dp),dimension(4,no_roots) :: neoclassics_calc_nu - real(dp) :: t,erfn,phixmgx,expxk,xk, lnlambda,x,v - real(dp),dimension(4) :: temp, mass,density,z - - integer :: jj,ii,kk - - temp = temperatures - density = densities - - ! e D T a (He) - mass = (/me_,mp_*2.0d0,mp_*3.0d0,mp_*4.0d0/) - z = (/-1.0d0,1.0d0,1.0d0,2.0d0/) * e_ - - ! transform the temperature back in eV - ! Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962 - lnlambda = 32.2d0 - 1.15d0*log10(density(1)) + 2.3d0*log10(temp(1)/e_) - - neoclassics_calc_nu(:,:) = 0.0 - - do jj = 1, 4 - do ii = 1, no_roots - x = roots(ii) - do kk = 1,4 - xk = (mass(kk)/mass(jj))*(temp(jj)/temp(kk))*x - expxk = exp(-xk) - t = 1.0d0/(1.0d0+0.3275911d0*sqrt(xk)) - erfn = 1.0d0-t*(.254829592d0 + t*(-.284496736d0 + t*(1.421413741d0 & - + t*(-1.453152027d0 +t*1.061405429d0))))*expxk - phixmgx = (1.0d0-0.5d0/xk)*erfn + expxk/sqrt(pi*xk) - v = sqrt(2.0d0*x*temp(jj)/mass(jj)) - neoclassics_calc_nu(jj,ii) = neoclassics_calc_nu(jj,ii) + density(kk)*(z(jj)*z(kk))**2 & - *lnlambda *phixmgx/(4.0*pi*eps0_**2*mass(jj)**2*v**3) - enddo - enddo - enddo - - end function neoclassics_calc_nu - - subroutine init_profile_values_from_PROCESS(rho, densities, temperatures, dr_densities, dr_temperatures) - !! Initializes the profile_values object from PROCESS' parabolic profiles - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use physics_variables, only: ne0,te0,alphan,& - alphat,ti0,ni0,f_deuterium, dnalp, rminor - use const_and_precisions, only: keV_ - - real(dp), intent(in) :: rho - - real(dp),dimension(4) :: dens,temp, dr_dens, dr_temp - real(dp) :: dense, densD,densT,densa, & - tempD,tempT,tempa,tempe, & - dr_tempe, dr_tempT, dr_tempD, dr_tempa,& - dr_dense, dr_densT, dr_densD, dr_densa, r - real(dp), dimension(4), intent(out):: densities - real(dp), dimension(4), intent(out):: temperatures - real(dp), dimension(4), intent(out):: dr_densities - real(dp), dimension(4), intent(out):: dr_temperatures - - r = rho * rminor - - tempe = te0 * (1-rho**2)**alphat * keV_ ! To SI units bc.. convenience I guess? - tempT = ti0 * (1-rho**2)**alphat * keV_ - tempD = ti0 * (1-rho**2)**alphat * keV_ - tempa = ti0 * (1-rho**2)**alphat * keV_ - - dense = ne0 * (1-rho**2)**alphan - densT = (1-f_deuterium) * ni0 * (1-rho**2)**alphan - densD = f_deuterium *ni0 * (1-rho**2)**alphan - densa = dnalp*(1+alphan) * (1-rho**2)**alphan - - ! Derivatives in real space - dr_tempe = -2.0d0 * 1.0d0/rminor * te0 * rho * (1.0d0-rho**2)**(alphat-1.0d0) * alphat * keV_ - dr_tempT = -2.0d0 * 1.0d0/rminor * ti0 * rho * (1.0d0-rho**2)**(alphat-1.0d0) * alphat * keV_ - dr_tempD = -2.0d0 * 1.0d0/rminor * ti0 * rho * (1.0d0-rho**2)**(alphat-1.0d0) * alphat * keV_ - dr_tempa = -2.0d0 * 1.0d0/rminor * ti0 * rho * (1.0d0-rho**2)**(alphat-1.0d0) * alphat * keV_ - - dr_dense = -2.0d0 * 1.0d0/rminor * rho * ne0 * (1.0d0-rho**2)**(alphan-1.0d0) * alphan - dr_densT = -2.0d0 * 1.0d0/rminor * rho * (1-f_deuterium) * ni0 * (1.0d0-rho**2)**(alphan-1.0d0) * alphan - dr_densD = -2.0d0 * 1.0d0/rminor * rho * f_deuterium *ni0 * (1.0d0-rho**2)**(alphan-1.0d0) * alphan - dr_densa = -2.0d0 * 1.0d0/rminor * rho * dnalp*(1+alphan)* (1.0d0-rho**2)**(alphan-1.0d0) * alphan - - dens(1) = dense - dens(2) = densD - dens(3) = densT - dens(4) = densa - - temp(1) = tempe - temp(2) = tempD - temp(3) = tempT - temp(4) = tempa - - dr_dens(1) = dr_dense - dr_dens(2) = dr_densD - dr_dens(3) = dr_densT - dr_dens(4) = dr_densa - - dr_temp(1) = dr_tempe - dr_temp(2) = dr_tempD - dr_temp(3) = dr_tempT - dr_temp(4) = dr_tempa - - densities = dens - temperatures = temp - dr_densities = dr_dens - dr_temperatures = dr_temp - - end subroutine init_profile_values_from_PROCESS - - subroutine gauss_laguerre_30_roots(roots) - - !! Sets the gauss Laguerre roots and weights for 30 - !! discretization points. Used for integration in this module. - !! roots - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp), dimension(no_roots), intent(out):: roots - roots= (/4.740718054080526184d-02,& - 2.499239167531593919d-01,& - 6.148334543927683749d-01,& - 1.143195825666101451d+00,& - 1.836454554622572344d+00,& - 2.696521874557216147d+00,& - 3.725814507779509288d+00,& - 4.927293765849881879d+00,& - 6.304515590965073635d+00,& - 7.861693293370260349d+00,& - 9.603775985479263255d+00,& - 1.153654659795613924d+01,& - 1.366674469306423489d+01,& - 1.600222118898106771d+01,& - 1.855213484014315029d+01,& - 2.132720432178312819d+01,& - 2.434003576453269346d+01,& - 2.760555479678096091d+01,& - 3.114158670111123683d+01,& - 3.496965200824907072d+01,& - 3.911608494906788991d+01,& - 4.361365290848483056d+01,& - 4.850398616380419980d+01,& - 5.384138540650750571d+01,& - 5.969912185923549686d+01,& - 6.618061779443848991d+01,& - 7.344123859555988076d+01,& - 8.173681050672767867d+01,& - 9.155646652253683726d+01,& - 1.041575244310588886d+02/) - - end subroutine gauss_laguerre_30_roots - - subroutine gauss_laguerre_30_weights(weights) - - real(dp), dimension(no_roots), intent(out):: weights - weights = (/1.160440860204388913d-01,& - 2.208511247506771413d-01,& - 2.413998275878537214d-01,& - 1.946367684464170855d-01,& - 1.237284159668764899d-01,& - 6.367878036898660943d-02,& - 2.686047527337972682d-02,& - 9.338070881603925677d-03,& - 2.680696891336819664d-03,& - 6.351291219408556439d-04,& - 1.239074599068830081d-04,& - 1.982878843895233056d-05,& - 2.589350929131392509d-06,& - 2.740942840536013206d-07,& - 2.332831165025738197d-08,& - 1.580745574778327984d-09,& - 8.427479123056716393d-11,& - 3.485161234907855443d-12,& - 1.099018059753451500d-13,& - 2.588312664959080167d-15,& - 4.437838059840028968d-17,& - 5.365918308212045344d-19,& - 4.393946892291604451d-21,& - 2.311409794388543236d-23,& - 7.274588498292248063d-26,& - 1.239149701448267877d-28,& - 9.832375083105887477d-32,& - 2.842323553402700938d-35,& - 1.878608031749515392d-39,& - 8.745980440465011553d-45/) - - - end subroutine gauss_laguerre_30_weights - - end module neoclassics_module diff --git a/tests/unit/test_neoclassics.py b/tests/unit/test_neoclassics.py new file mode 100644 index 00000000..7bc762ff --- /dev/null +++ b/tests/unit/test_neoclassics.py @@ -0,0 +1,1053 @@ +import pytest +import numpy +from typing import NamedTuple, Any + +from process.fortran import neoclassics_module, physics_variables +from process.stellarator import Neoclassics + + +@pytest.fixture +def neoclassics(): + """Provides Neoclassics object for testing. + + :returns: initialised Neoclassics object + :rtype: process.stellerator.Neoclassics + """ + return Neoclassics() + + +class InitNeoclassicsParam(NamedTuple): + ne0: Any = None + te0: Any = None + alphan: Any = None + alphat: Any = None + ti0: Any = None + ni0: Any = None + fdeut: Any = None + dnalp: Any = None + rminor: Any = None + rmajor: Any = None + bt: Any = None + te: Any = None + ti: Any = None + dene: Any = None + deni: Any = None + densities: Any = None + temperatures: Any = None + dr_densities: Any = None + dr_temperatures: Any = None + roots: Any = None + weights: Any = None + nu: Any = None + nu_star: Any = None + nu_star_averaged: Any = None + vd: Any = None + iota: Any = None + q_flux: Any = None + eps_eff: Any = None + r_eff: Any = None + r_effin: Any = None + eps_effin: Any = None + iotain: Any = None + expected_densities: Any = None + expected_temperatures: Any = None + expected_dr_densities: Any = None + expected_dr_temperatures: Any = None + expected_roots: Any = None + expected_weights: Any = None + expected_nu: Any = None + expected_nu_star: Any = None + expected_nu_star_averaged: Any = None + expected_vd: Any = None + expected_q_flux: Any = None + + +@pytest.mark.parametrize( + "initneoclassicsparam", + ( + InitNeoclassicsParam( + ne0=2.7956610000000002e20, + te0=13.241800000000001, + alphan=0.35000000000000003, + alphat=1.2, + ti0=12.579710000000002, + ni0=2.3930858160000005e20, + fdeut=0.5, + dnalp=2.9820384000000004e19, + rminor=1.7993820274145451, + rmajor=22.16, + bt=5.2400000000000002, + te=6.0190000000000001, + ti=5.7180500000000007, + dene=2.07086e20, + deni=1.47415411616e20, + densities=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + temperatures=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + dr_densities=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + dr_temperatures=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + roots=numpy.array( + numpy.array( + ( + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ), + order="F", + ), + order="F", + ).transpose(), + weights=numpy.array( + numpy.array( + ( + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ), + order="F", + ), + order="F", + ).transpose(), + nu=numpy.array( + ( + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + ), + order="F", + ).transpose(), + nu_star=numpy.array( + ( + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + ), + order="F", + ).transpose(), + nu_star_averaged=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + vd=numpy.array( + ( + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + (0, 0, 0, 0), + ), + order="F", + ).transpose(), + iota=1, + q_flux=numpy.array( + numpy.array((0, 0, 0, 0), order="F"), order="F" + ).transpose(), + eps_eff=1.0000000000000001e-05, + r_eff=0, + r_effin=0.59999999999999998, + eps_effin=0.01464553, + iotain=0.90000000000000002, + expected_densities=( + 2.391373976836772e20, + 1.0235080620861384e20, + 1.0235080620861384e20, + 3.4435785266449519e19, + ), + expected_temperatures=( + 1.2416608937807636e-15, + 1.1795778490917256e-15, + 1.1795778490917256e-15, + 1.1795778490917256e-15, + ), + expected_dr_densities=( + -8.7215452215783653e19, + -3.7328213548355412e19, + -3.7328213548355412e19, + -1.2559025119072848e19, + ), + expected_dr_temperatures=( + -1.5526091560561597e-15, + -1.4749786982533515e-15, + -1.4749786982533515e-15, + -1.4749786982533515e-15, + ), + expected_roots=numpy.array( + numpy.array( + ( + 0.047407180540805262, + 0.24992391675315939, + 0.61483345439276837, + 1.1431958256661015, + 1.8364545546225723, + 2.6965218745572161, + 3.7258145077795093, + 4.9272937658498819, + 6.3045155909650736, + 7.8616932933702603, + 9.6037759854792633, + 11.536546597956139, + 13.666744693064235, + 16.002221188981068, + 18.55213484014315, + 21.327204321783128, + 24.340035764532693, + 27.605554796780961, + 31.141586701111237, + 34.969652008249071, + 39.11608494906789, + 43.613652908484831, + 48.5039861638042, + 53.841385406507506, + 59.699121859235497, + 66.18061779443849, + 73.441238595559881, + 81.736810506727679, + 91.556466522536837, + 104.15752443105889, + ), + order="F", + ), + order="F", + ).transpose(), + expected_weights=( + 0.11604408602043889, + 0.22085112475067714, + 0.24139982758785372, + 0.19463676844641709, + 0.12372841596687649, + 0.063678780368986609, + 0.026860475273379727, + 0.0093380708816039257, + 0.0026806968913368197, + 0.00063512912194085564, + 0.00012390745990688301, + 1.9828788438952331e-05, + 2.5893509291313925e-06, + 2.7409428405360132e-07, + 2.3328311650257382e-08, + 1.580745574778328e-09, + 8.4274791230567164e-11, + 3.4851612349078554e-12, + 1.0990180597534515e-13, + 2.5883126649590802e-15, + 4.437838059840029e-17, + 5.3659183082120453e-19, + 4.3939468922916045e-21, + 2.3114097943885432e-23, + 7.2745884982922481e-26, + 1.2391497014482679e-28, + 9.8323750831058875e-32, + 2.8423235534027009e-35, + 1.8786080317495154e-39, + 8.7459804404650116e-45, + ), + expected_nu=numpy.array( + ( + ( + 3695464.938733798, + 11862.660224787611, + 7933.188679624609, + 23814.418394017364, + ), + ( + 343316.72274464089, + 2123.530294833155, + 1450.1975058356002, + 4403.93523774, + ), + ( + 97325.601943625719, + 782.90996572377162, + 550.46150951687332, + 1698.5476568456329, + ), + ( + 40967.943564057088, + 373.3569294227749, + 270.83070812853452, + 851.34936648572614, + ), + ( + 21050.27896218509, + 204.31685155512426, + 152.39427303068328, + 487.91926000009528, + ), + ( + 12196.141744781391, + 122.4627549323195, + 93.37151291439298, + 303.82292065987531, + ), + ( + 7662.855304124023, + 78.497186099853351, + 60.822352763230512, + 200.532111745934, + ), + ( + 5107.6181124745644, + 53.008588809359139, + 41.543171451498225, + 138.36568682166154, + ), + ( + 3562.1896935110958, + 37.316304271771997, + 29.478466134989471, + 98.931218624509825, + ), + ( + 2575.1493019063892, + 27.17114569450117, + 21.584157832183529, + 72.845381503203299, + ), + ( + 1916.5626815978521, + 20.341237798241099, + 16.222662407259669, + 54.977961633660357, + ), + ( + 1461.0373149373809, + 15.58437516643423, + 12.464365949260028, + 42.371363838150558, + ), + ( + 1136.3293377214109, + 12.174519399214853, + 9.7573846405961202, + 33.245772325916597, + ), + ( + 898.87370477934678, + 9.6691588841578611, + 7.7612618185478333, + 26.490662053885323, + ), + ( + 721.36428722004121, + 7.7886363224542317, + 6.2588248584263582, + 21.391037801763538, + ), + ( + 586.10784951621235, + 6.3505475050965412, + 5.1074008781905587, + 17.473656809416916, + ), + ( + 481.30424894598372, + 5.2325813464582778, + 4.2107797561372902, + 14.417477347614785, + ), + ( + 398.88189106951904, + 4.3506921903280844, + 3.5025544153548962, + 11.999840321777146, + ), + ( + 333.19454764267971, + 3.6458582076791037, + 2.9359214172811088, + 10.063206020290355, + ), + ( + 280.21278073497348, + 3.0758139209321369, + 2.4772660955148873, + 8.494068483233054, + ), + ( + 237.00858514975135, + 2.6097529846718648, + 2.1020279060352189, + 7.2092769760436468, + ), + ( + 201.41883814132996, + 2.2248531750213925, + 1.7919748195144374, + 6.1469708670487027, + ), + ( + 171.82051852815496, + 1.9039461312634718, + 1.5333681490810633, + 5.260448904005937, + ), + ( + 146.97724533175426, + 1.6339205441745102, + 1.3157006442816281, + 4.5139376970122962, + ), + ( + 125.93193667489498, + 1.4046010159515565, + 1.130807807982801, + 3.8796018446013152, + ), + ( + 107.92910284917983, + 1.2079326736732123, + 0.97221975350348488, + 3.335356496099501, + ), + ( + 92.354617547302752, + 1.0373444948341899, + 0.83465363159517136, + 2.8631484655280879, + ), + ( + 78.68026674668576, + 0.88715494590298993, + 0.71353811705010739, + 2.447338816503998, + ), + ( + 66.385967521689778, + 0.75172062458023658, + 0.60432964829100722, + 2.0723651470862698, + ), + ( + 54.725799270711548, + 0.62283558834560981, + 0.50041870619678697, + 1.7155597347573166, + ), + ), + order="F", + ).transpose(), + expected_nu_star=numpy.array( + ( + ( + 10.191074937923535, + 2.0331966316122507, + 1.6652931278769849, + 5.7723445649189209, + ), + ( + 0.41298040929536733, + 0.15851635531154254, + 0.13258314893305531, + 0.4649127239401013, + ), + ( + 0.074848077541108754, + 0.037260863337660155, + 0.03208584666735572, + 0.11432298738955823, + ), + ( + 0.023197077266323505, + 0.01303119334609255, + 0.011577190579052517, + 0.042022540114017584, + ), + ( + 0.0094525677652312987, + 0.0056264565943291555, + 0.0051397801137085886, + 0.019001722106245209, + ), + ( + 0.0045482132226884351, + 0.0027830670017342913, + 0.0025988361335270451, + 0.0097645901733738206, + ), + ( + 0.0024492474025462241, + 0.001517630205259042, + 0.0014401895657356347, + 0.0054828859338448806, + ), + ( + 0.0014317922467676544, + 0.00089118019395046364, + 0.00085538801619877453, + 0.003289732501027202, + ), + ( + 0.00089132611440290003, + 0.00055462245458156732, + 0.00053659613981289387, + 0.0020794309035960076, + ), + ( + 0.00058320422633617944, + 0.00036163954478676596, + 0.00035184101734703275, + 0.0013711394734907242, + ), + ( + 0.0003973238458688693, + 0.00024495361972729277, + 0.00023926078180568864, + 0.00093628060138880047, + ), + ( + 0.00027986715624797473, + 0.00017123030253161586, + 0.00016772730703600486, + 0.00065837564428247015, + ), + ( + 0.00020271700094666678, + 0.0001228995355882499, + 0.00012063507255191103, + 0.00047461779631969821, + ), + ( + 0.00015035093166913791, + 9.0205186378467152e-05, + 8.8678069715659373e-05, + 0.0003494972842656118, + ), + ( + 0.00011379147025090292, + 6.7483810604850391e-05, + 6.6415788403404255e-05, + 0.00026210609671657975, + ), + ( + 8.7635514886206191e-05, + 5.1319413760932046e-05, + 5.0548705442571374e-05, + 0.00019969180051899852, + ), + ( + 6.8517128585695404e-05, + 3.9581790112080437e-05, + 3.9010425484079875e-05, + 0.00015423159537601957, + ), + ( + 5.4275274472921955e-05, + 3.0903164171056251e-05, + 3.0469636469663725e-05, + 0.00012053781020124761, + ), + ( + 4.3485101647524803e-05, + 2.4382332471219346e-05, + 2.4046766615879637e-05, + 9.5172958985582652e-05, + ), + ( + 3.5184659036602886e-05, + 1.9411697131159546e-05, + 1.914748805113672e-05, + 7.5808745485372899e-05, + ), + ( + 2.8710573684483027e-05, + 1.5573074040919529e-05, + 1.5362006532190414e-05, + 6.0836608331866169e-05, + ), + ( + 2.3596485999114507e-05, + 1.2573222214622194e-05, + 1.2402532788994784e-05, + 4.9125028134176529e-05, + ), + ( + 1.9508685715836919e-05, + 1.0202969147411138e-05, + 1.0063526905956131e-05, + 3.9864751628194759e-05, + ), + ( + 1.6204365827687939e-05, + 8.310707531476441e-06, + 8.1958593174412469e-06, + 3.246793683964068e-05, + ), + ( + 1.3503921608329314e-05, + 6.7848331921788916e-06, + 6.6896542854973157e-06, + 2.6501040497028117e-05, + ), + ( + 1.1272082960471908e-05, + 5.5418275702450897e-06, + 5.462627168620374e-06, + 2.1639107879962109e-05, + ), + ( + 9.4045207292483327e-06, + 4.5178819909572341e-06, + 4.4518747990969652e-06, + 1.7633538292126709e-05, + ), + ( + 7.8173435227433508e-06, + 3.6625163533856096e-06, + 3.6076109486006602e-06, + 1.4287438740418093e-05, + ), + ( + 6.4359614189096637e-06, + 2.93230577098992e-06, + 2.8869969639588791e-06, + 1.1431285631297126e-05, + ), + ( + 5.1694249624199883e-06, + 2.2779077252622501e-06, + 2.2413635447347678e-06, + 8.8723608152712553e-06, + ), + ), + order="F", + ).transpose(), + expected_nu_star_averaged=numpy.array( + numpy.array( + ( + 0.030792808467814459, + 0.01661476041777428, + 0.014667916022065946, + 0.053033247562414364, + ), + order="F", + ), + order="F", + ).transpose(), + expected_vd=numpy.array( + ( + ( + 3.1645071193768142, + 3.0062817634079737, + 3.0062817634079737, + 1.5031408817039869, + ), + ( + 16.682831690173259, + 15.848690105664597, + 15.848690105664597, + 7.9243450528322983, + ), + ( + 41.04114232193708, + 38.989085205840233, + 38.989085205840233, + 19.494542602920117, + ), + ( + 76.310197904479367, + 72.494688009255412, + 72.494688009255412, + 36.247344004627706, + ), + ( + 122.58635603762467, + 116.45703823574345, + 116.45703823574345, + 58.228519117871727, + ), + ( + 179.99726143272298, + 170.99739836108685, + 170.99739836108685, + 85.498699180543426, + ), + ( + 248.70423427095045, + 236.26902255740296, + 236.26902255740296, + 118.13451127870148, + ), + ( + 328.90494695992078, + 312.45969961192475, + 312.45969961192475, + 156.22984980596237, + ), + ( + 420.8367644782997, + 399.79492625438479, + 399.79492625438479, + 199.8974631271924, + ), + ( + 524.78093220104745, + 498.54188559099515, + 498.54188559099515, + 249.27094279549758, + ), + ( + 641.06781150569088, + 609.01442093040646, + 609.01442093040646, + 304.50721046520323, + ), + ( + 770.08342250666237, + 731.57925138132941, + 731.57925138132941, + 365.78962569066471, + ), + ( + 912.27764204794528, + 866.66375994554812, + 866.66375994554812, + 433.33187997277406, + ), + ( + 1068.17453180507, + 1014.7658052148167, + 1014.7658052148167, + 507.38290260740837, + ), + ( + 1238.3854536706519, + 1176.4661809871195, + 1176.4661809871195, + 588.23309049355976, + ), + ( + 1423.6258968110376, + 1352.4446019704858, + 1352.4446019704858, + 676.2223009852429, + ), + ( + 1624.7373411386941, + 1543.5004740817596, + 1543.5004740817596, + 771.75023704087982, + ), + ( + 1842.7160968488222, + 1750.5802920063811, + 1750.5802920063811, + 875.29014600319056, + ), + ( + 2078.7520308138296, + 1974.8144292731388, + 1974.8144292731388, + 987.40721463656939, + ), + ( + 2334.2816737853227, + 2217.5675900960568, + 2217.5675900960568, + 1108.7837950480284, + ), + ( + 2611.0628789014031, + 2480.5097349563334, + 2480.5097349563334, + 1240.2548674781667, + ), + ( + 2911.2829228925325, + 2765.7187767479063, + 2765.7187767479063, + 1382.8593883739532, + ), + ( + 3237.7206951043381, + 3075.8346603491214, + 3075.8346603491214, + 1537.9173301745607, + ), + ( + 3594.0008558270988, + 3414.3008130357439, + 3414.3008130357439, + 1707.1504065178719, + ), + ( + 3985.0143794458554, + 3785.763660473563, + 3785.763660473563, + 1892.8818302367815, + ), + ( + 4417.6648724129982, + 4196.7816287923488, + 4196.7816287923488, + 2098.3908143961744, + ), + ( + 4902.3232291036566, + 4657.2070676484745, + 4657.2070676484745, + 2328.6035338242373, + ), + ( + 5456.0662712488684, + 5183.2629576864265, + 5183.2629576864265, + 2591.6314788432132, + ), + ( + 6111.5444291433778, + 5805.9672076862098, + 5805.9672076862098, + 2902.9836038431049, + ), + ( + 6952.6857290119615, + 6605.0514425613637, + 6605.0514425613637, + 3302.5257212806819, + ), + ), + order="F", + ).transpose(), + expected_q_flux=numpy.array( + numpy.array( + ( + 13499.211058929142, + 472549.72072198224, + 598376.09017838724, + 14996.485377861178, + ), + order="F", + ), + order="F", + ).transpose(), + ), + ), +) +def test_init_neoclassics(initneoclassicsparam, monkeypatch, neoclassics): + """ + Automatically generated Regression Unit Test for init_neoclassics. + + This test was generated using data from stellarator/stellarator.IN.DAT. + + :param initneoclassicsparam: the data used to mock and assert in this test. + :type initneoclassicsparam: initneoclassicsparam + + :param monkeypatch: pytest fixture used to mock module/class variables + :type monkeypatch: _pytest.monkeypatch.monkeypatch + """ + monkeypatch.setattr(physics_variables, "ne0", initneoclassicsparam.ne0) + monkeypatch.setattr(physics_variables, "te0", initneoclassicsparam.te0) + monkeypatch.setattr(physics_variables, "alphan", initneoclassicsparam.alphan) + monkeypatch.setattr(physics_variables, "alphat", initneoclassicsparam.alphat) + monkeypatch.setattr(physics_variables, "ti0", initneoclassicsparam.ti0) + monkeypatch.setattr(physics_variables, "ni0", initneoclassicsparam.ni0) + monkeypatch.setattr(physics_variables, "fdeut", initneoclassicsparam.fdeut) + monkeypatch.setattr(physics_variables, "dnalp", initneoclassicsparam.dnalp) + monkeypatch.setattr(physics_variables, "rminor", initneoclassicsparam.rminor) + monkeypatch.setattr(physics_variables, "rmajor", initneoclassicsparam.rmajor) + monkeypatch.setattr(physics_variables, "bt", initneoclassicsparam.bt) + monkeypatch.setattr(physics_variables, "te", initneoclassicsparam.te) + monkeypatch.setattr(physics_variables, "ti", initneoclassicsparam.ti) + monkeypatch.setattr(physics_variables, "dene", initneoclassicsparam.dene) + monkeypatch.setattr(physics_variables, "deni", initneoclassicsparam.deni) + monkeypatch.setattr(neoclassics_module, "densities", initneoclassicsparam.densities) + monkeypatch.setattr( + neoclassics_module, "temperatures", initneoclassicsparam.temperatures + ) + monkeypatch.setattr( + neoclassics_module, "dr_densities", initneoclassicsparam.dr_densities + ) + monkeypatch.setattr( + neoclassics_module, "dr_temperatures", initneoclassicsparam.dr_temperatures + ) + monkeypatch.setattr(neoclassics_module, "roots", initneoclassicsparam.roots) + monkeypatch.setattr(neoclassics_module, "weights", initneoclassicsparam.weights) + monkeypatch.setattr(neoclassics_module, "nu", initneoclassicsparam.nu) + monkeypatch.setattr(neoclassics_module, "nu_star", initneoclassicsparam.nu_star) + monkeypatch.setattr( + neoclassics_module, "nu_star_averaged", initneoclassicsparam.nu_star_averaged + ) + monkeypatch.setattr(neoclassics_module, "vd", initneoclassicsparam.vd) + monkeypatch.setattr(neoclassics_module, "iota", initneoclassicsparam.iota) + monkeypatch.setattr(neoclassics_module, "q_flux", initneoclassicsparam.q_flux) + monkeypatch.setattr(neoclassics_module, "eps_eff", initneoclassicsparam.eps_eff) + monkeypatch.setattr(neoclassics_module, "r_eff", initneoclassicsparam.r_eff) + + neoclassics.init_neoclassics( + r_effin=initneoclassicsparam.r_effin, + eps_effin=initneoclassicsparam.eps_effin, + iotain=initneoclassicsparam.iotain, + ) + + assert neoclassics_module.densities == pytest.approx( + initneoclassicsparam.expected_densities, rel=0.001 + ) + assert neoclassics_module.temperatures == pytest.approx( + initneoclassicsparam.expected_temperatures, rel=0.001 + ) + assert neoclassics_module.dr_densities == pytest.approx( + initneoclassicsparam.expected_dr_densities, rel=0.001 + ) + assert neoclassics_module.dr_temperatures == pytest.approx( + initneoclassicsparam.expected_dr_temperatures, rel=0.001 + ) + assert neoclassics_module.roots == pytest.approx( + initneoclassicsparam.expected_roots, rel=0.001 + ) + assert neoclassics_module.weights == pytest.approx( + initneoclassicsparam.expected_weights, rel=0.001 + ) + assert neoclassics_module.nu == pytest.approx( + initneoclassicsparam.expected_nu, rel=0.001 + ) + assert neoclassics_module.nu_star == pytest.approx( + initneoclassicsparam.expected_nu_star, rel=0.001 + ) + assert neoclassics_module.nu_star_averaged == pytest.approx( + initneoclassicsparam.expected_nu_star_averaged, rel=0.001 + ) + assert neoclassics_module.vd == pytest.approx( + initneoclassicsparam.expected_vd, rel=0.001 + ) + assert neoclassics_module.q_flux == pytest.approx( + initneoclassicsparam.expected_q_flux, rel=0.001 + ) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 2c644e59..a6cf8c62 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -14,7 +14,7 @@ impurity_radiation_module, ) from process.power import Power -from process.stellarator import Stellarator +from process.stellarator import Stellarator, Neoclassics from process.vacuum import Vacuum from process.availability import Availability from process.buildings import Buildings @@ -44,6 +44,7 @@ def stellarator(): CCFE_HCPB(BlanketLibrary(Fw())), CurrentDrive(PlasmaProfile()), Physics(PlasmaProfile(), CurrentDrive(PlasmaProfile())), + Neoclassics(), )