src/gray_core.f90: make integrator configurable
This change makes the integration method of the raytracing equations configurable and significantly simplifies the integrator subroutine by moving the implementation details outside.
This commit is contained in:
Michele Guerini Rocco
parent
1b814fcb8a
commit
c9c20198a7
GPG Key ID:
BFBAF4C975F76450
@ -26,6 +26,14 @@ nstep = 12000
|
||||
; 2: real part of the eikonal (S_r)
|
||||
idst = 0
|
||||
|
||||
; Choice of the integration method
|
||||
; 0: Explicit Euler (1⁰ order)
|
||||
; 1: Semi-implicit Euler (1⁰ order, symplectic)
|
||||
; 2: Velocity Verlet (2⁰ order, symplectic)
|
||||
; 3: 2-stage Runge-Kutta (2⁰ order)
|
||||
; 4: 4-stage Runge-Kutta (4⁰ order)
|
||||
integrator = 4
|
||||
|
||||
|
||||
[ecrh_cd]
|
||||
; Choice of the power absorption model
|
||||
|
||||
@ -3,14 +3,14 @@ module beamdata
|
||||
implicit none
|
||||
|
||||
integer, save :: nray,nrayr,nrayth,nstep,jkray1
|
||||
real(wp_), save :: dst,h,hh,h6,rwmax,twodr2
|
||||
real(wp_), save :: dst,rwmax,twodr2
|
||||
|
||||
contains
|
||||
|
||||
subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
|
||||
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
|
||||
use gray_params, only : raytracing_parameters
|
||||
use const_and_precisions, only : half,two
|
||||
use const_and_precisions, only : two
|
||||
implicit none
|
||||
type(raytracing_parameters), intent(inout) :: rtrparam
|
||||
real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, &
|
||||
@ -24,10 +24,6 @@ contains
|
||||
integer :: jray1
|
||||
|
||||
dst = rtrparam%dst
|
||||
h = dst
|
||||
hh = h*half
|
||||
h6 = h/6.0_wp_
|
||||
|
||||
nrayr = rtrparam%nrayr
|
||||
nrayth = rtrparam%nrayth
|
||||
rwmax = rtrparam%rwmax
|
||||
|
||||
@ -5,6 +5,15 @@ module gray_core
|
||||
|
||||
implicit none
|
||||
|
||||
abstract interface
|
||||
function rhs_function(y) result(f)
|
||||
! Function passed to the integrator subroutine
|
||||
use const_and_precisions, only : wp_
|
||||
real(wp_), intent(in) :: y(6)
|
||||
real(wp_) :: f(6)
|
||||
end function
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine gray_main(params, data, results, error, rhout)
|
||||
@ -265,7 +274,9 @@ 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(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk),igrad_b)
|
||||
call integrator(yw(:, jk), ypw(:, jk), rhs, &
|
||||
params%raytracing%dst, &
|
||||
params%raytracing%integrator)
|
||||
end do
|
||||
! update position and grad
|
||||
if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
|
||||
@ -613,6 +624,39 @@ contains
|
||||
stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
|
||||
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
|
||||
! =========== free memory END ===========
|
||||
|
||||
contains
|
||||
|
||||
! Functions that needs the scope of gray_main
|
||||
|
||||
function rhs(y) result(dery)
|
||||
! 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̅)
|
||||
! result
|
||||
real(wp_), dimension(6) :: dery
|
||||
|
||||
! local variables
|
||||
real(wp_) :: xg, yg
|
||||
real(wp_), dimension(3) :: xv, anv, bv, derxg, deryg
|
||||
real(wp_), dimension(3, 3) :: derbv
|
||||
|
||||
xv = y(1:3) ! x̅ position
|
||||
anv = y(4:6) ! N̅ refractive index
|
||||
|
||||
! computes derivatives of plasma quantities: B̅, ∇X, ∇Y, ∇B̅
|
||||
! (these are needed for the next part)
|
||||
call plas_deriv(xv, bres, xgcn, dens, btot, bv, derbv, &
|
||||
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)
|
||||
end function rhs
|
||||
|
||||
end subroutine gray_main
|
||||
|
||||
|
||||
@ -966,55 +1010,65 @@ contains
|
||||
|
||||
|
||||
|
||||
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
|
||||
! Runge-Kutta integrator
|
||||
! use gray_params, only : igrad
|
||||
use beamdata, only : h,hh,h6
|
||||
subroutine integrator(y, yp, f, h, method)
|
||||
! Integrator of the raytracing equations
|
||||
implicit none
|
||||
real(wp_), intent(in) :: bres,xgcn
|
||||
real(wp_), dimension(6), intent(inout) :: y
|
||||
real(wp_), dimension(6), intent(in) :: yp
|
||||
real(wp_), dimension(3), intent(in) :: dgr
|
||||
real(wp_), dimension(3,3), intent(in) :: ddgr
|
||||
integer, intent(in) :: igrad,sox
|
||||
|
||||
real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
|
||||
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
|
||||
|
||||
fk1 = yp
|
||||
yy = y + fk1*hh
|
||||
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk2,igrad)
|
||||
yy = y + fk2*hh
|
||||
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk3,igrad)
|
||||
yy = y + fk3*h
|
||||
call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk4,igrad)
|
||||
! local variables
|
||||
real(wp_), dimension(6) :: yy, k1, k2, k3, k4
|
||||
|
||||
y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
|
||||
end subroutine rkstep
|
||||
select case (method)
|
||||
case (0)
|
||||
! Explicit Euler (1⁰ order)
|
||||
y = 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
|
||||
|
||||
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
|
||||
|
||||
subroutine rhs(sox,bres,xgcn,y,dgr,ddgr,dery,igrad)
|
||||
! Compute right-hand side terms of the ray equations (dery)
|
||||
! used in R-K integrator
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), dimension(6), intent(in) :: y
|
||||
real(wp_), intent(in) :: bres,xgcn
|
||||
real(wp_), dimension(3), intent(in) :: dgr
|
||||
real(wp_), dimension(3,3), intent(in) :: ddgr
|
||||
real(wp_), dimension(6), intent(out) :: dery
|
||||
integer, intent(in) :: igrad,sox
|
||||
! local variables
|
||||
real(wp_) :: dens,btot,xg,yg
|
||||
real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg
|
||||
real(wp_), dimension(3,3) :: derbv
|
||||
|
||||
xv = y(1:3)
|
||||
call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg)
|
||||
anv = y(4:6)
|
||||
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad,dery)
|
||||
end subroutine rhs
|
||||
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
|
||||
|
||||
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
|
||||
end select
|
||||
end subroutine integrator
|
||||
|
||||
|
||||
subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
|
||||
|
||||
@ -5,6 +5,9 @@ module gray_params
|
||||
integer, parameter :: lenfnm = 256
|
||||
integer, parameter :: headw = 132, headl = 21
|
||||
|
||||
! Note: when adding/removing a parameter remember to keep
|
||||
! the gray_params.sh script in sync with this file
|
||||
|
||||
! Antenna/wave launcher parameters
|
||||
type antenna_parameters
|
||||
! From gray_params.data:
|
||||
@ -55,6 +58,7 @@ module gray_params
|
||||
integer :: igrad ! Complex eikonal switch
|
||||
integer :: nstep ! Max number of integration steps
|
||||
integer :: idst ! Choice of the integration variable
|
||||
integer :: integrator ! Choice of the integration method
|
||||
integer :: ipass ! Number of plasma passes
|
||||
integer :: ipol ! Whether to compute wave polarisation
|
||||
end type
|
||||
@ -365,6 +369,11 @@ contains
|
||||
read(u, *) params%output%istpr, params%output%istpl
|
||||
|
||||
close(u)
|
||||
|
||||
! Default values of parameters introduced after
|
||||
! gray_params.data has been deprecated
|
||||
params%raytracing%integrator = 4
|
||||
|
||||
end subroutine read_gray_params
|
||||
|
||||
|
||||
|
||||
@ -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'
|
||||
raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator'
|
||||
ecrh_cd='iwarm ilarm imx ieccd'
|
||||
output='ipec nrho istpr istpl'
|
||||
misc='rwall'
|
||||
|
||||
Loading…
Reference in New Issue
Block a user