From fdf5ef72fe9a52d81fc003c83cade15b74d3473c Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 17 Oct 2024 00:56:02 +0200 Subject: [PATCH] src/magsurf_data.f90: cleanup --- src/eccd.f90 | 17 +- src/gray_core.f90 | 5 +- src/gray_equil.f90 | 41 ++- src/gray_tables.f90 | 72 ++--- src/magsurf_data.f90 | 698 +++++++++++++++++++------------------------ src/main.f90 | 3 +- src/multipass.f90 | 41 ++- 7 files changed, 408 insertions(+), 469 deletions(-) diff --git a/src/eccd.f90 b/src/eccd.f90 index de585ad..dd94257 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -139,21 +139,24 @@ contains end subroutine setcdcoeff_cohen subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cst2,eccdpar) - use magsurf_data, only : ch,tjp,tlm,njpt,nlmt + use magsurf_data, only : cd_eff use dierckx, only : profil use logger, only : log_warning - integer, parameter :: ksp=3 + real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop real(wp_), intent(out) :: cst2 real(wp_), dimension(:), allocatable, intent(out) :: eccdpar - real(wp_), dimension(nlmt) :: chlm - integer :: nlm,ierr,npar + real(wp_), dimension(cd_eff%nknots_y) :: chlm + integer :: nlm, ierr, npar cst2=1.0_wp_-rbx if(cst2 0) then block character(256) :: msg @@ -167,7 +170,7 @@ contains eccdpar(1)=zeff eccdpar(2) = fc eccdpar(3) = rbx - eccdpar(4:3+nlm) = tlm + eccdpar(4:3+nlm) = cd_eff%knots_y eccdpar(4+nlm:npar) = chlm end subroutine setcdcoeff_ncl diff --git a/src/gray_core.f90 b/src/gray_core.f90 index abb3edb..0d0496c 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -21,7 +21,7 @@ contains print_err_raytracing, print_err_ecrh_cd use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight - use magsurf_data, only : compute_flux_averages, dealloc_surfvec + use magsurf_data, only : compute_flux_averages use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab use utils, only : vmaxmin @@ -660,7 +660,6 @@ contains ! ============ main loop END ============ ! Free memory - call dealloc_surfvec call dealloc_pec end block @@ -1799,7 +1798,7 @@ contains error = error + ierrcd if (allocated(eccdpar)) deallocate(eccdpar) - ! current drive efficiency [A⋅m/W] + ! current drive efficiency R* = / [A⋅m/W] effjcdav = rbavi*effjcd dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii) diff --git a/src/gray_equil.f90 b/src/gray_equil.f90 index b772e12..886b6cc 100644 --- a/src/gray_equil.f90 +++ b/src/gray_equil.f90 @@ -176,19 +176,22 @@ module gray_equil contains - pure subroutine b_field(self, R, z, B_R, B_z, B_phi) - ! Computes the magnetic field as a function of - ! (R, z) in cylindrical coordinates + pure subroutine b_field(self, R, z, phi, B_R, B_z, B_phi, B, B_tot) + ! Computes the magnetic field as a function of (R, z) + ! The outputs include cylindrical coordinates, cartesian vector and norm. ! ! Note: all output arguments are optional. ! subroutine arguments class(abstract_equil), intent(in) :: self real(wp_), intent(in) :: R, z + real(wp_), intent(in), optional :: phi real(wp_), intent(out), optional :: B_R, B_z, B_phi + real(wp_), intent(out), optional :: B(3), B_tot ! local variables real(wp_) :: psi_n, fpol, dpsidr, dpsidz + real(wp_) :: B_R_, B_z_, B_phi_ call self%pol_flux(R, z, psi_n, dpsidr, dpsidz) call self%pol_curr(psi_n, fpol) @@ -204,10 +207,34 @@ contains ! and carrying out the cross products: ! ! B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R + ! = F(ψ)φ^/R - ∂ψ/∂z R^/R + ∂ψ/∂R z^/R ! - if (present(B_R)) B_R = - 1/R * dpsidz * self%psi_a - if (present(B_z)) B_z = + 1/R * dpsidr * self%psi_a - if (present(B_phi)) B_phi = fpol / R + B_R_ = - 1/R * dpsidz * self%psi_a + B_z_ = + 1/R * dpsidr * self%psi_a + B_phi_ = fpol / R + + if (present(B_R)) B_R = B_R_ + if (present(B_z)) B_z = B_z_ + if (present(B_phi)) B_phi = B_phi_ + + ! Convert to Cartesian coordinates + ! + ! Since the unit vectors are given by: + ! + ! { R^ = x^ cosφ + y^ sinφ + ! { φ^ = -x^ sinφ + y^ cosφ, + ! + ! B = φ^ B_φ + R^ B_R + z^ B_z + ! = (B_R cosφ - B_φ sinφ) x^ + (B_R sinφ + B_φ cosφ) y^ + z^ B_Z + ! = x^ B_x + y^ B_y + B_Z z^ + ! + if (present(B)) then + B(1) = B_R_*cos(phi) - B_phi_*sin(phi) + B(2) = B_R_*sin(phi) + B_phi_*cos(phi) + B(3) = B_z_ + end if + + if (present(B_tot)) B_tot = norm2([B_R_, B_z_, B_phi_]) end subroutine b_field @@ -571,7 +598,7 @@ contains if (present(dfpol)) dfpol = 0 end if end subroutine numeric_pol_curr - + pure function numeric_safety(self, psi_n) result(q) ! function arguments diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index b84878c..d0de7ce 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -116,7 +116,6 @@ contains use gray_params, only : gray_parameters, EQ_VACUUM use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma - use magsurf_data, only : npsi, vajphiav ! function arguments type(gray_parameters), intent(in) :: params @@ -125,7 +124,8 @@ contains type(table) :: tbl ! local variables - integer :: i, ntail + integer, parameter :: n_grid = 100 + integer :: i, n_tail real(wp_) :: rho_t, J_phi, dens, ddens real(wp_) :: psi_n, rho_p, drho_p @@ -137,22 +137,17 @@ contains if (params%equilibrium%iequil == EQ_VACUUM) return ! Δρ_p for the uniform ρ_p grid - ! Note: we don't use ψ_n because J_phi has been - ! sampled uniformly in ρ_p (see flux_average) - drho_p = 1.0_wp_ / (npsi - 1) + drho_p = 1.0_wp_ / (n_grid - 1) ! extra points to reach ψ=ψ_bnd - ntail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p) + n_tail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p) - do i = 0, npsi + ntail + do i = 0, n_grid + n_tail rho_p = i * drho_p rho_t = equil%pol2tor(rho_p) psi_n = rho_p**2 - if (psi_n < 1) then - J_phi = vajphiav(i+1) * 1.e-6_wp_ - else - J_phi = 0 - end if + !call equil%flux_average(rho_p, J_phi=J_phi) + J_phi = 0 call plasma%density(psi_n, dens, ddens) call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), & equil%safety(psi_n), J_phi]) @@ -164,10 +159,9 @@ contains function ec_resonance(params, equil, B_res) result(tbl) ! Generates the EC resonance surface table - use const_and_precisions, only : comp_eps + use const_and_precisions, only : comp_eps, comp_huge use gray_params, only : gray_parameters, EQ_VACUUM use gray_equil, only : abstract_equil - use magsurf_data, only : npsi use types, only : item use cniteq, only : cniteq_f @@ -185,30 +179,30 @@ contains block ! local variables + integer, parameter :: n_grid=101, n_cont=2002 integer :: j, k, n integer :: n_conts, n_points(10), n_tot real(wp_) :: dr, dz, B_min, B_max - real(wp_) :: B, B_R, B_z, B_phi, B_tot(npsi, npsi) - real(wp_), dimension(2002) :: R_cont, z_cont - real(wp_), dimension(npsi) :: R, z + real(wp_) :: B, B_R, B_z, B_phi, B_tot(n_grid, n_grid) + real(wp_), dimension(n_cont) :: R_cont, z_cont + real(wp_), dimension(n_grid) :: R, z ! Build a regular (R, z) grid - dr = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1) - dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1) - do j=1,npsi + dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(n_grid - 1) + dz = (equil%z_range(2) - equil%z_range(1))/(n_grid - 1) + do j = 1, n_grid R(j) = comp_eps + equil%r_range(1) + dr*(j - 1) z(j) = equil%z_range(1) + dz*(j - 1) end do - ! B_tot on psi grid - B_max = -1.0e30_wp_ - B_min = +1.0e30_wp_ - do k = 1, npsi - do j = 1, npsi - call equil%b_field(R(j), z(k), B_R, B_z, B_phi) - B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2) - if(B_tot(j,k) >= B_max) B_max = B_tot(j,k) - if(B_tot(j,k) <= B_min) B_min = B_tot(j,k) + ! Compute magnetic field on the (R, z) grid + B_max = -comp_huge + B_min = +comp_huge + do k = 1, n_grid + do j = 1, n_grid + call equil%b_field(R(j), z(k), B_tot=B_tot(j,k)) + if (B_tot(j,k) >= B_max) B_max = B_tot(j,k) + if (B_tot(j,k) <= B_min) B_min = B_tot(j,k) end do end do @@ -237,7 +231,6 @@ contains use gray_params, only : gray_parameters, EQ_VACUUM use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma - use magsurf_data, only : npsi ! function arguments type(gray_parameters), intent(in) :: params @@ -248,11 +241,12 @@ contains type(table) :: tbl ! local variables + integer, parameter :: n_grid = 100 integer :: j, k real(wp_) :: R0, Npl0 real(wp_) :: dR, dz, B_R, B_phi, B_z, B real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl - real(wp_), dimension(npsi) :: R, z + real(wp_), dimension(n_grid) :: R, z call tbl%init(title='inputs-maps', id=72, active=is_active(params, 72), & labels=[character(64) :: 'R', 'z', 'ψ_n', 'B_R', 'B_z', 'B', & @@ -265,20 +259,19 @@ contains Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥ ! Build a regular (R, z) grid - dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1) - dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1) - do concurrent (j = 1:npsi) + dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(n_grid - 1) + dz = (equil%z_range(2) - equil%z_range(1))/(n_grid - 1) + do concurrent (j = 1:n_grid) R(j) = comp_eps + equil%r_range(1) + dR*(j - 1) z(j) = equil%z_range(1) + dz*(j - 1) end do - do j = 1, npsi + do j = 1, n_grid Npl = Npl0 * R0/r(j) - do k = 1, npsi + do k = 1, n_grid call equil%pol_flux(r(j), z(k), psi_n) - call equil%b_field(r(j), z(k), B_R, B_z, B_phi) + call equil%b_field(r(j), z(k), B_tot=B) 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 = plasma%temp(psi_n) @@ -296,7 +289,6 @@ contains use gray_params, only : gray_parameters use gray_equil, only : abstract_equil - use magsurf_data, only : npoints use logger, only : log_info use minpack, only : hybrj1 use types, only : table @@ -336,7 +328,7 @@ contains ! Compute and print the ψ_n(R,z) = ψ_n₀ contour block - real(wp_), dimension(npoints) :: R_cont, z_cont + real(wp_), dimension(101) :: R_cont, z_cont real(wp_) :: R_hi, R_lo, z_hi, z_lo ! Guesses for the high,low horizontal-tangent points diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index e38d1d2..849b02b 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -1,448 +1,372 @@ module magsurf_data use const_and_precisions, only : wp_ - use splines, only : spline_simple + use splines, only : spline_simple, spline_2d implicit none - - integer, save :: npsi, npoints !# sup mag, # punti per sup - integer, save :: njpt, nlmt - real(wp_), save :: rarea - - real(wp_), dimension(:), allocatable, save :: psicon,pstab,rhot_eq, & - rhotqv,bav,varea,vcurrp,vajphiav,qqv,ffc,vratja,vratjb - real(wp_), dimension(:), allocatable, save :: rpstab - real(wp_), dimension(:), allocatable, save :: vvol,rri,rbav,bmxpsi,bmnpsi - real(wp_), dimension(:), allocatable, save :: tjp,tlm,ch,ch01 - - real(wp_), dimension(:,:), allocatable, save :: rcon,zcon - - type(spline_simple), save :: cvol, crri, crbav, cbmx, cbmn, carea, cfc - type(spline_simple), save :: crhotq - type(spline_simple), save :: cratja, cratjb, cratjpl - type(spline_simple), save :: cdadrhot, cdvdrhot + type(spline_simple), save :: spline_volume, spline_R_J, spline_B_avg_min + type(spline_simple), save :: spline_B_max, spline_B_min, spline_area, spline_f_c + type(spline_simple), save :: spline_J_astra, spline_J_jintrac, spline_J_paral + type(spline_simple), save :: spline_dAdrho_t, spline_dVdrho_t + type(spline_2d), save :: cd_eff contains - subroutine alloc_cnt(ierr) - integer, intent(out) :: ierr - - if(npsi.le.0.or.npoints.le.0) then - ierr = -1 - return - end if - - call dealloc_cnt - allocate(psicon(npsi),rcon(npoints,npsi),zcon(npoints,npsi)) - end subroutine alloc_cnt - - subroutine dealloc_cnt - if(allocated(psicon)) deallocate(psicon) - if(allocated(rcon)) deallocate(rcon) - if(allocated(zcon)) deallocate(zcon) - end subroutine dealloc_cnt - - - subroutine alloc_surfvec(ierr) - integer, intent(out) :: ierr - - if(npsi.le.0.or.npoints.le.0) then - ierr = -1 - return - end if - - call dealloc_surfvec - call alloc_cnt(ierr) - allocate(pstab(npsi), & - rhot_eq(npsi),rhotqv(npsi),bav(npsi),bmxpsi(npsi),bmnpsi(npsi),varea(npsi), & - vvol(npsi),vcurrp(npsi),vajphiav(npsi),qqv(npsi),ffc(npsi),vratja(npsi), & - vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi)) - end subroutine alloc_surfvec - - subroutine dealloc_surfvec - call dealloc_cnt - if(allocated(pstab)) deallocate(pstab) - if(allocated(rhot_eq)) deallocate(rhot_eq) - if(allocated(rhotqv)) deallocate(rhotqv) - if(allocated(bav)) deallocate(bav) - if(allocated(bmxpsi)) deallocate(bmxpsi) - if(allocated(bmnpsi)) deallocate(bmnpsi) - if(allocated(varea)) deallocate(varea) - if(allocated(vvol)) deallocate(vvol) - if(allocated(vcurrp)) deallocate(vcurrp) - if(allocated(vajphiav)) deallocate(vajphiav) - if(allocated(qqv)) deallocate(qqv) - if(allocated(ffc)) deallocate(ffc) - if(allocated(vratja)) deallocate(vratja) - if(allocated(vratjb)) deallocate(vratjb) - if(allocated(rpstab)) deallocate(rpstab) - if(allocated(rri)) deallocate(rri) - if(allocated(rbav)) deallocate(rbav) - if(allocated(tjp)) deallocate(tjp,tlm,ch) - - call cvol%deinit - call crbav%deinit - call crri%deinit - call cbmx%deinit - call cbmn%deinit - call cratja%deinit - call cratjb%deinit - call cratjpl%deinit - call carea%deinit - call cfc%deinit - call cdadrhot%deinit - call cdvdrhot%deinit - end subroutine dealloc_surfvec + pure function cross(v, w) + ! z-component of the cross product of 2D vectors + real(wp_), intent(in) :: v(2), w(2) + real(wp_) :: cross + cross = v(1)*w(2) - v(2)*w(1) + end function cross subroutine compute_flux_averages(params, equil, tables) - use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use dierckx, only : regrid, coeff_parder - use types, only : table, wrap - use gray_params, only : gray_parameters, gray_tables - use gray_equil, only : abstract_equil + ! Computes the splines of the flux surface averages - ! subroutine arguments + use const_and_precisions, only : zero, one, comp_huge, pi, mu0inv + use types, only : table, wrap + use gray_params, only : gray_parameters, gray_tables + use gray_tables, only : store_flux_surface + use gray_equil, only : abstract_equil + + ! subroutine arguments type(gray_parameters), intent(in) :: params class(abstract_equil), intent(in) :: equil type(gray_tables), intent(inout), optional :: tables ! local constants - integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, & - njest=nnintp+ksp+1,nlest=nlam+ksp+1, & - lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, & - kwrk=nnintp+nlam+njest+nlest+3 + integer, parameter :: nsurf=101, npoints=101, nlam=101 ! local variables - integer :: ier,ierr,l,jp,inc,inc1,iopt,njp,nlm,ninpr - integer, dimension(kwrk) :: iwrk - real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, & - ratio_cdbtor,ratio_pltor,fc,height,r2iav,currp, & - area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, & - dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, & - shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, & - bphi,brr,bzz,riav,fp,psinjp,rhopjp,rhotjp,qq,rup,rlw,zup,zlw - real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl - real(wp_), dimension(2*ncnt) :: dlpv - real(wp_), dimension(2*ncnt+1) :: bv,bpv - real(wp_), dimension(nlam) :: alam,weights - real(wp_), dimension(nnintp,nlam) :: fhlam - real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam - real(wp_), dimension(lwrk) :: wrk - real(wp_), dimension(:), allocatable :: rctemp,zctemp - real(wp_) :: fpolv, ddpsidrr, ddpsidzz + integer :: i, j, l - npsi=nnintp - npoints = 2*ncnt+1 + ! Arrays of quantities evaluated on a uniform ρ_p grid + real(wp_) :: rho_p(nsurf) ! ρ_p ∈ [0, 1) + real(wp_) :: psi_n(nsurf) ! ψ_n = ρ_p² + real(wp_) :: B_avg(nsurf) ! ⟨B⟩ = 1/N ∮ B dl/B_p, average magnetic field + real(wp_) :: J_phi_avg(nsurf) ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p, average toroidal current + real(wp_) :: I_p(nsurf) ! I_p = 1/μ₀ ∫ B_p⋅dS_φ, plasma current + real(wp_) :: q(nsurf) ! q = 1/2π ∮ B_φ/B_p dl/R, safety factor + real(wp_) :: R_J(nsurf) ! R_J = ⟨B⟩ / (F(ψ) ⋅ ⟨1/R²⟩), effective J_cd radius + real(wp_) :: area(nsurf) ! poloidal area enclosed by the flux surface + real(wp_) :: volume(nsurf) ! volume V enclosed by the flux surface + real(wp_) :: dAdrho_t(nsurf) ! dA/dρ_t, where A is the poloidal area + real(wp_) :: dVdrho_t(nsurf) ! dV/dρ_t, where V is the enclosed volume + real(wp_) :: B_min(nsurf) ! min |B| on the flux surface + real(wp_) :: B_max(nsurf) ! max |B| on the flux surface + real(wp_) :: f_c(nsurf) ! fraction of circulating particles - call alloc_surfvec(ierr) - if(allocated(tjp)) deallocate(tjp) - if(allocated(tlm)) deallocate(tlm) - if(allocated(ch)) deallocate(ch) - allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), & - rctemp(npoints),zctemp(npoints),stat=ierr) - if (ierr.ne.0) return + ! Ratios of different J_cd definitions w.r.t. ⟨J_φ⟩ + real(wp_) :: ratio_astra_tor(nsurf) ! J_cd_astra = ⟨J_cd⋅B⟩/B₀ + real(wp_) :: ratio_jintrac_tor(nsurf) ! J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩ + real(wp_) :: ratio_paral_tor(nsurf) ! Jcd_∥ = ⟨J_cd⋅B/|B|⟩ -! computation of flux surface averaged quantities + ! Arrays of quantities evaluated on the poloidal flux contour + real(wp_) :: R(npoints), z(npoints) ! R, z coordinates + real(wp_) :: B_tots(npoints) ! magnetic field + real(wp_) :: B_pols(npoints) ! poloidal field + real(wp_) :: dls(npoints) ! line element - dlam=1.0_wp_/dble(nlam-1) - do l=1,nlam-1 - alam(l)=dble(l-1)*dlam - fhlam(1,l)=sqrt(1.0_wp_-alam(l)) - ffhlam(l)=fhlam(1,l) - dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) - weights(l)=1.0_wp_ + ! Intermediate results + real(wp_) :: norm ! N = ∮ dl/B_p, normalisation of the ⟨⋅⟩ integral + real(wp_) :: B2_avg ! ⟨B²⟩ = 1/N ∮ B² dl/B_p, average squared field + real(wp_) :: dAdpsi ! dA/dψ, where A is the poloidal area + real(wp_) :: dVdpsi ! dV/dψ, where V is the enclosed volume + real(wp_) :: R2_inv_avg ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p + real(wp_) :: R_inv_avg ! ⟨R⁻¹⟩ = 1/N ∮ R⁻¹ dl/B_p + real(wp_) :: fpol ! F(ψ), poloidal current function + + ! Neoclassical CD efficiency variables + real(wp_) :: dlam, lambda(nlam) ! uniform grid λ ∈ [0,1] + real(wp_) :: H(nlam, nsurf) ! H(ρ_p, λ) sampled on a uniform ρ_p×λ grid + + ! Miscellanea + real(wp_) :: B_R, B_z, B_phi, R_hi, R_lo, z_hi, z_lo + + ! On the magnetic axis, ρ_p = 0 + psi_n(1) = 0 + rho_p(1) = 0 + area(1) = 0 + volume(1) = 0 + I_p(1) = 0 + J_phi_avg(1) = equil%tor_curr(equil%axis(1), equil%axis(2)) + B_max(1) = abs(equil%b_axis) + B_min(1) = abs(equil%b_axis) + B_avg(1) = abs(equil%b_axis) + R_J(1) = equil%axis(1) + f_c(1) = 1 + dAdrho_t(1) = 0 + dVdrho_t(1) = 0 + + ratio_paral_tor(1) = 1 + ratio_astra_tor(1) = abs(equil%b_axis / equil%b_centre) + ratio_jintrac_tor(1) = 1 + + ! The safety factor at ρ=0 is given by + ! q₀ = B₀ / √(∂²ψ/∂R² ⋅ ∂²ψ/∂z²) + ! Source: Spheromaks, Bellan + block + real(wp_) :: a, b + call equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b) + q(1) = abs(equil%b_axis) / sqrt(a*b)/equil%psi_a + end block + + ! Compute the uniform λ grid and H(0, λ) + H = 0 + dlam = one/(nlam - 1) + do l = 1, nlam + lambda(l) = (l-1) * dlam ! λ + H(l,1) = sqrt(1 - lambda(l)) ! H(0, λ) = √(1-λ) end do - weights(1)=0.5_wp_ - weights(nlam)=0.5_wp_ - alam(nlam)=1.0_wp_ - fhlam(1,nlam)=0.0_wp_ - ffhlam(nlam)=0.0_wp_ - dffhlam(nlam)=-99999.0_wp_ - jp=1 - anorm=2.0_wp_*pi*equil%axis(1)/abs(equil%b_axis) - dvdpsi=2.0_wp_*pi*anorm - dadpsi=2.0_wp_*pi/abs(equil%b_axis) - b2av=equil%b_axis**2 - ratio_cdator=abs(equil%b_axis/equil%b_centre) - ratio_cdbtor=1.0_wp_ - ratio_pltor=1.0_wp_ - fc=1.0_wp_ - call equil%pol_flux(equil%axis(1), equil%axis(2), & - ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) - qq=abs(equil%b_axis)/sqrt(ddpsidrr*ddpsidzz)/equil%psi_a - ajphiav=-ccj*(ddpsidrr+ddpsidzz)*equil%psi_a/equil%axis(1) - - psicon(1)=0.0_wp_ - rcon(:,1)=equil%axis(1) - zcon(:,1)=equil%axis(2) - pstab(1)=0.0_wp_ - rpstab(1)=0.0_wp_ - vcurrp(1)=0.0_wp_ - vajphiav(1)=ajphiav - bmxpsi(1)=abs(equil%b_axis) - bmnpsi(1)=abs(equil%b_axis) - bav(1)=abs(equil%b_axis) - rbav(1)=1.0_wp_ - rri(1)=equil%axis(1) - varea(1)=0.0_wp_ - vvol(1)=0.0_wp_ - vratjpl(1)=ratio_pltor - vratja(1)=ratio_cdator - vratjb(1)=ratio_cdbtor - ffc(1)=fc - qqv(1)=qq - dadrhotv(1)=0.0_wp_ - dvdrhotv(1)=0.0_wp_ + ! Guess for first contour reconstruction + R_hi = equil%axis(1) + R_lo = equil%axis(1) + z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10 + z_lo = equil%axis(2) - (equil%axis(2) - equil%z_boundary(1))/10 - rup=equil%axis(1) - rlw=equil%axis(1) - zup=equil%axis(2)+(equil%z_boundary(2)-equil%axis(2))/10.0_wp_ - zlw=equil%axis(2)-(equil%axis(2)-equil%z_boundary(1))/10.0_wp_ + ! For each surface with ρ_p > 0 + surfaces_loop: do i = 2, nsurf + ! Note: keep ρ_p < 1 to avoid the X point + rho_p(i) = merge((i - 1) * one/(nsurf - 1), 0.9999_wp_, i < nsurf) + psi_n(i) = rho_p(i)**2 - do jp=2,npsi - height=dble(jp-1)/dble(npsi-1) - if(jp.eq.npsi) height=0.9999_wp_ - rhopjp=height - psinjp=height*height - rhotjp=equil%pol2tor(rhopjp) - if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!' - psicon(jp)=height + call equil%flux_contour(psi_n(i), params%misc%rwall, & + R, z, R_hi, z_hi, R_lo, z_lo) - call equil%flux_contour(psinjp, params%misc%rwall, & - rctemp, zctemp, rup, zup, rlw, zlw) - rcon(:,jp) = rctemp - zcon(:,jp) = zctemp + if ((mod(i, 10) == 0 .or. i == nsurf) .and. present(tables)) then + call store_flux_surface(tables%flux_surfaces, psi_n(i), R, z) + end if - r2iav=0.0_wp_ - anorm=0.0_wp_ - dadpsi=0.0_wp_ - currp=0.0_wp_ - b2av=0.0_wp_ - area=0.0_wp_ - volume=0.0_wp_ - ajphiav=0.0_wp_ - bbav=0.0_wp_ - bmmx=-1.0e+30_wp_ - bmmn=1.0e+30_wp_ + block + ! Variables for the previous and current contour point + real(wp_) :: B_tot0, B_pol0, B_tot1, B_pol1, J_phi0, J_phi1 + real(wp_) :: dl ! line element - ajphi0 = equil%tor_curr(rctemp(1), zctemp(1)) - call equil%b_field(rctemp(1), zctemp(1), brr, bzz, bphi) - fpolv=bphi*rctemp(1) - btot0=sqrt(bphi**2+brr**2+bzz**2) - bpoloid0=sqrt(brr**2+bzz**2) - bv(1)=btot0 - bpv(1)=bpoloid0 - rpsim0=rctemp(1) + ! Initialise all integrals + area(i) = 0 + volume(i) = 0 + norm = 0 + dAdpsi = 0 + I_p(i) = 0 + B_avg(i) = 0 + B2_avg = 0 + R2_inv_avg = 0 + J_phi_avg(i) = 0 + B_max(i) = -comp_huge + B_min(i) = +comp_huge - do inc=1,npoints-1 - inc1=inc+1 - dla=sqrt((rctemp(inc)-equil%axis(1))**2+(zctemp(inc)-equil%axis(2))**2) - dlb=sqrt((rctemp(inc1)-equil%axis(1))**2+(zctemp(inc1)-equil%axis(2))**2) - dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2) - drc=(rctemp(inc1)-rctemp(inc)) + ! Evaluate integrands at the first "previous" point + J_phi0 = equil%tor_curr(R(1), z(1)) + call equil%b_field(R(1), z(1), B_R=B_R, B_z=B_z, & + B_phi=B_phi, B_tot=B_tot0) + fpol = B_phi*R(1) + B_pol0 = sqrt(B_R**2 + B_z**2) + B_tots(1) = B_tot0 + B_pols(1) = B_pol0 -! compute length, area and volume defined by psi=psinjp=height^2 - ph=0.5_wp_*(dla+dlb+dlp) - area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) - area=area+sqrt(area2) - rzp=rctemp(inc1)*zctemp(inc1) - rz=rctemp(inc)*zctemp(inc) - volume=pi*(rzp+rz)*drc+volume + contour_loop: do j = 1, npoints-1 + block + ! Compute the area by decomposing the contour into triangles + ! formed by the magnetic axis and two adjacent points, whose + ! area is given by the cross-product of the vertices, ½|v×w|. + real(wp_) :: O(2), P(2), Q(2), tri_area + O = equil%axis + P = [R(j), z(j)] + Q = [R(j+1), z(j+1)] + tri_area = abs(cross(P-O, Q-O))/2 + area(i) = area(i) + tri_area -! compute line integrals on the contour psi=psinjp=height^2 - rpsim=rctemp(inc1) - zpsim=zctemp(inc1) - call equil%b_field(rpsim, zpsim, brr, bzz) - ajphi = equil%tor_curr(rpsim, zpsim) - bphi=fpolv/rpsim - btot=sqrt(bphi**2+brr**2+bzz**2) - bpoloid=sqrt(brr**2+bzz**2) - dlpv(inc)=dlp - bv(inc1)=btot - bpv(inc1)=bpoloid - - dlph=0.5_wp_*dlp - anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) - dadpsi=dadpsi+dlph*(1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) - currp=currp+dlph*(bpoloid+bpoloid0) - b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) - bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) - r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2)+1.0_wp_/(bpoloid0*rpsim0**2)) - ajphiav=ajphiav+dlph*(ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) - - ajphi0=ajphi - rpsim0=rpsim - bpoloid0=bpoloid - btot0=btot + ! Compute the volume using the centroid theorem (V = 2πR⋅A) with + ! the centroid R of the contour as the mean of the triangles + ! centroids weighted by their area. Note: the total area cancels. + volume(i) = volume(i) + 2*pi * (O(1) + P(1)+ Q(1))/3 * tri_area -! computation maximum/minimum B values on given flux surface - if(btot.le.bmmn) bmmn=btot - if(btot.ge.bmmx) bmmx=btot - end do + ! The line differential is difference of the two vertices + dl = norm2(P - Q) + end block -! bav= [T] , b2av= [T^2] , rbav=/b_min -! anorm = int d l_p/B_p = dV/dpsi/(2pi) -! r2iav=<1/R^2> [m^-2] , -! riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), -! rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] -! currp = plasma current within psi=const + ! Evaluate integrands at the current point for line integrals ∮ dl + call equil%b_field(R(j+1), z(j+1), B_R=B_R, B_z=B_z, B_tot=B_tot1) + J_phi1 = equil%tor_curr(R(j+1), z(j+1)) + B_pol1 = sqrt(B_R**2 + B_z**2) + dls(j) = dl + B_tots(j+1) = B_tot1 + B_pols(j+1) = B_pol1 - bbav=bbav/anorm - r2iav=r2iav/anorm - dvdpsi=2.0_wp_*pi*anorm - riav=dadpsi/anorm - b2av=b2av/anorm - vcurrp(jp)=ccj*currp - vajphiav(jp)=ajphiav/dadpsi + ! Do one step of the trapezoid method + ! N = ∮ dl/B_p + ! dA/dψ = ∮ (1/R) dl/B_p + ! I_p = 1/μ₀ ∫ B_p⋅dS_φ = 1/μ₀ ∮ B_p dl + ! ⟨B⟩ = 1/N ∮ B dl/B_p + ! ⟨B²⟩ = 1/N ∮ B² dl/B_p + ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p + ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p + norm = norm + dl*(1/B_pol1 + 1/B_pol0)/2 + dAdpsi = dAdpsi + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2 + I_p(i) = I_p(i) + dl*(B_pol1 + B_pol0)/2 * mu0inv + B_avg(i) = B_avg(i) + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2 + B2_avg = B2_avg + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2 + R2_inv_avg = R2_inv_avg + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2 + J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2 -! area == varea, volume == vvol -! flux surface minor radius == (area/pi)^1/2 -! ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 -! ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / -! ratio_pltor = Jcd_||/J_phi Jcd_|| = - pstab(jp)=psinjp - rpstab(jp)=rhopjp - vvol(jp)=abs(volume) - varea(jp)=area - bav(jp)=bbav - rbav(jp)=bbav/bmmn - bmxpsi(jp)=bmmx - bmnpsi(jp)=bmmn - rri(jp)=bav(jp)/abs(fpolv*r2iav) - ratio_cdator=abs(b2av*riav/(fpolv*r2iav*equil%b_centre)) - ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) - ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) - vratjpl(jp)=ratio_pltor - vratja(jp)=ratio_cdator - vratjb(jp)=ratio_cdbtor - qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) - qqv(jp)=qq - dadrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dadpsi/pi - dvdrhotv(jp)=equil%phi_a*rhotjp/equil%safety(psinjp)*dvdpsi/pi - -! computation of fraction of circulating/trapped fraction fc, ft -! and of function H(lambda,rhop) -! ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ - fc=0.0_wp_ - shlam=0.0_wp_ - do l=nlam,1,-1 - lam=alam(l) - srl=0.0_wp_ - rl2=1.0_wp_-lam*bv(1)/bmmx - rl0=0.0_wp_ - if(rl2.gt.0) rl0=sqrt(rl2) - do inc=1,npoints-1 - rl2=1.0_wp_-lam*bv(inc+1)/bmmx - rl=0.0_wp_ - if(rl2.gt.0) rl=sqrt(rl2) - srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) - rl0=rl + ! Update the previous values + J_phi0 = J_phi1 + B_pol0 = B_pol1 + B_tot0 = B_tot1 + + ! Compute max,min magnetic field + if(B_tot1 <= B_min(i)) B_min(i) = B_tot1 + if(B_tot1 >= B_max(i)) B_max(i) = B_tot1 + end do contour_loop + end block + + ! Normalise integrals to obtain average value + B_avg(i) = B_avg(i) / norm + B2_avg = B2_avg / norm + R2_inv_avg = R2_inv_avg / norm + R_inv_avg = dAdpsi / norm + J_phi_avg(i) = J_phi_avg(i)/dAdpsi + dVdpsi = 2*pi * norm + + ! J_cd_astra/⟨J_φ⟩ where J_cd_astra = ⟨J_cd⋅B⟩/B₀ + ! J_cd_jintrac/⟨J_φ⟩ where J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩ + ! J_cd_∥/⟨J_φ⟩ where Jcd_∥ = ⟨J_cd⋅B/|B|⟩ + ratio_astra_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*equil%b_centre)) + ratio_jintrac_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*B_avg(i))) + ratio_paral_tor(i) = abs(B_avg(i)*R_inv_avg/(fpol*R2_inv_avg)) + + R_J(i) = B_avg(i)/abs(fpol*R2_inv_avg) + + ! By definition q(γ) = 1/2π ∮_γ dφ, where γ is a field line: + ! + ! { γ̅(0) = γ̅₀ + ! { dγ̅/dt = B̅(γ(t)) + ! + ! Using the poloidal angle θ, we rewrite it as: + ! + ! q(γ) = 1/2π ∮_γ dφ/dθ dθ = 1/2π ∮_γ dφ/dt / dθ/dt dθ + ! + ! By definition of directional derivative, d/dt = B̅⋅∇, so + ! + ! q(γ) = 1/2π ∮_γ B̅⋅∇φ / B̅⋅∇θ dθ = 1/2π ∮_γ B_φ/B_p ρ/R dθ + ! + ! Finally, since ρdθ=dl in the poloidal plane and B_φ=F/R: + ! + ! q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩ + ! + q(i) = abs(fpol/(2*pi) * norm * R2_inv_avg) + + ! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have + ! + ! ∂ψ/∂ρ_t = ∂ψ/∂Φ ∂Φ/∂ρ_t = 1/(2πq) Φ_a 2ρ_t + ! = Φ_a ρ_t / (πq) + block + real(wp_) :: rho_t, dpsidrho_t + + rho_t = equil%pol2tor(rho_p(i)) + dpsidrho_t = equil%phi_a * rho_t / (equil%safety(psi_n(i)) * pi) + dAdrho_t(i) = dAdpsi * dpsidrho_t + dVdrho_t(i) = dVdpsi * dpsidrho_t + end block + + ! Compute the fraction of circulating particles: + ! + ! f_c(ρ_p) = ¾ ∫₀¹ dλ λ/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩ + ! + ! and the function: + ! + ! H(ρ_p, λ) = ½ B_max/B_min/f_c ∫_λ^1 dλ'/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩ + ! + f_c(i) = 0 + + block + real(wp_) :: srl, rl0, rl1, dHdlam1, dHdlam0 + + do l = nlam, 1, -1 + ! Compute srl = ⟨√[1-λ B(ρ_p)/B_max]⟩ with the trapezoid method + srl = 0 + rl0 = sqrt(max(1 - lambda(l)*B_tots(1)/B_max(i), zero)) + do j = 1, npoints-1 + rl1 = sqrt(max(1 - lambda(l)*B_tots(j+1)/B_max(i), zero)) + srl = srl + dls(j) * (rl1/B_pols(j+1) + rl0/B_pols(j))/2 + rl0 = rl1 + end do + srl = srl / norm + + ! Do one step of the Riemann sum of f_c + f_c(i) = f_c(i) + 3/4.0_wp_ * B2_avg/B_max(i)**2 & + * dlam * lambda(l)/srl * merge(1.0_wp_, 0.5_wp_, l > 1 .and. l < nlam) + + ! Do one step of the trapezoid method for ∫_λ^1 dλ'/srl + dHdlam1 = 1 / srl + if (l /= nlam) H(l, i) = H(l+1, i) + dlam*(dHdlam1 + dHdlam0)/2 + dHdlam0 = dHdlam1 end do - srl=srl/anorm - dhlam=0.5_wp_/srl - fc=fc+lam/srl*weights(l) - if(l.eq.nlam) then - fhlam(jp,l)=0.0_wp_ - ffhlam(nlam*(jp-1)+l)=0.0_wp_ - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - else - shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam - fhlam(jp,l)=shlam - dffhlam(nlam*(jp-1)+l)=-dhlam - dhlam0=dhlam - end if - end do - fc=0.75_wp_*b2av/bmmx**2*fc*dlam - ffc(jp)=fc - - ccfh=bmmn/bmmx/fc - do l=1,nlam - ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) - dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) - end do - end do - rpstab(npsi)=1.0_wp_ - pstab(npsi)=1.0_wp_ + end block - ! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs - ! used for computations of dP/dV and J_cd - call cvol%init(rpstab, vvol, npsi) - call crbav%init(rpstab, rbav, npsi) - call crri%init(rpstab, rri, npsi) - call cbmx%init(rpstab, bmxpsi, npsi) - call cbmn%init(rpstab, bmnpsi, npsi) - call cratja%init(rpstab, vratja, npsi) - call cratjb%init(rpstab, vratjb, npsi) - call cratjpl%init(rpstab, vratjpl, npsi) - call carea%init(rpstab, varea, npsi) - call cfc%init(rpstab, ffc, npsi) - call cdadrhot%init(rpstab, dadrhotv, npsi) - call cdvdrhot%init(rpstab, dvdrhotv, npsi) + ! Add leading factors + H(:, i) = H(:, i)/2 * B_min(i)/B_max(i)/f_c(i) - ! spline interpolation of H(lambda,rhop) and dH/dlambda - iopt=0 - s=0.0_wp_ - call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam,zero,one,zero,one, & - ksp,ksp,s,njest,nlest,njp,tjp,nlm,tlm,ch,fp, & - wrk,lwrk,iwrk,kwrk,ier) - njpt=njp - nlmt=nlm + end do surfaces_loop - deallocate(rctemp, zctemp) + ! Fake ρ_p = 1 on the last surface (actually slightly < 1) + rho_p(nsurf) = 1 + psi_n(nsurf) = 1 + + ! Compute splines + call spline_volume%init(rho_p, volume, nsurf) + call spline_B_avg_min%init(rho_p, B_avg/B_min, nsurf) + call spline_R_J%init(rho_p, R_J, nsurf) + call spline_B_max%init(rho_p, B_max, nsurf) + call spline_B_min%init(rho_p, B_min, nsurf) + call spline_J_astra%init(rho_p, ratio_astra_tor, nsurf) + call spline_J_jintrac%init(rho_p, ratio_jintrac_tor, nsurf) + call spline_J_paral%init(rho_p, ratio_paral_tor, nsurf) + call spline_area%init(rho_p, area, nsurf) + call spline_f_c%init(rho_p, f_c, nsurf) + call spline_dAdrho_t%init(rho_p, dAdrho_t, nsurf) + call spline_dVdrho_t%init(rho_p, dVdrho_t, nsurf) + call cd_eff%init(rho_p, lambda, reshape(H, [nsurf*nlam]), & + nsurf, nlam, range=[zero, one, zero, one]) if (.not. present(tables)) return - - do jp=1,npsi - if (.not. tables%flux_averages%active) exit + if (.not. tables%flux_averages%active) return + do i = 1, nsurf call tables%flux_averages%append([ & - rpstab(jp), equil%pol2tor(rpstab(jp)), bav(jp), bmxpsi(jp), & - bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), & - ffc(jp), vratja(jp), vratjb(jp), qqv(jp)]) - end do - - ninpr = (npsi-1)/10 - do jp = ninpr+1, npsi, ninpr - ! Note: can't use store_flux_surface to avoid a circular dependency - if (.not. tables%flux_surfaces%active) exit - do l = 1, npoints - call tables%flux_surfaces%append( & - [wrap(l), wrap(psicon(jp)), wrap(rcon(l,jp)), wrap(zcon(l,jp))]) - end do + rho_p(i), equil%pol2tor(rho_p(i)), B_avg(i), B_max(i), & + B_min(i), area(i), volume(i), I_p(i), J_phi_avg(i), & + f_c(i), ratio_astra_tor(i), ratio_jintrac_tor(i), q(i)]) end do - end subroutine compute_flux_averages - subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, & - rri,rbav,bmn,bmx,fc,ratja,ratjb,ratjpl) + subroutine fluxval(rho_p, area, vol, dervol, dadrhot, dvdrhot, & + rri, rbav, bmn, bmx, fc, ratja, ratjb, ratjpl) use const_and_precisions, only : wp_ ! subroutine arguments - real(wp_), intent(in) :: rhop + real(wp_), intent(in) :: rho_p real(wp_), intent(out), optional :: & vol, area, rri, rbav, dervol, bmn, bmx, fc, & ratja, ratjb, ratjpl, dadrhot, dvdrhot - if (present(area)) area = carea%eval(rhop) - if (present(vol)) vol = cvol%eval(rhop) + if (present(area)) area = spline_area%eval(rho_p) + if (present(vol)) vol = spline_volume%eval(rho_p) - if (present(dervol)) dervol = cvol%deriv(rhop) - if (present(dadrhot)) dadrhot = cdadrhot%eval(rhop) - if (present(dvdrhot)) dvdrhot = cdvdrhot%eval(rhop) + if (present(dervol)) dervol = spline_volume%deriv(rho_p) + if (present(dadrhot)) dAdrhot = spline_dAdrho_t%eval(rho_p) + if (present(dvdrhot)) dVdrhot = spline_dVdrho_t%eval(rho_p) - if (present(rri)) rri = crri%eval(rhop) - if (present(rbav)) rbav = crbav%eval(rhop) - if (present(bmn)) bmn = cbmn%eval(rhop) - if (present(bmx)) bmx = cbmx%eval(rhop) - if (present(fc)) fc = cfc%eval(rhop) - - if (present(ratja)) ratja = cratja%eval(rhop) - if (present(ratjb)) ratjb = cratjb%eval(rhop) - if (present(ratjpl)) ratjpl = cratjpl%eval(rhop) + if (present(rri)) rri = spline_R_J%eval(rho_p) + if (present(rbav)) rbav = spline_B_avg_min%eval(rho_p) + if (present(bmn)) bmn = spline_B_min%eval(rho_p) + if (present(bmx)) bmx = spline_B_max%eval(rho_p) + if (present(fc)) fc = spline_f_c%eval(rho_p) + if (present(ratja)) ratja = spline_J_astra%eval(rho_p) + if (present(ratjb)) ratjb = spline_J_jintrac%eval(rho_p) + if (present(ratjpl)) ratjpl = spline_J_paral%eval(rho_p) end subroutine fluxval end module magsurf_data diff --git a/src/main.f90 b/src/main.f90 index 08771ca..358d1c6 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -238,7 +238,7 @@ contains ! (of different beams) and recomputes the summary table use gray_params, only : gray_parameters use gray_equil, only : abstract_equil - use magsurf_data, only : compute_flux_averages, dealloc_surfvec + use magsurf_data, only : compute_flux_averages use pec, only : pec_init, postproc_profiles, dealloc_pec, & rhot_tab @@ -383,7 +383,6 @@ contains end do ! Free memory - call dealloc_surfvec call dealloc_pec end subroutine sum_profiles diff --git a/src/multipass.f90 b/src/multipass.f90 index 22264da..a69d34f 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -26,7 +26,7 @@ contains logical, intent(in) :: perfect ! whether to assume perfect coupling ! local variables - real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z + real(wp_) :: R, z, phi real(wp_) :: B(3) real(wp_) :: c complex(wp_) :: e_mode(2), e_ray(2) @@ -34,15 +34,11 @@ contains ! Update the inside/outside flag iop(i) = iop(i) + 1 - ! Compute B in cartesian coordinates + ! Compute magnetic field R = norm2(x(1:2)) * cm z = x(3) * cm - cosphi = x(1)/R * cm - sinphi = x(2)/R * cm - call equil%b_field(R, z, B_R, B_z, B_phi) - B(1) = B_R*cosphi - B_phi*sinphi - B(2) = B_R*sinphi + B_phi*cosphi - B(3) = B_z + phi = atan2(x(2), x(1)) + call equil%b_field(R, z, phi, B=B) ! Get the polarisation vector of the given mode call pol_limit(N, B, Bres, sox, e_mode(1), e_mode(2)) @@ -72,33 +68,32 @@ contains end subroutine plasma_in - subroutine plasma_out(i, equil, xv, anv, bres, sox, iop, ext, eyt) + subroutine plasma_out(i, equil, x, N, Bres, sox, iop, ext, eyt) ! Ray exits plasma + use const_and_precisions, only : cm ! subroutine arguments integer, intent(in) :: i ! ray index class(abstract_equil), intent(in) :: equil ! MHD equilibrium - real(wp_), intent(in) :: xv(3), anv(3) - real(wp_), intent(in) :: bres + real(wp_), intent(in) :: x(3), N(3) + real(wp_), intent(in) :: Bres integer, intent(in) :: sox integer, intent(inout) :: iop(:) ! in/out plasma flag complex(wp_), intent(out) :: ext(:), eyt(:) ! local variables - real(wp_) :: rm, csphi, snphi, bphi, br, bz - real(wp_), dimension(3) :: bv, xmv + real(wp_) :: R, z, phi, B(3) - iop(i)=iop(i)+1 ! in->out + iop(i) = iop(i)+1 ! in->out - xmv=xv*0.01_wp_ ! convert from cm to m - rm=sqrt(xmv(1)**2+xmv(2)**2) - csphi=xmv(1)/rm - snphi=xmv(2)/rm - call equil%b_field(rm,xmv(3),br,bz,bphi) - bv(1)=br*csphi-bphi*snphi - bv(2)=br*snphi+bphi*csphi - bv(3)=bz - call pol_limit(anv,bv,bres,sox,ext(i),eyt(i)) ! polarization at plasma exit + ! Compute magnetic field + R = norm2(x(1:2)) * cm + z = x(3) * cm + phi = atan2(x(2), x(1)) + call equil%b_field(R, z, phi, B=B) + + ! Compute polarisation on the boundary + call pol_limit(N, B, Bres, sox, ext(i), eyt(i)) end subroutine plasma_out