From 92b3ad9bd3805975c7656596b1531e452f0a0b4f Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 19 Jul 2022 18:22:50 +0100 Subject: [PATCH] src/gray_core.f90: implement adaptive step control MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This implements a method to control the integrator step size based on the integration error and resonance conditions. The main advantages are that: - the ray trajectories have a bounded error; - the initial step size can be large as to quickly traverse the vacuum; - the results no longer depend on the choice of the step size. The error is estimated from the real part of the dispersion relation Λ(x̅, N̅), which if solved exactly should be zero. The error bound is set to a strict value when crossing the plasma boundary to ensure a correct coupling and is relaxed afterwards. Finally, when the ray is approaching a resonance the controller ensures the step size is small compared to the absorption profile. --- input/gray.ini | 5 + src/dispersion.f90 | 83 ++++++++++++- src/gray_core.f90 | 280 ++++++++++++++++++++++++++++++++++---------- src/gray_params.f90 | 2 + src/gray_params.sh | 2 +- src/multipass.f90 | 175 ++++++++++++++++++--------- 6 files changed, 421 insertions(+), 126 deletions(-) diff --git a/input/gray.ini b/input/gray.ini index 20929c4..18798d1 100644 --- a/input/gray.ini +++ b/input/gray.ini @@ -34,6 +34,11 @@ idst = 0 ; 4: 4-stage Runge-Kutta (4⁰ order) integrator = 4 +; Whether to automatically adjust the integration step +; size based on the local error. If true `dst` will set +; the initial step size. +adaptive_step = false + [ecrh_cd] ; Choice of the power absorption model diff --git a/src/dispersion.f90 b/src/dispersion.f90 index a2c0eda..05654c9 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -22,12 +22,16 @@ module dispersion real(wp_), parameter :: tmax = 5.0_wp_ real(wp_), parameter :: dtex = 2*tmax/dble(npts) + ! Maximum value for the argument of the electron distribution, + ! μ⋅(γ-1), above which the function is considered to be 0. + real(wp_), parameter :: expcr = 16.0_wp_ + ! global variables real(wp_), dimension(npts+1), save :: ttv, extv private public expinit, colddisp, warmdisp - public zetac, harmnumber + public zetac, harmnumber, resonance_width contains @@ -131,10 +135,6 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) ! local constants - ! Maximum value for μ⋅(γ-1) above which the - ! distribution function is considered to be 0. - real(wp_), parameter :: expcr = 16.0_wp_ - nhmin = 0 nhmax = 0 @@ -218,6 +218,79 @@ pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) end subroutine harmnumber +function resonance_width(Y, mu, Npl2, R) + ! Estimates the width, as the extent in major radius, of the + ! resonating plasma layer at the smallest possible harmonic + ! + ! This computes, following the same logic of `harmnumber`, the + ! smallest possible harmonic, then estimates the width as the + ! value ΔR of major radius that causes a change Δγ ≈ dγ/dR ΔR + ! in the argument of the electron distribution equals to 1/μ. + implicit none + + ! subroutine arguments + + ! CMA Y variable: Y=ω_c/ω + real(wp_), intent(in) :: Y + ! squared parallel refractive index: N∥²=N²cos²θ + real(wp_), intent(in) :: Npl2 + ! reciprocal of adimensional temperature: μ=mc²/kT + real(wp_), intent(in) :: mu + ! major radius R=√(x² + y²) + real(wp_), intent(in) :: R + ! width of the resonance + real(wp_) :: resonance_width + + ! local variables + integer :: nh, nhc, nhmin + real(wp_) :: Yc, Yn, gamma, rdu2, argexp + + ! Derivatives of the argument A ≡ μ⋅γ(Y, N∥, n) + real(wp_) :: dAdR, dAdY + + ! local constants + + nhmin = 0 ! Minimum harmonic number + resonance_width = 0 ! The result, ΔR + + Yc = sqrt(max(1 - npl2, zero)) ! Critical Y value + nhc = ceiling(Yc/Y) ! Critical harmonic number + + ! Check a few numbers starting with nhc + do nh = nhc, nhc + 10 + Yn = Y*nh + + ! γ = [Yn - √(N∥²(Yn²-Yc²))]/Yc² + rdu2 = Yn**2 - Yc**2 + gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2) + argexp = mu*(gamma - one) + + if (argexp <= expcr) then + ! The are enough electrons with this energy + nhmin = nh + exit + end if + end do + + ! No harmonics possible + if (nhmin == 0) return + + ! The derivative of the minimum γ = γ(N∥, Y, n) is + ! + ! dγ/dY = n (Yn - γ)/(Yn - γ⋅Yc²) + ! + dAdY = mu * nhmin * (Yn - gamma)/(Yn - gamma*Yc**2) + + ! dAdR = dAdY⋅dY/dR, assuming the magnetic field varies + ! solely as B(R) ≈ B₀⋅R₀/R, then dY/dR = -Y/R + dAdR = dAdY * (-Y/R) + + ! Take ΔR = |ΔA / A/dR|, where ΔA = μΔγ = 1 + resonance_width = 1 / abs(dAdR) + +end function resonance_width + + subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & error, Npr, e, & model, nlarmor, max_iters, fallback) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index d8a1040..68b3cb9 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -6,11 +6,13 @@ module gray_core implicit none abstract interface - function rhs_function(y) result(f) + function rhs_function(y, e) result(f) ! Function passed to the integrator subroutine + ! This represent the right-hand side of the ODE: dy/ds = f(y) use const_and_precisions, only : wp_ - real(wp_), intent(in) :: y(6) - real(wp_) :: f(6) + real(wp_), intent(in) :: y(6) ! variable y + real(wp_), intent(inout), optional :: e ! error estimator + real(wp_) :: f(6) ! f(y) end function end interface @@ -85,7 +87,7 @@ contains real(wp_), dimension(:), pointer :: p0jk=>null() real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), & currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), & - psipv=>null(),chipv=>null() + psipv=>null(),chipv=>null(),dst=>null() complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() integer, dimension(:), pointer :: iiv=>null(),iop=>null(),iow=>null() logical, dimension(:), pointer :: iwait=>null() @@ -119,13 +121,11 @@ contains call pec_init(params%output%ipec, rhout) nnd = size(rhop_tab) ! number of radial profile points - call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & - stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & - pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) - - ! Allocate memory for the results... - allocate(results%dpdv(params%output%nrho)) - allocate(results%jcd(params%output%nrho)) + ! Initialise multipass module + call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, & + yw0, ypw0, stnext, stv, dst, p0ray, taus, tau1, & + etau1, cpls, cpl1, lgcpl1, jphi_beam, pins_beam, & + currins_beam, dpdv_beam, jcd_beam, psipv, chipv) ! ...and initialise them results%pabs = zero @@ -208,8 +208,8 @@ contains index_rt = index_rt +1 iO = 2*index_rt +1 ! * index_rt of O-mode derived ray (iX=iO+1) - call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam) + call initbeam(params, index_rt, iroff, iboff, iwait, stv, dst, & + jphi_beam, pins_beam, currins_beam, dpdv_beam, jcd_beam) write(msg, '(" beam: ",g0," (",a1," mode)")') index_rt, mode(iox) call log_info(msg, mod='gray_core', proc='gray_main') @@ -273,10 +273,11 @@ contains ! advance one step with "frozen" grad(S_I) 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 integrator(yw(:, jk), ypw(:, jk), rhs, & - params%raytracing%dst, & - params%raytracing%integrator) + call step_controller( & + y=yw(:, jk), yp=ypw(:, jk), f=rhs, & + h=dst(jk), method=params%raytracing%integrator, & + adaptive=params%raytracing%adaptive_step, Bres=Bres) + stv(jk) = stv(jk) + dst(jk) ! current ray step end do ! update position and grad if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) @@ -346,7 +347,7 @@ contains if(ip < params%raytracing%ipass) then ! + not last pass yynext(:,jk,index_rt) = yw0(:,jk) ! . copy starting coordinates yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step - stnext(jk,index_rt) = stv(jk) - params%raytracing%dst ! . starting step for next pass = last step + stnext(jk,index_rt) = stv(jk) - dst(jk) ! . starting step for next pass = last step if(cpl(1) < etaucr) then ! . low coupled power for O-mode => de-activate derived rays call turnoffray(jk,ip+1,2*ib-1,iroff) @@ -468,15 +469,15 @@ contains ! Computation of the ray τ, dP/ds, P(s), dI/ds, I(s) - ! optical depth: τ = ∫α(s)ds using the trapezoid rule - tau = tau0(jk) + 0.5_wp_*(alphaabs0(jk) + alpha) * dersdst * params%raytracing%dst + ! optical depth: τ = ∫α(s)ds using the trapezoidal rule + tau = tau0(jk) + 0.5_wp_*(alphaabs0(jk) + alpha) * dersdst * dst(jk) pow = p0ray(jk) * exp(-tau) ! residual power: P = P₀exp(-τ) ppabs(jk,i) = p0ray(jk) - pow ! absorbed power: P_abs = P₀ - P dids = didp * pow * alpha ! current driven: dI/ds = dI/dP⋅dP/ds = dI/dP⋅P⋅α - ! current: I = ∫dI/ds⋅ds using the trapezoid rule - ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * params%raytracing%dst + ! current: I = ∫dI/ds⋅ds using the trapezoidal rule + ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * dst(jk) tau0(jk) = tau alphaabs0(jk) = alpha @@ -621,7 +622,7 @@ contains alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) call dealloc_pec call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, & - stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & + stnext,stv,dst,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) ! =========== free memory END =========== @@ -629,15 +630,16 @@ contains ! Functions that needs the scope of gray_main - function rhs(y) result(dery) - ! Computes the right-hand-side terms of the ray equations + function rhs(y, e) result(f) + ! Computes the right-hand side terms of the ray equations ! To be passed to the integrator subroutine implicit none ! function arguments - real(wp_), dimension(6), intent(in) :: y ! (x̅, N̅) + real(wp_), intent(in) :: y(6) ! (x̅, N̅) + real(wp_), intent(inout), optional :: e ! |Λ(x̅, N̅)| as an error ! result - real(wp_), dimension(6) :: dery + real(wp_) :: f(6) ! (dx̅/ds, dN̅/ds) ! local variables real(wp_) :: xg, yg @@ -653,8 +655,15 @@ contains xg, yg, derxg, deryg) ! computes derivatives of dispersion relation: ∂Λ/∂x̅, ∂Λ/∂N̅ - call disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, & - gri(:, jk), ggri(:, :, jk), igrad_b, dery) + call disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, & + gri(:, jk), ggri(:, :, jk), igrad_b, dery=f, & + ddr=e) + + ! make the error positive and correct it for an unknown bias: + ! on the correct trajectory |Λ(x̅, N̅)| ≈ -kX, instead of zero. + if (present(e)) then + e = abs(e + 4.15e-4 * xg) + end if end function rhs end subroutine gray_main @@ -1009,66 +1018,209 @@ contains end subroutine ic_gb - - subroutine integrator(y, yp, f, h, method) + function integrator(y, yp, f, h, method) result(y1) ! Integrator of the raytracing equations + implicit none - real(wp_), dimension(6), intent(inout) :: y ! y̅ = (x̅, N̅) - real(wp_), dimension(6), intent(in) :: yp ! y̅˙ = (dx̅/dσ, dN̅/dσ) - procedure(rhs_function) :: f ! dy̅/dσ = f̅(y̅) - real(wp_), intent(in) :: h ! step size - integer, intent(in) :: method + real(wp_), dimension(6), intent(in) :: y ! y̅ = (x̅, N̅) + real(wp_), dimension(6), intent(in) :: yp ! y̅˙ = f(y̅) + procedure(rhs_function) :: f ! dy̅/dσ = f̅(y̅) + real(wp_), intent(in) :: h ! step size + integer, intent(in) :: method ! kind of integrator + real(wp_), dimension(6) :: y1 ! the new y̅ ! local variables - real(wp_), dimension(6) :: yy, k1, k2, k3, k4 + real(wp_), dimension(6) :: k1, k2, k3, k4 select case (method) case (0) ! Explicit Euler (1⁰ order) - y = y + yp*h + y1 = y + yp*h case (1) ! Semi-implicit Euler (1⁰ order, symplectic) - ! P = p - ∂H/∂q(q, p) - ! Q = q + ∂H/∂p(q, P) - k1 = yp - y(4:6) = y(4:6) + k1(4:6)*h - k2 = f(y) - y(1:3) = y(1:3) + k2(1:3)*h + ! P = p - ∂H/∂q(q, p)⋅h + ! Q = q + ∂H/∂p(q, P)⋅h + k1 = h*yp + y1(1:3) = y(1:3) + y1(4:6) = y(4:6) + k1(4:6) + k2 = h*f(y1) + y1(1:3) = y1(1:3) + k2(1:3) case (2) ! Velocity Verlet (2⁰ order, symplectic) ! p(n+½) = p(n) - ∂H/∂q(q(n), p(n))⋅h/2 ! q(n+1) = q(n) + ∂H/∂p(q(n), p(n+½))⋅h ! p(n+1) = p(n+½) - ∂H/∂q(q(n+1), p(n+½))⋅h/2 - k1 = yp - y(4:6) = y(4:6) + k1(4:6)*h/2 - k2 = f(y) - y(1:3) = y(1:3) + k2(1:3)*h - k3 = f(y) - y(4:6) = y(4:6) + k3(4:6)*h/2 + k1 = h*yp + y1(1:3) = y(1:3) + y1(4:6) = y(4:6) + k1(4:6)/2 + k2 = h*f(y1) + y1(1:3) = y(1:3) + k2(1:3) + k3 = h*f(y1) + y1(4:6) = y1(4:6) + k3(4:6)/2 case (3) ! 2-stage Runge-Kutta (2⁰ order) - k1 = yp - yy = y + k1*h - k2 = f(yy) - yy = y + k2*h - y = y + (k1 + k2)*h/2 + k1 = h*yp + k2 = h*f(y + k1/2) + y1 = y + k2 case default ! 4-stage Runge-Kutta (4⁰ order) - k1 = yp - yy = y + k1*h/2 - k2 = f(yy) - yy = y + k2*h/2 - k3 = f(yy) - yy = y + k3*h - k4 = f(yy) - y = y + (k1 + 2*k2 + 2*k3 + k4)*h/6 + k1 = h*yp + k2 = h*f(y + k1/2) + k3 = h*f(y + k2/2) + k4 = h*f(y + k3) + y1 = y + k1/6 + k2/3 + k3/3 + k4/6 end select - end subroutine integrator + end function integrator + + + subroutine step_controller(y, yp, f, h, method, adaptive, Bres) + ! Advances the integration of dy/dσ = f(y) by one step while + ! controlling the error, the latter estimated by the dispersion relation: + ! |Λ(x̅, N̅)| ≠ 0. + use logger, only : log_debug, log_warning + + implicit none + + ! subroutine arguments + real(wp_), dimension(6), intent(inout) :: y ! y̅ = (x̅, N̅) + real(wp_), dimension(6), intent(in) :: yp ! y̅˙ = (dx̅/dσ, dN̅/dσ) + procedure(rhs_function) :: f ! dy̅/dσ = f̅(y̅) + real(wp_), intent(inout) :: h ! step size + integer, intent(in) :: method ! kind of integrator + logical, intent(in) :: adaptive ! whether to change the step + + ! arguments for adaptive control + real(wp_), optional, intent(in) :: Bres ! resonant magnetic field + + ! local variables + real(wp_), dimension(6) :: y1, dummy ! new position, dummy variable + real(wp_) :: e ! error at new position + real(wp_) :: h_max ! max step size + character(256) :: msg + + ! local constants + real(wp_), parameter :: h_min = 1e-2_wp_ ! min step size (cm) + real(wp_), parameter :: e_min = 1e-4_wp_ ! min error + real(wp_), parameter :: e_max = 1e-3_wp_ ! max error + + ! Compute the max step size: this is 1/10 the width of the + ! resonating plasma layer or 5cm otherwise + h_max = max_plasma_step(y(1:3), y(4:6), Bres) + + ! Check if the step could miss the resonance + if (h > h_max) then + if (.not. adaptive) then + write(msg, '("step size is too large: h=", g0.2, " h_max=", g0.2)') h, h_max + call log_warning(msg, mod='gray_core', proc='step_controller') + else + h = h_max + end if + end if + + do + ! Advance by one step + y1 = integrator(y, yp, f, h, method) + + if (.not. adaptive) exit + + ! Compute the error + dummy = f(y1, e) + + ! Try to keep the error bounded to e_max + if (e > e_max .and. h/2 > h_min) then + h = h/2 + write(msg, '("e=", 1pe8.2, ": decreasing step to h=", g0.2)') e, h + call log_debug(msg, mod='gray_core', proc='step_controller') + cycle + end if + + if (e < e_min .and. h*2 < h_max) then + h = h*2 + write(msg, '("e=", 1pe8.2, ": increasing step to h=", g0.2)') e, h + call log_debug(msg, mod='gray_core', proc='step_controller') + end if + + exit + end do + + ! Update the position + y = y1 + + end subroutine step_controller + + + function max_plasma_step(x, N, Bres) result(ds) + ! Takes the position x̅, refractive index N̅, resonant magnetic field + ! Bres and returns the maximum integration step `ds` that can be + ! taken inside the plasma such that it's still possible to resolve + ! the resonance profile well enough. + use equilibrium, only : bfield, equinum_psi + use dispersion, only : resonance_width + use coreprofiles, only : temp + use const_and_precisions, only : mc2=>mc2_ + + implicit none + + ! function arguments + real(wp_), intent(in) :: x(3), N(3) ! position, refractive index + real(wp_), intent(in) :: Bres ! resonant magnetic field + real(wp_) :: ds ! maximum step size + + ! local variables + real(wp_) :: R, dR, z ! cylindrical coordinates + real(wp_) :: B(3), BR, Bphi ! magnetic field components + real(wp_) :: Te, psi ! temperature, poloidal flux + real(wp_) :: Npl ! parallel component of N̅ + + ! To convert from CGS (internal) to SI (equilibrium data) + real(wp_), parameter :: cm = 1e-2_wp_ + + ! Initially assume we are in a vacuum outside the plasma + ! and return ds=5cm as a fallback value + ds = 5 + + ! Compute the local flux and temperature to check + ! whether we are inside the plasma + R = norm2(x(1:2)) + z = x(3) + call equinum_psi(R*cm, z*cm, psi) + + ! No flux data ⇒ outside plasma + if (psi <= 0) return + + ! No temperature data ⇒ outside plasma + Te = temp(psi) + if (Te <= 0) return + + ! Inside the plasma, check for possible harmonics + + ! Compute magnetic field in Cartesian coordinates + call bfield(R*cm, z*cm, Bphi=Bphi, BR=BR, Bz=B(3)) + B(1) = (BR*x(1) - Bphi*x(2)) / R + B(2) = (BR*x(2) + Bphi*x(1)) / R + + Npl = dot_product(N, B)/norm2(B) ! N∥ = N̅⋅b̅ + + ! Compute the extent ΔR of the resonating plasma + ! layer in the major radius R + dR = resonance_width(Y=norm2(B)/Bres, mu=mc2/Te, Npl2=Npl**2, R=R) + + ! No resonance, return the vacuum step + if (dR == 0) return + + ! Compute the extent in the ray direction: + ! + ! ΔR = Δs⋅sinθ ⇒ Δs = ΔR/sinθ = ΔR/√[1 - (N∥/N)²] + ! + ! and divide by a safety factor of 10: + ds = dR / sqrt(1 - (Npl/norm2(N))**2) / 10 + + end function max_plasma_step subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 790d483..c29f1f5 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -59,6 +59,7 @@ module gray_params integer :: nstep ! Max number of integration steps integer :: idst ! Choice of the integration variable integer :: integrator ! Choice of the integration method + logical :: adaptive_step ! Allow variable step sizes integer :: ipass ! Number of plasma passes integer :: ipol ! Whether to compute wave polarisation end type @@ -373,6 +374,7 @@ contains ! Default values of parameters introduced after ! gray_params.data has been deprecated params%raytracing%integrator = 4 + params%raytracing%adaptive_step = .false. end subroutine read_gray_params diff --git a/src/gray_params.sh b/src/gray_params.sh index 3351119..76da5dc 100644 --- a/src/gray_params.sh +++ b/src/gray_params.sh @@ -11,7 +11,7 @@ sets='antenna equilibrium profiles raytracing ecrh_cd output misc' antenna='alpha beta power psi chi iox ibeam filenm fghz pos w ri phi' equilibrium='ssplps ssplf factb sgnb sgni ixp iequil icocos ipsinorm idesc ifreefmt filenm' profiles='psnbnd sspld factne factte iscal irho iprof filenm' -raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator' +raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator adaptive_step' ecrh_cd='iwarm ilarm imx ieccd' output='ipec nrho istpr istpl' misc='rwall' diff --git a/src/multipass.f90 b/src/multipass.f90 index 7df3349..3cba3d5 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -133,24 +133,34 @@ contains end if end subroutine wall_out -! ------------------------------ - subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam) ! initialization at beam propagation start - use logger, only : log_info, log_warning + + + subroutine initbeam(params, i, iroff, iboff, iwait, stv, dst, jphi_beam, & + pins_beam, currins_beam, dpdv_beam, jcd_beam) + ! Initialises the beam variables at the start of the beam propagation + + use gray_params, only : gray_parameters + use logger, only : log_warning implicit none -! arguments - integer, intent(in) :: i ! beam index - logical, dimension(:,:), intent(in), pointer :: iroff ! global ray status (F = active, T = inactive) - logical, intent(out) :: iboff - logical, dimension(:), intent(out), pointer :: iwait - real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam, & - currins_beam,dpdv_beam,jcd_beam,stv - character(256) :: msg ! buffer for formatting log messages - - iboff = .false. ! beam status (F = active, T = inactive) - iwait = iroff(:,i) ! copy ray status for current beam from global ray status - if(all(iwait)) then ! no rays active => stop beam + + ! subroutine arguments + type(gray_parameters), intent(in) :: params + integer, intent(in) :: i ! beam index + logical, pointer, intent(in) :: iroff(:,:) ! global ray status (F = active, T = inactive) + logical, intent(out) :: iboff ! beam status (F = active, T = inactive) + logical, pointer, intent(out) :: iwait(:) + real(wp_), pointer, intent(out), dimension(:) :: & + jphi_beam, pins_beam, currins_beam, dpdv_beam, jcd_beam, stv, dst + + ! local variables + character(256) :: msg ! buffer for formatting log messages + + iboff = .false. + iwait = iroff(:,i) ! copy current beam status from the global one + + if (all(iwait)) then + ! no rays active => stop beam iboff = .true. else if (any(iwait)) then ! only some rays active @@ -158,15 +168,18 @@ contains call log_warning(msg, mod='multipass', proc='initbeam') end if - stv = zero ! starting step - - jphi_beam = zero ! 1D beam profiles + stv = zero ! starting ray parameter (s, c⋅t, S_R) + dst = params%raytracing%dst ! starting step size (ds, c⋅dt, dS_R) + + ! 1D beam profiles + jphi_beam = zero pins_beam = zero currins_beam = zero dpdv_beam = zero jcd_beam = zero end subroutine initbeam -! ------------------------------ + + subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop implicit none @@ -219,71 +232,121 @@ contains end if end subroutine turnoffray -! ------------------------------ - subroutine alloc_multipass(dim,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & - stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) + + + subroutine alloc_multipass(& + dim, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & + stv, dst, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & + pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) + implicit none - integer :: dim - logical, dimension(:), intent(out), pointer :: iwait - logical, dimension(:,:), intent(out), pointer :: iroff - integer, dimension(:), intent(out), pointer :: iop,iow - real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam,stv,tau1,etau1,cpl1,lgcpl1,p0ray,psipv,chipv - real(wp_), dimension(:,:), intent(out), pointer :: taus,cpls,stnext,yw0,ypw0 - real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext + + ! subroutine arguments + integer, intent(in) :: dim + logical, pointer, dimension(:), intent(out) :: iwait + logical, pointer, dimension(:,:), intent(out) :: iroff + integer, pointer, dimension(:), intent(out) :: iop, iow + real(wp_), pointer, dimension(:), intent(out) :: & + jphi_beam, pins_beam, currins_beam, dpdv_beam, & + jcd_beam, stv, dst, tau1, etau1, cpl1, lgcpl1, & + p0ray, psipv, chipv + real(wp_), pointer, dimension(:, :), intent(out) :: taus, cpls, stnext + real(wp_), pointer, dimension(:, :), intent(out) :: yw0, ypw0 + real(wp_), pointer, dimension(:, :, :), intent(out) :: yynext, yypnext - call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext,stv, & - p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam,psipv,chipv) + call dealloc_multipass(& + iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & + stv, dst, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & + pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) + - nbeam_max = 2**ipass ! max n of beams active at a time - nbeam_tot = 2**(ipass+1)-2 ! total n of beams + nbeam_max = 2**ipass ! max n of beams active at a time + nbeam_tot = 2**(ipass+1)-2 ! total n of beams - allocate(iwait(nray),iroff(nray,nbeam_tot),iop(nray),iow(nray), & - yynext(6,nray,nbeam_max-2),yypnext(6,nray,nbeam_max-2), & - yw0(6,nray),ypw0(6,nray),stnext(nray,nbeam_tot),stv(nray), & - p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), & - cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), & - pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_beam(dim), & - psipv(nbeam_tot),chipv(nbeam_tot)) - + allocate(iwait(nray)) + allocate(iroff(nray, nbeam_tot)) + allocate(iop(nray)) + allocate(iow(nray)) + + allocate(yynext(6, nray, nbeam_max-2)) + allocate(yypnext(6, nray, nbeam_max-2)) + + allocate(yw0(6, nray)) + allocate(ypw0(6, nray)) + + allocate(stnext(nray, nbeam_tot)) + allocate(stv(nray)) + allocate(dst(nray)) + + allocate(p0ray(nray)) + allocate(taus(nray, nbeam_tot)) + allocate(tau1(nray)) + allocate(etau1(nray)) + + allocate(cpls(nray, nbeam_tot)) + allocate(cpl1(nray)) + allocate(lgcpl1(nray)) + allocate(jphi_beam(dim)) + + allocate(pins_beam(dim)) + allocate(currins_beam(dim)) + allocate(dpdv_beam(dim)) + allocate(jcd_beam(dim)) + + allocate(psipv(nbeam_tot)) + allocate(chipv(nbeam_tot)) + end subroutine alloc_multipass -! ------------------------------ - subroutine dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & - stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) + + + subroutine dealloc_multipass(& + iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & + stv, dst, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & + pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) + implicit none - logical, dimension(:), intent(out), pointer :: iwait - logical, dimension(:,:), intent(out), pointer :: iroff - integer, dimension(:), intent(out), pointer :: iop,iow - real(wp_), dimension(:), intent(out), pointer :: stv,p0ray,tau1,etau1,cpl1,lgcpl1, & - jphi_beam,pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv - real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls - real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext + + logical, pointer, dimension(:), intent(out) :: iwait + logical, pointer, dimension(:,:), intent(out) :: iroff + integer, pointer, dimension(:), intent(out) :: iop, iow + real(wp_), pointer, dimension(:), intent(out) :: & + jphi_beam, pins_beam, currins_beam, dpdv_beam, & + jcd_beam, stv, dst, tau1, etau1, cpl1, lgcpl1, & + p0ray, psipv, chipv + real(wp_), pointer, dimension(:, :), intent(out) :: taus, cpls, stnext + real(wp_), pointer, dimension(:, :), intent(out) :: yw0, ypw0 + real(wp_), pointer, dimension(:, :, :), intent(out) :: yynext, yypnext if (associated(iwait)) deallocate(iwait) if (associated(iroff)) deallocate(iroff) if (associated(iop)) deallocate(iop) if (associated(iow)) deallocate(iow) + if (associated(yynext)) deallocate(yynext) if (associated(yypnext)) deallocate(yypnext) + if (associated(yw0)) deallocate(yw0) if (associated(ypw0)) deallocate(ypw0) + if (associated(stnext)) deallocate(stnext) if (associated(stv)) deallocate(stv) + if (associated(dst)) deallocate(dst) + if (associated(p0ray)) deallocate(p0ray) if (associated(taus)) deallocate(taus) if (associated(tau1)) deallocate(tau1) if (associated(etau1)) deallocate(etau1) + if (associated(cpls)) deallocate(cpls) if (associated(cpl1)) deallocate(cpl1) if (associated(lgcpl1)) deallocate(lgcpl1) if (associated(jphi_beam)) deallocate(jphi_beam) + if (associated(pins_beam)) deallocate(pins_beam) if (associated(currins_beam)) deallocate(currins_beam) if (associated(dpdv_beam)) deallocate(dpdv_beam) if (associated(jcd_beam)) deallocate(jcd_beam) + if (associated(psipv)) deallocate(psipv) if (associated(chipv)) deallocate(chipv) end subroutine dealloc_multipass