From a6cc20e1e7a4f231eb3e9fdc4b24a6f786c90b2c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 21 Nov 2023 10:18:26 +0000 Subject: [PATCH 01/12] docs update --- .../proc-pages/eng-models/heating-and-current-drive.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/documentation/proc-pages/eng-models/heating-and-current-drive.md b/documentation/proc-pages/eng-models/heating-and-current-drive.md index 26effb3a43..cea6e5c553 100644 --- a/documentation/proc-pages/eng-models/heating-and-current-drive.md +++ b/documentation/proc-pages/eng-models/heating-and-current-drive.md @@ -17,7 +17,8 @@ The fraction of the required plasma current to be produced by non-inductive mean - `iefrf` = 9: Oscillating Field current drive (RFPs only - OBSOLETE-REMOVED), - `iefrf` = 10: ECRH user input gamma, - `iefrf` = 11: ECRH "HARE" model [^3], -- `iefrf` = 12: EBW user scaling input. Scaling (S. Freethy) +- `iefrf` = 12: EBW user scaling input[^4], +- `iefrf` = 13: ECRH O-mode cutoff with Zeff and Te [^4], (Note that, at present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile.) @@ -63,4 +64,6 @@ Switch `ignite` can be used to denote whether the plasma is ignited, i.e. fully [^2]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) -[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 \ No newline at end of file +[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 + +[^4]: Laqua, H & Maassberg, H & Marushchenko, Nikolai & Volpe, Francesco & Weller, A & Kasparek, W. (2003). *"Electron-Bernstein-Wave Current Drive in an Overdense Plasma at the Wendelstein 7-AS Stellarator"*, Physical review letters. 90. 075003. 10.1103/PhysRevLett.90.075003. \ No newline at end of file From 22c4adaa8b4d3be32bb48976436c3a9eacf8856a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 24 Nov 2023 08:53:11 +0000 Subject: [PATCH 02/12] missing cutoff in ebw --- process/current_drive.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/process/current_drive.py b/process/current_drive.py index 34958c0be4..bf43cb18e3 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -403,11 +403,39 @@ def cudriv(self, output: bool): dene20 * physics_variables.rmajor ) current_drive_variables.effcd = effrfss - # EBWs can only couple to plasma if cyclotron harmonic is above plasma density cut-off; # this behaviour is captured in the following function (ref issue #1262): # current_drive_variables.harnum = cyclotron harmonic number (fundamental used as default) # contant 'a' controls sharpness of transition + a = 0.1e0 + + fc = ( + 1.0e0 + / (2.0e0 * np.pi) + * current_drive_variables.harnum + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1.0e0 + / (2.0e0 * np.pi) + * np.sqrt( + physics_variables.dene + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + + density_factor = 0.5e0 * ( + 1.0e0 + np.tanh((2.0e0 / a) * ((fp - fc) / fp - a)) + ) + + current_drive_variables.effcd = ( + current_drive_variables.effcd * density_factor + ) + effrfss = effrfss * density_factor + else: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrf = }" From 4f5044c042b1ed5d8c532b52242e4fedbd664237 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 24 Nov 2023 12:21:06 +0000 Subject: [PATCH 03/12] initial model commit on .py --- process/current_drive.py | 96 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 95 insertions(+), 1 deletion(-) diff --git a/process/current_drive.py b/process/current_drive.py index bf43cb18e3..2d962fae82 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -194,8 +194,54 @@ def cudriv(self, output: bool): ) effcdfix = effcdfix * density_factor - effrfssfix = effrfssfix * density_factor + elif current_drive_variables.iefrffix == 13: + # ECCD model for O-mode cut-off with added Te and Zeff dependance + # Scaling author: Simon Freethy + # Ref : PROCESS issue #2994 + + fc = ( + 1 + / (2 * np.pi) + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1 + / (2 * np.pi) + * np.sqrt( + ( + (physics_variables.dene / 1.0e19) + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + ) + + xi_CD = 0.18e0 # This is tuned to the results of a GRAY study + xi_CD = xi_CD * ( + 4.8e0 / (2 + physics_variables.zeff) + ) # Zeff correction + effrfssfix = ( + xi_CD + * physics_variables.te + / ( + 3.27e0 + * physics_variables.rmajor + * (physics_variables.dene / 1.0e19) + ) + ) + + # O-mode case + f_cutoff = fp + + # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic + a = 0.1 # This controls how sharply the transition is reached + cutoff_factor = 0.5 * ( + 1 + np.tanh((2 / (a)) * ((2 * fc - f_cutoff) / fp - a)) + ) + effcdfix = effrfssfix * cutoff_factor elif current_drive_variables.iefrffix != 0: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrffix = }" @@ -436,6 +482,54 @@ def cudriv(self, output: bool): ) effrfss = effrfss * density_factor + elif current_drive_variables.iefrf == 13: + # ECCD model for O-mode cut-off with added Te and Zeff dependance + # Scaling author: Simon Freethy + # Ref : PROCESS issue #2994 + + fc = ( + 1 + / (2 * np.pi) + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1 + / (2 * np.pi) + * np.sqrt( + ( + (physics_variables.dene / 1.0e19) + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + ) + + xi_CD = 0.18e0 # This is tuned to the results of a GRAY study + xi_CD = xi_CD * ( + 4.8e0 / (2 + physics_variables.zeff) + ) # Zeff correction + effrfss = ( + xi_CD + * physics_variables.te + / ( + 3.27e0 + * physics_variables.rmajor + * (physics_variables.dene / 1.0e19) + ) + ) + + # O-mode case + f_cutoff = fp + + # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic + a = 0.1 # This controls how sharply the transition is reached + cutoff_factor = 0.5 * ( + 1 + np.tanh((2 / (a)) * ((2 * fc - f_cutoff) / fp - a)) + ) + current_drive_variables.effcd = effrfss * cutoff_factor + else: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrf = }" From b8ce1d40358bd47452759e7f148349a2f2c72c57 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 09:28:57 +0000 Subject: [PATCH 04/12] input.f90 limits extended --- source/fortran/input.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 22b23caf11..a86511b01e 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -1166,10 +1166,10 @@ subroutine parse_input_file(in_file,out_file,show_changes) call parse_real_variable('xi_ebw', xi_ebw, 0.0D0, 1.0D0, & 'User input EBW scaling for Plasma Heating') case ('iefrf') - call parse_int_variable('iefrf', iefrf, 1, 12, & + call parse_int_variable('iefrf', iefrf, 1, 13, & 'Switch for curr drive efficiency model') case ('iefrffix') - call parse_int_variable('iefrffix', iefrffix, 0, 12, & + call parse_int_variable('iefrffix', iefrffix, 0, 13, & 'Switch for 2nd curr drive efficiency model') case ('irfcd') call parse_int_variable('irfcd', irfcd, 0, 1, & From 5221bf047ec46507b8b7050a0636142a8f65dd98 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 09:35:21 +0000 Subject: [PATCH 05/12] added case 13 to output array --- process/current_drive.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index 2d962fae82..a1bb08ad58 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -272,7 +272,7 @@ def cudriv(self, output: bool): * 1.0e6 ) faccdfix = auxiliary_cdfix / physics_variables.plascur - elif current_drive_variables.iefrffix in [3, 7, 10, 11, 12]: + elif current_drive_variables.iefrffix in [3, 7, 10, 11, 12, 13]: # Injected power pinjemwfix = current_drive_variables.pinjfixmw @@ -564,7 +564,7 @@ def cudriv(self, output: bool): gamrf = effrfss * (dene20 * physics_variables.rmajor) current_drive_variables.gamcd = gamrf # ECCD - elif current_drive_variables.iefrf in [3, 7, 10, 11, 12]: + elif current_drive_variables.iefrf in [3, 7, 10, 11, 12, 13]: # Injected power (set to close to close the Steady-state current equilibrium) current_drive_variables.echpwr = ( 1.0e-6 From 0889df6f77c18e40aa464afdb5ea066733f3ac20 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 09:38:51 +0000 Subject: [PATCH 06/12] output in mfile --- process/current_drive.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/process/current_drive.py b/process/current_drive.py index a1bb08ad58..a13cd03a3d 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -710,6 +710,11 @@ def cudriv(self, output: bool): ) elif current_drive_variables.iefrf == 12: po.ocmmnt(self.outfile, "EBW current drive") + elif current_drive_variables.iefrf == 13: + po.ocmmnt( + self.outfile, + "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", + ) po.ovarin( self.outfile, @@ -732,6 +737,11 @@ def cudriv(self, output: bool): ) elif current_drive_variables.iefrffix == 12: po.ocmmnt(self.outfile, "EBW current drive") + elif current_drive_variables.iefrffix == 13: + po.ocmmnt( + self.outfile, + "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", + ) po.ovarin( self.outfile, From 8413b5012fddcfa33212a66eaf6c55beb75b9f17 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 09:47:03 +0000 Subject: [PATCH 07/12] added harmonic number dependance --- process/current_drive.py | 12 ++++++++++-- source/fortran/current_drive_variables.f90 | 2 +- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index a13cd03a3d..6ec149b644 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -239,7 +239,11 @@ def cudriv(self, output: bool): # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic a = 0.1 # This controls how sharply the transition is reached cutoff_factor = 0.5 * ( - 1 + np.tanh((2 / (a)) * ((2 * fc - f_cutoff) / fp - a)) + 1 + + np.tanh( + (2 / (a)) + * ((current_drive_variables.harnum * fc - f_cutoff) / fp - a) + ) ) effcdfix = effrfssfix * cutoff_factor elif current_drive_variables.iefrffix != 0: @@ -526,7 +530,11 @@ def cudriv(self, output: bool): # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic a = 0.1 # This controls how sharply the transition is reached cutoff_factor = 0.5 * ( - 1 + np.tanh((2 / (a)) * ((2 * fc - f_cutoff) / fp - a)) + 1 + + np.tanh( + (2 / (a)) + * ((current_drive_variables.harnum * fc - f_cutoff) / fp - a) + ) ) current_drive_variables.effcd = effrfss * cutoff_factor diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index adbe82f375..5821dce1ac 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -246,7 +246,7 @@ subroutine init_current_drive_variables echpwr = 0.0D0 echwpow = 0.0D0 effcd = 0.0D0 - harnum = 1.0D0 + harnum = 2.0 enbeam = 1.0D3 etacd = 0.0D0 etacdfix = 0.0D0 From b40b73f2938507ec5f01442fb78c4fc8ddfdee4c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 10:06:55 +0000 Subject: [PATCH 08/12] created wave_mode variable --- source/fortran/current_drive_variables.f90 | 9 ++++++++- source/fortran/input.f90 | 7 +++++-- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index 5821dce1ac..4d30d4e6c6 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -68,8 +68,14 @@ module current_drive_variables !! current drive efficiency (A/W) real(dp) :: harnum - !! cyclotron harmonic frequency number, used in EBW cut-off + !! cyclotron harmonic frequency number, used in cut-off function + integer :: wave_mode + !! Switch for ECRH wave mode : + !! + !! - =0 O-mode + !! - =1 X-mode + real(dp) :: enbeam !! neutral beam energy (keV) (`iteration variable 19`) @@ -247,6 +253,7 @@ subroutine init_current_drive_variables echwpow = 0.0D0 effcd = 0.0D0 harnum = 2.0 + wave_mode = 0 enbeam = 1.0D3 etacd = 0.0D0 etacdfix = 0.0D0 diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index a86511b01e..78501cb1d8 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -260,7 +260,7 @@ subroutine parse_input_file(in_file,out_file,show_changes) use current_drive_variables, only: pinjfixmw, etaech, pinjalw, etanbi, & ftritbm, gamma_ecrh, pheat, rho_ecrh, beamwd, enbeam, pheatfix, bscfmax, & forbitloss, nbshield, tbeamin, feffcd, iefrf, iefrffix, irfcd, cboot, & - etalh, frbeam, harnum, xi_ebw + etalh, frbeam, harnum, xi_ebw, wave_mode use divertor_variables, only: fdfs, anginc, divdens, divclfr, c4div, & c5div, ksic, fififi, flux_exp, divplt, delld, c2div, beta_div, betao, divdum, tdiv, c6div, & omegan, prn1, fgamp, frrp, xpertin, c1div, betai, bpsout, xparain, fdiva, & @@ -1158,7 +1158,10 @@ subroutine parse_input_file(in_file,out_file,show_changes) 'User input ECRH gamma_CD') case ('harnum') call parse_real_variable('harnum', harnum, 1.0D0, 10.0D0, & - 'cyclotron harmonic frequency number') + 'Cyclotron harmonic frequency number') + case ('wave_mode') + call parse_int_variable('wave_mode', wave_mode, 0, 1, & + 'Cyclotron wave mode switch') case ('rho_ecrh') call parse_real_variable('rho_ecrh', rho_ecrh, 0.0D0, 1.0D0, & 'normalised minor radius at which electron cyclotron current drive is maximum') From de5fdcc37bc4c22520c9d71b7c04db085261925f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 10:07:29 +0000 Subject: [PATCH 09/12] added X and O mode cutoff function --- process/current_drive.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index 6ec149b644..efb29914cf 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -234,7 +234,17 @@ def cudriv(self, output: bool): ) # O-mode case - f_cutoff = fp + if current_drive_variables.wave_mode == 0: + f_cutoff = fp + + # X-mode case + elif current_drive_variables.wave_mode == 1: + f_cutoff = 0.5 * ( + fc + + np.sqrt( + current_drive_variables.harnum * fc**2 + 4 * fp**2 + ) + ) # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic a = 0.1 # This controls how sharply the transition is reached @@ -525,7 +535,17 @@ def cudriv(self, output: bool): ) # O-mode case - f_cutoff = fp + if current_drive_variables.wave_mode == 0: + f_cutoff = fp + + # X-mode case + elif current_drive_variables.wave_mode == 1: + f_cutoff = 0.5 * ( + fc + + np.sqrt( + current_drive_variables.harnum * fc**2 + 4 * fp**2 + ) + ) # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic a = 0.1 # This controls how sharply the transition is reached From 1b58cfc5b4cb22b33aa826617ee6ef1f487d51f0 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 10:36:04 +0000 Subject: [PATCH 10/12] output wave_mode and harnum --- process/current_drive.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index efb29914cf..d63aaf3808 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -557,7 +557,6 @@ def cudriv(self, output: bool): ) ) current_drive_variables.effcd = effrfss * cutoff_factor - else: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrf = }" @@ -881,13 +880,27 @@ def cudriv(self, output: bool): "(gamma_ecrh)", current_drive_variables.gamma_ecrh, ) - elif current_drive_variables.iefrf == 12: + if current_drive_variables.iefrf == 12: po.ovarre( self.outfile, "EBW plasma heating efficiency", "(xi_ebw)", current_drive_variables.xi_ebw, ) + if current_drive_variables.iefrf in [12, 13]: + po.ovarre( + self.outfile, + "EC harmonic number", + "(harnum)", + current_drive_variables.harnum, + ) + if current_drive_variables.iefrf == 13: + po.ovarin( + self.outfile, + "EC cutoff wave mode switch", + "(wave_mode)", + current_drive_variables.wave_mode, + ) if current_drive_variables.iefrffix != 0: po.ovarre( From 5a2595ce8ed6a43b3b3033bbfbac6ff7f513114f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 28 Nov 2023 10:49:24 +0000 Subject: [PATCH 11/12] remove HARE model remnants --- process/current_drive.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/process/current_drive.py b/process/current_drive.py index d63aaf3808..ff897f43bf 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -286,7 +286,7 @@ def cudriv(self, output: bool): * 1.0e6 ) faccdfix = auxiliary_cdfix / physics_variables.plascur - elif current_drive_variables.iefrffix in [3, 7, 10, 11, 12, 13]: + elif current_drive_variables.iefrffix in [3, 7, 10, 12, 13]: # Injected power pinjemwfix = current_drive_variables.pinjfixmw @@ -591,7 +591,7 @@ def cudriv(self, output: bool): gamrf = effrfss * (dene20 * physics_variables.rmajor) current_drive_variables.gamcd = gamrf # ECCD - elif current_drive_variables.iefrf in [3, 7, 10, 11, 12, 13]: + elif current_drive_variables.iefrf in [3, 7, 10, 12, 13]: # Injected power (set to close to close the Steady-state current equilibrium) current_drive_variables.echpwr = ( 1.0e-6 From f319f5ea22d4a0708f20521b647df3d413ad9750 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 11 Dec 2023 13:19:03 +0000 Subject: [PATCH 12/12] ensure new model capture in plot_proc.py --- process/io/plot_proc.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 06d07af343..9610947753 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2283,7 +2283,7 @@ def plot_current_drive_info(axis, mfile_data, scan): if (iefrf == 5) or (iefrf == 8): nbi = True axis.text(-0.05, 1, "Neutral Beam Current Drive:", ha="left", va="center") - if (iefrf == 3) or (iefrf == 7) or (iefrf == 10) or (iefrf == 11): + if (iefrf == 3) or (iefrf == 7) or (iefrf == 10) or (iefrf == 11) or (iefrf == 13): ecrh = True axis.text(-0.05, 1, "Electron Cyclotron Current Drive:", ha="left", va="center") if iefrf == 12: @@ -2301,7 +2301,13 @@ def plot_current_drive_info(axis, mfile_data, scan): if (iefrffix == 5) or (iefrffix == 8): secondary_heating = "NBI" - if (iefrffix == 3) or (iefrffix == 7) or (iefrffix == 10) or (iefrffix == 11): + if ( + (iefrffix == 3) + or (iefrffix == 7) + or (iefrffix == 10) + or (iefrffix == 11) + or (iefrffix == 13) + ): secondary_heating = "ECH" if iefrffix == 12: secondary_heating = "EBW"