From 25b872a9608b53aad5daef1c63f15d8c1d6d2050 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 14 Aug 2024 12:35:16 +0000 Subject: [PATCH] Convert physics functions to Python --- CMakeLists.txt | 1 - process/physics.py | 198 ++++- process/physics_functions.py | 541 +++++++++++-- process/stellarator.py | 42 +- source/fortran/init_module.f90 | 2 - source/fortran/physics_functions.f90 | 1088 -------------------------- tests/unit/test_physics.py | 14 + tests/unit/test_physics_functions.py | 268 ++----- 8 files changed, 782 insertions(+), 1372 deletions(-) delete mode 100644 source/fortran/physics_functions.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 908764827a..9988cdd44b 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,7 +85,6 @@ LIST(APPEND PROCESS_SRCS divertor_variables.f90 fwbs_variables.f90 physics.f90 - physics_functions.f90 physics_variables.f90 tfcoil_variables.f90 times_variables.f90 diff --git a/process/physics.py b/process/physics.py index cff586ad89..fb794f3541 100644 --- a/process/physics.py +++ b/process/physics.py @@ -14,7 +14,6 @@ reinke_module, impurity_radiation_module, constants, - physics_functions_module, physics_variables, physics_module, pulse_variables, @@ -1543,10 +1542,14 @@ def physics(self): self.plasma_profile.run() # Calculate total magnetic field [T] - physics_variables.btot = physics_functions_module.total_mag_field() + physics_variables.btot = np.sqrt( + physics_variables.bt**2 + physics_variables.bp**2 + ) # Calculate physics_variables.beta poloidal [-] - physics_variables.betap = physics_functions_module.beta_poloidal() + physics_variables.betap = beta_poloidal( + physics_variables.btot, physics_variables.bp, physics_variables.beta + ) # Set PF coil ramp times if pulse_variables.lpulse != 1: @@ -1799,7 +1802,6 @@ def physics(self): # Calculate fusion power + components - # physics_funcs.palph(self.plasma_profile) fusion_rate = physics_funcs.FusionReactionRate(self.plasma_profile) fusion_rate.calculate_fusion_rates() fusion_rate.set_physics_variables() @@ -1819,7 +1821,7 @@ def physics(self): physics_variables.betanb, physics_variables.dnbeam2, physics_variables.palpnb, - ) = physics_functions_module.beamfus( + ) = physics_funcs.beamfus( physics_variables.beamfus0, physics_variables.betbm0, physics_variables.bp, @@ -1858,30 +1860,32 @@ def physics(self): # Create some derived values and add beam contribution to fusion power ( + physics_variables.pneutpv, physics_variables.palpmw, physics_variables.pneutmw, physics_variables.pchargemw, physics_variables.betaft, + physics_variables.palppv, physics_variables.palpipv, physics_variables.palpepv, physics_variables.pfuscmw, physics_variables.powfmw, - ) = physics_functions_module.palph2( - physics_variables.bt, + ) = physics_funcs.palph2( physics_variables.bp, + physics_variables.bt, physics_variables.dene, physics_variables.deni, physics_variables.dnitot, physics_variables.falpe, physics_variables.falpi, physics_variables.palpnb, - physics_variables.ifalphap, physics_variables.pchargepv, physics_variables.pneutpv, physics_variables.ten, physics_variables.tin, physics_variables.vol, physics_variables.palppv, + physics_variables.ifalphap, ) # Nominal mean neutron wall load on entire first wall area including divertor and beam holes @@ -1954,20 +1958,20 @@ def physics(self): ) # Calculate L- to H-mode power threshold for different scalings - - physics_variables.pthrmw = physics_functions_module.pthresh( + physics_variables.pthrmw = pthresh( physics_variables.dene, physics_variables.dnla, physics_variables.bt, physics_variables.rmajor, + physics_variables.rminor, physics_variables.kappa, physics_variables.sarea, physics_variables.aion, physics_variables.aspect, + physics_variables.plasma_current, ) # Enforced L-H power threshold value (if constraint 15 is turned on) - physics_variables.plhthresh = physics_variables.pthrmw[ physics_variables.ilhthresh - 1 ] @@ -2005,7 +2009,9 @@ def physics(self): ) # Resistive diffusion time = current penetration time ~ mu0.a^2/resistivity - physics_variables.res_time = physics_functions_module.res_diff_time() + physics_variables.res_time = res_diff_time( + physics_variables.rmajor, physics_variables.rplas, physics_variables.kappa95 + ) # Power transported to the first wall by escaped alpha particles physics_variables.palpfwmw = physics_variables.palpmw * ( @@ -6607,3 +6613,171 @@ def pcond( taueff = (ratio + 1.0e0) / (ratio / tauei + 1.0e0 / tauee) return kappaa, ptrepv, ptripv, tauee, tauei, taueff, powerht + + +def beta_poloidal(btot, bp, beta): + """Calculates total poloidal beta + + Author: James Morris (UKAEA) + + J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) + Page 270 ISBN 0521851076 + + :param btot: sum of the toroidal and poloidal fields (T) + :param bp: poloidal field (T) + :param beta: total plasma beta + """ + return beta * (btot / bp) ** 2 + + +def res_diff_time(rmajor, rplas, kappa95): + """Calculates resistive diffusion time + + Author: James Morris (UKAEA) + + :param rmajor: plasma major radius (m) + :param rplas: plasma resistivity (Ohms) + :param kappa95: plasma elongation at 95% flux surface + """ + + return 2 * constants.rmu0 * rmajor / (rplas * kappa95) + + +def pthresh(dene, dnla, bt, rmajor, rminor, kappa, sarea, aion, aspect, plascur): + """L-mode to H-mode power threshold calculation + + Author: P J Knight, CCFE, Culham Science Centre + + - ITER Physics Design Description Document, p.2-2 + - ITER-FDR Plasma Performance Assessments, p.III-9 + - Snipes, 24th EPS Conference, Berchtesgaden 1997, p.961 + - Martin et al, 11th IAEA Tech. Meeting on H-mode Physics and + Transport Barriers, Journal of Physics: Conference Series + 123 (2008) 012033 + - J A Snipes and the International H-mode Threshold Database + Working Group, 2000, Plasma Phys. Control. Fusion, 42, A299 + + :param dene: volume-averaged electron density (/m3) + :param dnla: line-averaged electron density (/m3) + :param bt: toroidal field on axis (T) + :param rmajor: plasma major radius (m) + :param rminor: plasma minor radius (m) + :param kappa: plasma elongation + :param sarea: plasma surface area (m2) + :param aion: average mass of all ions (amu) + :param aspect: aspect ratio + :param plascur: plasma current (A) + + :returns: array of power thresholds (18 different scalings) + """ + + dene20 = 1e-20 * dene + dnla20 = 1e-20 * dnla + + # ITER-DDD, D.Boucher + # Fit to 1996 H-mode power threshold database: nominal + iterdd = 0.45 * dene20**0.75 * bt * rmajor**2 + + # Fit to 1996 H-mode power threshold database: upper bound + iterdd_ub = 0.37 * dene20 * bt * rmajor**2.5 + + # Fit to 1996 H-mode power threshold database: lower bound + iterdd_lb = 0.54 * dene20**0.5 * bt * rmajor**1.5 + + # J. A. Snipes, ITER H-mode Threshold Database Working Group, + # Controlled Fusion and Plasma Physics, 24th EPS Conference, + # Berchtesgaden, June 1997, vol.21A, part III, p.961 + snipes_1997 = 0.65 * dnla20**0.93 * bt**0.86 * rmajor**2.15 + + pthrmw5 = 0.42 * dnla20**0.80 * bt**0.90 * rmajor**1.99 * kappa**0.76 + + # Martin et al (2008) for recent ITER scaling, with mass correction + # and 95% confidence limits + martin = 0.0488 * dnla20**0.717 * bt**0.803 * sarea**0.941 * (2.0 / aion) + martin_error = ( + np.sqrt( + 0.057**2 + + (0.035 * np.log(dnla20)) ** 2 + + (0.032 * np.log(bt)) ** 2 + + (0.019 * np.log(sarea)) ** 2 + ) + * martin + ) + martin_ub = martin + 2 * martin_error + martin_lb = martin - 2 * martin_error + + # Snipes et al (2000) scaling with mass correction + # Nominal, upper and lower + snipes_2000 = ( + 1.42 * dnla20**0.58 * bt**0.82 * rmajor * rminor**0.81 * (2.0 / aion) + ) + snipes_2000_ub = ( + 1.547 + * dnla20**0.615 + * bt**0.851 + * rmajor**1.089 + * rminor**0.876 + * (2.0 / aion) + ) + snipes_2000_lb = ( + 1.293 + * dnla20**0.545 + * bt**0.789 + * rmajor**0.911 + * rminor**0.744 + * (2.0 / aion) + ) + + # Snipes et al (2000) scaling (closed divertor) with mass correction + # Nominal, upper and lower + + snipes_2000_cd = 0.8 * dnla20**0.5 * bt**0.53 * rmajor**1.51 * (2.0 / aion) + snipes_2000_cd_ub = ( + 0.867 * dnla20**0.561 * bt**0.588 * rmajor**1.587 * (2.0 / aion) + ) + snipes_2000_cd_lb = ( + 0.733 * dnla20**0.439 * bt**0.472 * rmajor**1.433 * (2.0 / aion) + ) + + # Hubbard et al. 2012 L-I threshold scaling + hubbard_2012 = 2.11 * (plascur / 1e6) ** 0.94 * dnla20**0.65 + hubbard_2012_lb = 2.11 * (plascur / 1e6) ** 0.70 * dnla20**0.47 + hubbard_2012_ub = 2.11 * (plascur / 1e6) ** 1.18 * dnla20**0.83 + + # Hubbard et al. 2017 L-I threshold scaling + hubbard_2017 = 0.162 * dnla20 * sarea * (bt) ** 0.26 + + pthrmw = [ + iterdd, + iterdd_ub, + iterdd_lb, + snipes_1997, + pthrmw5, + martin, + martin_ub, + martin_lb, + snipes_2000, + snipes_2000_ub, + snipes_2000_lb, + snipes_2000_cd, + snipes_2000_cd_ub, + snipes_2000_cd_lb, + hubbard_2012, + hubbard_2012_lb, + hubbard_2012_ub, + hubbard_2017, + ] + + # Aspect ratio corrected Martin et al (2008) + # Correction: Takizuka 2004, Plasma Phys. Control Fusion 46 A227 + if aspect <= 2.7: + takizuka_correction = 0.098 * aspect / (1.0 - (2.0 / (1.0 + aspect)) ** 0.5) + pthrmw += [ + martin * takizuka_correction, + martin_ub * takizuka_correction, + martin_lb * takizuka_correction, + ] + return pthrmw + + pthrmw += [martin, martin_ub, martin_lb] + return pthrmw diff --git a/process/physics_functions.py b/process/physics_functions.py index 806e3147c0..a832f7a9f3 100644 --- a/process/physics_functions.py +++ b/process/physics_functions.py @@ -1,15 +1,64 @@ import logging import numpy from scipy import integrate -from process.fortran import constants from dataclasses import dataclass -from process.fortran import physics_variables -from process.fortran import physics_module +from process.fortran import physics_variables, physics_module, constants import process.impurity_radiation as impurity logger = logging.getLogger(__name__) +ATOMIC_MASS_DEUTERIUM = 2.0 +ATOMIC_MASS_TRITIUM = 3.0 + +REACTION_CONSTANTS_DT = dict( + bg=34.3827, + mrc2=1.124656e6, + cc1=1.17302e-9, + cc2=1.51361e-2, + cc3=7.51886e-2, + cc4=4.60643e-3, + cc5=1.35000e-2, + cc6=-1.06750e-4, + cc7=1.36600e-5, +) + +REACTION_CONSTANTS_DHE3 = dict( + bg=68.7508, + mrc2=1.124572e6, + cc1=5.51036e-10, + cc2=6.41918e-3, + cc3=-2.02896e-3, + cc4=-1.91080e-5, + cc5=1.35776e-4, + cc6=0.0, + cc7=0.0, +) + +REACTION_CONSTANTS_DD1 = dict( + bg=31.3970, + mrc2=0.937814e6, + cc1=5.43360e-12, + cc2=5.85778e-3, + cc3=7.68222e-3, + cc4=0.0, + cc5=-2.96400e-6, + cc6=0.0, + cc7=0.0, +) + +REACTION_CONSTANTS_DD2 = dict( + bg=31.3970, + mrc2=0.937814e6, + cc1=5.65718e-12, + cc2=3.41267e-3, + cc3=1.99167e-3, + cc4=0.0, + cc5=1.05060e-5, + cc6=0.0, + cc7=0.0, +) + class FusionReactionRate: """Calculate the fusion reaction rate for each reaction case (DT,DHE3,DD1,DD2). @@ -37,17 +86,7 @@ def __init__(self, plasma_profile): def dt(self): """D + T --> 4He + n reaction""" - dt = BoschHaleConstants( - bg=34.3827, - mrc2=1.124656e6, - cc1=1.17302e-9, - cc2=1.51361e-2, - cc3=7.51886e-2, - cc4=4.60643e-3, - cc5=1.35000e-2, - cc6=-1.06750e-4, - cc7=1.36600e-5, - ) + dt = BoschHaleConstants(**REACTION_CONSTANTS_DT) sigmav = integrate.simpson( fint(self.plasma_profile, dt), x=self.plasma_profile.neprofile.profile_x, @@ -75,17 +114,7 @@ def dt(self): def dhe3(self): """D + 3He --> 4He + p reaction""" - dhe3 = BoschHaleConstants( - bg=68.7508, - mrc2=1.124572e6, - cc1=5.51036e-10, - cc2=6.41918e-3, - cc3=-2.02896e-3, - cc4=-1.91080e-5, - cc5=1.35776e-4, - cc6=0.0, - cc7=0.0, - ) + dhe3 = BoschHaleConstants(**REACTION_CONSTANTS_DHE3) sigmav = integrate.simpson( fint(self.plasma_profile, dhe3), x=self.plasma_profile.neprofile.profile_x, @@ -112,17 +141,7 @@ def dhe3(self): def dd1(self): """D + D --> 3He + n reaction""" - dd1 = BoschHaleConstants( - bg=31.3970, - mrc2=0.937814e6, - cc1=5.43360e-12, - cc2=5.85778e-3, - cc3=7.68222e-3, - cc4=0.0, - cc5=-2.96400e-6, - cc6=0.0, - cc7=0.0, - ) + dd1 = BoschHaleConstants(**REACTION_CONSTANTS_DD1) sigmav = integrate.simpson( fint(self.plasma_profile, dd1), x=self.plasma_profile.neprofile.profile_x, @@ -150,17 +169,7 @@ def dd1(self): def dd2(self): """D + D --> T + p reaction""" - dd2 = BoschHaleConstants( - bg=31.3970, - mrc2=0.937814e6, - cc1=5.65718e-12, - cc2=3.41267e-3, - cc3=1.99167e-3, - cc4=0.0, - cc5=1.05060e-5, - cc6=0.0, - cc7=0.0, - ) + dd2 = BoschHaleConstants(**REACTION_CONSTANTS_DD2) sigmav = integrate.simpson( fint(self.plasma_profile, dd2), x=self.plasma_profile.neprofile.profile_x, @@ -246,13 +255,13 @@ def radpwr(plasma_profile): pedgeradpv = imp_rad.radtot - imp_rad.radcore - # Synchrotron radiation power/volume; assumed to be from core only. + # Synchrotron radiation power/volume; assumed to be from core only. psyncpv = psync_albajar_fidone() - # Total core radiation power/volume. + # Total core radiation power/volume. pcoreradpv = imp_rad.radcore + psyncpv - # Total radiation power/volume. + # Total radiation power/volume. pradpv = imp_rad.radtot + psyncpv return RadpwrData(psyncpv, pcoreradpv, pedgeradpv, pradpv) @@ -273,8 +282,8 @@ def psync_albajar_fidone(): """ tbet = 2.0e0 - # rpow is the (1-Rsyn) power dependence based on plasma shape - # (see Fidone) + # rpow is the (1-Rsyn) power dependence based on plasma shape + # (see Fidone) rpow = 0.62e0 @@ -291,8 +300,8 @@ def psync_albajar_fidone(): 2.0e0 * numpy.pi**2 * physics_variables.rmajor * physics_variables.rminor**2 ) - # No account is taken of pedestal profiles here, other than use of - # the correct physics_variables.ne0 and physics_variables.te0... + # No account is taken of pedestal profiles here, other than use of + # the correct physics_variables.ne0 and physics_variables.te0... de2o = 1.0e-20 * physics_variables.ne0 pao = 6.04e3 * (physics_variables.rminor * de2o) / physics_variables.bt @@ -315,7 +324,7 @@ def psync_albajar_fidone(): * (1.0e0 - physics_variables.ssync) ** 0.41e0 ) - # Very high T modification, from Fidone + # Very high T modification, from Fidone dum = dum ** (-1.51e0) @@ -335,7 +344,7 @@ def psync_albajar_fidone(): * kfun ) - # psyncpv should be per unit volume; Albajar gives it as total + # psyncpv should be per unit volume; Albajar gives it as total psyncpv = psync / physics_variables.vol @@ -463,3 +472,423 @@ class RadpwrData: pcoreradpv: float pedgeradpv: float pradpv: float + + +def palph2( + bp, + bt, + dene, + deni, + dnitot, + falpe, + falpi, + palpnb, + pchargepv, + pneutpv, + ten, + tin, + vol, + palppv, + ifalphap, +): + """ + Fusion power and fast alpha pressure calculations. + ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + D J Ward, UKAEA Fusion: F/PL/PJK/PROCESS/CODE/050 + + :param bp: poloidal field (T) + :param bt: totoidal field on axis (T) + :param dene: electron density (m^-3) + :param deni: fuel ion density (m^-3) + :param dnitot: total ion density (m^-3) + :param falpe: fraction of alpha energy to electrons + :param falpi: fraction of alpha energy to ions + :param palpnb: alpha power from hot neutral beam ions (MW) + :param pchargepv: other charged particle fusion power/volume (MW/m3) + :param pneutpv: neutron fusion power per volume (MW/m3) + :param ten: density-weighted electron temperature (keV) + :param tin: density-weighted ion temperature (keV) + :param vol: plasma volume (m3) + :param palppv: alpha power per volume (MW/m3) + :param ifalphap: switch for fast alpha pressure method + + :return: neutron fusion power per volume (MW/m3), alpha power (MW), + neutron fusion power (MW), other charged particle fusion power (MW), + fast alpha beta component, alpha power per volume (MW/m3), + alpha power per volume to electrons (MW/m3), alpha power per volume to ions (MW/m3), + charged particle fusion power (MW), fusion power (MW) + """ + + # Add neutral beam alpha power / volume + palppv_out = palppv + palpnb / vol + + # Add extra neutron power + pneutpv_out = pneutpv + 4.0 * palpnb / vol + + # Total alpha power + palpmw = palppv_out * vol + + # Total non-alpha charged particle power + pchargemw = pchargepv * vol + + # Total neutron power + pneutmw = pneutpv_out * vol + + # Total fusion power + powfmw = palpmw + pneutmw + pchargemw + + # Charged particle fusion power + pfuscmw = palpmw + pchargemw + + # Alpha power to electrons and ions (used with electron + # and ion power balance equations only) + # No consideration of pchargepv here... + palpipv = physics_variables.falpha * palppv_out * falpi + palpepv = physics_variables.falpha * palppv_out * falpe + + # Determine average fast alpha density + if physics_variables.fdeut < 1.0: + + betath = ( + 2.0e3 + * constants.rmu0 + * constants.echarge + * (dene * ten + dnitot * tin) + / (bt**2 + bp**2) + ) + + # jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) + # IPDG89 fast alpha scaling + if ifalphap == 0: + fact = min(0.3, 0.29 * (deni / dene) ** 2 * ((ten + tin) / 20.0 - 0.37)) + + # Modified scaling, D J Ward + else: + fact = min( + 0.30, + 0.26 + * (deni / dene) ** 2 + * numpy.sqrt(max(0.0, ((ten + tin) / 20.0 - 0.65))), + ) + + fact = max(fact, 0.0) + fact2 = palppv_out / palppv + betaft = betath * fact * fact2 + + else: # negligible alpha production, palppv = palpnb = 0 + betaft = 0.0 + + return ( + pneutpv_out, + palpmw, + pneutmw, + pchargemw, + betaft, + palppv_out, + palpepv, + palpipv, + pfuscmw, + powfmw, + ) + + +def beamfus( + beamfus0, + betbm0, + bp, + bt, + cnbeam, + dene, + deni, + dlamie, + ealphadt, + enbeam, + fdeut, + ftrit, + ftritbm, + sigvdt, + ten, + tin, + vol, + zeffai, +): + """Routine to calculate beam slowing down properties + author: P J Knight, CCFE, Culham Science Centre + + :param beamfus0: multiplier for beam-background fusion calculation + :param betbm0: leading coefficient for neutral beam beta fraction + :param bp: poloidal field (T) + :param bt: toroidal field on axis (T) + :param cnbeam: neutral beam current (A) + :param dene: electron density (m^-3) + :param deni: fuel ion density (m^-3) + :param dlamie: ion-electron coulomb logarithm + :param ealphadt: alpha particle birth energy (D-T) (keV) + :param enbeam: neutral beam energy (keV) + :param fdeut: deuterium fraction of main plasma + :param ftrit: tritium fraction of main plasma + :param ftritbm: tritium fraction of neutral beam + :param sigvdt: profile averaged for D-T (m3/s) + :param ten: density-weighted electron temperature (keV) + :param tin: density-weighted ion temperature (keV) + :param vol: plasma volume (m3) + :param zeffai: mass weighted plasma effective charge + + :returns: neutral beam beta component, hot beam ion density (m^-3), + alpha power from hot neutral beam ions (MW) + """ + + tausl = ( + 1.99e19 + * (2.0 * (1.0 - ftritbm) + (3.0 * ftritbm)) + * (ten**1.5 / dene) + / dlamie + ) + + # Critical energy for electron/ion slowing down of the beam ion + # (deuterium and tritium neutral beams, respectively) (keV) + + ecritd = 14.8 * ten * 2.0 * zeffai**0.6666 * (dlamie + 4.0) / dlamie + ecritt = ecritd * 1.5 + + # Deuterium and tritium ion densities + + denid = deni * fdeut + denit = deni * ftrit + + palpdb, palptb, dnbeam2, ehotnb = beamcalc( + denid, + denit, + ealphadt, + enbeam, + ecritd, + ecritt, + tausl, + ftritbm, + cnbeam, + tin, + vol, + sigvdt, + ) + + # Neutral beam alpha power + + palpnb = beamfus0 * (palpdb + palptb) + + # Neutral beam beta + + betanb = betbm0 * 4.03e-22 * 0.66666 * dnbeam2 * ehotnb / (bt**2 + bp**2) + + return betanb, dnbeam2, palpnb + + +def beamcalc( + nd, nt, ealphadt, ebeam, ecritd, ecritt, tausbme, ftritbm, ibeam, ti, vol, svdt +): + """Neutral beam alpha power and ion energy + author: P J Knight, CCFE, Culham Science Centre + + :param nd: thermal deuterium density (/m3) + :param nt: thermal tritium density (/m3) + :param ealphadt: alpha particle birth energy (D-T) (keV) + :param ebeam: beam energy (keV) + :param ecritd: critical energy for electron/ion slowing down of + the beam ion (deuterium neutral beam) (keV) + :param ecritt: critical energy for beam slowing down (tritium neutral beam) (keV) + :param ftritbm: beam tritium fraction (0.0 = deuterium beam) + :param ibeam: beam current (A) + :param svdt: profile averaged for D-T (m3/s) + :param tausbme: beam ion slowing down time on electrons (s) + :param ti: thermal ion temperature (keV) + :param vol: plasma volume (m3) (95% flux surface) + + :returns: alpha power from deut. beam-background fusion (MW), + alpha power from trit. beam-background fusion (MW), hot beam ion density (m^-3), + average hot beam ion energy (keV) + """ + + # D and T beam current fractions + ifbmd = ibeam * (1.0 - ftritbm) + ifbmt = ibeam * ftritbm + + ebmratd = ebeam / ecritd + vcritd = numpy.sqrt( + 2.0 + * constants.echarge + * 1000.0 + * ecritd + / (constants.mproton * ATOMIC_MASS_DEUTERIUM) + ) + tauseffd = tausbme / 3.0 * numpy.log(1.0 + (ebmratd) ** 1.5) + nhotmsd = (1.0 - ftritbm) * ibeam * tauseffd / (constants.echarge * vol) + + ebmratt = ebeam / ecritt + vcritt = numpy.sqrt( + 2.0 + * constants.echarge + * 1000.0 + * ecritt + / (constants.mproton * ATOMIC_MASS_TRITIUM) + ) + tausefft = tausbme / 3.0 * numpy.log(1.0 + (ebmratt) ** 1.5) + nhotmst = ftritbm * ibeam * tausefft / (constants.echarge * vol) + + nhot = nhotmsd + nhotmst + ndhot = nhotmsd + nthot = nhotmst + + # Average hot ion energy from Deng & Emmert, UWFDM-718, Jan 87 + vcds = 2.0 * ecritd * constants.echarge * 1000.0 / (2.0 * constants.mproton) + vcts = 2.0 * ecritt * constants.echarge * 1000.0 / (3.0 * constants.mproton) + + s0d = ifbmd / (constants.echarge * vol) + s0t = ifbmt / (constants.echarge * vol) + + xcoefd = ( + ATOMIC_MASS_DEUTERIUM + * constants.mproton + * tausbme + * vcds + * s0d + / (constants.echarge * 1000.0 * 3.0) + ) + xcoeft = ( + ATOMIC_MASS_TRITIUM + * constants.mproton + * tausbme + * vcts + * s0t + / (constants.echarge * 1000.0 * 3.0) + ) + + presd = xcoefd * xbrak(ebeam, ecritd) + prest = xcoeft * xbrak(ebeam, ecritt) + + ehotd = 1.5 * presd / ndhot + ehott = 1.5 * prest / nthot + ehot = (ndhot * ehotd + nthot * ehott) / nhot + + iabm = 2 + svdhotn = 1e-4 * sgvhot(iabm, vcritd, ebeam) + iabm = 3 + svthotn = 1e-4 * sgvhot(iabm, vcritt, ebeam) + + palfdb = palphabm(ealphadt, ndhot, nt, svdhotn, vol, ti, svdt) + palftb = palphabm(ealphadt, nthot, nd, svthotn, vol, ti, svdt) + + return palfdb, palftb, nhot, ehot + + +def xbrak(e0, ec): + """Hot ion energy parameter + author: P J Knight, CCFE, Culham Science Centre + + :param e0: neutral beam energy (keV) + :param ec: critical energy for electron/ion slowing down of the beam ion (keV) + """ + xcs = e0 / ec + xc = numpy.sqrt(xcs) + + t1 = xcs / 2.0 + t2 = (numpy.log((xcs + 2.0 * xc + 1.0) / (xcs - xc + 1.0))) / 6.0 + + xarg = (2.0 * xc - 1.0) / numpy.sqrt(3.0) + t3 = (numpy.arctan(xarg)) / numpy.sqrt(3.0) + t4 = 0.3022999 + + return t1 + t2 - t3 - t4 + + +def palphabm(ealphadt, nbm, nblk, sigv, vol, ti, svdt): + """Alpha power from beam-background fusion + author: P J Knight, CCFE, Culham Science Centre + + :param ealphadt: alpha particle birth energy (D-T) (keV) + :param nblk: thermal ion density (/m3) + :param nbm: hot beam ion density (/m3) + :param sigv: hot beam fusion reaction rate (m3/s) + :param vol: plasma volume (m3) + :param ti: thermal ion temperature (keV) + :param svdt: profile averaged for D-T (m3/s) + """ + + # [ti] because bosch_hale expects temperature profile + # so we pass it a profile of length 1 + ratio = svdt / bosch_hale( + numpy.array([ti]), BoschHaleConstants(**REACTION_CONSTANTS_DT) + ) + return ( + constants.echarge / 1000.0 * nbm * nblk * sigv * ealphadt * vol * ratio.item() + ) + + +def sgvhot(rmass_ion, vcrx, ebeam): + """Hot beam fusion reaction rate + author: P J Knight, CCFE, Culham Science Centre + + :param rmass_ion: relative atomic mass of the ion (of D or T) + :param vcrx: critical velocity for electron/ion slowing down of the beam ion (m/s) + :param ebeam: neutral beam energy (keV) + """ + # Beam velocity + + vbeams = ebeam * constants.echarge * 1000.0 * 2.0 / (rmass_ion * constants.mproton) + vbeam = numpy.sqrt(vbeams) + + xv = vbeam / vcrx + t1 = 3.0 * vcrx / numpy.log(1.0 + (xv**3)) + + svint = integrate.quad( + _hot_beam_fusion_reaction_rate_integrand, 0.0, xv, args=(vcrx,) + )[0] + + return t1 * svint + + +def _hot_beam_fusion_reaction_rate_integrand(u, vcritx): + """Integrand function for the hot beam fusion reaction rate + author: P J Knight, CCFE, Culham Science Centre + + :param u: ratio of beam velocity to the critical velocity + """ + t1 = (u**3) / (1.0 + u**3) + + # vcritx : critical velocity for electron/ion slowing down of beam ion (m/s) + xvc = vcritx * u + xvcs = xvc * xvc * constants.mproton / (constants.echarge * 1000.0) + t2 = _sigbmfus(xvcs) + + return t1 * t2 + + +def _sigbmfus(vrelsq): + """Fusion reaction cross-section + author: P J Knight, CCFE, Culham Science Centre + + The functional form of the cross-section is in terms of the equivalent + deuterium energy, i.e. for a tritium beam at 500 keV the energy + used in the cross-section function is 333 keV. + + :param vrelsq: square of the speed of the beam ion (keV/amu) + """ + a1 = 45.95 + a2 = 5.02e4 + a3 = 1.368e-2 + a4 = 1.076 + a5 = 4.09e2 + + # Beam kinetic energy + + ebm = 0.5 * ATOMIC_MASS_DEUTERIUM * vrelsq + + # Set limits on cross-section at low and high beam energies + + if ebm < 10.0: + return 1.0e-27 + elif ebm > 1.0e4: + return 8.0e-26 + else: + t1 = a2 / (1.0e0 + (a3 * ebm - a4) ** 2) + a5 + t2 = ebm * (numpy.exp(a1 / numpy.sqrt(ebm)) - 1.0) + return 1.0e-24 * t1 / t2 diff --git a/process/stellarator.py b/process/stellarator.py index 62ff4a483b..34c233ab48 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -24,7 +24,6 @@ constraint_variables, rebco_variables, maths_library, - physics_functions_module, neoclassics_module, impurity_radiation_module, sctfcoil_module, @@ -3934,30 +3933,13 @@ def stphys(self, output): # Calculate fusion power - ( - physics_variables.palppv, - physics_variables.pchargepv, - physics_variables.pneutpv, - sigvdt, - physics_variables.fusionrate, - physics_variables.alpharate, - physics_variables.protonrate, - pdtpv, - pdhe3pv, - pddpv, - ) = physics_functions_module.palph( - physics_variables.alphan, - physics_variables.alphat, - physics_variables.deni, - physics_variables.fdeut, - physics_variables.fhe3, - physics_variables.ftrit, - physics_variables.ti, - ) + fusion_rate = physics_funcs.FusionReactionRate(self.plasma_profile) + fusion_rate.calculate_fusion_rates() + fusion_rate.set_physics_variables() - physics_variables.pdt = pdtpv * physics_variables.vol - physics_variables.pdhe3 = pdhe3pv * physics_variables.vol - physics_variables.pdd = pddpv * physics_variables.vol + physics_variables.pdt = physics_module.pdtpv * physics_variables.vol + physics_variables.pdhe3 = physics_module.pdhe3pv * physics_variables.vol + physics_variables.pdd = physics_module.pddpv * physics_variables.vol # Calculate neutral beam slowing down effects # If ignited, then ignore beam fusion effects @@ -3969,7 +3951,7 @@ def stphys(self, output): physics_variables.betanb, physics_variables.dnbeam2, physics_variables.palpnb, - ) = physics_functions_module.beamfus( + ) = physics_funcs.beamfus( physics_variables.beamfus0, physics_variables.betbm0, physics_variables.bp, @@ -3983,7 +3965,7 @@ def stphys(self, output): physics_variables.fdeut, physics_variables.ftrit, current_drive_variables.ftritbm, - sigvdt, + physics_module.sigvdt, physics_variables.ten, physics_variables.tin, physics_variables.vol, @@ -4007,30 +3989,32 @@ def stphys(self, output): physics_variables.pdt = physics_variables.pdt + 5.0e0 * physics_variables.palpnb ( + physics_variables.pneutpv, physics_variables.palpmw, physics_variables.pneutmw, physics_variables.pchargemw, physics_variables.betaft, + physics_variables.palppv, physics_variables.palpipv, physics_variables.palpepv, physics_variables.pfuscmw, physics_variables.powfmw, - ) = physics_functions_module.palph2( - physics_variables.bt, + ) = physics_funcs.palph2( physics_variables.bp, + physics_variables.bt, physics_variables.dene, physics_variables.deni, physics_variables.dnitot, physics_variables.falpe, physics_variables.falpi, physics_variables.palpnb, - physics_variables.ifalphap, physics_variables.pchargepv, physics_variables.pneutpv, physics_variables.ten, physics_variables.tin, physics_variables.vol, physics_variables.palppv, + physics_variables.ifalphap, ) # Neutron wall load diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90 index 18d80364be..b18ca215d8 100644 --- a/source/fortran/init_module.f90 +++ b/source/fortran/init_module.f90 @@ -49,7 +49,6 @@ subroutine init_all_module_vars use rebco_variables, only: init_rebco_variables use reinke_variables, only: init_reinke_variables use define_iteration_variables, only: init_define_iteration_variables - use physics_functions_module, only: init_physics_functions use reinke_module, only: init_reinke_module use water_usage_variables, only: init_watuse_variables use CS_fatigue_variables, only: init_CS_fatigue_variables @@ -91,7 +90,6 @@ subroutine init_all_module_vars call init_rebco_variables call init_reinke_variables call init_define_iteration_variables - call init_physics_functions call init_reinke_module call init_watuse_variables call init_CS_fatigue_variables diff --git a/source/fortran/physics_functions.f90 b/source/fortran/physics_functions.f90 deleted file mode 100644 index f0b498e693..0000000000 --- a/source/fortran/physics_functions.f90 +++ /dev/null @@ -1,1088 +0,0 @@ -module physics_functions_module - - !! Module containing physics subfunctions - !! author: K Ellis, CCFE, Culham Science Centre - !! N/A - !! This module contains physics routines which can be called by physics or - !! other modules - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - implicit none - - public :: beamfus, palph, palph2 - - ! Module-level variables - real(dp) :: vcritx - -contains - - subroutine init_physics_functions - !! Initialise module variables - implicit none - - vcritx = 0.0D0 - end subroutine init_physics_functions - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine pthresh(dene,dnla,bt,rmajor,kappa,sarea,aion,aspect,pthrmw) - !! author: P J Knight, CCFE, Culham Science Centre - !! L-mode to H-mode power threshold calculation - !! dene : input real : volume-averaged electron density (/m3) - !! dnla : input real : line-averaged electron density (/m3) - !! bt : input real : toroidal field on axis (T) - !! rmajor : input real : plasma major radius (m) - !! kappa : input real : plasma elongation - !! sarea : input real : plasma surface area (m**2) - !! aion : input real : average mass of all ions (amu) - !! aspect : input real : aspect ratio - !! pthrmw(17) : output real array : power threshold (different scalings) - !! This routine calculates the power threshold for the L-mode to - !! H-mode transition. - !! ITER Physics Design Description Document, p.2-2 - !! ITER-FDR Plasma Performance Assessments, p.III-9 - !! Snipes, 24th EPS Conference, Berchtesgaden 1997, p.961 - !! Martin et al, 11th IAEA Tech. Meeting on H-mode Physics and - !! Transport Barriers, Journal of Physics: Conference Series - !! 123 (2008) 012033 - !! J A Snipes and the International H-mode Threshold Database - !! Working Group, 2000, Plasma Phys. Control. Fusion, 42, A299 - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use physics_variables, only: rminor, plasma_current - implicit none - - ! Arguments - - real(dp), intent(in) :: dene,dnla,bt,rmajor,kappa,sarea,aion,aspect - real(dp), dimension(21), intent(out) :: pthrmw - - ! Local variables - - real(dp) :: dene20,dnla20,marterr - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - dene20 = 1.0D-20*dene - dnla20 = 1.0D-20*dnla - - ! ITER-DDD, D.Boucher - ! Fit to 1996 H-mode power threshold database: nominal - - pthrmw(1) = 0.45D0 * dene20**0.75D0 * bt * rmajor**2 - - ! Fit to 1996 H-mode power threshold database: upper bound - - pthrmw(2) = 0.37D0 * dene20 * bt * rmajor**2.5D0 - - ! Fit to 1996 H-mode power threshold database: lower bound - - pthrmw(3) = 0.54D0 * dene20**0.5D0 * bt * rmajor**1.5D0 - - ! J. A. Snipes, ITER H-mode Threshold Database Working Group, - ! Controlled Fusion and Plasma Physics, 24th EPS Conference, - ! Berchtesgaden, June 1997, vol.21A, part III, p.961 - - pthrmw(4) = 0.65D0 * dnla20**0.93D0 * bt**0.86D0 * rmajor**2.15D0 - - pthrmw(5) = 0.42D0 * dnla20**0.80D0 * bt**0.90D0 * rmajor**1.99D0 & - * kappa**0.76D0 - - ! Martin et al (2008) for recent ITER scaling, with mass correction - ! and 95% confidence limits - - pthrmw(6) = 0.0488D0 * dnla20**0.717D0 * bt**0.803D0 & - * sarea**0.941D0 * (2.0D0/aion) - - marterr = 0.057D0**2 + (0.035D0 * log(dnla20))**2 & - + (0.032D0 * log(bt))**2 + (0.019D0 * log(sarea))**2 - marterr = sqrt(marterr) * pthrmw(6) - - pthrmw(7) = pthrmw(6) + 2.0D0*marterr - pthrmw(8) = pthrmw(6) - 2.0D0*marterr - - ! Snipes et al (2000) scaling with mass correction - ! Nominal, upper and lower - - pthrmw(9) = 1.42D0 * dnla20**0.58D0 * bt**0.82D0 * rmajor & - * rminor**0.81D0 * (2.0D0/aion) - - pthrmw(10) = 1.547D0 * dnla20**0.615D0 * bt**0.851D0 & - * rmajor**1.089D0 * rminor**0.876D0 * (2.0D0/aion) - - pthrmw(11) = 1.293D0 * dnla20**0.545D0 * bt**0.789D0 & - * rmajor**0.911D0 * rminor**0.744D0 * (2.0D0/aion) - - ! Snipes et al (2000) scaling (closed divertor) with mass correction - ! Nominal, upper and lower - - pthrmw(12) = 0.8D0 * dnla20**0.5D0 * bt**0.53D0 * rmajor**1.51D0 & - * (2.0D0/aion) - - pthrmw(13) = 0.867D0 * dnla20**0.561D0 * bt**0.588D0 * rmajor**1.587D0 & - * (2.0D0/aion) - - pthrmw(14) = 0.733D0 * dnla20**0.439D0 * bt**0.472D0 * rmajor**1.433D0 & - * (2.0D0/aion) - - ! Hubbard et al. 2012 L-I threshold scaling - - ! Nominal - pthrmw(15) = 2.11 * (plasma_current/1.0D6)**0.94 * dnla20**0.65 - - ! Lower bound - pthrmw(16) = 2.11 * (plasma_current/1.0D6)**0.70 * dnla20**0.47 - - ! Upper bound - pthrmw(17) = 2.11 * (plasma_current/1.0D6)**1.18 * dnla20**0.83 - - ! Hubbard et al. 2017 L-I threshold scaling - pthrmw(18) = 0.162 * dnla20 * sarea * (bt)**0.26 - - ! Aspect ratio corrected Martin et al (2008) - ! Correction: Takizuka 2004, Plasma Phys. Control Fusion 46 A227 - if (aspect.le.2.7D0) then - pthrmw(19) = pthrmw(6) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0)) - pthrmw(20) = pthrmw(7) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0)) - pthrmw(21) = pthrmw(8) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0)) - else - pthrmw(19) = pthrmw(6) - pthrmw(20) = pthrmw(7) - pthrmw(21) = pthrmw(8) - end if - - end subroutine pthresh - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine palph(alphan,alphat,deni,fdeut,fhe3,ftrit,ti, & - palppv,pchargepv,pneutpv,sigvdt,fusionrate,alpharate,protonrate, & - pdtpv,pdhe3pv,pddpv) - - !! (Initial part of) fusion power and fast alpha pressure calculations - !! author: P J Knight, CCFE, Culham Science Centre - !! alphan : input real : density profile index - !! alphat : input real : temperature profile index - !! deni : input real : fuel ion density (/m3) - !! fdeut : input real : deuterium fuel fraction - !! fhe3 : input real : helium-3 fuel fraction - !! ftrit : input real : tritium fuel fraction - !! ti : input real : ion temperature (keV) - !! palppv : output real : alpha particle fusion power per volume (MW/m3) - !! pchargepv : output real : other charged particle fusion power/volume (MW/m3) - !! pneutpv : output real : neutron fusion power per volume (MW/m3) - !! sigvdt : output real : profile averaged (m3/s) - !! fusionrate : output real : fusion reaction rate (reactions/m3/s) - !! alpharate : output real : alpha particle production rate (/m3/s) - !! protonrate : output real : proton production rate (/m3/s) - !! pdtpv : output real : D-T fusion power (MW/m3) - !! pdhe3pv : output real : D-He3 fusion power (MW/m3) - !! pddpv : output real : D-D fusion power (MW/m3) - !! This subroutine numerically integrates over plasma cross-section to - !! find the core plasma fusion power. - !! T&M/PKNIGHT/LOGBOOK24, p.6 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: echarge - use maths_library, only: quanc8 - implicit none - - ! Arguments - - real(dp), intent(in) :: alphan, alphat, deni, fdeut, & - fhe3, ftrit, ti - real(dp), intent(out) :: palppv, pchargepv, pneutpv, sigvdt, & - fusionrate, alpharate, protonrate, pdtpv, pdhe3pv, pddpv - - ! Local variables - - integer, parameter :: DT=1, DHE3=2, DD1=3, DD2=4 - integer :: ireaction,nofun - real(dp) :: alow,arate,bhigh,epsq8,errest,etot,flag, & - fpow,frate,pa,pc,pn,prate,sigmav - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Initialise local quantities - alow = 0.0D0 - bhigh = 1.0D0 - epsq8 = 1.0D-9 - - ! Find fusion power - ! Integrate over plasma profiles to obtain fusion reaction rate - palppv = 0.0D0 - pchargepv = 0.0D0 - pneutpv = 0.0D0 - fusionrate = 0.0D0 - alpharate = 0.0D0 - protonrate = 0.0D0 - pddpv = 0.0D0 - - do ireaction = 1,4 - ! Fusion reaction rate (m3/s) is calculated in fint for each ireaction - ! sigmav is the volume-averaged fusion reaction rate (m3/s) - ! = integral(2 rho sigv(rho).ni(rho)^2 drho) / (deni**2) - - call quanc8(fint,alow,bhigh,epsq8,epsq8,sigmav,errest,nofun,flag) - if (ireaction == DT) sigvdt = sigmav - - select case (ireaction) - - case (DT) ! D + T --> 4He + n reaction - - etot = 17.59D0 * echarge ! MJ - fpow = 1.0D0 * sigmav * etot * fdeut*ftrit * deni*deni ! MW/m3 - pa = 0.2D0 * fpow - pc = 0.0D0 - pn = 0.8D0 * fpow - frate = fpow/etot ! reactions/m3/second - arate = frate - prate = 0.0D0 - pdtpv = fpow - - case (DHE3) ! D + 3He --> 4He + p reaction - - etot = 18.35D0 * echarge ! MJ - fpow = 1.0D0 * sigmav * etot * fdeut*fhe3 * deni*deni ! MW/m3 - pa = 0.2D0 * fpow - pc = 0.8D0 * fpow - pn = 0.0D0 - frate = fpow/etot ! reactions/m3/second - arate = frate - prate = frate ! proton production /m3/second - pdhe3pv = fpow - - case (DD1) ! D + D --> 3He + n reaction - - ! The 0.5 branching ratio is assumed to be included in sigmav - etot = 3.27D0 * echarge ! MJ - fpow = 1.0D0 * sigmav * etot * 0.5D0*fdeut*fdeut * deni*deni ! MW/m3 - pa = 0.0D0 - pc = 0.25D0 * fpow - pn = 0.75D0 * fpow - frate = fpow/etot ! reactions/m3/second - arate = 0.0D0 - prate = 0.0D0 ! Issue #557: No proton production - pddpv = pddpv + fpow - - case (DD2) ! D + D --> T + p reaction - - ! The 0.5 branching ratio is assumed to be included in sigmav - etot = 4.03D0 * echarge ! MJ - fpow = 1.0D0 * sigmav * etot * 0.5D0*fdeut*fdeut * deni*deni ! MW/m3 - pa = 0.0D0 - pc = fpow - pn = 0.0D0 - frate = fpow/etot ! reactions/m3/second - arate = 0.0D0 - prate = frate ! proton production /m3/second - pddpv = pddpv + fpow - - end select - - palppv = palppv + pa - pchargepv = pchargepv + pc - pneutpv = pneutpv + pn - fusionrate = fusionrate + frate - alpharate = alpharate + arate - protonrate = protonrate + prate - - end do - - contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function fint(rho) - - !! Integrand for fusion power integration - !! author: P J Knight, CCFE, Culham Science Centre - !! rho : input real : Abscissa of the integration, = normalised - !! plasma minor radius (0.0 <= rho < 1.0) - !! This function evaluates the integrand for the fusion power - !! integration, performed using routine - !! QUANC8 - !! in routine PALPH. - !! The fusion reaction assumed is controlled by flag - !! ireaction set in PALPH. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use physics_variables, only: te, rhopedt, te0, teped, tesep, tbeta, & - dene, rhopedn, ne0, neped, nesep - use profiles_module, only: tprofile, nprofile - implicit none - - real(dp) :: fint - - ! Arguments - - real(dp), intent(in) :: rho - - ! Local variables - - real(dp) :: nprof, nprofsq, sigv, tiofr - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Local ion temperature (keV) at r/a = rho - - tiofr = ti/te * tprofile(rho,rhopedt,te0,teped,tesep,alphat,tbeta) - - ! Fusion reaction rate (m3/s) - - sigv = bosch_hale(tiofr,ireaction) - - ! Integrand for the volume averaged fusion reaction rate sigmav: - ! sigmav = integral(2 rho (sigv(rho) ni(rho)^2) drho), - ! divided by the square of the volume-averaged ion density - ! to retain the dimensions m3/s (this is multiplied back in later) - - nprof = 1.0D0/dene * nprofile(rho,rhopedn,ne0,neped,nesep,alphan) - nprofsq = nprof*nprof - - fint = 2.0D0 * rho * sigv * nprofsq - - end function fint - - end subroutine palph - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine palph2(bt,bp,dene,deni,dnitot,falpe,falpi,palpnb, & - ifalphap,pchargepv,pneutpv,ten,tin,vol,palpmw,pneutmw,pchargemw, & - betaft,palppv,palpipv,palpepv,pfuscmw,powfmw) - - !! (Concluding part of) fusion power and fast alpha pressure - !! calculations - !! author: P J Knight, CCFE, Culham Science Centre - !! bp : input real : poloidal field (T) - !! bt : input real : toroidal field on axis (T) - !! dene : input real : electron density (/m3) - !! deni : input real : fuel ion density (/m3) - !! dnitot : input real : total ion density (/m3) - !! falpe : input real : fraction of alpha energy to electrons - !! falpi : input real : fraction of alpha energy to ions - !! ifalphap : input integer : switch for fast alpha pressure method - !! palpnb : input real : alpha power from hot neutral beam ions (MW) - !! pchargepv : input real : other charged particle fusion power/volume (MW/m3) - !! pneutpv : input/output real : neutron fusion power per volume (MW/m3) - !! ten : input real : density-weighted electron temperature (keV) - !! tin : input real : density-weighted ion temperature (keV) - !! vol : input real : plasma volume (m3) - !! palpmw : output real : alpha power (MW) - !! pneutmw : output real : neutron fusion power (MW) - !! pchargemw : output real : other charged particle fusion power (MW) - !! betaft : output real : fast alpha beta component - !! palppv : input/output real : alpha power per volume (MW/m3) - !! palpepv : output real : alpha power per volume to electrons (MW/m3) - !! palpipv : output real : alpha power per volume to ions (MW/m3) - !! pfuscmw : output real : charged particle fusion power (MW) - !! powfmw : output real : fusion power (MW) - !! This subroutine completes the calculation of the fusion power - !! fast alpha pressure, and determines other alpha particle quantities. - !! ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - !! ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - !! D J Ward, UKAEA Fusion: F/PL/PJK/PROCESS/CODE/050 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: echarge, rmu0 - use physics_variables, only: falpha, fdeut - implicit none - - ! Arguments - - integer, intent(in) :: ifalphap - real(dp), intent(in) :: bp, bt, dene, deni, dnitot, falpe, & - falpi, palpnb, pchargepv, ten, tin, vol - real(dp), intent(inout) :: palppv, pneutpv - real(dp), intent(out) :: palpmw, pneutmw, pchargemw, betaft, palpepv, & - palpipv, pfuscmw, powfmw - - ! Local variables - - real(dp) :: betath, fact, fact2, palppv_no_nb - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Store the palppv value without the NB alpha power - palppv_no_nb = palppv - - ! Add neutral beam alpha power / volume - palppv = palppv + palpnb/vol - - ! Add extra neutron power - pneutpv = pneutpv + 4.0D0*palpnb/vol - - ! Total alpha power - palpmw = palppv*vol - - ! Total non-alpha charged particle power - pchargemw = pchargepv*vol - - ! Total neutron power - pneutmw = pneutpv*vol - - ! Total fusion power - powfmw = palpmw + pneutmw + pchargemw - - ! Charged particle fusion power - pfuscmw = palpmw + pchargemw - - ! Alpha power to electrons and ions (used with electron - ! and ion power balance equations only) - ! No consideration of pchargepv here... - palpipv = falpha * palppv*falpi - palpepv = falpha * palppv*falpe - - ! Determine average fast alpha density - if (fdeut < 1.0D0) then - - betath = 2.0D3*rmu0*echarge * (dene*ten + dnitot*tin)/(bt**2 + bp**2) - - ! jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) - ! IPDG89 fast alpha scaling - if (ifalphap == 0) then - fact = min( 0.30D0, & - 0.29D0*(deni/dene)**2 * ( (ten+tin)/20.0D0 - 0.37D0) ) - - ! Modified scaling, D J Ward - else - fact = min( 0.30D0, & - 0.26D0*(deni/dene)**2 * & - sqrt( max(0.0D0, ((ten+tin)/20.0D0 - 0.65D0)) ) ) - end if - - fact = max(fact,0.0D0) - fact2 = palppv / palppv_no_nb - betaft = betath * fact*fact2 - - else ! negligible alpha production, palppv = palpnb = 0 - betaft = 0.0D0 - end if - - end subroutine palph2 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function bosch_hale(t,reaction) - - !! Routine to calculate the fusion reaction rate - !! author: R Kemp, CCFE, Culham Science Centre - !! author: P J Knight, CCFE, Culham Science Centre - !! t : input real : Maxwellian density-weighted ion temperature (keV) - !! reaction : input integer : flag for fusion reaction to use: - !! 1 : D-T reaction - !! 2 : D-3He reaction - !! 3 : D-D 1st reaction (50% probability) - !! 4 : D-D 2nd reaction (50% probability) - !! This routine calculates the volumetric fusion reaction rate - !! <sigma v> in m3/s for one of four nuclear reactions, - !! using the Bosch-Hale parametrization. - !!

The valid range of the fit is 0.2 keV < t < 100 keV - !! Bosch and Hale, Nuclear Fusion 32 (1992) 611-631 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - real(dp) :: bosch_hale - - ! Arguments - - real(dp), intent(in) :: t - integer, intent(in) :: reaction - - ! Local variables - - integer, parameter :: DT=1, DHE3=2, DD1=3, DD2=4 - real(dp) :: theta1, theta, xi - real(dp), dimension(4) :: bg, mrc2 - real(dp), dimension(4,7) :: cc - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if (t == 0.0D0) then - bosch_hale = 0.0D0 - return - end if - - ! Gamov constant, BG - - bg(DT) = 34.3827D0 ! D + T --> 4He + n reaction - bg(DHE3) = 68.7508D0 ! D + 3He --> 4He + p reaction - bg(DD1) = 31.3970D0 ! D + D --> 3He + n reaction - bg(DD2) = 31.3970D0 ! D + D --> T + p reaction - - ! Reduced mass of the particles, keV - - mrc2(DT) = 1.124656D6 - mrc2(DHE3) = 1.124572D6 - mrc2(DD1) = 0.937814D6 - mrc2(DD2) = 0.937814D6 - - ! Parametrization coefficients - - cc(DT,1) = 1.17302D-9 - cc(DT,2) = 1.51361D-2 - cc(DT,3) = 7.51886D-2 - cc(DT,4) = 4.60643D-3 - cc(DT,5) = 1.35000D-2 - cc(DT,6) = -1.06750D-4 - cc(DT,7) = 1.36600D-5 - - cc(DHE3,1) = 5.51036D-10 - cc(DHE3,2) = 6.41918D-3 - cc(DHE3,3) = -2.02896D-3 - cc(DHE3,4) = -1.91080D-5 - cc(DHE3,5) = 1.35776D-4 - cc(DHE3,6) = 0.00000D0 - cc(DHE3,7) = 0.00000D0 - - cc(DD1,1) = 5.43360D-12 - cc(DD1,2) = 5.85778D-3 - cc(DD1,3) = 7.68222D-3 - cc(DD1,4) = 0.00000D0 - cc(DD1,5) = -2.96400D-6 - cc(DD1,6) = 0.00000D0 - cc(DD1,7) = 0.00000D0 - - cc(DD2,1) = 5.65718D-12 - cc(DD2,2) = 3.41267D-3 - cc(DD2,3) = 1.99167D-3 - cc(DD2,4) = 0.00000D0 - cc(DD2,5) = 1.05060D-5 - cc(DD2,6) = 0.00000D0 - cc(DD2,7) = 0.00000D0 - - theta1 = t*(cc(reaction,2) + t*(cc(reaction,4) + t*cc(reaction,6))) / & - (1.0D0 + t*(cc(reaction,3) + t*(cc(reaction,5) + t*cc(reaction,7)))) - theta = t/(1.0D0 - theta1) - - xi = ((bg(reaction)**2)/(4.0D0*theta))**0.3333333333D0 - - ! Volumetric reaction rate (m3/s) - - bosch_hale = 1.0D-6 * cc(reaction,1) * theta * & - sqrt( xi/(mrc2(reaction)*t**3) ) * exp(-3.0D0*xi) - - end function bosch_hale - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine beamfus(beamfus0,betbm0,bp,bt,cnbeam,dene,deni,dlamie, & - ealphadt,enbeam,fdeut,ftrit,ftritbm,sigvdt,ten,tin,vol,zeffai, & - betanb,dnbeam2,palpnb) - - !! Routine to calculate beam slowing down properties - !! author: P J Knight, CCFE, Culham Science Centre - !! beamfus0: input real : multiplier for beam-background fusion calculation - !! betbm0 : input real : leading coefficient for neutral beam beta fraction - !! bp : input real : poloidal field (T) - !! bt : input real : toroidal field on axis (T) - !! cnbeam : input real : neutral beam current (A) - !! dene : input real : electron density (/m3) - !! deni : input real : fuel ion density (/m3) - !! dlamie : input real : ion-electron coulomb logarithm - !! ealphadt : input real : alpha particle birth energy (D-T) (keV) - !! enbeam : input real : neutral beam energy (keV) - !! fdeut : input real : deuterium fraction of main plasma - !! ftrit : input real : tritium fraction of main plasma - !! ftritbm: input real : tritium fraction of neutral beam - !! sigvdt : input real : profile averaged for D-T (m3/s) - !! ten : input real : density weighted average electron temperature (keV) - !! tin : input real : density weighted average ion temperature (keV) - !! vol : input real : plasma volume (m3) - !! zeffai : input real : mass weighted plasma effective charge - !! betanb : output real : neutral beam beta component - !! dnbeam2: output real : hot beam ion density (/m3) - !! palpnb : output real : alpha power from hot neutral beam ions (MW) - !! This routine calculates the beam slowing down properties. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Arguments - - real(dp), intent(in) :: beamfus0, betbm0, bp, bt, cnbeam, & - dene, deni, dlamie, ealphadt, enbeam, fdeut, ftrit, ftritbm, & - sigvdt, ten, tin, vol, zeffai - real(dp), intent(out) :: betanb, dnbeam2, palpnb - - ! Local variables - - real(dp) :: denid,denit,ecritd,ecritt,ehotnb,palpdb, & - palptb,tausl - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Velocity slowing down time - - tausl = 1.99D19 * (2.0D0*(1.0D0-ftritbm) + 3.0D0*ftritbm) * & - ten**1.5D0 / dene / dlamie - - ! Critical energy for electron/ion slowing down of the beam ion - ! (deuterium and tritium neutral beams, respectively) (keV) - - ecritd = 14.8D0 * ten * 2.0D0 * zeffai**0.6666D0 * & - (dlamie+4.0D0)/dlamie - ecritt = ecritd * 1.5D0 - - ! Deuterium and tritium ion densities - - denid = deni * fdeut - denit = deni * ftrit - - ! Perform beam calculations - - call beamcalc(denid,denit,ealphadt,enbeam,ecritd,ecritt,tausl, & - ftritbm,cnbeam,tin,vol,sigvdt,palpdb,palptb,dnbeam2,ehotnb) - - ! Neutral beam alpha power - - palpnb = beamfus0 * (palpdb + palptb) - - ! Neutral beam beta - - betanb = betbm0 * 4.03D-22 * 0.66666D0 * dnbeam2 * ehotnb / & - (bt**2 + bp**2) - - end subroutine beamfus - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine beamcalc(nd,nt,ealphadt,ebeam,ecritd,ecritt,tausbme, & - ftritbm,ibeam,ti,vol,svdt,palfdb,palftb,nhot,ehot) - - !! Neutral beam alpha power and ion energy - !! author: P J Knight, CCFE, Culham Science Centre - !! ealphadt : input real : alpha particle birth energy (D-T) (keV) - !! ebeam : input real : beam energy (keV) - !! ecritd : input real : critical energy for electron/ion slowing down of - !! the beam ion (deuterium neutral beam) (keV) - !! ecritt : input real : critical energy for beam slowing down - !! (tritium neutral beam) (keV) - !! ftritbm: input real : beam tritium fraction (0.0 = deuterium beam) - !! ibeam : input real : beam current (A) - !! nd : input real : thermal deuterium density (/m3) - !! nt : input real : thermal tritium density (/m3) - !! svdt : input real : profile averaged for D-T (m3/s) - !! tausbme: input real : beam ion slowing down time on electrons (s) - !! ti : input real : thermal ion temperature (keV) - !! vol : input real : plasma volume (m3) (95% flux surface) - !! ehot : output real : average hot beam ion energy (keV) - !! nhot : output real : hot beam ion density (/m3) - !! palfdb : output real : alpha power from deut. beam-background fusion (MW) - !! palftb : output real : alpha power from trit. beam-background fusion (MW) - !! This routine calculates the neutral beam alpha power and ion energy. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: echarge, mproton - - implicit none - - ! Arguments - real(kind(1.0D0)), intent(in) :: ealphadt, ebeam, ecritd, ecritt, & - ftritbm, ibeam, nd, nt, svdt, tausbme, ti, vol - real(kind(1.0D0)), intent(out) :: ehot, nhot, palfdb, palftb - - ! Local variables - integer :: iabm - real(kind(1.0D0)) :: ebmratd, ebmratt, ehotd, ehott, ifbmd, ifbmt, & - ndhot, nhotmsd, nhotmst, nthot, presd, prest, s0d, s0t, svdhotn, & - svthotn, tauseffd, tausefft, vcds, vcritd, vcritt, vcts, xcoefd, & - xcoeft - real(kind(1.0D0)) :: atmd, atmt, epsabs, epsrel - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Initialise shared variables - atmd = 2.0D0 ! atomic mass of deuterium - atmt = 3.0D0 ! atomic mass of tritium - epsabs = 1.0D-7 ! absolute error - epsrel = 1.0D-7 ! relative error - - ! D and T beam current fractions - ifbmd = ibeam * (1.0D0 - ftritbm) - ifbmt = ibeam * ftritbm - - ebmratd = ebeam/ecritd - vcritd = sqrt(2.0D0*echarge*1000.0D0*ecritd/(mproton*atmd)) - tauseffd = tausbme/3.0D0 * log(1.0D0+(ebmratd)**1.5D0) - nhotmsd = (1.0D0-ftritbm) * ibeam * tauseffd/(echarge * vol) - - ebmratt = ebeam/ecritt - vcritt = sqrt(2.0D0*echarge*1000.0D0*ecritt/(mproton*atmt)) - tausefft = tausbme/3.0D0 * log(1.0D0+(ebmratt)**1.5D0) - nhotmst = ftritbm * ibeam * tausefft/(echarge * vol) - - nhot = nhotmsd + nhotmst - ndhot = nhotmsd - nthot = nhotmst - - ! Average hot ion energy from Deng & Emmert, UWFDM-718, Jan 87 - vcds = 2.0D0 * ecritd * echarge * 1000.0D0/(2.0D0 * mproton) - vcts = 2.0D0 * ecritt * echarge * 1000.0D0/(3.0D0 * mproton) - - s0d = ifbmd/(echarge * vol) - s0t = ifbmt/(echarge * vol) - - xcoefd = atmd * mproton * tausbme * vcds * s0d / & - (echarge * 1000.0D0 * 3.0D0) - xcoeft = atmt * mproton * tausbme * vcts * s0t / & - (echarge * 1000.0D0 * 3.0D0) - - presd = xcoefd * xbrak(ebeam,ecritd) - prest = xcoeft * xbrak(ebeam,ecritt) - - ehotd = 1.5D0 * presd/ndhot - ehott = 1.5D0 * prest/nthot - ehot = (ndhot*ehotd + nthot*ehott)/nhot - - iabm = 2 ; svdhotn = 1.0D-4 * sgvhot(iabm,vcritd,ebeam) - iabm = 3 ; svthotn = 1.0D-4 * sgvhot(iabm,vcritt,ebeam) - - palfdb = palphabm(ealphadt,ndhot,nt,svdhotn,vol,ti,svdt) - palftb = palphabm(ealphadt,nthot,nd,svthotn,vol,ti,svdt) - - contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function xbrak(e0,ec) - !! Hot ion energy parameter - !! author: P J Knight, CCFE, Culham Science Centre - !! e0 : input real : neutral beam energy (keV) - !! ec : input real : critical energy for electron/ion slowing down of - !! the beam ion (keV) - !! This routine calculates something to do with the hot ion energy... - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - real(dp) :: xbrak - - ! Arguments - real(dp), intent(in) :: e0, ec - - real(dp) :: ans,t1,t2,t3,t4,xarg,xc,xcs - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - xcs = e0/ec - xc = sqrt(xcs) - - t1 = xcs/2.0D0 - t2 = (log((xcs + 2.0D0*xc + 1.0D0)/(xcs - xc + 1.0D0)))/6.0D0 - - xarg = (2.0D0*xc -1.0D0)/sqrt(3.0D0) - t3 = (atan(xarg))/sqrt(3.0D0) - t4 = 0.3022999D0 - - ans = t1 + t2 - t3 - t4 - xbrak = ans - - end function xbrak - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function palphabm(ealphadt,nbm,nblk,sigv,vol,ti,svdt) - - !! Alpha power from beam-background fusion - !! author: P J Knight, CCFE, Culham Science Centre - !! ealphadt : input real : alpha particle birth energy (D-T) (keV) - !! nblk : input real : thermal ion density (/m3) - !! nbm : input real : hot beam ion density (/m3) - !! sigv : input real : hot beam fusion reaction rate (m3/s) - !! svdt : input real : profile averaged for D-T (m3/s) - !! ti : input real : thermal ion temperature (keV) - !! vol : input real : plasma volume (m3) - !! This routine calculates the alpha power from - !! beam-background fusion. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - real(dp) :: palphabm - - ! Arguments - - real(dp), intent(in) :: ealphadt,nblk,nbm,sigv,svdt,ti,vol - - ! Local variables - - real(dp) :: ratio - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ratio = svdt / bosch_hale(ti,1) - - palphabm = echarge/1000.0D0 * nbm * nblk * sigv * ealphadt * vol * ratio - - end function palphabm - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function sgvhot(iabm,vcrx,ebeam) - - !! Hot beam fusion reaction rate - !! author: P J Knight, CCFE, Culham Science Centre - !! ebeam : input real : neutral beam energy (keV) - !! iabm : input integer : switch denoting type of ion (2=D,3=T) - !! vcrx : input real : critical velocity for electron/ion slowing down of - !! the beam ion (m/s) - !! This routine calculates the hot beam fusion reaction rate in m3/s. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: idiags, report_error - use maths_library, only: quanc8 - - implicit none - - real(dp) :: sgvhot - - ! Arguments - integer, intent(in) :: iabm - real(dp), intent(in) :: ebeam, vcrx - - ! Local variables - - integer :: nofun - real(dp) :: abm,abserr,epsabs1,flag,svint,t1,t2, & - vbeam,vbeams,xv - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - epsabs1 = 1.0D-33 - - if (iabm == 2) then - abm = atmd - else if (iabm == 3) then - abm = atmt - else - idiags(1) = iabm ; call report_error(84) - end if - - ! Initialise global variables - - vcritx = vcrx - - ! Beam velocity - - vbeams = ebeam * echarge * 1000.0D0 * 2.0D0/(abm * mproton) - vbeam = sqrt(vbeams) - - xv = vbeam/vcrx - t1 = 3.0D0 * vcrx/log(1.0D0+(xv**3)) - - call quanc8(fsv,0.0D0,xv,epsabs1,epsrel,svint,abserr,nofun,flag) - t2 = svint - - sgvhot = t1 * t2 - - end function sgvhot - - end subroutine beamcalc - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function fsv(u) - - !! Integrand function for the hot beam fusion reaction rate - !! author: P J Knight, CCFE, Culham Science Centre - !! u : input real : abscissa of integration, = ratio of beam velocity - !! to the critical velocity - !! This is the integrand function for the hot beam fusion reaction rate. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: mproton, echarge - - implicit none - - real(dp) :: fsv - - ! Arguments - - real(dp), intent(in) :: u - - ! Local variables - - real(dp) :: t1,t2,xvc,xvcs - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - t1 = (u**3)/(1.0D0+u**3) - - ! vcritx : critical velocity for electron/ion slowing down of beam ion (m/s) - - xvc = vcritx*u - xvcs = xvc * xvc * mproton/(echarge * 1000.0D0) - t2 = sigbmfus(xvcs) - - fsv = t1 * t2 - - contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function sigbmfus(vrelsq) - - !! Fusion reaction cross-section - !! author: P J Knight, CCFE, Culham Science Centre - !! vrelsq : input real : square of the speed of the beam ion (keV/amu) - !! This function evaluates the fusion reaction cross-section as a - !! function of beam ion velocity (squared). - !! The functional form of the cross-section is in terms of the equivalent - !! deuterium energy, i.e. for a tritium beam at 500 keV the energy - !! used in the cross-section function is 333 keV. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - real(dp) :: sigbmfus - - ! Arguments - - real(dp), intent(in) :: vrelsq - - ! Local variables - - real(dp) :: a1,a2,a3,a4,a5,atmd,ebm,t1,t2 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - a1 = 45.95D0 - a2 = 5.02D4 - a3 = 1.368D-2 - a4 = 1.076D0 - a5 = 4.09D2 - - ! Deuterium atomic mass - - atmd = 2.0D0 - - ! Beam kinetic energy - - ebm = 0.5D0 * atmd * vrelsq - - ! Set limits on cross-section at low and high beam energies - - if (ebm < 10.0D0) then - sigbmfus = 1.0D-27 - else if (ebm > 1.0D4) then - sigbmfus = 8.0D-26 - else - t1 = a2/(1.0D0 + (a3 * ebm - a4)**2) + a5 - t2 = ebm * (exp (a1/sqrt(ebm)) - 1.0D0) - sigbmfus = 1.0D-24 * t1/t2 - end if - - end function sigbmfus - - end function fsv - - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - function plasma_elongation_IPB() - !! author: H Lux (UKAEA) - !! - !! Volume measure of plasma elongation using the IPB definition - !! - !! See Otto Kardaun et al 2008 Nucl. Fusion 48 099801 - - ! Module variables - use physics_variables, only : vol, rminor, rmajor - use constants, only : pi - - real(dp) :: plasma_elongation_IPB - !! Plasma elongation (IPB) - - plasma_elongation_IPB = vol / ( 2.0D0 * pi*pi * rminor*rminor * rmajor ) - !! \begin{equation} \kappa_{IPB} = \frac{V}{2\pi a^2 R_0} \end{equation} - !! - !! - \( V \) -- Plasma volume [m\(^3\)] - !! - \( a \) -- Plasma minor radius [m] - !! - \( R_0 \) -- Plasma major radius [m] - - end function plasma_elongation_IPB - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function total_mag_field() - !! author: J. Morris (UKAEA) - !! - !! Calculates the total magnetic field - - ! Module variables - use physics_variables, only : bt, bp - - ! Return value - real(dp) :: total_mag_field - - total_mag_field = sqrt(bt**2 + bp**2) - !! \begin{equation} B_{tot} = \sqrt{B_T^2 + B_p^2} \end{equation} - - end function total_mag_field - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function beta_poloidal() - !! author: J. Morris (UKAEA) - !! - !! Calculates total poloidal beta - - ! Module variables - use physics_variables, only : btot, bp, beta - - ! Return value - real(dp) :: beta_poloidal - - beta_poloidal = beta * ( btot/bp )**2 - !! \begin{equation} \beta_p = \beta \left( \frac{B_{tot}}{B_p} \right)^2 \end{equation} - !! See J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) - !! Page 270 ISBN 0521851076 - - end function beta_poloidal - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function res_diff_time() - !! author: J. Morris (UKAEA) - !! - !! Calculates resistive diffusion time - - ! Module variables - use physics_variables, only : rmajor, rplas, kappa95 - use constants, only : rmu0 - - ! Return value - real(dp) :: res_diff_time - - res_diff_time = 2.0D0*rmu0*rmajor / (rplas*kappa95) - !! Resistive diffusion time equals the current penetration time which is approximated by: - !! \begin{equation} t_{\text{res-diff}} \sim - !! \frac{2\mu_0.R_0}{\rho_{\text{plasma}}\kappa_{95}}\end{equation} - !! - !! * \( \mu_0 \) -- permittivity of free space [H/m] - !! * \( R_0 \) -- plasma major radius [m] - !! - \( \rho_{\text{plasma}} \) -- plasma resistivity [Ohms] - !! - \( \kappa_{95} \) -- plasma elongation at 95% flux surface - !! - !! #TODO Reference needed - - end function res_diff_time - -end module physics_functions_module diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 5644d3bc33..402a67322f 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -21,6 +21,8 @@ calculate_current_coefficient_hastie, vscalc, rether, + beta_poloidal, + res_diff_time, ) from process.plasma_profiles import PlasmaProfile from process.current_drive import CurrentDrive @@ -37,6 +39,18 @@ def physics(): return Physics(PlasmaProfile(), CurrentDrive()) +def test_beta_poloidal(): + """Test beta_poloidal()""" + betap = beta_poloidal(5.347, 0.852, 0.0307) + assert betap == pytest.approx(1.209, abs=0.001) + + +def test_res_diff_time(): + """Test res_diff_time()""" + res_time = res_diff_time(9.137, 2.909e-9, 1.65) + assert res_time == pytest.approx(4784.3, abs=0.1) + + def test_diamagnetic_fraction_hender(): """Test diamagnetic_fraction_hender().""" beta = 0.14 diff --git a/tests/unit/test_physics_functions.py b/tests/unit/test_physics_functions.py index 6198767bde..c7c384338a 100644 --- a/tests/unit/test_physics_functions.py +++ b/tests/unit/test_physics_functions.py @@ -1,60 +1,11 @@ """Unit tests for physics_functions.f90.""" + from typing import Any, NamedTuple -from process.fortran import physics_functions_module as pf from process.fortran import physics_variables as pv +from process import physics_functions +import numpy as np import pytest from pytest import approx -import numpy - - -def test_plasma_elongation_ipb(monkeypatch): - """Test plasma_elongation_IPB(). - :param monkeypatch: Mock fixture - :type monkeypatch: object - """ - monkeypatch.setattr(pv, "rmajor", 9.137) - monkeypatch.setattr(pv, "rminor", 2.947) - monkeypatch.setattr(pv, "vol", 2634.0) - kappaa_ipb = pf.plasma_elongation_ipb() - assert kappaa_ipb == approx(1.682, abs=0.001) - - -def test_total_mag_field(monkeypatch): - """Test total_mag_field(). - - :param monkeypatch: Mock fixture - :type monkeypatch: object - """ - monkeypatch.setattr(pv, "bt", 5.278) - monkeypatch.setattr(pv, "bp", 0.852) - btot = pf.total_mag_field() - assert btot == approx(5.347, abs=0.001) - - -def test_beta_poloidal(monkeypatch): - """Test beta_poloidal(). - - :param monkeypatch: Mock fixture - :type monkeypatch: object - """ - monkeypatch.setattr(pv, "btot", 5.347) - monkeypatch.setattr(pv, "beta", 0.0307) - monkeypatch.setattr(pv, "bp", 0.852) - betap = pf.beta_poloidal() - assert betap == approx(1.209, abs=0.001) - - -def test_res_diff_time(monkeypatch): - """Test res_diff_time(). - - :param monkeypatch: Mock fixture - :type monkeypatch: object - """ - monkeypatch.setattr(pv, "rmajor", 9.137) - monkeypatch.setattr(pv, "rplas", 2.909e-9) - monkeypatch.setattr(pv, "kappa95", 1.650) - res_time = pf.res_diff_time() - assert res_time == approx(4784.3, abs=0.1) class Palph2Param(NamedTuple): @@ -249,11 +200,18 @@ def test_palph2(palph2param, monkeypatch): monkeypatch.setattr(pv, "falpha", palph2param.falpha) monkeypatch.setattr(pv, "fdeut", palph2param.fdeut) - # allow inout params to be mutated - palppv = numpy.array(palph2param.palppv) - pneutpv = numpy.array(palph2param.pneutpv) - - (palpmw, pneutmw, pchargemw, betaft, palpipv, palpepv, pfuscmw, powfmw) = pf.palph2( + ( + pneutpv, + palpmw, + pneutmw, + pchargemw, + betaft, + palppv, + palpepv, + palpipv, + pfuscmw, + powfmw, + ) = physics_functions.palph2( ifalphap=palph2param.ifalphap, bp=palph2param.bp, bt=palph2param.bt, @@ -267,8 +225,8 @@ def test_palph2(palph2param, monkeypatch): ten=palph2param.ten, tin=palph2param.tin, vol=palph2param.vol, - palppv=palppv, - pneutpv=pneutpv, + palppv=palph2param.palppv, + pneutpv=palph2param.pneutpv, ) assert palppv == pytest.approx(palph2param.expected_palppv) @@ -286,13 +244,12 @@ def test_palph2(palph2param, monkeypatch): @pytest.mark.parametrize( "t, reaction, expected_bosch_hale", ( - (0.0, -1, 0.0), - (55.73, 1, 8.832857074192583e-22), - (55.73, 2, 7.067916724597656e-23), - (55.73, 3, 1.3127277533210717e-23), - (55.73, 4, 1.1329338540436287e-23), + (55.73, physics_functions.REACTION_CONSTANTS_DT, 8.832857074192583e-22), + (55.73, physics_functions.REACTION_CONSTANTS_DHE3, 7.067916724597656e-23), + (55.73, physics_functions.REACTION_CONSTANTS_DD1, 1.3127277533210717e-23), + (55.73, physics_functions.REACTION_CONSTANTS_DD2, 1.1329338540436287e-23), ), - ids=["t_0", "DT", "DHE3", "DD1", "DD2"], + ids=["DT", "DHE3", "DD1", "DD2"], ) def test_bosch_hale(t, reaction, expected_bosch_hale): """ @@ -305,134 +262,77 @@ def test_bosch_hale(t, reaction, expected_bosch_hale): :param expected_bosch_hale: expected return value from the bosch_hale function :type expected_bosch_hale: float """ - bosch_hale = pf.bosch_hale(t, reaction) + bosch_hale = physics_functions.bosch_hale( + np.array([t]), physics_functions.BoschHaleConstants(**reaction) + ) assert bosch_hale == approx(expected_bosch_hale, abs=1e-23) -class PhalphParams(NamedTuple): - alphan: float - alphat: float - deni: float - fdeut: float - fhe3: float - ftrit: float - ti: float - - ipedestal: int - te: float - rhopedt: float - te0: float - teped: float - tesep: float - tbeta: float - dene: float - rhopedn: float - ne0: float - neped: float - nesep: float - - expected_palppv: float - expected_pchargepv: float - expected_pneutpv: float - expected_sigvdt: float - expected_fusionrate: float - expected_alpharate: float - expected_protonrate: float - expected_pdtpv: float - expected_pdhe3pv: float - expected_pddpv: float +def test_beamfus(): + betanb, dnbeam2, palpnb = physics_functions.beamfus( + 1.0, + 1.5, + 0.85, + 5.3, + 130, + 7.8e19, + 6.6e19, + 17.8, + 3520.0, + 1000.0, + 0.5, + 0.5, + 1e-06, + 2.8e-22, + 13.5, + 13.5, + 1888.0, + 0.425, + ) + assert betanb == pytest.approx(0.002616169278788316) + assert dnbeam2 == pytest.approx(4.2028390908892986e17) + assert palpnb == pytest.approx(11.506114015489336) + + +def test_beamcalc(): + palfdb, palftb, nhot, ehot = physics_functions.beamcalc( + 3.3e19, + 3.3e19, + 3520.0, + 1000.0, + 276.7, + 415.0, + 1.42, + 1e-06, + 130, + 13.5, + 1888.0, + 2.8e-22, + ) -@pytest.mark.parametrize( - "phalphparams", - ( - PhalphParams( - 1.0, - 1.45, - 6.2262793637240177e19, - 0.5, - 0.0, - 0.5, - 13.84, - 1, - 12.33, - 0.94, - 28.089723663920328, - 3.7775374842470044, - 0.1, - 2.0, - 7.4321e19, - 0.94, - 9.7756974320342041e19, - 5.8300851381352219e19, - 3.4294618459618943e19, - 0.19030547335201128, - 0.000813064815368787, - 0.7616652065443165, - 3.48375533764882e-22, - 3.3979157349166214e17, - 3.3763297977777126e17, - 1030380967354505.6, - 0.9515273667600563, - 0.0, - 0.0012563779516400874, - ), - ), -) -def test_phalph(phalphparams, monkeypatch): - """ - Automatically generated Integration for palph. + assert palfdb == pytest.approx(11.489365278680932) + assert palftb == pytest.approx(1.0379265294979434e-05) + assert nhot == pytest.approx(4.1968331737565126e17) + assert ehot == pytest.approx(445.05787301616635) - This test was generated using data from tracking/baseline_2018/baseline_2018_IN.DAT. - :param palphparam: the data used to mock and assert in this test. - :type palphparam: palphparam +def test_xbrak(): + xbrak = physics_functions.xbrak(1000.0, 276.7) - :param monkeypatch: pytest fixture used to mock module/class variables - :type monkeypatch: _pytest.monkeypatch.monkeypatch - """ - monkeypatch.setattr(pv, "ipedestal", phalphparams.ipedestal) - monkeypatch.setattr(pv, "te", phalphparams.te) - monkeypatch.setattr(pv, "rhopedt", phalphparams.rhopedt) - monkeypatch.setattr(pv, "te0", phalphparams.te0) - monkeypatch.setattr(pv, "teped", phalphparams.teped) - monkeypatch.setattr(pv, "tesep", phalphparams.tesep) - monkeypatch.setattr(pv, "tbeta", phalphparams.tbeta) - monkeypatch.setattr(pv, "dene", phalphparams.dene) - monkeypatch.setattr(pv, "rhopedn", phalphparams.rhopedn) - monkeypatch.setattr(pv, "ne0", phalphparams.ne0) - monkeypatch.setattr(pv, "neped", phalphparams.neped) - monkeypatch.setattr(pv, "nesep", phalphparams.nesep) + assert xbrak == pytest.approx(1.1061397270783706) - ( - palppv, - pchargepv, - pneutpv, - sigvdt, - fusionrate, - alpharate, - protonrate, - pdtpv, - pdhe3pv, - pddpv, - ) = pf.palph( - phalphparams.alphan, - phalphparams.alphat, - phalphparams.deni, - phalphparams.fdeut, - phalphparams.fhe3, - phalphparams.ftrit, - phalphparams.ti, + +def test_palphabm(): + palphabm = physics_functions.palphabm( + 3520.0, 316000000000, 3.3e19, 7.5e-22, 1888.0, 13.5, 2.8e-22 ) - assert palppv == pytest.approx(phalphparams.expected_palppv) - assert pchargepv == pytest.approx(phalphparams.expected_pchargepv) - assert pneutpv == pytest.approx(phalphparams.expected_pneutpv) - assert sigvdt == pytest.approx(phalphparams.expected_sigvdt) - assert fusionrate == pytest.approx(phalphparams.expected_fusionrate) - assert alpharate == pytest.approx(phalphparams.expected_alpharate) - assert protonrate == pytest.approx(phalphparams.expected_protonrate) - assert pdtpv == pytest.approx(phalphparams.expected_pdtpv) - assert pdhe3pv == pytest.approx(phalphparams.expected_pdhe3pv) - assert pddpv == pytest.approx(phalphparams.expected_pddpv) + assert palphabm == pytest.approx(1.0413228502045627e-05) + + +def test_sgvhot(): + sgvhot = physics_functions.sgvhot(3, 5140000.0, 1000.0) + + assert sgvhot == pytest.approx(7.465047902975452e-18)