From eb648039136cdfa85c5e71f005d9ba5944797a5f Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 14 Aug 2024 10:47:39 +0200 Subject: [PATCH] replace coreprofiles module with an object This change replaces the `coreprofiles` module with a new `gray_plasma` module providing the same functionality without using global variables. - `read_profiles`, `read_profiles_an` are replaced by a single `load_plasma` routines that handles both profiles kind (numerical, analytical). - `scale_profiles` is merged into `load_plasma`, which besides reading the profiles from file peforms the rescaling and interpolation based on the `gray_parameters` settings. - `set_profiles_spline`, `set_profiles_an`, `unset_profiles_spline` are completely removed as the module no longer has any internal state. - `density`, `ftemp`, `fzeff` are replaced by the `abstract_plasma` type which provides the `dens`, `temp` and `zeff` methods for either `numeric_plasma` or `analytic_plasma` subtypes. --- src/coreprofiles.f90 | 519 -------------------------------------- src/gray_core.f90 | 95 ++++--- src/gray_jetto1beam.f90 | 25 +- src/gray_plasma.f90 | 536 ++++++++++++++++++++++++++++++++++++++++ src/gray_tables.f90 | 30 +-- src/main.f90 | 81 +----- 6 files changed, 622 insertions(+), 664 deletions(-) delete mode 100644 src/coreprofiles.f90 create mode 100644 src/gray_plasma.f90 diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 deleted file mode 100644 index 55a0a85..0000000 --- a/src/coreprofiles.f90 +++ /dev/null @@ -1,519 +0,0 @@ -! This module handles the loading, interpolation and evaluation of the -! plasma profiles (density, temperature, effective charge) -! -! Two kinds of profiles are supported: analytical (suffix `_an` in the -! subroutine names) or numerical (suffix `_spline`). For the latter, the -! the data is interpolated using splines. -module coreprofiles - use const_and_precisions, only : wp_, zero, one - use splines, only : spline_simple, spline_1d - - implicit none - - ! Parameters of the C² polynomial tail of the density spline - type density_tail - real(wp_) :: start ! ψ₀, start of the tail - real(wp_) :: end ! ψ₁, end of the end - real(wp_) :: value ! s(ψ₀), value at the start - real(wp_) :: deriv1 ! s'(ψ₀), first derivative at the start - real(wp_) :: deriv2 ! s"(ψ₀), second derivative at the start - end type - - ! Parameters of the analytical profiles model - type analytic_model - real(wp_) :: dens0 ! Density scaling factor - real(wp_) :: n1, n2 ! Density exponents - real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1 - real(wp_) :: t1, t2 ! Temperature exponents - real(wp_) :: zeff ! Effective charge - end type - - ! Global variables storing the state of the module - type(spline_1d), save :: dens_spline - type(spline_simple), save :: temp_spline - type(spline_simple), save :: zeff_spline - type(density_tail), save :: tail - type(analytic_model), save :: model - - integer, save :: iprof - - private - public read_profiles, read_profiles_an ! Reading data files - public scale_profiles ! Applying rescaling - public density, temp, fzeff ! Accessing interpolated data - public set_profiles_spline, set_profiles_an ! Initialising internal state - public unset_profiles_spline ! Deinitialising internal state - -contains - - subroutine density(psin, dens, ddens) - ! Computes the density its first derivative as a function of - ! normalised poloidal flux. - ! - ! Note: density has units of 10¹⁹ m⁻³. - use logger, only : log_error - - ! subroutine arguments - real(wp_), intent(in) :: psin ! normalised poloidal flux - real(wp_), intent(out) :: dens, ddens ! density and first derivative - - ! local variables - character(256) :: msg ! for log messages formatting - - ! Initialise both to zero - dens = 0 - ddens = 0 - - ! Bail out when ψ is not available - if (psin < 0) return - - if (iprof == 0) then - ! Use the analytical model - ! - ! n(ψ) = dens0⋅(1 - ψ^aln1)^aln2 - ! - if (psin > 1) return - dens = model%dens0 * (1 - psin**model%n1)**model%n2 - ddens = -model%dens0 * model%n1*model%n2 * psin**(model%n1 - 1) & - * (1 - psin**model%n1)**(model%n2 - 1) - else - ! Use the numerical data - - ! Past the tail end both density and its - ! derivatives are identically zero - if (psin >= tail%end) return - - if (psin < tail%start) then - ! Use the interpolating spline when in range - - ! Evaluate the spline - dens = dens_spline%eval(psin) - ddens = dens_spline%deriv(psin) - - ! Evaluate the spline 1st derivative - if (abs(dens) < 1.0e-10_wp_) dens = zero - else - ! Use a C² polynomial extension outside (ψ > ψ₀) - - ! The tail consists of the product p(ψ)⋅h(t), where: - ! - ! - p(ψ) is the 2nd order Taylor polynomial of the spline, - ! centered at ψ₀. See set_profiles_spline for details. - ! - ! - h(t) is a "smoothing" polynomial in the variable - ! t = (ψ - ψ₀)/(ψ₁ - ψ₀), defined as: - ! - ! h(t) = (1 - t)³⋅(1 + 3t + 6t²) - ! - ! with the following properties: - ! - ! h(0) = 1 h'(0)=0 h"(0)=0 - ! h(1) = 0 h'(1)=0 h"(1)=0 - block - real(wp_) :: dpsi, t, p, dp, h, dh - dpsi = psin - tail%start ! Δψ = (ψ - ψ₀) - - ! Taylor polynomial p(ψ) and its derivative - p = tail%value + dpsi*tail%deriv1 + dpsi**2*tail%deriv2/2 - dp = tail%deriv1 + dpsi*tail%deriv2 - - ! Smoothing polynomial h(t) and its derivative - t = dpsi/(tail%end - tail%start) - h = (1 - t)**3 * (1 + 3*t + 6*t**2) - dh = -30*(1 - t)**2 * t**2 / (tail%end - tail%start) - - dens = p*h - ddens = dp*h + p*dh - end block - end if - - if (dens < 0) then - write (msg, '("negative density:", 2(x,a,"=",g0.3))') & - 'ne', dens, 'ψ', psin - call log_error(msg, mod='coreprofiles', proc='density') - end if - - end if - end subroutine density - - - function temp(psin) - ! Computes the temperature as a function of the - ! normalised poloidal flux. - ! - ! Note: temperature has units of keV. - - ! subroutine arguments - real(wp_), intent(in) :: psin - real(wp_) :: temp - - ! local variables - real(wp_) :: proft - - temp = 0 - if (psin >= 1 .or. psin < 0) return - if (iprof == 0) then - ! Use the analytical model - ! - ! T(ψ) = (te0 - te1)⋅(1 - ψ^t1)^t2 + te1 - ! - proft = (1 - psin**model%t1)**model%t2 - temp = (model%te0 - model%te1)*proft + model%te1 - else - ! Use the interpolated numerical data - temp = temp_spline%eval(psin) - endif - end function temp - - - function fzeff(psin) - ! Computes the effective charge Zeff as a - ! function of the normalised poloidal flux. - - ! subroutine arguments - real(wp_), intent(in) :: psin - real(wp_) :: fzeff - - fzeff = one - if (psin >= 1 .or. psin < 0) return - if (iprof == 0) then - ! Use the analytical model (just a constant) - fzeff = model%zeff - else - ! Use the interpolated numerical data - fzeff = zeff_spline%eval(psin) - endif - end function fzeff - - - subroutine read_profiles(filenm, data, err, unit) - ! Reads the radial plasma profiles from `file` and store them - ! into `data`. If given, the file is opened in the `unit` number. - ! Format notes: - ! 1. The file is formatted as a table with the following columns: - ! radial coordinate, temperature, density, effective charge. - ! 2. The first line is a header specifying the number of rows. - use utils, only : get_free_unit - use gray_params, only : profiles_data, PROF_NUMERIC - use logger, only : log_error - - ! subroutine arguments - character(len=*), intent(in) :: filenm - type(profiles_data), intent(out) :: data - integer, intent(out) :: err - integer, optional, intent(in) :: unit - - ! local variables - integer :: u, i, nrows - - iprof = PROF_NUMERIC - - ! Free the arrays when already allocated - if (allocated(data%psrad)) deallocate(data%psrad) - if (allocated(data%terad)) deallocate(data%terad) - if (allocated(data%derad)) deallocate(data%derad) - if (allocated(data%zfc)) deallocate(data%zfc) - - u = get_free_unit(unit) - - ! Read number of rows and allocate the arrays - open(file=filenm, status='old', action='read', unit=u, iostat=err) - if (err /= 0) then - call log_error('opening profiles file ('//trim(filenm)//') failed!', & - mod='coreprofiles', proc="read_profiles") - err = 1 - return - end if - - read(u, *) nrows - allocate(data%psrad(nrows), data%terad(nrows), & - data%derad(nrows), data%zfc(nrows)) - - ! Read the table rows - do i=1,nrows - read(u, *) data%psrad(i), data%terad(i), data%derad(i), data%zfc(i) - end do - close(u) - - ! ?? - data%psrad(1) = max(data%psrad(1), zero) - - end subroutine read_profiles - - - subroutine read_profiles_an(filenm, data, err, unit) - ! Reads the plasma profiles `data` in the analytical format - ! from params%filenm. - ! If given, the file is opened in the `unit` number. - ! - ! The file should be formatted as follows: - ! - ! 1 dens0 n1 n2 - ! 2 te0 te1 t1 t2 - ! 3 zeff - ! - ! See `density`, `temp`, `fzeff` subroutines for the meaning - ! of the parameters (i.e. the formulae for n,T,Zeff). - use gray_params, only : profiles_data, PROF_ANALYTIC - use utils, only : get_free_unit - use logger, only : log_error - - ! subroutine arguments - character(len=*), intent(in) :: filenm - type(profiles_data), intent(out) :: data - integer, intent(out) :: err - integer, optional, intent(in) :: unit - - ! local variables - integer :: u - - u = get_free_unit(unit) - - iprof = PROF_ANALYTIC - - if (allocated(data%terad)) deallocate(data%terad) - if (allocated(data%derad)) deallocate(data%derad) - if (allocated(data%zfc)) deallocate(data%zfc) - allocate(data%terad(4), data%derad(3), data%zfc(1)) - - open(file=filenm, status='old', action='read', unit=u, iostat=err) - if (err /= 0) then - call log_error('opening profiles file ('//trim(filenm)//') failed!', & - mod='coreprofiles', proc='read_profiles_an') - err = 1 - return - end if - - read (u,*) data%derad(1:3) ! dens0, n1, n2 - read (u,*) data%terad(1:4) ! te0, te1, t1, t2 - read (u,*) data%zfc(1) ! zeff - close(u) - end subroutine read_profiles_an - - - subroutine scale_profiles(params, factb, data) - ! Rescales the temperature and density profiles by `params%factte` - ! and `params%factne`, respectively, according to the model - ! specified by `params%iscal`. - ! See the GRAY user manual for the explanation. - use gray_params, only : profiles_parameters, profiles_data - - ! subroutine arguments - type(profiles_parameters), intent(in) :: params - real(wp_), intent(in) :: factb - type(profiles_data), intent(inout) :: data - - ! local variables - real(wp_) :: aat, aan, ffact - integer :: last_te, last_ne - - if (params%iscal==0) then - aat = 2.0_wp_/3.0_wp_ - aan = 4.0_wp_/3.0_wp_ - else - aat = one - aan = one - end if - - if (params%iscal==2) then - ffact = one - else - ffact = factb - end if - - if (params%iprof==0) then - last_te = 2 - last_ne = 1 - else - last_te = size(data%terad) - last_ne = size(data%derad) - end if - - data%terad(1:last_te) = data%terad(1:last_te) * ffact**aat * params%factte - data%derad(1:last_ne) = data%derad(1:last_ne) * ffact**aan * params%factne - end subroutine scale_profiles - - - subroutine set_profiles_spline(params, data, err, launch_pos) - ! Computes splines for the plasma profiles data and stores them - ! in their respective global variables, see the top of this file. - ! - ! When `launch_pos` (cartesian launch coordinates in cm) is present, - ! the subroutine will also check that the wave launcher is strictly - ! outside the reconstructed plasma density boundary. - ! - ! `err` is 1 if I/O errors occured, 2 if other initialisation failed. - use gray_params, only : profiles_parameters, profiles_data - use logger, only : log_debug, log_info, log_warning, log_error - - ! subroutine arguments - type(profiles_parameters), intent(inout) :: params - type(profiles_data), intent(inout) :: data - integer, intent(out) :: err - real(wp_), optional, intent(in) :: launch_pos(3) - - ! local variables - integer :: n - - ! for log messages formatting - character(256) :: msg - - n = size(data%psrad) - - ! Spline interpolation of temperature and effective charge - call temp_spline%init(data%psrad, data%terad, n) - call zeff_spline%init(data%psrad, data%zfc, n) - - ! Spline interpolation of density (smooth spline) - call dens_spline%init(data%psrad, data%derad, n, range=[zero, data%psrad(n)], & - tension=params%sspld, err=err) - - ! Tension is very small, re-fit with an interpolating spline - if (err == -1) then - call log_warning('re-fitting with density curve with zero tension', & - mod='coreprofiles', proc='density') - call dens_spline%init(data%psrad, data%derad, n, & - range=[zero, data%psrad(n)], tension=zero) - else if (err > 0) then - write (msg, '(a, g0)') 'density fit failed with error ', err - call log_error(msg, mod='coreprofiles', proc='density') - err = 2 - return - else - err = 0 - end if - - ! Computation of the polynomial tail parameters - ! - ! Note: The density is the only quantity that needs to be evaluated - ! at the edge. The spline thus has to be extended to transition - ! smoothly from the last profile point to 0 outside the plasma. - block - real(wp_) :: s0, s1, s2 ! spline, 1st, 2nd derivative - real(wp_) :: delta4 ! discriminant Δ/4 of q(x) - real(wp_) :: x0, x1 ! vertex of q(x), solution - - ! Compute the coefficients of a 2nd order Taylor polinomial to - ! extend the spline beyond the last point: - ! - ! p(ψ) = s(ψ₀) + (ψ - ψ₀)s'(ψ₀) + ½(ψ - ψ₀)²s"(ψ₀) - ! - ! where s(ψ) is the spline and ψ₀ the last point. - ! - s0 = dens_spline%eval(data%psrad(n)) - s1 = dens_spline%deriv(data%psrad(n), order=1) - s2 = dens_spline%deriv(data%psrad(n), order=2) - - ! Determine where to end the tail (to ensure the density remains - ! positive) from the zeros of the Taylor polynomial p(ψ) - ! - ! Define x=(ψ - ψ₀), then p(ψ)=0 is rewritten as - ! - ! q(x) = x² + 2s'/s" x + 2s/s" = 0 - ! - ! The discriminant is Δ/4 = (s'/s")² - 2(s/s") and - ! the solutions are x = x₀ ± √(Δ/4), with x₀ = -s'/s". - ! - x0 = -s1 / s2 ! vertex of parabola y=q(x) - delta4 = (s1 / s2)**2 - 2*s0/s2 ! Δ/4 of q(x) - - if (delta4 > 0) then - ! Pick the smallest positive zero (implying >ψ₀) - x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0) - else - ! There are no zeros, use the parabola vertex - x1 = x0 - call log_debug('spline extension has no zeros', & - mod='coreprofiles', proc='set_profiles_spline') - end if - - ! Store the tail parameters - tail%start = data%psrad(n) - tail%end = tail%start + x1 - tail%value = s0 - tail%deriv1 = s1 - tail%deriv2 = s2 - end block - - ! Make sure the wave launcher does not fall inside the tail - ! Note: if it does, the initial wave conditions become - ! invalid as they are given assuming a vacuum (N=1) - if (present(launch_pos)) then - block - use equilibrium, only : pol_flux - use const_and_precisions, only : cm - real(wp_) :: R, Z, psi - - ! Convert from cartesian to cylindrical coordinates - R = hypot(launch_pos(1), launch_pos(2)) * cm ! R = √(x²+y²) - z = launch_pos(3) * cm - - ! Get the poloidal flux at the launcher - ! Note: this returns -1 when the data is not available - call pol_flux(R, z, psi) - - if (psi > tail%start .and. psi < tail%end) then - ! Fall back to the midpoint of ψ₀ and the launcher ψ - call log_warning('downscaled tail to not reach the wave launcher', & - mod='coreprofiles', proc='set_profiles_spline') - write (msg, '(a, g0.3, a, g0.3, " → ", g0.3)') & - 'launcher: ψ=', psi, ' boundary: ψ=', tail%end, (tail%start + psi)/2 - call log_debug(msg, mod='coreprofiles', proc='set_profiles_spline') - tail%end = (tail%start + psi)/2 - end if - - if (psi > 0 .and. psi < tail%start) then - ! This must be a user error, stop here - write (msg, '(a, a, g0.3, a, g0.3)') & - 'wave launcher is inside the plasma! ', & - 'launcher: ψ=', psi, ' boundary: ψ=', tail%end - call log_error(msg, mod='coreprofiles', proc='set_profiles_spline') - err = 2 - return - end if - end block - end if - - ! Set the density boundary ψ - ! Note: this is used to detect entrance in the plasma - params%psnbnd = tail%end - - write (msg, '(a,g0.4)') 'density boundary: ψ=', tail%end - call log_info(msg, mod='coreprofiles', proc='set_profiles_spline') - end subroutine set_profiles_spline - - - subroutine unset_profiles_spline - ! Unsets the splines global variables, see the top of this file. - - call dens_spline%deinit - call temp_spline%deinit - call zeff_spline%deinit - end subroutine unset_profiles_spline - - - subroutine set_profiles_an(params, data) - ! Stores the analytical profiles data in their respective - ! global variables, see the top of this file. - use gray_params, only : profiles_parameters, profiles_data - - ! subroutine arguments - type(profiles_parameters), intent(inout) :: params - type(profiles_data), intent(in) :: data - - model%te0 = data%terad(1) - model%te1 = data%terad(2) - model%t1 = data%terad(3) - model%t2 = data%terad(4) - model%dens0 = data%derad(1) - model%n1 = data%derad(2) - model%n2 = data%derad(3) - model%zeff = data%zfc(1) - - ! Define the plasma boundary to be exactly ψ=1 - ! Note: this is used to detect entrance in the plasma - params%psnbnd = one - - end subroutine set_profiles_an - -end module coreprofiles diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 59eab95..dd244d8 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -7,12 +7,12 @@ module gray_core contains - subroutine gray_main(params, data, results, error, rhout) + subroutine gray_main(params, data, plasma, results, error, rhout) use const_and_precisions, only : zero, one, comp_tiny - use coreprofiles, only : temp, fzeff use polarization, only : ellipse_to_field use types, only : table, wrap use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM + use gray_plasma, only : abstract_plasma use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & store_beam_shape, find_flux_surfaces, & kinetic_profiles, ec_resonance, inputs_maps @@ -28,9 +28,10 @@ contains use logger, only : log_info, log_debug ! subroutine arguments - type(gray_parameters), intent(inout) :: params - type(gray_data), intent(in) :: data - type(gray_results), intent(out) :: results + type(gray_parameters), intent(inout) :: params + type(gray_data), intent(in) :: data + class(abstract_plasma), intent(in) :: plasma + type(gray_results), intent(out) :: results ! Predefined grid for the output profiles (optional) real(wp_), dimension(:), intent(in), optional :: rhout @@ -132,9 +133,9 @@ contains ! ========= set environment END ========= ! Pre-determinted tables - results%tables%kinetic_profiles = kinetic_profiles(params) + results%tables%kinetic_profiles = kinetic_profiles(params, plasma) results%tables%ec_resonance = ec_resonance(params, bres) - results%tables%inputs_maps = inputs_maps(params, bres, xgcn) + results%tables%inputs_maps = inputs_maps(params, plasma, bres, xgcn) ! print ψ surface for q=3/2 and q=2/1 call find_flux_surfaces( & @@ -257,8 +258,8 @@ contains do jk=1,params%raytracing%nray if(iwait(jk)) cycle ! jk ray is waiting for next pass stv(jk) = stv(jk) + params%raytracing%dst ! current ray step - call rkstep(params, sox, bres, xgcn, yw(:,jk), ypw(:,jk), & - gri(:,jk), ggri(:,:,jk), igrad_b) + call rkstep(params, plasma, sox, bres, xgcn, yw(:,jk), & + ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b) end do ! update position and grad if(igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, du1, gri, ggri) @@ -273,9 +274,10 @@ contains ! compute derivatives with updated gradient and local plasma values xv = yw(1:3,jk) anv = yw(4:6,jk) - call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & - xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, derxg, & - anpl, anpr, ddr, ddi, dersdst, derdnm, ierrn, igrad_b) + call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & + sox, bres, xgcn,ypw(:,jk), psinv, dens, btot, bv, & + xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & + derdnm, ierrn, igrad_b) ! update global error code and print message if(ierrn/=0) then error = ior(error,ierrn) @@ -372,9 +374,9 @@ contains yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) igrad_b = 0 ! * switch to ray-tracing - call ywppla_upd(params, xv, anv, gri(:,jk), ggri(:,:,jk), sox, & - bres, xgcn, ypw(:,jk), psinv, dens, btot, bv, & - xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & + call ywppla_upd(params, plasma, xv, anv, gri(:,jk), ggri(:,:,jk), & + sox, bres, xgcn, ypw(:,jk), psinv, dens, btot, & + bv, xg, yg, derxg, anpl, anpr, ddr, ddi, dersdst, & derdnm, ierrn, igrad_b) ! * update derivatives after reflection if(ierrn/=0) then ! * update global error code and print message error = ior(error,ierrn) @@ -429,12 +431,13 @@ contains if(ierrn==0 .and. params%ecrh_cd%iwarm>0 .and. ins_pl .and. & (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check - tekev=temp(psinv) + tekev = plasma%temp(psinv) block complex(wp_) :: Npr_warm - call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, & - ak0, bres, derdnm, anpl, anpr, sox, Npr_warm, & - alpha, didp, nharm, nhf, iokhawa, ierrhcd) + call alpha_effj(params%ecrh_cd, plasma, psinv, xg, yg, dens, & + tekev, ak0, bres, derdnm, anpl, anpr, sox, & + Npr_warm, alpha, didp, nharm, nhf, iokhawa, & + ierrhcd) anprre = Npr_warm%re anprim = Npr_warm%im if(ierrhcd /= 0) then @@ -995,18 +998,21 @@ contains - subroutine rkstep(params, sox, bres, xgcn, y, yp, dgr, ddgr, igrad) + subroutine rkstep(params, plasma, sox, bres, xgcn, y, yp, dgr, ddgr, igrad) ! Runge-Kutta integrator use gray_params, only : gray_parameters + use gray_plasma, only : abstract_plasma ! subroutine arguments - type(gray_parameters), intent(in) :: params - real(wp_), intent(in) :: bres, xgcn - real(wp_), intent(inout) :: y(6) - real(wp_), intent(in) :: yp(6) - real(wp_), intent(in) :: dgr(3) - real(wp_), intent(in) :: ddgr(3,3) - integer, intent(in) :: igrad,sox + type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma + + real(wp_), intent(in) :: bres, xgcn + real(wp_), intent(inout) :: y(6) + real(wp_), intent(in) :: yp(6) + real(wp_), intent(in) :: dgr(3) + real(wp_), intent(in) :: ddgr(3,3) + integer, intent(in) :: igrad,sox ! local variables real(wp_), dimension(6) :: k1, k2, k3, k4 @@ -1024,18 +1030,20 @@ contains function f(y) real(wp_), intent(in) :: y(6) real(wp_) :: f(6) - call rhs(params, sox, bres, xgcn, y, dgr, ddgr, f, igrad) + call rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) end function end subroutine rkstep - subroutine rhs(params, sox, bres, xgcn, y, dgr, ddgr, dery, igrad) + subroutine rhs(params, plasma, sox, bres, xgcn, y, dgr, ddgr, dery, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator use gray_params, only : gray_parameters + use gray_plasma, only : abstract_plasma ! subroutine arguments type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: dgr(3) @@ -1049,7 +1057,7 @@ contains real(wp_), dimension(3,3) :: derbv xv = y(1:3) - call plas_deriv(params, xv, bres, xgcn, dens, btot, & + call plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg) anv = y(4:6) call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & @@ -1057,16 +1065,18 @@ contains end subroutine rhs - subroutine ywppla_upd(params, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & + subroutine ywppla_upd(params, plasma, xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & psinv, dens, btot, bv, xg, yg, derxg, anpl, anpr, & ddr, ddi, dersdst, derdnm, error, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use gray_errors, only : raise_error, large_npl use gray_params, only : gray_parameters + use gray_plasma, only : abstract_plasma ! subroutine rguments type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: xv(3), anv(3) real(wp_), intent(in) :: dgr(3) @@ -1086,7 +1096,7 @@ contains real(wp_), dimension(3,3) :: derbv real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ - call plas_deriv(params, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv) + call plas_deriv(params, plasma, xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv) call disp_deriv(params, anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, & dery,anpl,anpr,ddr,ddi,dersdst,derdnm) @@ -1305,15 +1315,16 @@ contains end subroutine solg3 - subroutine plas_deriv(params, xv, bres, xgcn, dens, btot, bv, & + subroutine plas_deriv(params, plasma, xv, bres, xgcn, dens, btot, bv, & derbv, xg, yg, derxg, deryg, psinv_opt) use const_and_precisions, only : zero, cm use gray_params, only : gray_parameters, EQ_VACUUM + use gray_plasma, only : abstract_plasma use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi - use coreprofiles, only : density ! subroutine arguments type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma real(wp_), dimension(3), intent(in) :: xv real(wp_), intent(in) :: xgcn, bres real(wp_), intent(out) :: dens, btot, xg, yg @@ -1385,7 +1396,7 @@ contains end if ! compute xg and derivative - call density(psinv, dens, ddensdpsi) + call plasma%density(psinv, dens, ddensdpsi) xg = xgcn*dens dxgdpsi = xgcn*ddensdpsi @@ -1735,14 +1746,14 @@ contains - subroutine alpha_effj(params, psinv, X, Y, density, temperature, & - k0, Bres, derdnm, Npl, Npr_cold, sox, & - Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) + subroutine alpha_effj(params, plasma, psinv, X, Y, density, temperature, & + k0, Bres, derdnm, Npl, Npr_cold, sox, Npr, & + alpha, dIdp, nhmin, nhmax, iokhawa, error) ! Computes the absorption coefficient and effective current use const_and_precisions, only : pi, mc2=>mc2_ use gray_params, only : ecrh_cd_parameters - use coreprofiles, only : fzeff + use gray_plasma, only : abstract_plasma use equilibrium, only : sgnbphi use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl @@ -1755,7 +1766,9 @@ contains ! Inputs ! ECRH & CD parameters - type(ecrh_cd_parameters) :: params + type(ecrh_cd_parameters), intent(in) :: params + ! plasma object + class(abstract_plasma), intent(in) :: plasma ! poloidal flux ψ real(wp_), intent(in) :: psinv ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω @@ -1855,7 +1868,7 @@ contains ! Current drive computation if (params%ieccd <= 0) return - zeff = fzeff(psinv) + zeff = plasma%zeff(psinv) ithn = 1 if (nlarmor > nhmin) ithn = 2 rhop = sqrt(psinv) diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index e96abc2..8b94cee 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -4,7 +4,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & use const_and_precisions, only: wp_ use gray_params, only: gray_parameters, gray_data, gray_results use gray_core, only: gray_main + use gray_plasma, only: numeric_plasma + implicit none ! subroutine arguments integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam @@ -21,6 +23,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! local variables type(gray_parameters) :: params type(gray_data) :: data + type(numeric_plasma) :: plasma type(gray_results) :: res logical, save :: firstcall = .true. @@ -86,19 +89,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & if (err /= 0) return end block init_equilibrium - ! Set plasma kinetic profiles - init_profiles: block - use coreprofiles, only : set_profiles_spline - - ! Copy argument arrays - data%profiles%derad = dne - data%profiles%terad = te - data%profiles%zfc = zeff - - ! Compute splines - call set_profiles_spline(params%profiles, data%profiles, err) - if (err /= 0) return - end block init_profiles + ! Compute splines of kinetic profiles + call plasma%init(params, psrad, te, dne, zeff, err) + if (err /= 0) return ! Set wave launcher parameters init_antenna: block @@ -115,7 +108,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & end block init_antenna ! Call main subroutine for the ibeam-th beam - call gray_main(params, data, res, err, rhout=sqrt(psrad)) + call gray_main(params, data, plasma, res, err, rhout=sqrt(psrad)) write_debug_files: block integer :: i, err @@ -130,13 +123,9 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Free memory free_memory: block use equilibrium, only : unset_equil_spline - use coreprofiles, only : unset_profiles_spline ! Unset global variables of the `equilibrium` module call unset_equil_spline - - ! Unset global variables of the `coreprofiles` module - call unset_profiles_spline end block free_memory ! Copy over the results diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90 new file mode 100644 index 0000000..a3f79e4 --- /dev/null +++ b/src/gray_plasma.f90 @@ -0,0 +1,536 @@ +! This module handles the loading, interpolation and evaluation of the +! plasma profiles (density, temperature, effective charge) +! +! The type `abstract_plasma` provides an abstract interface for +! accessing the profiles data; these implementations are available: +! - `numeric_plasma`, numerical data interpolated using splines +! - `analytic_plasma`, simple analytic model +! +module gray_plasma + use const_and_precisions, only : wp_, zero + use splines, only : spline_simple, spline_1d + + implicit none + + ! Generic plasma interface + type, abstract :: abstract_plasma + contains + procedure(density_sub), deferred :: density + procedure(temp_fun), deferred :: temp + procedure(zeff_fun), deferred :: zeff + end type + + ! Analytical plasma description + type, extends(abstract_plasma) :: analytic_plasma + private + real(wp_) :: dens0 ! Density scaling factor + real(wp_) :: n1, n2 ! Density exponents + real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1 + real(wp_) :: t1, t2 ! Temperature exponents + real(wp_) :: z ! Effective charge + contains + procedure :: density => analytic_density + procedure :: temp => analytic_temp + procedure :: zeff => analytic_zeff + end type + + ! Parameters of the C² polynomial tail of the density spline + type density_tail + real(wp_) :: start ! ψ₀, start of the tail + real(wp_) :: end ! ψ₁, end of the end + real(wp_) :: value ! s(ψ₀), value at the start + real(wp_) :: deriv1 ! s'(ψ₀), first derivative at the start + real(wp_) :: deriv2 ! s"(ψ₀), second derivative at the start + end type + + ! Numerical plasma description + type, extends(abstract_plasma) :: numeric_plasma + private + type(spline_1d) :: dens_spline + type(spline_simple) :: temp_spline + type(spline_simple) :: zeff_spline + type(density_tail) :: tail + contains + procedure :: init => init_splines + procedure :: density => numeric_density + procedure :: temp => numeric_temp + procedure :: zeff => numeric_zeff + end type + + abstract interface + subroutine density_sub(self, psin, dens, ddens) + ! Computes the density its first derivative as a function of + ! normalised poloidal flux. + ! + ! Note: density has units of 10¹⁹ m⁻³. + import :: abstract_plasma, wp_ + class(abstract_plasma), intent(in) :: self + real(wp_), intent(in) :: psin ! normalised poloidal flux + real(wp_), intent(out) :: dens, ddens ! density and first derivative + end subroutine density_sub + + function temp_fun(self, psin) result(temp) + ! Computes the temperature as a function of the + ! normalised poloidal flux. + ! + ! Note: temperature has units of keV. + import :: abstract_plasma, wp_ + class(abstract_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: temp + end function temp_fun + + function zeff_fun(self, psin) result(zeff) + ! Computes the effective charge Z_eff as a + ! function of the normalised poloidal flux. + import :: abstract_plasma, wp_ + class(abstract_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: zeff + end function zeff_fun + end interface + + ! Parameters of the analytical profiles model + type analytic_model + real(wp_) :: dens0 ! Density scaling factor + real(wp_) :: n1, n2 ! Density exponents + real(wp_) :: te0, te1 ! Temperature at ψ=0, ψ=1 + real(wp_) :: t1, t2 ! Temperature exponents + real(wp_) :: zeff ! Effective charge + end type + + private + public abstract_plasma ! The abstract plasma object + public analytic_plasma, numeric_plasma ! Implementations + public load_plasma ! To load plasma from file + +contains + + subroutine analytic_density(self, psin, dens, ddens) + ! subroutine arguments + class(analytic_plasma), intent(in) :: self + real(wp_), intent(in) :: psin ! normalised poloidal flux + real(wp_), intent(out) :: dens, ddens ! density and first derivative + + if (psin > 1) then + ! Outside: the density and every derivative is identically zero + dens = 0 + ddens = 0 + else + ! Inside: the density is given by n(ψ) = dens0⋅(1 - ψ^n1)^n2 + dens = self%dens0 * (1 - psin**self%n1)**self%n2 + ddens = -self%dens0 * self%n1*self%n2 * psin**(self%n1 - 1) & + * (1 - psin**self%n1)**(self%n2 - 1) + end if + end subroutine analytic_density + + + function analytic_temp(self, psin) result(temp) + ! function arguments + class(analytic_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: temp + + ! local variables + real(wp_) :: proft + + if (psin >= 1 .or. psin < 0) then + ! Outside: temperature is undefined, return 0 + temp = 0 + else + ! Inside: the temperature is given by + ! + ! T(ψ) = (te0 - te1)⋅(1 - ψ^t1)^t2 + te1 + ! + proft = (1 - psin**self%t1)**self%t2 + temp = (self%te0 - self%te1)*proft + self%te1 + end if + end function analytic_temp + + + function analytic_zeff(self, psin) result(zeff) + ! function arguments + class(analytic_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: zeff + + if (psin >= 1 .or. psin < 0) then + ! Outside: Zeff is undefined, return 0 + zeff = 0 + else + ! Inside: just a constant + zeff = self%z + end if + end function analytic_zeff + + + subroutine numeric_density(self, psin, dens, ddens) + use logger, only : log_error + + ! subroutine arguments + class(numeric_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_), intent(out) :: dens, ddens + + ! local variables + character(256) :: msg ! for log messages formatting + + ! Initialise both to zero + dens = 0 + ddens = 0 + + ! Bail out when ψ is not available + if (psin < 0) return + + ! Past the tail end both density and its + ! derivatives are identically zero + if (psin >= self%tail%end) return + + if (psin < self%tail%start) then + ! Use the interpolating spline when in range + + ! Evaluate the spline + dens = self%dens_spline%eval(psin) + ddens = self%dens_spline%deriv(psin) + + ! Evaluate the spline 1st derivative + if (abs(dens) < 1.0e-10_wp_) dens = zero + else + ! Use a C² polynomial extension outside (ψ > ψ₀) + + ! The tail consists of the product p(ψ)⋅h(t), where: + ! + ! - p(ψ) is the 2nd order Taylor polynomial of the spline, + ! centered at ψ₀. See set_profiles_spline for details. + ! + ! - h(t) is a "smoothing" polynomial in the variable + ! t = (ψ - ψ₀)/(ψ₁ - ψ₀), defined as: + ! + ! h(t) = (1 - t)³⋅(1 + 3t + 6t²) + ! + ! with the following properties: + ! + ! h(0) = 1 h'(0)=0 h"(0)=0 + ! h(1) = 0 h'(1)=0 h"(1)=0 + block + real(wp_) :: dpsi, t, p, dp, h, dh + dpsi = psin - self%tail%start ! Δψ = (ψ - ψ₀) + + ! Taylor polynomial p(ψ) and its derivative + p = self%tail%value + dpsi*self%tail%deriv1 + dpsi**2*self%tail%deriv2/2 + dp = self%tail%deriv1 + dpsi*self%tail%deriv2 + + ! Smoothing polynomial h(t) and its derivative + t = dpsi/(self%tail%end - self%tail%start) + h = (1 - t)**3 * (1 + 3*t + 6*t**2) + dh = -30*(1 - t)**2 * t**2 / (self%tail%end - self%tail%start) + + dens = p*h + ddens = dp*h + p*dh + end block + end if + + if (dens < 0) then + write (msg, '("negative density:", 2(x,a,"=",g0.3))') & + 'ne', dens, 'ψ', psin + call log_error(msg, mod='coreprofiles', proc='density') + end if + + end subroutine numeric_density + + + function numeric_temp(self, psin) result(temp) + ! function arguments + class(numeric_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: temp + + temp = 0 + + ! Outside: temperature is undefined, return 0 + if (psin >= 1 .or. psin < 0) return + + ! Inside: use the interpolated numerical data + temp = self%temp_spline%eval(psin) + end function numeric_temp + + + function numeric_zeff(self, psin) result(zeff) + ! function arguments + class(numeric_plasma), intent(in) :: self + real(wp_), intent(in) :: psin + real(wp_) :: zeff + + zeff = 1 + + ! Outside: Zeff is undefined, return 1 + if (psin >= 1 .or. psin < 0) return + + ! Inside: use the interpolated numerical data + zeff = self%zeff_spline%eval(psin) + end function numeric_zeff + + + subroutine load_plasma(params, plasma, err) + use gray_params, only : gray_parameters + use gray_params, only : PROF_ANALYTIC, PROF_NUMERIC + use gray_params, only : RHO_TOR, RHO_POL, RHO_PSI + use gray_params, only : SCALE_COLLISION, SCALE_GREENWALD, SCALE_OFF + + use logger, only : log_error, log_debug + use equilibrium, only : frhopol + + ! subroutine arguments + type(gray_parameters), intent(inout) :: params + class(abstract_plasma), allocatable, intent(out) :: plasma + integer, intent(out) :: err + + ! local variables + integer :: u ! unit number + real(wp_) :: at, an, fact ! rescaling exponents and factor + type(numeric_plasma) :: np + type(analytic_plasma) :: ap + + open(file=params%profiles%filenm, status='old', & + action='read', newunit=u, iostat=err) + if (err /= 0) then + call log_error( & + 'opening profiles file ('//trim(params%profiles%filenm)//') failed!', & + mod='coreprofiles', proc='load_plasma') + err = 1 + return + end if + + ! Set appropriate rescaling factors + select case(params%profiles%iscal) + case (SCALE_COLLISION) + call log_debug('rescaling preserving ν*', mod='coreprofiles', proc='load_plasma') + at = 2/3.0_wp_ + an = 4/3.0_wp_ + fact = params%equilibrium%factb + case (SCALE_GREENWALD) + call log_debug('rescaling preserving n_G', mod='coreprofiles', proc='load_plasma') + at = 1 + an = 1 + fact = params%equilibrium%factb + case (SCALE_OFF) + at = 1 + an = 1 + fact = 1 + end select + + select case(params%profiles%iprof) + case (PROF_ANALYTIC) + ! The analytical file format is: + ! + ! 1 dens0 n1 n2 + ! 2 te0 te1 t1 t2 + ! 3 zeff + call log_debug('loading analytical file', mod='coreprofiles', proc='load_plasma') + + ! Read model parameters + read (u, *) ap%dens0, ap%n1, ap%n2 + read (u, *) ap%te0, ap%te1, ap%t1, ap%t2 + read (u, *) ap%z + + ! Apply rescaling factors + ap%dens0 = ap%dens0 * fact**an * params%profiles%factne + ap%te0 = ap%te0 * fact**at * params%profiles%factte + ap%te1 = ap%te1 * fact**at * params%profiles%factte + + allocate(plasma, source=ap) + + case (PROF_NUMERIC) + ! The numerical file format is: + ! + ! 1 nrows + ! 2 ψ_n T_e n_e Z_eff + ! 3 ... + call log_debug('loading numerical file', mod='coreprofiles', proc='load_plasma') + + block + integer :: i, nrows + real(wp_), allocatable, dimension(:) :: psi, temp, dens, zeff + + ! Read header and allocate arrays + read(u, *) nrows + allocate(psi(nrows), temp(nrows), dens(nrows), zeff(nrows)) + + ! Read the table rows + do i=1, nrows + read(u, *) psi(i), temp(i), dens(i), zeff(i) + end do + + ! ?? + psi(1) = max(psi(1), zero) + + ! Apply rescaling factors + dens = dens * fact**an * params%profiles%factne + temp = temp * fact**at * params%profiles%factte + + ! Convert psi to ψ_n (normalised poloidal flux) + select case (params%profiles%irho) + case (RHO_TOR) ! psi is ρ_t = √Φ_n (toroidal flux) + psi = [(frhopol(psi(i))**2, i=1,nrows)] + case (RHO_POL) ! psi is ρ_p = √ψ_n (poloidal flux) + psi = psi**2 + case (RHO_PSI) ! psi is already ψ_n + end select + + ! Interpolate + call np%init(params, psi, temp, dens, zeff, err) + end block + + allocate(plasma, source=np) + + end select + + close(u) + + end subroutine load_plasma + + + subroutine init_splines(self, params, psi, temp, dens, zeff, err) + ! Computes splines for the plasma profiles data and stores them + ! in their respective global variables, see the top of this file. + ! + ! `err` is 1 if I/O errors occured, 2 if other initialisation failed. + use gray_params, only : gray_parameters + use logger, only : log_debug, log_info, log_warning, log_error + + ! subroutine arguments + class(numeric_plasma), intent(out) :: self + type(gray_parameters), intent(inout) :: params + real(wp_), dimension(:), intent(in) :: psi, temp, dens, zeff + integer, intent(out) :: err + + ! local variables + integer :: n + + ! for log messages formatting + character(256) :: msg + + n = size(psi) + + ! Spline interpolation of temperature and effective charge + call self%temp_spline%init(psi, temp, n) + call self%zeff_spline%init(psi, zeff, n) + + ! Spline interpolation of density (smooth spline) + call self%dens_spline%init(psi, dens, n, range=[zero, psi(n)], & + tension=params%profiles%sspld, err=err) + + ! Tension is very small, re-fit with an interpolating spline + if (err == -1) then + call log_warning('re-fitting with density curve with zero tension', & + mod='coreprofiles', proc='init_splines') + call self%dens_spline%init(psi, dens, n, range=[zero, psi(n)], tension=zero) + else if (err > 0) then + write (msg, '(a, g0)') 'density fit failed with error ', err + call log_error(msg, mod='coreprofiles', proc='init_splines') + err = 2 + return + else + err = 0 + end if + + ! Computation of the polynomial tail parameters + ! + ! Note: The density is the only quantity that needs to be evaluated + ! at the edge. The spline thus has to be extended to transition + ! smoothly from the last profile point to 0 outside the plasma. + block + real(wp_) :: s0, s1, s2 ! spline, 1st, 2nd derivative + real(wp_) :: delta4 ! discriminant Δ/4 of q(x) + real(wp_) :: x0, x1 ! vertex of q(x), solution + + ! Compute the coefficients of a 2nd order Taylor polinomial to + ! extend the spline beyond the last point: + ! + ! p(ψ) = s(ψ₀) + (ψ - ψ₀)s'(ψ₀) + ½(ψ - ψ₀)²s"(ψ₀) + ! + ! where s(ψ) is the spline and ψ₀ the last point. + ! + s0 = self%dens_spline%eval(psi(n)) + s1 = self%dens_spline%deriv(psi(n), order=1) + s2 = self%dens_spline%deriv(psi(n), order=2) + + ! Determine where to end the tail (to ensure the density remains + ! positive) from the zeros of the Taylor polynomial p(ψ) + ! + ! Define x=(ψ - ψ₀), then p(ψ)=0 is rewritten as + ! + ! q(x) = x² + 2s'/s" x + 2s/s" = 0 + ! + ! The discriminant is Δ/4 = (s'/s")² - 2(s/s") and + ! the solutions are x = x₀ ± √(Δ/4), with x₀ = -s'/s". + ! + x0 = -s1 / s2 ! vertex of parabola y=q(x) + delta4 = (s1 / s2)**2 - 2*s0/s2 ! Δ/4 of q(x) + + if (delta4 > 0) then + ! Pick the smallest positive zero (implying >ψ₀) + x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0) + else + ! There are no zeros, use the parabola vertex + x1 = x0 + call log_debug('spline extension has no zeros', & + mod='coreprofiles', proc='init_splines') + end if + + ! Store the tail parameters + self%tail%start = psi(n) + self%tail%end = self%tail%start + x1 + self%tail%value = s0 + self%tail%deriv1 = s1 + self%tail%deriv2 = s2 + end block + + ! Make sure the wave launcher does not fall inside the tail + ! Note: if it does, the initial wave conditions become + ! invalid as they are given assuming a vacuum (N=1) + block + use equilibrium, only : pol_flux + use const_and_precisions, only : cm + real(wp_) :: R, Z, psi + + ! Convert from cartesian to cylindrical coordinates + associate (launch_pos => params%antenna%pos) + R = hypot(launch_pos(1), launch_pos(2)) * cm ! R = √(x²+y²) + z = launch_pos(3) * cm + end associate + + ! Get the poloidal flux at the launcher + ! Note: this returns -1 when the data is not available + call pol_flux(R, z, psi) + + if (psi > self%tail%start .and. psi < self%tail%end) then + ! Fall back to the midpoint of ψ₀ and the launcher ψ + call log_warning('downscaled tail to not reach the wave launcher', & + mod='coreprofiles', proc='init_splines') + write (msg, '(a, g0.3, a, g0.3, " → ", g0.3)') & + 'launcher: ψ=', psi, ' boundary: ψ=', self%tail%end, (self%tail%start + psi)/2 + call log_debug(msg, mod='coreprofiles', proc='init_splines') + self%tail%end = (self%tail%start + psi)/2 + end if + + if (psi > 0 .and. psi < self%tail%start) then + ! This must be a user error, stop here + write (msg, '(a, a, g0.3, a, g0.3)') & + 'wave launcher is inside the plasma! ', & + 'launcher: ψ=', psi, ' boundary: ψ=', self%tail%end + call log_error(msg, mod='coreprofiles', proc='init_splines') + err = 2 + return + end if + end block + + ! Set the density boundary ψ + ! Note: this is used to detect entrance in the plasma + params%profiles%psnbnd = self%tail%end + + write (msg, '(a,g0.4)') 'density boundary: ψ=', self%tail%end + call log_info(msg, mod='coreprofiles', proc='init_splines') + end subroutine init_splines + +end module gray_plasma diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index 7c9fc75..22d6f56 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -110,17 +110,18 @@ contains end subroutine init_tables - function kinetic_profiles(params) result(tbl) + function kinetic_profiles(params, plasma) result(tbl) ! Generates the plasma kinetic profiles use gray_params, only : gray_parameters, EQ_VACUUM + use gray_plasma, only : abstract_plasma use equilibrium, only : fq, frhotor - use coreprofiles, only : density, temp use magsurf_data, only : npsi, vajphiav ! function arguments - type(gray_parameters), intent(in) :: params - type(table) :: tbl + type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma + type(table) :: tbl ! local variables integer :: i, ntail @@ -151,8 +152,8 @@ contains else J_phi = 0 end if - call density(psi_n, dens, ddens) - call tbl%append([psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi]) + call plasma%density(psi_n, dens, ddens) + call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), fq(psi_n), J_phi]) end do end function kinetic_profiles @@ -226,20 +227,21 @@ contains end function ec_resonance - function inputs_maps(params, B_res, X_norm) result(tbl) + function inputs_maps(params, plasma, B_res, X_norm) result(tbl) ! Generates 2D maps of several input quantities use const_and_precisions, only : comp_eps, cm, degree use gray_params, only : gray_parameters, EQ_VACUUM + use gray_plasma, only : abstract_plasma use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield - use coreprofiles, only : density, temp use magsurf_data, only : npsi ! function arguments - type(gray_parameters), intent(in) :: params - real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω - real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω² - type(table) :: tbl + type(gray_parameters), intent(in) :: params + class(abstract_plasma), intent(in) :: plasma + real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω + real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω² + type(table) :: tbl ! local variables integer :: j, k @@ -271,11 +273,11 @@ contains do k = 1, npsi call pol_flux(r(j), z(k), psi_n) call bfield(r(j), z(k), B_R, B_z, B_phi) - call density(psi_n, ne, dne) + call plasma%density(psi_n, ne, dne) B = sqrt(B_R**2 + B_phi**2 + B_z**2) X = X_norm*ne Y = B/B_res - Te = temp(psi_n) + Te = plasma%temp(psi_n) call tbl%append([R(j), z(k), psi_n, B_R, B_phi, & B_z, B, ne, Te, X, Y, Npl]) end do diff --git a/src/main.f90 b/src/main.f90 index 98e0b96..897fcd2 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -5,6 +5,7 @@ program main use gray_cli, only : cli_options, parse_cli_options, & parse_param_overrides use gray_core, only : gray_main + use gray_plasma, only : abstract_plasma, load_plasma use gray_params, only : gray_parameters, gray_data, gray_results, & read_gray_params, read_gray_config, & make_gray_header @@ -16,10 +17,11 @@ program main type(cli_options) :: opts ! gray_main subroutine arguments - type(gray_parameters) :: params ! Inputs - type(gray_data) :: data ! - type(gray_results) :: results ! Outputs - integer :: err ! Exit code + type(gray_parameters) :: params ! Inputs + type(gray_data) :: data ! + class(abstract_plasma), allocatable :: plasma ! + type(gray_results) :: results ! Outputs + integer :: err ! Exit code ! Store the original working directory character(len=256) :: cwd @@ -81,8 +83,8 @@ program main call init_equilibrium(params, data, err) if (err /= 0) call exit(err) - call init_profiles(params%profiles, params%equilibrium%factb, & - params%antenna%pos, data%profiles, err) + ! Read and initialise the plasma state object + call load_plasma(params, plasma, err) if (err /= 0) call exit(err) call init_misc(params, data) @@ -98,7 +100,7 @@ program main params%misc%active_tables = opts%tables ! Finally run the simulation - call gray_main(params, data, results, err) + call gray_main(params, data, plasma, results, err) end if print_res: block @@ -136,10 +138,7 @@ program main ! Free memory free_mem: block - use coreprofiles, only : unset_profiles_spline use equilibrium, only : unset_equil_spline - - call unset_profiles_spline call unset_equil_spline end block free_mem @@ -199,68 +198,6 @@ contains end subroutine init_equilibrium - subroutine init_profiles(params, factb, launch_pos, data, err) - ! Reads the plasma kinetic profiles file (containing the elecron - ! temperature, density and plasma effective charge) and initialises - ! the respective GRAY data structure. - use gray_params, only : profiles_parameters, profiles_data, & - RHO_TOR, RHO_POL, RHO_PSI, & - PROF_ANALYTIC, PROF_NUMERIC - use equilibrium, only : frhopol - use coreprofiles, only : read_profiles_an, read_profiles, & - scale_profiles, set_profiles_an, & - set_profiles_spline - use logger, only : log_debug - - ! subroutine arguments - type(profiles_parameters), intent(inout) :: params - real(wp_), intent(in) :: factb - real(wp_), intent(in) :: launch_pos(3) - type(profiles_data), intent(out) :: data - integer, intent(out) :: err - - ! local variables - integer :: i - - select case (params%iprof) - case (PROF_ANALYTIC) - call log_debug('loading analytical file', & - mod='main', proc='init_profiles') - call read_profiles_an(params%filenm, data, err) - if (err /= 0) return - - case (PROF_NUMERIC) - call log_debug('loading numerical file', & - mod='main', proc='init_profiles') - call read_profiles(params%filenm, data, err) - if (err /= 0) return - - ! Convert psrad to ψ - select case (params%irho) - case (RHO_TOR) ! psrad is ρ_t = √Φ (toroidal flux) - data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))] - case (RHO_POL) ! psrad is ρ_p = √ψ (poloidal flux) - data%psrad = data%psrad**2 - case (RHO_PSI) ! psrad is already ψ - end select - end select - - ! Rescale input data - call scale_profiles(params, factb, data) - - ! Set global variables - select case (params%iprof) - case (PROF_ANALYTIC) - call set_profiles_an(params, data) - case (PROF_NUMERIC) - call log_debug('computing splines...', mod='main', proc='init_profiles') - call set_profiles_spline(params, data, err, launch_pos) - call log_debug('splines computed', mod='main', proc='init_profiles') - end select - - end subroutine init_profiles - - subroutine init_antenna(params, err) ! Reads the wave launcher file (containing the wave frequency, launcher ! position, direction and beam description) and initialises the respective