From 15fc89110aa61d0146545f0d3c68f0b97e15f140 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 2 May 2024 00:47:59 +0200 Subject: [PATCH] abstract the outputs units This change replaces the output files (Fortran units) with a derived type called table, that hold the data in memory until further processing. The data stored in a table can be dumped to a file, as before, or processed in other ways, for example converted to other derived type. --- src/gray_cli.f90 | 32 +- src/gray_core.f90 | 813 ++++------------------------------------ src/gray_jetto1beam.f90 | 24 +- src/gray_params.f90 | 135 ++++--- src/gray_tables.f90 | 612 ++++++++++++++++++++++++++++++ src/magsurf_data.f90 | 89 ++--- src/main.f90 | 278 +++++++------- src/types.f90 | 354 +++++++++++++++++ src/units.f90 | 70 ---- src/utils.f90 | 16 + src/vendor/cniteq.f90 | 218 +++++++++++ 11 files changed, 1526 insertions(+), 1115 deletions(-) create mode 100644 src/gray_tables.f90 create mode 100644 src/types.f90 delete mode 100644 src/units.f90 create mode 100644 src/vendor/cniteq.f90 diff --git a/src/gray_cli.f90 b/src/gray_cli.f90 index 7288cfb..497d1c0 100644 --- a/src/gray_cli.f90 +++ b/src/gray_cli.f90 @@ -19,7 +19,7 @@ module gray_cli character(len=:), allocatable :: sum_filelist ! others integer :: verbose - integer, allocatable :: units(:) + integer, allocatable :: tables(:) end type private @@ -54,8 +54,8 @@ contains print '(a)', ' -c, --config-file FILE set the (new) GRAY config file' print '(a)', ' (default: gray.ini)' print '(a)', ' -s, --sum FILE sum the output profiles from a list of files' - print '(a)', ' -u, --units ID[,ID...] select which units to output (default: 4, 7);' - print '(a)', ' use `all` to enable all, or `none` for no units.' + print '(a)', ' -t, --tables ID[,ID...] select which tables to write (default: 4, 7);' + print '(a)', ' use `all` to enable all, or `none` for no tables.' print '(a)', ' see the manual for all unit IDs.' print '(a)', ' -g, --gray-param ID=VAL set a GRAY parameter, overriding the value' print '(a)', ' specified via --params-file/--config-file;' @@ -94,7 +94,7 @@ contains print '(a,a)' , ' param-file: ' , opts%params_file print '(a,a)' , ' sum: ' , opts%sum_filelist print '(a)' , 'others:' - print '(a,20i3)' , ' - units: ' , opts%units + print '(a,20i3)' , ' - tables: ' , opts%tables print '(a,l)' , ' - verbose: ' , opts%verbose end subroutine @@ -102,7 +102,6 @@ contains subroutine parse_cli_options(opts) ! Parse the CLI arguments and initialise the options - use units, only : ucenr, usumm, all_enabled use logger, only : WARNING ! subroutine arguments @@ -117,7 +116,7 @@ contains opts%verbose = WARNING opts%quiet = .false. opts%params_file = 'gray_params.data' - opts%units = [ucenr, usumm] + opts%tables = [4, 7] nargs = command_argument_count() i = 1 @@ -156,24 +155,25 @@ contains case ('-s', '--sum') call get_next_command(i, opts%sum_filelist) - case ('-u', '--units') + case ('-t', '--tables') call get_next_command(i, temp) if (temp == 'none') then - ! disable all output units - deallocate(opts%units) - allocate(opts%units(0)) + ! disable all output tables + deallocate(opts%tables) + allocate(opts%tables(0)) elseif (temp == 'all') then - ! enable all output units - all_enabled = .true. + ! enable all output tables + deallocate(opts%tables) + opts%tables = [-1] else ! resize the array commas = count([(temp(i:i) == ',', i = 1, len(temp))]) - deallocate(opts%units) - allocate(opts%units(commas + 1)) + deallocate(opts%tables) + allocate(opts%tables(commas + 1)) ! read the list of table IDs - read (temp, *, iostat=error) opts%units + read (temp, *, iostat=error) opts%tables if (error > 0) then print '(a,a)', 'invalid table IDs: ', temp deallocate(argument) @@ -296,7 +296,7 @@ contains if (allocated(opts%params_file)) deallocate(opts%params_file) if (allocated(opts%config_file)) deallocate(opts%config_file) if (allocated(opts%sum_filelist)) deallocate(opts%sum_filelist) - if (allocated(opts%units)) deallocate(opts%units) + if (allocated(opts%tables)) deallocate(opts%tables) end subroutine deinit_cli_options end module gray_cli diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 7989a32..7620a63 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -8,16 +8,19 @@ module gray_core contains subroutine gray_main(params, data, results, error, rhout) - use const_and_precisions, only : zero, one, degree, comp_tiny + use const_and_precisions, only : zero, one, comp_tiny use coreprofiles, only : temp, fzeff use dispersion, only : expinit use polarization, only : ellipse_to_field - use gray_params, only : gray_parameters, gray_data, gray_results, & - print_parameters, EQ_VACUUM + use types, only : table, wrap + use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM + use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & + store_beam_shape, find_flux_surfaces, & + kinetic_profiles, ec_resonance, inputs_maps + use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight, rayi2jk - use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd - use magsurf_data, only : flux_average, dealloc_surfvec + use magsurf_data, only : compute_flux_averages, dealloc_surfvec use beamdata, only : init_btr use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab @@ -63,12 +66,12 @@ contains real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg - associate ( nray => params%raytracing%nray, & - nrayr => params%raytracing%nrayr, & - nrayth => params%raytracing%nrayth, & - nstep => params%raytracing%nstep, & - nbeam_tot => 2**(params%raytracing%ipass+1)-2, & - nbeam_max => 2**(params%raytracing%ipass)) + associate (nray => params%raytracing%nray, & + nrayr => params%raytracing%nrayr, & + nrayth => params%raytracing%nrayth, & + nstep => params%raytracing%nstep, & + nbeam_tot => 2**(params%raytracing%ipass+1)-2, & + nbeam_max => 2**(params%raytracing%ipass)) block ! ray variables @@ -106,6 +109,9 @@ contains ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) + ! Initialise the output tables + call init_tables(params, results%tables) + ! Compute the initial cartesian wavevector (anv0) ! from launch angles α,β and the position call launchangles2n(params%antenna, anv0) @@ -118,7 +124,7 @@ contains if (params%equilibrium%iequil /= EQ_VACUUM) then ! Initialise the magsurf_data module - call flux_average ! requires frhotor for dadrhot,dvdrhot + call compute_flux_averages(results%tables) ! Initialise the output profiles call pec_init(params%output%ipec, rhout) @@ -136,11 +142,13 @@ contains results%jcd = zero ! ========= set environment END ========= - ! ======== pre-proc prints BEGIN ======== - call print_headers(params) + ! Pre-determinted tables + results%tables%kinetic_profiles = kinetic_profiles(params) + results%tables%ec_resonance = ec_resonance(params, bres) + results%tables%inputs_maps = inputs_maps(params, bres, xgcn) - ! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop - call print_surfq([1.5_wp_, 2.0_wp_]) + ! print ψ surface for q=3/2 and q=2/1 + call find_flux_surfaces([1.5_wp_, 2.0_wp_], results%tables%flux_surfaces) ! print initial position write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos @@ -150,17 +158,6 @@ contains 'α', params%antenna%alpha, 'β', params%antenna%beta call log_info(msg, mod='gray_core', proc='gray_main') - ! print Btot=Bres - ! print ne, Te, q, Jphi versus psi, rhop, rhot - if (params%equilibrium%iequil /= EQ_VACUUM) then - call print_bres(bres) - call print_prof(params%profiles) - call print_maps(bres, xgcn, & - norm2(params%antenna%pos(1:2)) * 0.01_wp_, & - sin(params%antenna%beta*degree)) - end if - ! ========= pre-proc prints END ========= - ! =========== main loop BEGIN =========== call initmultipass(params%raytracing%ipol, params%antenna%iox, & iroff,yynext,yypnext,yw0,ypw0, & @@ -229,7 +226,8 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri, index_rt) ! * initial conditions + call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, & + gri, ggri, index_rt, results%tables) ! * initial conditions do jk=1,params%raytracing%nray zzm = yw(3,jk)*0.01_wp_ @@ -259,7 +257,7 @@ contains iop = 0 ! start propagation outside plasma if(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8 - call print_projxyzt(stv,yw,0) + call store_beam_shape(params%raytracing, results%tables, stv, yw) ! ======= propagation loop BEGIN ======= call log_debug(' propagation loop start', mod='gray_core', proc='gray_main') @@ -487,9 +485,10 @@ contains ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1) psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1) else - call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, & - btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & - nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) + call store_ray_data(results%tables, & + i, jk, stv(jk), p0ray(jk), xv, psinv, btot, bv, ak0, & + anpl, anpr, anv, anprim, dens, tekev, alpha, tau, dids, & + nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) end if end do rays_loop @@ -504,7 +503,7 @@ contains ! print ray positions for j=nrayr in local reference system if(mod(i,params%output%istpr) == 0) then if(params%raytracing%nray > 1 .and. all(.not.iwait)) & - call print_projxyzt(stv,yw,0) + call store_beam_shape(params%raytracing, results%tables, stv, yw) end if ! test whether further trajectory integration is unnecessary @@ -522,7 +521,9 @@ contains ! ======== propagation loop END ======== ! print all ray positions in local reference system - if(params%raytracing%nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1) + if(params%raytracing%nray > 1 .and. all(.not.iwait)) & + call store_beam_shape(params%raytracing, results%tables, & + stv, yw, full=.true.) ! =========== post-proc BEGIN =========== ! compute total absorbed power and driven current for current beam @@ -581,17 +582,24 @@ contains end if if (params%equilibrium%iequil /= EQ_VACUUM) then - call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam, & - dpdv_beam,currins_beam,pins_beam,ip) ! *print power and current density profiles for current beam + call store_ec_profiles( & + results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, & + jcd_beam, dpdv_beam, currins_beam, pins_beam, ip) call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, & jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam - call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav, & - rhotjava,drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj, & - drhotp,drhotj,ratjamx,ratjbmx,stv(1),psipv(index_rt), & - chipv(index_rt),index_rt,sum(p0ray),cpl_beam1,cpl_beam2) ! *print 0D results for current beam + if (results%tables%summary%active) then + call results%tables%summary%append([ & + wrap(icd_beam), wrap(pabs_beam), wrap(jphip), wrap(dpdvp), & + wrap(rhotj), wrap(rhotjava), wrap(rhotp), wrap(rhotpav), & + wrap(drhotjava), wrap(drhotpav), wrap(ratjamx), wrap(ratjbmx), & + wrap(stv(1)), wrap(psipv(index_rt)), wrap(chipv(index_rt)), & + wrap(index_rt), wrap(jphimx), wrap(dpdvmx), wrap(drhotj), & + wrap(drhotp), wrap(sum(p0ray)), wrap(cpl_beam1), wrap(cpl_beam2)]) ! *print 0D results for current beam + end if + end if ! ============ post-proc END ============ @@ -621,10 +629,9 @@ contains call log_debug('pass loop end', mod='gray_core', proc='gray_main') ! ============ main loop END ============ - ! ========== free memory BEGIN ========== + ! Free memory call dealloc_surfvec call dealloc_pec - ! =========== free memory END =========== end block end associate @@ -663,12 +670,15 @@ contains subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, & - stv, xc0, du10, gri, ggri, index_rt) + stv, xc0, du10, gri, ggri, index_rt, & + tables) ! beam tracing initial conditions igrad=1 ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im - use gray_params, only : gray_parameters + use gray_params, only : gray_parameters, gray_tables use beamdata, only : nray,nrayr,nrayth,rwmax + use types, only : table + use gray_tables, only : store_ray_data ! subroutine arguments type(gray_parameters), intent(in) :: params @@ -681,6 +691,9 @@ contains real(wp_), dimension(3, 3, nray), intent(out) :: ggri integer, intent(in) :: index_rt + ! tables for storing initial rays data + type(gray_tables), intent(inout), optional :: tables + ! local variables real(wp_), dimension(3) :: xv0c real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir @@ -973,10 +986,18 @@ contains ddr = anx**2 + any**2 + anz**2 - an20 ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) - ! i=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0 - call print_output(0,jk,stv(jk),one,xc0(:,k,j),-one,zero,[zero,zero,zero], & - ak0,zero,zero,[zero,zero,zero],zero,zero,zero,zero,zero,zero, & - 0,0,0,index_rt,ddr,ddi,zero,zero,[zero,zero,zero]) + + ! save step "zero" data + if (present(tables)) & + call store_ray_data(tables, & + i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), & + psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & + N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, & + n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, & + nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, & + lambda_r=ddr, lambda_i=ddi, X=zero, Y=zero, & + grad_X=[zero,zero,zero]) + end do end subroutine ic_gb @@ -1848,700 +1869,4 @@ contains end subroutine alpha_effj - - subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon) -! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. -! (based on an older code) -! arguments - integer, intent(in) :: nr,nz - real(wp_), dimension(nr), intent(in) :: rqgrid - real(wp_), dimension(nz), intent(in) :: zqgrid - real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid - real(wp_), intent(in) :: h - integer, intent(inout) :: ncon, icount - integer, dimension(ncon), intent(out) :: npts - real(wp_), dimension(icount), intent(out) :: rcon,zcon -! local variables - integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb - integer :: jabs,jnb,kx,ikx,itm,inext,in - integer, dimension(3,2) :: ja - integer, dimension(icount/2-1) :: lx - real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y - real(wp_), dimension(nr*nz) :: a - logical :: flag1 - - px = 0.5_wp_ - - a = reshape(matr2dgrid, [nr*nz]) - - rcon = 0.0_wp_ - zcon = 0.0_wp_ - - nrqmax = nr - drgrd = rqgrid(2) - rqgrid(1) - dzgrd = zqgrid(2) - zqgrid(1) - - ncon = 0 - - npts = 0 - - iclast = 0 - icount = 0 - mpl = 0 - ix = 0 - mxr = nrqmax * (nz - 1) - n1 = nr - 1 - - do jx=2,n1 - do jm=jx,mxr,nrqmax - j = jm + nrqmax - ah=a(j)-h - if (ah <= 0.0_wp_ .and. a(jm) > h .or. & - ah > 0.0_wp_ .and. a(jm) <= h) then - ix=ix+1 - lx(ix)=-j - end if - if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & - ah > 0.0_wp_ .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - end do - end do - - do jm=nr,mxr,nrqmax - j = jm + nrqmax - ah=a(j)-h - if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & - ah > 0.0_wp_ .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - if (ah <= 0.0_wp_ .and. a(jm) > h .or. & - ah > 0.0_wp_ .and. a(jm) <= h) then - ix=ix+1 - lx(ix)=-j - end if - end do - - do jm=1,mxr,nrqmax - j = jm + nrqmax - if (a(j) <= h .and. a(jm) > h .or. & - a(j) > h .and. a(jm) <= h) then - ix=ix+1 - lx(ix) =-j - end if - end do - - do j=2,nr - if (a(j) <= h .and. a(j-1) > h .or. & - a(j) > h .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - end do - - if(ix<=0) return - -bb: do - in=ix - jx=lx(in) - jfor=0 - lda=1 - ldb=2 - - do - if(jx<0) then - jabs=-jx - jnb = jabs - nrqmax - else - jabs=jx - jnb=jabs-1 - end if - - adn=a(jabs)-a(jnb) - if(adn/=0) px=(a(jabs)-h)/adn - kx = (jabs - 1) / nrqmax - ikx = jabs - nrqmax * kx - 1 - - if(jx<0) then - x = drgrd * ikx - y = dzgrd * (kx - px) - else - x = drgrd * (ikx - px) - y = dzgrd * kx - end if - - icount = icount + 1 - rcon(icount) = x + rqgrid(1) - zcon(icount) = y + zqgrid(1) - mpl= icount - itm = 1 - ja(1,1) = jabs + nrqmax - j=1 - - if(jx<=0) then - ja(1,1) = -jabs-1 - j=2 - end if - - ja(2,1) = -ja(1,1) - ja(3,1) = -jx + 1 - nrqmax - ja(3,2) = -jx - ja(j,2) = jabs - nrqmax - k= 3-j - ja(k,2) = 1-jabs - - if (kx<=0 .or. ikx<=0) then - lda=1 - ldb=lda - else if (ikx + 1 - nr >= 0 .and. jx <= 0) then - lda=2 - ldb=lda - else if(jfor/=0) then - lda=2 - do i=1,3 - if(jfor==ja(i,2)) then - lda=1 - exit - end if - end do - ldb=lda - end if - - flag1=.false. - aa: do k=1,3 - do l=lda,ldb - do i=1,ix - if(lx(i)==ja(k,l)) then - itm=itm+1 - inext= i - if(jfor/=0) exit aa - if(itm > 3) then - flag1=.true. - exit aa - end if - end if - end do - end do - end do aa - - if(.not.flag1) then - lx(in)=0 - if(itm == 1) exit - end if - - jfor=jx - jx=lx(inext) - in = inext - end do - - do - if(lx(ix)/=0) then - if(mpl>=4) then - ncon = ncon + 1 - npts(ncon) = icount - iclast - iclast = icount - end if - exit - end if - ix= ix-1 - if(ix<=0) exit bb - end do - - end do bb - - if(mpl >= 4) then - ncon = ncon + 1 - npts(ncon) = icount - iclast - iclast = icount - end if - end subroutine cniteq - - - - subroutine print_headers(params) - ! Prints the headers of all gray_main output tables - - use units, only : uprj0, uwbm, udisp, ucenr, uoutr, upec, usumm, & - unit_active, active_units - use gray_params, only : gray_parameters, headw, headl, print_parameters - - ! subroutine arguments - type(gray_parameters), intent(in) :: params - - ! local variables - integer :: i, j - integer :: main_units(8) - character(256) :: main_headers(8) - character(len=headw), dimension(headl) :: header - - main_units = [uprj0, uprj0+1, uwbm, udisp, ucenr, uoutr, upec, usumm] - - main_headers(1) = '#sst j k xt yt zt rt' - main_headers(2) = '#sst j k xt yt zt rt' - main_headers(3) = '#sst w1 w2' - main_headers(4) = '#sst Dr_Nr Di_Nr' - main_headers(5) = '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp ' & - // 'Npl Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax ' & - // 'iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz' - main_headers(6) = '#i k sst x y R z psin tau Npl alpha index_rt' - main_headers(7) = '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt' - main_headers(8) = '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' & - // 'drhotjava drhotpav ratjamx ratjbmx stmx psipol ' & - // 'chipol index_rt Jphimx dPdVmx drhotj drhotp P0 ' & - // 'cplO cplX' - - if (all(active_units == 0)) return - - call print_parameters(params, header) - - do i=1,size(main_units) - if (unit_active(main_units(i))) then - do j=1,size(header) - write (main_units(i), '(1x,a)') header(j) - end do - write (main_units(i), '(1x,a)') trim(main_headers(i)) - end if - end do - - end subroutine print_headers - - - subroutine print_prof(params) - ! Prints the (input) plasma kinetic profiles (unit 55) - - use gray_params, only : profiles_parameters - use equilibrium, only : fq, frhotor - use coreprofiles, only : density, temp - use units, only : uprfin, unit_active - use magsurf_data, only : npsi, vajphiav - - ! suborutine arguments - type(profiles_parameters), intent(in) :: params - - ! local variables - integer :: i, ntail - real(wp_) :: rho_t, J_phi, dens, ddens - real(wp_) :: psi_n, rho_p, drho_p - - if (.not. unit_active(uprfin)) 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) - - ! extra points to reach ψ=ψ_bnd - ntail = ceiling((sqrt(params%psnbnd) - 1) / drho_p) - - write (uprfin, *) '#psi rhot ne Te q Jphi' - do i = 0, npsi + ntail - rho_p = i * drho_p - rho_t = frhotor(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 density(psi_n, dens, ddens) - write (uprfin, '(12(1x,e12.5))') & - psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi - end do - end subroutine print_prof - - - subroutine print_bres(bres) - ! Prints the EC resonance surface table (unit 70) - - use const_and_precisions, only : comp_eps - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield - use units, only : ubres, unit_active - use magsurf_data, only : npsi - - ! subroutine arguments - real(wp_), intent(in) :: bres - - ! local constants - integer, parameter :: icmx = 2002 - - ! local variables - integer :: j,k,n,nconts,nctot - integer, dimension(10) :: ncpts - real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb - real(wp_), dimension(icmx) :: rrcb,zzcb - real(wp_) :: rv(npsi), zv(npsi) - real(wp_), dimension(npsi,npsi) :: btotal - - if (.not. unit_active(ubres)) return - - ! Build a regular (R, z) grid - dr = (rmxm - rmnm - comp_eps)/(npsi - 1) - dz = (zmxm - zmnm)/(npsi - 1) - do j=1,npsi - rv(j) = comp_eps + rmnm + dr*(j - 1) - zv(j) = zmnm + dz*(j - 1) - end do - - ! Btotal on psi grid - btmx = -1.0e30_wp_ - btmn = 1.0e30_wp_ - do k = 1, npsi - zzk = zv(k) - do j = 1, npsi - rrj = rv(j) - call bfield(rrj, zzk, bbr, bbz, bbphi) - btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2) - if(btotal(j,k) >= btmx) btmx = btotal(j,k) - if(btotal(j,k) <= btmn) btmn = btotal(j,k) - end do - end do - - ! compute Btot=Bres/n with n=1,5 - write (ubres, *) '#i Btot R z' - do n = 1, 5 - bbb = bres/n - if (bbb >= btmn .and. bbb <= btmx) then - nconts = size(ncpts) - nctot = size(rrcb) - call cniteq(rv, zv, btotal, npsi, npsi, bbb, & - nconts, ncpts, nctot, rrcb, zzcb) - do j = 1, nctot - write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j) - end do - end if - write (ubres, *) - end do - - end subroutine print_bres - - - subroutine print_maps(B_res, xgcn, R0, Npl0) - ! Prints several input quantities on a regular grid (unit 72) - - use const_and_precisions, only : comp_eps - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield - use coreprofiles, only : density, temp - use units, only : umaps, unit_active - use magsurf_data, only : npsi - - ! subroutine arguments - real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω - real(wp_), intent(in) :: xgcn ! X normalisation, e²/ε₀m_eω² - real(wp_), intent(in) :: R0 ! initial value of R - real(wp_), intent(in) :: Npl0 ! initial value of N∥ - - ! local variables - integer :: j, k - 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 - - if (.not. unit_active(umaps)) return - - ! Build a regular (R, z) grid - dR = (rmxm - rmnm - comp_eps)/(npsi - 1) - dz = (zmxm - zmnm)/(npsi - 1) - do j = 1, npsi - R(j) = comp_eps + rmnm + dR*(j - 1) - z(j) = zmnm + dz*(j - 1) - end do - - write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl' - do j = 1, npsi - Npl = Npl0 * R0/r(j) - 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) - B = sqrt(B_R**2 + B_phi**2 + B_z**2) - X = xgcn*ne - Y = B/B_res - Te = temp(psi_n) - write (umaps,'(12(x,e12.5))') & - R(j), z(k), psi_n, B_R, B_phi, B_z, B, ne, Te, X, Y, Npl - end do - write (umaps,*) - end do - - end subroutine print_maps - - - subroutine print_surfq(qvals) - ! Prints the constant ψ surfaces for a set of values of q - - use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf - use magsurf_data, only : contours_psi, npoints, print_contour - use logger, only : log_info - use minpack, only : hybrj1 - - ! subroutine arguments - real(wp_), intent(in) :: qvals(:) - - ! local variables - integer :: i - real(wp_) :: rho_t, rho_p, psi_n - character(256) :: msg ! for log messages formatting - - call log_info('constant ψ surfaces for:', & - mod='gray_core', proc='print_surfq') - do i = 1, size(qvals) - ! Find value of ψ_n for the current q - block - real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7) - integer :: info - - ! Note: since we are mostly interested in q=3/2,2 we start - ! searching near the boundary in case q is not monotonic. - sol = [0.8_wp_] ! first guess - - ! Solve fq(ψ_n) = qvals(i) for ψ_n - call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, & - ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7) - ! Solution - psi_n = sol(1) - - ! Handle spurious or no solutions - if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle - end block - - ! Compute and print the ψ_n(R,z) = ψ_n₀ contour - block - real(wp_), dimension(npoints) :: R_cont, z_cont - real(wp_) :: R_hi, R_lo, z_hi, z_lo - - ! Guesses for the high,low horzizontal-tangent points - R_hi = rmaxis; - z_hi = (zbsup + zmaxis)/2 - R_lo = rmaxis - z_lo = (zbinf + zmaxis)/2 - - call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) - call print_contour(psi_n, R_cont, z_cont) - end block - - ! Log the values found for ψ_n, ρ_p, ρ_t - rho_p = sqrt(psi_n) - rho_t = frhotor(rho_p) - write (msg, '(4(x,a,"=",g0.3))') & - 'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t - call log_info(msg, mod='gray_core', proc='print_surfq') - end do - - contains - - subroutine equation(n, x, f, df, ldf, flag) - ! The equation to solve: f(x) = q(x) - q₀ = 0 - use const_and_precisions, only : comp_eps - - ! subroutine arguments - integer, intent(in) :: n, ldf, flag - real(wp_), intent(in) :: x(n) - real(wp_), intent(inout) :: f(n), df(ldf,n) - - ! optimal step size - real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) - - if (flag == 1) then - ! return f(x) - f(1) = fq(x(1)) - qvals(i) - else - ! return f'(x), computed numerically - df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e) - end if - end subroutine - - end subroutine print_surfq - - - subroutine print_projxyzt(stv, ywrk, iproj) - ! Prints the beam shape tables (unit 8, 12) - - use const_and_precisions, only : comp_huge, zero, one - use beamdata, only : nray, nrayr, nrayth, rayi2jk - use units, only : uprj0, uwbm, unit_active - - ! subroutine arguments - real(wp_), dimension(:), intent(in) :: stv - real(wp_), dimension(:,:), intent(in) :: ywrk - integer, intent(in) :: iproj - - ! local variables - integer :: jk, jkz, uprj - integer, dimension(2) ::jkv - real(wp_), dimension(3) :: xv1, dir, dxv - real(wp_) :: rtimn, rtimx, csth1, snth1, & - csps1, snps1, xti, yti, zti, rti - - if (.not.(unit_active(uprj0) .or. unit_active(uwbm))) return - - ! common/external functions/variables - - uprj = uprj0 + iproj - - xv1 = ywrk(1:3, 1) - dir = ywrk(4:6, 1) - dir = dir/norm2(dir) - csth1 = dir(3) - snth1 = sqrt(one - csth1**2) - if(snth1 > zero) then - csps1 = dir(2)/snth1 - snps1 = dir(1)/snth1 - else - csps1 = one - snps1 = zero - end if - - if(iproj==0) then - jkz = nray - nrayth + 1 - else - jkz = 1 - end if - - rtimn = comp_huge - rtimx = zero - do jk = jkz, nray - dxv = ywrk(1:3,jk) - xv1 - xti = dxv(1)*csps1 - dxv(2)*snps1 - yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1 - zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1 - rti = hypot(xti, yti) - - jkv=rayi2jk(jk) - if(.not.(iproj==0 .and. jk==1)) & - write(uprj, '(1x,e16.8e3,2i5,4(1x,e16.8e3))') & - stv(jk), jkv, xti, yti, zti, rti - if(iproj==1 .and. jkv(2)==nrayth) write (uprj, *) - - if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti - if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti - end do - write (uprj, *) - write (uwbm, '(3(1x,e16.8e3))') stv(1), rtimn, rtimx - - end subroutine print_projxyzt - - - subroutine print_output(i, jk, st, qj, xv, psinv, btot, bv, ak0, & - anpl, anpr, anv, anprim, dens, tekev, & - alpha, tau, dids, nhm, nhf, iokhawa, & - index_rt, ddr, ddi, xg, yg, derxg) - ! Prints the ray variables and Hamiltonian error (units 4, 33, 17) - - use const_and_precisions, only : degree,zero,one - use equilibrium, only : frhotor - use gray_params, only : istpl0 - use beamdata, only : nray,nrayth,jkray1 - use units, only : ucenr, uoutr, udisp, unit_active - - ! subroutine arguments - integer, intent(in) :: i, jk, nhm, nhf, iokhawa, index_rt - real(wp_), dimension(3), intent(in) :: xv, bv, anv - real(wp_), intent(in) :: st, qj, psinv, btot, ak0, anpl, anpr, anprim - real(wp_), intent(in) :: dens, tekev, alpha, tau, dids, ddr, ddi - real(wp_), intent(in) :: xg, yg - real(wp_), dimension(3), intent(in) :: derxg - - ! local variables - real(wp_) :: stm, xxm, yym, zzm, rrm, phideg, rhot, akim, pt, didsn - integer :: k - - stm = st * 1.0e-2_wp_ - xxm = xv(1) * 1.0e-2_wp_ - yym = xv(2) * 1.0e-2_wp_ - zzm = xv(3) * 1.0e-2_wp_ - rrm = hypot(xxm, yym) - - ! print central ray trajectory. dIds in A/m/W, ki in m⁻¹ - if(jk == 1) then - phideg = atan2(yym, xxm)/degree - if(psinv>=zero .and. psinv<=one) then - rhot = frhotor(sqrt(psinv)) - else - rhot = 1.0_wp_ - end if - akim = anprim * ak0 * 1.0e2_wp_ - pt = qj*exp(-tau) - didsn = dids * 1.0e2_wp_/qj - - ! print central ray data - if (unit_active(ucenr)) then - write (ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') & - stm, rrm, zzm, phideg, psinv, rhot, dens, tekev, & - btot, bv, anpr, anpl, anv, akim, alpha, tau, pt, & - didsn, nhm, nhf, iokhawa, index_rt, ddr, xg, yg, & - derxg - end if - end if - - ! print conservation of dispersion relation (Hamiltonian error) - if(unit_active(udisp) .and. jk==nray) then - write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi - end if - - ! print outer ray trajectories - if(mod(i, istpl0) == 0) then - k = jk - jkray1 + 1 - if(unit_active(uoutr) .and. k>0 .and. k<=nrayth) then - write(uoutr, '(2i5,9(1x,e16.8e3),i5)') & - i, k, stm, xxm, yym, rrm, zzm, & - psinv, tau, anpl, alpha, index_rt - end if - end if - end subroutine print_output - - - subroutine print_pec(rhop_tab, rhot_tab, jphi, jcd, & - dpdv, currins, pins, index_rt) - ! Prints the ECRH and CD profiles (unit 48) - - use units, only : upec, unit_active - - ! subroutine arguments - real(wp_), dimension(:), intent(in) :: rhop_tab, rhot_tab, jphi, jcd, & - dpdv, currins, pins - integer, intent(in) :: index_rt - - ! local variables - integer :: i - - if (.not. unit_active(upec)) return - - do i=1, size(rhop_tab) - write(upec, '(7(1x, e16.8e3), i5)') rhop_tab(i), rhot_tab(i), & - jphi(i), jcd(i), dpdv(i), currins(i), pins(i), index_rt - end do - write(upec, *) '' - end subroutine print_pec - - - subroutine print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, & - drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, & - drhotj, ratjamx, ratjbmx, stmx, psipol, chipol, & - index_rt, p0, cpl1, cpl2) - ! Prints global data and final results (unit 7) - - use units, only : ucenr, usumm, unit_active - - ! subroutine arguments - real(wp_), intent(in) :: pabs, icd, dpdvp, jphip, rhotpav, rhotjava, & - drhotpav, drhotjava, dpdvmx, jphimx, rhotp, & - rhotj, drhotp, drhotj, ratjamx, ratjbmx, & - stmx, psipol, chipol, p0, cpl1, cpl2 - integer, intent(in) :: index_rt - - if (unit_active(ucenr)) then - write(ucenr, *) '' - end if - - if (unit_active(usumm)) then - write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') & - icd, pabs, jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & - drhotjava, drhotpav, ratjamx, ratjbmx, stmx, psipol, & - chipol, index_rt, jphimx, dpdvmx, drhotj, drhotp, & - p0, cpl1, cpl2 - end if - end subroutine print_finals - end module gray_core diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index 5fec21a..682aa06 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -2,7 +2,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, & p0mw, alphain, betain, dpdv, jcd, pabs, icd, err) use const_and_precisions, only: wp_ - use units, only: set_active_units, close_units use gray_params, only: gray_parameters, gray_data, gray_results use gray_core, only: gray_main @@ -29,16 +28,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & first_call: if (firstcall) then firstcall = .false. - ! Activate some debug units - set_units: block - use units - call set_active_units([ & - ucenr, usumm, uprj0, uprj0+1, & - uwbm, udisp, ubres, ucnt, & - uoutr, umaps, uprfin, uflx, & - upec]) - end block set_units - ! Read parameters from external file init_params: block use gray_params, only : read_gray_params, set_globals @@ -47,6 +36,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Override some parameters params%misc%rwall = r(1) + params%misc%active_tables = [4, 7, 8, 9, 12, 17, 33, 70, 71, 55, 48] params%antenna%filenm = 'graybeam.data' params%equilibrium%filenm = 'JETTO' params%equilibrium%iequil = ijetto + 1 @@ -134,6 +124,16 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Call main subroutine for the ibeam-th beam call gray_main(params, data, res, err, rhout=sqrt(psrad)) + write_debug_files: block + integer :: i, err + ! Write results to file + do i = 1, size(res%tables%iter) + associate (tbl => res%tables%iter(i)%ptr) + call tbl%save(err, filepath='gray_'//trim(tbl%title)//'.txt') + end associate + end do + end block write_debug_files + ! Free memory free_memory: block use equilibrium, only : unset_equil_spline @@ -152,6 +152,4 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & dpdv = res%dpdv jcd = res%jcd - call close_units - end subroutine gray_jetto1beam diff --git a/src/gray_params.f90 b/src/gray_params.f90 index ea76ff6..349c47f 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -1,6 +1,7 @@ module gray_params use const_and_precisions, only : wp_ + use types, only : table implicit none @@ -157,13 +158,14 @@ module gray_params type output_parameters integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles integer :: nrho = 501 ! Number of grid steps for EC profiles + 1 - integer :: istpr = 5 ! Subsampling factor for beam cross-section (units 8, 12) - integer :: istpl = 5 ! Subsampling factor for outer rays (unit 33) + integer :: istpr = 5 ! Subsampling factor for beam cross-section (tables 8, 12) + integer :: istpl = 5 ! Subsampling factor for outer rays (table 33) end type ! Other parameters type misc_parameters - real(wp_) :: rwall ! R of the limiter (fallback) + real(wp_) :: rwall ! R of the limiter (fallback) + integer, allocatable :: active_tables(:) ! IDs of output tables to fill in end type ! MHD equilibrium data @@ -207,19 +209,19 @@ module gray_params integer, parameter :: max_name_size = 40 character(len=*), parameter :: mandatory(*) = [ & character(len=max_name_size) :: & - "antenna.alpha", & - "antenna.beta", & - "antenna.iox", & - "antenna.ibeam", & - "antenna.filenm", & - "equilibrium.filenm", & - "profiles.irho", & - "profiles.iprof", & - "profiles.filenm", & - "raytracing.dst", & - "raytracing.nrayr", & - "raytracing.nrayth", & - "misc.rwall" & + "antenna.alpha", & + "antenna.beta", & + "antenna.iox", & + "antenna.ibeam", & + "antenna.filenm", & + "equilibrium.filenm", & + "profiles.irho", & + "profiles.iprof", & + "profiles.filenm", & + "raytracing.dst", & + "raytracing.nrayr", & + "raytracing.nrayth", & + "misc.rwall" & ] ! All GRAY input data @@ -228,14 +230,33 @@ module gray_params type(profiles_data) :: profiles end type - ! GRAY final results (only some) - type gray_results - real(wp_) :: pabs ! Total absorbed power - real(wp_) :: icd ! Total driven current - real(wp_), allocatable :: dpdv(:) ! Absorbed power density - real(wp_), allocatable :: jcd(:) ! Driven current density + ! Wrapper type for array of pointers + type table_ptr + type(table), pointer :: ptr => null() end type + ! Extra outputs tables + type gray_tables + type(table) :: flux_averages, flux_surfaces + type(table) :: summary, ec_profiles + type(table) :: central_ray, outer_rays, dispersion + type(table) :: beam_shape, beam_final, beam_size + type(table) :: kinetic_profiles, ec_resonance, inputs_maps + + ! for iterating over the tables + type(table_ptr) :: iter(13) + end type + + ! GRAY final results + type gray_results + real(wp_) :: pabs ! Total absorbed power + real(wp_) :: icd ! Total driven current + real(wp_), allocatable :: dpdv(:) ! Absorbed power density + real(wp_), allocatable :: jcd(:) ! Driven current density + type(gray_tables) :: tables ! Extra outputs + end type + + ! Global variables exposed for gray_core integer, save :: iequil, iprof integer, save :: iwarm, ilarm, imx, ieccd @@ -245,11 +266,12 @@ module gray_params contains - subroutine print_parameters(params, strout) + subroutine write_gray_header(params, file) + ! Writes the Gray output header to a file ! subroutine arguments type(gray_parameters), intent(in) :: params - character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21) + integer :: file ! local variables character(len=8) :: date @@ -260,91 +282,90 @@ contains ! Date and time call date_and_time(date, time) - write(strout(1), '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') & + write(file, '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') & date(1:4), date(5:6), date(7:8), & time(1:2), time(3:4), time(5:10) ! Git revision - write(strout(2), '("# GRAY Git revision: ",a)') REVISION + write(file, '("# GRAY Git revision: ",a)') REVISION ! MHD equilibrium parameters if (params%equilibrium%iequil > 0) then - write(strout(3), '("# EQL input: ",a)') trim(params%equilibrium%filenm(1:95)) + write(file, '("# EQL input: ",a)') trim(params%equilibrium%filenm) ! TODO add missing values - write(strout(7), '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') & - 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_ + write(file, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') & + 0, 0, 0, 0, 0 else - write(strout(3), '("# EQL input: N/A (vacuum)")') - write(strout(7), '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') + write(file, '("# EQL input: N/A (vacuum)")') + write(file, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') end if - write(strout(4), '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') & + write(file, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') & params%equilibrium%iequil, params%equilibrium%sgnb, & params%equilibrium%sgni, params%equilibrium%factb if (params%equilibrium%iequil > 1) then - write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') & + write(file, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') & params%equilibrium%icocos, params%equilibrium%ipsinorm, & params%equilibrium%idesc, params%equilibrium%ifreefmt - write(strout(6), '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') & + write(file, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') & params%equilibrium%ssplps, params%equilibrium%ssplf, & params%equilibrium%ixp else - write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")') - write(strout(6), '("# EQL ssplps ssplf ixp: N/A (analytical)")') + write(file, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")') + write(file, '("# EQL ssplps ssplf ixp: N/A (analytical)")') end if ! Profiles parameters if (params%equilibrium%iequil > 0) then - write(strout(8), '("# PRF input: ",a)') trim(params%profiles%filenm(1:95)) - write(strout(9), '("# PRF iprof iscal factne factte:",4(1x,g0.4))') & + write(file, '("# PRF input: ",a)') trim(params%profiles%filenm) + write(file, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') & params%profiles%iprof, params%profiles%iscal, & params%profiles%factne, params%profiles%factte if (params%profiles%iprof > 0) then - write(strout(10), '("# PRF irho psnbnd sspld:",3(1x,g0.3))') & + write(file, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') & params%profiles%irho, params%profiles%psnbnd, params%profiles%sspld else - write(strout(10), '("# PRF irho psnbnd sspld: N/A (analytical)")') + write(file, '("# PRF irho psnbnd sspld: N/A (analytical)")') end if ! TODO: add missing values - write(strout(11), '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') & - 0._wp_, 0._wp_, 0._wp_ + write(file, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0 else - write(strout(8), '("# PRF input: N/A (vacuum)")') - write(strout(9), '("# PRF iprof iscal factne factte: N/A (vacuum)")') - write(strout(10), '("# PRF irho psnbnd sspld: N/A (vacuum)")') - write(strout(11), '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') + write(file, '("# PRF input: N/A (vacuum)")') + write(file, '("# PRF iprof iscal factne factte: N/A (vacuum)")') + write(file, '("# PRF irho psnbnd sspld: N/A (vacuum)")') + write(file, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') end if ! Antenna parameters - write(strout(12), '("# ANT input: ",a)') trim(params%antenna%filenm(1:95)) - write(strout(13), '("# ANT ibeam iox psi chi:",4(1x,g0.4))') & + write(file, '("# ANT input: ",a)') trim(params%antenna%filenm) + write(file, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') & params%antenna%ibeam, params%antenna%iox, & params%antenna%psi, params%antenna%chi - write(strout(14), '("# ANT alpha beta power fghz:",4(1x,g0.3))') & + write(file, '("# ANT alpha beta power fghz:",4(1x,g0.3))') & params%antenna%alpha, params%antenna%beta, & params%antenna%power, params%antenna%fghz - write(strout(15), '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos - write(strout(16), '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') & + write(file, '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos + write(file, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') & params%antenna%w, params%antenna%ri, params%antenna%phi ! Other parameters - write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall + write(file, '("# RFL rwall:",1x,e12.5)') params%misc%rwall ! code parameters - write(strout(18), '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') & + write(file, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') & params%raytracing%igrad, params%raytracing%idst, & params%raytracing%ipass, params%raytracing%ipol - write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') & + write(file, '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') & params%raytracing%nrayr, params%raytracing%nrayth, & params%raytracing%nstep, params%raytracing%rwmax, & params%raytracing%dst - write(strout(20), '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') & + write(file, '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') & params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, & params%ecrh_cd%imx, params%ecrh_cd%ieccd - write(strout(21), '("# COD ipec nrho istpr istpl:",4(1x,i4))') & + write(file, '("# COD ipec nrho istpr istpl:",4(1x,i4))') & params%output%ipec, params%output%nrho, & params%output%istpr, params%output%istpl - end subroutine print_parameters + end subroutine write_gray_header function update_parameter(params, name, value) result(err) diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 new file mode 100644 index 0000000..4e1046f --- /dev/null +++ b/src/gray_tables.f90 @@ -0,0 +1,612 @@ +! This module contains the gray ouput tables and routines to operate on them +module gray_tables + + use const_and_precisions, only : wp_ + use types, only : queue, table, wrap + + implicit none + + private + public init_tables ! Initialising main tables + public find_flux_surfaces, store_flux_surface ! Filling the main tables + public store_ray_data, store_ec_profiles ! + public store_beam_shape ! + public kinetic_profiles, ec_resonance, inputs_maps ! Extra tables + +contains + + pure function is_active(params, id) + ! Helper to check whether a table is active + use gray_params, only : gray_parameters + + ! function arguments + type(gray_parameters), intent(in) :: params + integer, intent(in) :: id + logical :: is_active + + associate (ids => params%misc%active_tables) + is_active = any(ids == id) + if (size(ids) == 1) then + ! special case, all enabled + is_active = is_active .or. (ids(1) == -1) + end if + end associate + end function + + + subroutine init_tables(params, tbls) + ! Initialises the gray_main output tables + + use gray_params, only : gray_parameters, gray_tables, ptr=>table_ptr + + ! subroutine arguments + type(gray_parameters), intent(in) :: params + type(gray_tables), target, intent(out) :: tbls + + call tbls%flux_averages%init( & + title='flux-averages', id=56, active=is_active(params, 56), & + labels=[character(64) :: 'ρ_p', 'ρ_t', 'B_avg', & + 'B_max', 'B_min', 'area', 'vol', 'I_pl', & + 'J_φ_avg', 'fc', 'ratio_Ja', 'ratio_Jb', 'q']) + + call tbls%flux_surfaces%init( & + title='flux-surfaces', id=71, active=is_active(params, 71), & + labels=[character(64) :: 'i', 'ψ_n', 'R', 'z']) + + call tbls%ec_profiles%init( & + title='ec-profiles', id=48, active=is_active(params, 48), & + labels=[character(64) :: 'ρ_p', 'ρ_t', 'J_φ', & + 'J_cdb', 'dPdV', 'I_cd_inside', & + 'P_inside', 'index_rt']) + + call tbls%summary%init( & + title='summary', id=7, active=is_active(params, 7), & + labels=[character(64) :: 'I_cd', 'P_abs', 'J_φ_peak', & + 'dPdV_peak', 'ρ_max_J', 'ρ_avg_J', 'ρ_max_P', & + 'ρ_avg_P', 'Δρ_avg_J', 'Δρ_avg_P', & + 'ratio_Ja_max', 'ratio_Jb_max', 's_max', & + 'ψ', 'χ', 'index_rt', 'J_φ_max', 'dPdV_max', & + 'Δρ_J', 'Δρ_P', 'P0', 'cpl1', 'cpl2']) + + call tbls%central_ray%init( & + title='central-ray', id=4, active=is_active(params, 4), & + labels=[character(64) :: 's', 'R', 'z', 'φ', & + 'ψ_n', 'ρ_t', 'n_e', 'T_e', 'B', 'b_x', & + 'b_y', 'b_z', 'N_⊥', 'N_∥', 'N_x', 'N_y', & + 'N_z', 'k_im', 'α', 'τ', 'P_t', 'dI/ds', & + 'n_harm_min', 'n_harm_max', 'i_okhawa', & + 'index_rt', 'Λ_r', 'X', 'Y', & + '∇X_x', '∇X_y', '∇X_z']) + + call tbls%outer_rays%init( & + title='outer-rays', id=33, active=is_active(params, 33), & + labels=[character(64) :: 'i', 'k', 's', 'x', 'y', & + 'R', 'z', 'ψ_n', 'τ', 'N_∥', 'α', 'index_rt']) + + call tbls%dispersion%init( & + title='dispersion', id=17, active=is_active(params, 17), & + labels=[character(64) :: 's', 'Λ_r', 'Λ_i']) + + call tbls%beam_shape%init( & + title='beam-shape', id=8, active=is_active(params, 8), & + labels=['s', 'j', 'k', 'x', 'y', 'z', 'r']) + + call tbls%beam_final%init( & + title='beam-shape-final', id=9, active=is_active(params, 9), & + labels=['s', 'j', 'k', 'x', 'y', 'z', 'r']) + + call tbls%beam_size%init( & + title='beam-size', id=12, active=is_active(params, 12), & + labels=[character(64) :: 's', 'r_min', 'r_max']) + + ! List of pointers to the tables for iterating over + tbls%iter = [ & + ptr(tbls%flux_averages), ptr(tbls%flux_surfaces), & + ptr(tbls%summary), ptr(tbls%ec_profiles), & + ptr(tbls%central_ray), ptr(tbls%outer_rays), ptr(tbls%dispersion), & + ptr(tbls%beam_shape), ptr(tbls%beam_final), ptr(tbls%beam_size), & + ptr(tbls%kinetic_profiles), ptr(tbls%ec_resonance), ptr(tbls%inputs_maps)] + + end subroutine init_tables + + + function kinetic_profiles(params) result(tbl) + ! Generates the plasma kinetic profiles + + use gray_params, only : gray_parameters, EQ_VACUUM + 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 + + ! local variables + integer :: i, ntail + + real(wp_) :: rho_t, J_phi, dens, ddens + real(wp_) :: psi_n, rho_p, drho_p + + call tbl%init(title='kinetic-profiles', id=55, active=is_active(params, 25), & + labels=[character(64) :: 'psi_n', 'rho_t', 'n_e', 'T_e', 'q', 'J_φ']) + + if (.not. tbl%active) return + 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) + + ! extra points to reach ψ=ψ_bnd + ntail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p) + + do i = 0, npsi + ntail + rho_p = i * drho_p + rho_t = frhotor(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 density(psi_n, dens, ddens) + call tbl%append([psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi]) + end do + + end function kinetic_profiles + + + function ec_resonance(params, B_res) result(tbl) + ! Generates the EC resonance surface table + + use const_and_precisions, only : comp_eps + use gray_params, only : gray_parameters, EQ_VACUUM + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield + use magsurf_data, only : npsi + use types, only : item + use cniteq, only : cniteq_f + + ! function arguments + type(gray_parameters), intent(in) :: params + real(wp_), intent(in) :: B_res + type(table) :: tbl + + call tbl%init(title='ec-resonance', id=70, active=is_active(params, 70), & + labels=[character(64) :: 'i', 'B_tot', 'R', 'z']) + + if (.not. tbl%active) return + if (params%equilibrium%iequil == EQ_VACUUM) return + + block + ! local variables + 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 + + ! Build a regular (R, z) grid + dr = (rmxm - rmnm - comp_eps)/(npsi - 1) + dz = (zmxm - zmnm)/(npsi - 1) + do j=1,npsi + R(j) = comp_eps + rmnm + dr*(j - 1) + z(j) = zmnm + 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 bfield(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) + end do + end do + + ! compute Btot=Bres/n with n=1,5 + do n = 1, 5 + B = B_res/n + if (B >= B_min .and. B <= B_max) then + n_conts = size(n_points) + n_tot = size(R_cont) + call cniteq_f(R, z, B_tot, size(R), size(z), B, & + n_conts, n_points, n_tot, R_cont, z_cont) + do j = 1, size(R_cont) + call tbl%append([wrap(j), wrap(B), wrap(R_cont(j)), wrap(z_cont(j))]) + end do + end if + end do + end block + + end function ec_resonance + + + function inputs_maps(params, 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 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 + + ! local variables + 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 + + call tbl%init(title='inputs-maps', id=72, active=is_active(params, 72), & + labels=[character(64) :: 'R', 'z', 'ψ_n', 'B_R', 'B_z', 'B', & + 'n_e', 'T_e', 'X', 'Y', 'N_∥']) + + if (.not. tbl%active) return + if (params%equilibrium%iequil == EQ_VACUUM) return + + R0 = norm2(params%antenna%pos(1:2)) * cm ! initial value of R + Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥ + + ! Build a regular (R, z) grid + dR = (rmxm - rmnm - comp_eps)/(npsi - 1) + dz = (zmxm - zmnm)/(npsi - 1) + do j = 1, npsi + R(j) = comp_eps + rmnm + dR*(j - 1) + z(j) = zmnm + dz*(j - 1) + end do + + do j = 1, npsi + Npl = Npl0 * R0/r(j) + 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) + B = sqrt(B_R**2 + B_phi**2 + B_z**2) + X = X_norm*ne + Y = B/B_res + Te = 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 + end do + + end function inputs_maps + + + subroutine find_flux_surfaces(qvals, tbl) + ! Finds the ψ for a set of values of q and stores the + ! associated surfaces to the flux surface table + + use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf + use magsurf_data, only : contours_psi, npoints + use logger, only : log_info + use minpack, only : hybrj1 + use types, only : table + + ! subroutine arguments + real(wp_), intent(in) :: qvals(:) + type(table), intent(inout) :: tbl + + ! local variables + integer :: i + real(wp_) :: rho_t, rho_p, psi_n + character(256) :: msg ! for log messages formatting + + call log_info('constant ψ surfaces for:', & + mod='gray_tables', proc='find_flux_surfaces') + do i = 1, size(qvals) + ! Find value of ψ_n for the current q + block + real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7) + integer :: info + + ! Note: since we are mostly interested in q=3/2,2 we start + ! searching near the boundary in case q is not monotonic. + sol = [0.8_wp_] ! first guess + + ! Solve fq(ψ_n) = qvals(i) for ψ_n + call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, & + ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7) + ! Solution + psi_n = sol(1) + + ! Handle spurious or no solutions + if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle + end block + + ! Compute and print the ψ_n(R,z) = ψ_n₀ contour + block + real(wp_), dimension(npoints) :: R_cont, z_cont + real(wp_) :: R_hi, R_lo, z_hi, z_lo + + ! Guesses for the high,low horzizontal-tangent points + R_hi = rmaxis; + z_hi = (zbsup + zmaxis)/2 + R_lo = rmaxis + z_lo = (zbinf + zmaxis)/2 + + call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo) + call store_flux_surface(tbl, psi_n, R_cont, z_cont) + end block + + ! Log the values found for ψ_n, ρ_p, ρ_t + rho_p = sqrt(psi_n) + rho_t = frhotor(rho_p) + write (msg, '(4(x,a,"=",g0.3))') & + 'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t + call log_info(msg, mod='gray_tables', proc='find_flux_surfaces') + end do + + contains + + subroutine equation(n, x, f, df, ldf, flag) + ! The equation to solve: f(x) = q(x) - q₀ = 0 + use const_and_precisions, only : comp_eps + + ! subroutine arguments + integer, intent(in) :: n, ldf, flag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: f(n), df(ldf,n) + + ! optimal step size + real(wp_), parameter :: e = comp_eps**(1/3.0_wp_) + + if (flag == 1) then + ! return f(x) + f(1) = fq(x(1)) - qvals(i) + else + ! return f'(x), computed numerically + df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e) + end if + end subroutine + + end subroutine find_flux_surfaces + + + subroutine store_flux_surface(tbl, psi_n, R, z) + ! Helper to add a contour to the flux surfaces table + + use const_and_precisions, only : wp_, comp_tiny + use types, only : table, wrap + + ! subroutine arguments + real(wp_), intent(in) :: psi_n + real(wp_), intent(in) :: R(:) + real(wp_), intent(in) :: z(:) + type(table), intent(inout) :: tbl + + ! local variables + integer :: i + + if (.not. tbl%active) return + if (psi_n < comp_tiny) return + + do i = 1, size(R) + call tbl%append([wrap(i), wrap(psi_n), wrap(R(i)), wrap(z(i))]) + end do + end subroutine store_flux_surface + + + subroutine store_ray_data(tables, & + i, jk, s, P0, pos, psi_n, B, b_n, k0, & + N_pl, N_pr, N, N_pr_im, n_e, T_e, & + alpha, tau, dIds, nhm, nhf, iokhawa, & + index_rt, lambda_r, lambda_i, X, Y, grad_X) + ! Stores some ray variables and local quantities + + use const_and_precisions, only : degree, cm + use equilibrium, only : frhotor + use gray_params, only : gray_tables, istpl0 + use beamdata, only : nray,nrayth,jkray1 + + ! subroutine arguments + + ! tables where to store the data + type(gray_tables), intent(inout) :: tables + + integer, intent(in) :: i, jk ! step, ray index + real(wp_), intent(in) :: s ! arclength + real(wp_), intent(in) :: P0 ! initial power + real(wp_), intent(in) :: pos(3) ! position vector + real(wp_), intent(in) :: B, b_n(3) ! magnetic field and unit vector + real(wp_), intent(in) :: N(3) ! refractive index vector, N + real(wp_), intent(in) :: psi_n ! normalised poloidal flux, ψ_n + real(wp_), intent(in) :: N_pl, N_pr ! N_∥, N_⊥ + real(wp_), intent(in) :: N_pr_im ! Im(N_⊥) from warm dispersion + real(wp_), intent(in) :: k0 ! vacuum wavenumber (k₀=ω/c) + real(wp_), intent(in) :: n_e, T_e ! electron density and temperature + real(wp_), intent(in) :: alpha, tau ! absorption coeff., optical depth + real(wp_), intent(in) :: dIds ! ECCD current density + integer, intent(in) :: nhm, nhf ! EC harmonic numbers + integer, intent(in) :: iokhawa ! + integer, intent(in) :: index_rt ! + real(wp_), intent(in) :: lambda_r, lambda_i ! quasioptical dispersion relation + real(wp_), intent(in) :: X, Y ! X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: grad_X(3) ! gradient of X + + ! local variables + real(wp_) :: s_m, pos_m(3), R_m, phi_deg ! coordinates in SI units + real(wp_) :: rho_t, k_im, Pt, dIds_n + integer :: k + + s_m = s * cm + pos_m = pos * cm + R_m = hypot(pos_m(1), pos_m(2)) + + ! Store central ray data + if (jk == 1 .and. tables%central_ray%active) then + phi_deg = atan2(pos_m(2), pos_m(1)) / degree + if(psi_n >= 0 .and. psi_n <= 1) then + rho_t = frhotor(sqrt(psi_n)) + else + rho_t = 1.0_wp_ + end if + k_im = N_pr_im * k0 / cm + Pt = P0 * exp(-tau) + dIds_n = (1/P0 * dIds) / cm + + call tables%central_ray%append([ & + wrap(s_m), wrap(R_m), wrap(pos_m(3)), wrap(phi_deg), wrap(psi_n), & + wrap(rho_t), wrap(n_e), wrap(T_e), wrap(B), wrap(b_n(1)), & + wrap(b_n(2)), wrap(b_n(3)), wrap(N_pr), wrap(N_pl), wrap(N(1)), & + wrap(N(2)), wrap(N(3)), wrap(k_im), wrap(alpha), wrap(tau), wrap(Pt), & + wrap(dIds_n), wrap(nhm), wrap(nhf), wrap(iokhawa), wrap(index_rt), & + wrap(lambda_r), wrap(X), wrap(Y), wrap(grad_X(1)), wrap(grad_X(2)), & + wrap(grad_X(3))]) + end if + + ! Store outer rays data + if (mod(i, istpl0) == 0) then + k = jk - jkray1 + 1 + if(k > 0 .and. k <= nrayth .and. tables%outer_rays%active) then + call tables%outer_rays%append([ & + wrap(i), wrap(k), wrap(s_m), wrap(pos_m(1)), wrap(pos_m(2)), & + wrap(R_m), wrap(pos_m(3)), wrap(psi_n), wrap(tau), wrap(N_pl), & + wrap(alpha), wrap(index_rt)]) + end if + end if + + ! Store dispersion relation data + if (jk == nray .and. tables%dispersion%active) then + call tables%dispersion%append([s, lambda_r, lambda_i]) + end if + + end subroutine store_ray_data + + + subroutine store_ec_profiles(tbl, rho_p, rho_t, J_phi, J_cd, & + dPdV, I_ins, P_ins, index_rt) + ! Helper to add data to the ECRH&CD profiles table + + use types, only : table, wrap + + ! subroutine arguments + type(table), intent(inout) :: tbl + real(wp_), dimension(:), intent(in) :: rho_p, rho_t, J_phi, J_cd, & + dPdV, I_ins, P_ins + integer, intent(in) :: index_rt + + integer :: i + + if (.not. tbl%active) return + + do i = 1, size(rho_p) + call tbl%append([wrap(rho_p(i)), wrap(rho_t(i)), wrap(J_phi(i)), & + wrap(J_cd(i)), wrap(dPdV(i)), wrap(I_ins(i)), & + wrap(P_ins(i)), wrap(index_rt)]) + end do + end subroutine store_ec_profiles + + + subroutine store_beam_shape(params, tables, s, y, full) + ! Computes the beam shape tables + + use const_and_precisions, only : zero, one, comp_huge + use types, only : table, wrap + use gray_params, only : raytracing_parameters, gray_tables + use beamdata, only : rayi2jk + + ! subroutine arguments + type(raytracing_parameters), intent(in) :: params + type(gray_tables), target, intent(inout) :: tables + real(wp_), intent(in) :: s(:) ! arclength + real(wp_), intent(in) :: y(:, :) ! (x̅, N̅) vector + + ! whether to show all rays or only the outer ring + logical, intent(in), optional :: full + + ! local variables + logical :: full_ + integer :: jk, jk0, jkv(2) + real(wp_) :: n(3), M(3,3), dx(3), dx_loc(3), c + real(wp_) :: r, r_min, r_max + type(table), pointer :: beam_shape + + if (.not. (tables%beam_shape%active .or. tables%beam_size%active)) return + + full_ = .false. + if (present(full)) full_ = full + + if (full_) then + beam_shape => tables%beam_final + else + beam_shape => tables%beam_shape + end if + + ! Convert to the local beam basis + ! + ! If n̅ is the normalised direction of the central ray, + ! and the {e̅₁,e̅₂,e̅₃} the standard basis, the unit vectors of + ! the beam basis are given by: + ! + ! x̅ = n̅×e̅₃ / |n̅×e̅₃| = (n₂/c)e̅₁ + (-n₁/c)e̅₂ + ! y̅ = n̅×(n̅×e̅₃) / |n̅×e̅₃| = (n₁n₃/c)e̅₁ + (n₂n₃/c)e̅₂ - ce̅₃ + ! z̅ = n̅ = n₁e̅₁ + n₂e̅₂ + n₃e̅₃ + ! + ! where c = |n̅×e̅₃| = √(n₁² + n₂²) + ! + ! Or in the case where n̅ ∥ e̅₃ by: + ! + ! x̅ = e̅₁ + ! y̅ = (n₃e̅₃)×e̅₁ = n₃e̅₂ + ! z̅ = n₃e̅₃ + ! + ! So, the change of basis matrix is + ! + ! [n₂/c -n₁/c 0] + ! M = [n₁n₃/c n₂n₃/c -c], c > 0 + ! [n₁ n₂ n₃] + ! + ! or M = diag(1, n₃, n₃) if c = 0. + ! + ! By definition then we have: x̅_loc = M x̅. + ! + n = y(4:6, 1) / norm2(y(4:6, 1)) + c = norm2(n(1:2)) + if (c > 0) then + M = reshape([ n(2)/c, n(1)*n(3)/c, n(1), & + -n(1)/c, n(2)*n(3)/c, n(2), & + zero, -c, n(3) ], [3, 3]) + else + M = reshape([ one, zero, zero, & + zero, n(3), zero, & + zero, zero, n(3) ], [3, 3]) + end if + + ! start from the central ray or the last ring + jk0 = 1 + merge(0, params%nray - params%nrayth, full_) + + r_min = comp_huge + r_max = 0 + do jk = jk0, params%nray + dx = y(1:3, jk) - y(1:3, 1) ! distance from the central ray + dx_loc = matmul(M, dx) ! local ray position + r = norm2(dx_loc(1:2)) ! cross-section radius + jkv = rayi2jk(jk) ! [j, k] vector + + if (full_ .or. jk /= 1) & + call beam_shape%append([ & + wrap(s(jk)), wrap(jkv(1)), wrap(jkv(2)), & + wrap(dx_loc(1)), wrap(dx_loc(2)), wrap(dx_loc(3)), wrap(r)]) + + ! Compute min/max radius + if (r>=r_max .and. jkv(1)==params%nrayr) r_max = r + if (r<=r_min .and. jkv(1)==params%nrayr) r_min = r + end do + + call tables%beam_size%append([s(1), r_min, r_max]) + + end subroutine store_beam_shape + +end module gray_tables diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 0725a97..2c392b3 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -95,11 +95,16 @@ contains end subroutine dealloc_surfvec - subroutine flux_average + subroutine compute_flux_averages(tables) use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & - bfield,frhotor,fq,tor_curr,psia,pol_flux - use dierckx, only : regrid,coeff_parder + use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & + bfield,frhotor,fq,tor_curr,psia,pol_flux + use dierckx, only : regrid,coeff_parder + use types, only : table, wrap + use gray_params, only : gray_tables + + ! subroutine arguments + type(gray_tables), intent(inout), optional :: tables ! local constants integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, & @@ -384,21 +389,29 @@ contains njpt=njp nlmt=nlm + deallocate(rctemp, zctemp) + + if (.not. present(tables)) return do jp=1,npsi - call print_fluxav(rpstab(jp),frhotor(rpstab(jp)),bav(jp),bmxpsi(jp), & - bmnpsi(jp),varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp), & - ffc(jp),vratja(jp),vratjb(jp),qqv(jp)) + if (.not. tables%flux_averages%active) exit + call tables%flux_averages%append([ & + rpstab(jp), frhotor(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 - call print_contour(psicon(jp),rcon(:,jp),zcon(:,jp)) + 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 end do - deallocate(rctemp, zctemp) - end subroutine flux_average - + end subroutine compute_flux_averages subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, & @@ -514,54 +527,4 @@ contains end if end subroutine contours_psi - - subroutine print_contour(psin, rc, zc) - ! Prints the flux surface contours (unit 71) - - use const_and_precisions, only : wp_, comp_tiny - use units, only : ucnt, unit_active - - ! subroutine arguments - real(wp_), intent(in) :: psin - real(wp_), dimension(:), intent(in) :: rc, zc - - ! local variables - integer :: npoints, ic - - if (.not. unit_active(ucnt)) return - - if (psin < comp_tiny) then - write(ucnt, *)' #i psin R z' - else - npoints = size(rc) - do ic=1,npoints - write(ucnt, '(i6,12(1x,e12.5))') ic, psin, rc(ic), zc(ic) - end do - write(ucnt, *) - write(ucnt, *) - end if - end subroutine print_contour - - - subroutine print_fluxav(psin, rhot, bav, bmx, bmn, area, vol, & - currp, ajphiav, ffc, ratja, ratjb, qq) - ! Prints the flux-averaged quantities table (unit 56) - - use const_and_precisions, only : wp_, comp_tiny - use units, only : uflx, unit_active - - ! subroutine arguments - real(wp_), intent(in) :: psin, rhot, bav, bmx, bmn, area, vol, & - currp, ajphiav, ffc, ratja, ratjb, qq - - if (.not. unit_active(uflx)) return - - if (psin < comp_tiny) & - write(uflx, *)' #rhop rhot || |Bmx| |Bmn| Area Vol' // & - ' |I_pl| fc ratJa ratJb qq' - write(uflx, '(20(1x,e12.5))') & - psin, rhot, bav, bmx, bmn, area, vol, currp, ajphiav, & - ffc, ratja, ratjb, qq - end subroutine print_fluxav - end module magsurf_data diff --git a/src/main.f90 b/src/main.f90 index 5a03e2b..5dd4e2f 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,13 +1,12 @@ program main - use const_and_precisions, only : wp_, zero + use const_and_precisions, only : wp_ use logger, only : INFO, ERROR, set_log_level, log_message - use units, only : set_active_units, close_units use utils, only : dirname use gray_cli, only : cli_options, parse_cli_options, & deinit_cli_options, parse_param_overrides use gray_core, only : gray_main use gray_params, only : gray_parameters, gray_data, gray_results, & - read_gray_params, read_gray_config, & + read_gray_params, read_gray_config, & params_set_globals => set_globals implicit none @@ -31,9 +30,6 @@ program main if (opts%quiet) opts%verbose = ERROR call set_log_level(opts%verbose) - ! Activate the given output units - call set_active_units(opts%units) - ! Load the parameters from file and move to its directory ! (all other filepaths are assumed relative to it) if (allocated(opts%config_file)) then @@ -79,6 +75,9 @@ program main call init_misc(params, data) + ! Activate the given output tables + params%misc%active_tables = opts%tables + ! Get back to the initial directory err = chdir(cwd) @@ -93,93 +92,7 @@ program main if (allocated(opts%sum_filelist)) then call log_message(level=INFO, mod='main', msg='summing profiles') - - sum: block - real(wp_), dimension(:), allocatable :: jphi - real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin - real(wp_), dimension(:), allocatable :: extracol - integer :: i, j, k, l, n, ngam, nextracol, irt - character(len=255) :: filename, headerline, fmtstr - real(wp_), dimension(5) :: f48v - real(wp_) :: jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - allocate(jphi(params%output%nrho), currins(params%output%nrho), & - pins(params%output%nrho), rtin(params%output%nrho), & - rpin(params%output%nrho)) - - open(100, file=opts%sum_filelist, action='read', status='old', iostat=err) - if (err /= 0) then - call log_message(level=ERROR, mod='main', & - msg='opening file list ('//opts%sum_filelist//') failed!') - call exit(1) - end if - - open(unit=100 -1, file='f48sum.txt', action='write', status='unknown') - open(unit=100 -2, file='f7sum.txt', action='write', status='unknown') - - read(100, *) n, ngam, nextracol - allocate(extracol(nextracol)) - do i=1,n - read(100, *) filename - open(100 + i, file=filename, action='read', status='old', iostat=err) - if (err /= 0) then - call log_message(level=ERROR, mod='main', & - msg='opening summand file ('//trim(filename)//') failed!') - call exit(1) - end if - - do j=1,22 - read(100 + i, '(a255)') headerline - if (i == 1) then - write(100 -1, '(a)') trim(headerline) - write(100 -2, '(a)') trim(headerline) - end if - end do - end do - - allocate(results%jcd(params%output%nrho), results%dpdv(params%output%nrho)) - do k=1,ngam - jphi = zero - results%jcd = zero - results%dpdv = zero - currins = zero - pins = zero - do j=1,params%output%nrho - do i=1,n - read(100+i, *) (extracol(l), l=1,nextracol), & - rpin(j), rtin(j), f48v(1:5), irt - jphi(j) = f48v(1) + jphi(j) - results%jcd(j) = f48v(2) + results%jcd(j) - results%dpdv(j) = f48v(3) + results%dpdv(j) - currins(j) = f48v(4) + currins(j) - pins(j) = f48v(5) + pins(j) - end do - write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7 - write(100 -1,trim(fmtstr)) & - (extracol(l), l=1,nextracol), rpin(j), rtin(j), & - jphi(j), results%jcd(j), results%dpdv(j), & - currins(j), pins(j), irt - end do - results%pabs = pins(params%output%nrho) - results%icd = currins(params%output%nrho) - write(100 -1, *) - call sum_profiles(params, jphi, results%jcd, results%dpdv, & - currins, pins, results%pabs, results%icd, & - jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & - drhotjava, drhotpav, ratjamx, ratjbmx) - write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12 - write(100 -2, trim(fmtstr)) & - (extracol(l), l=1,nextracol), results%icd, results%pabs, & - jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & - drhotjava, drhotpav, ratjamx, ratjbmx - end do - do i=-2,n - close(100 + i) - end do - deallocate(jphi, currins, pins, rtin, rpin) - deallocate(extracol) - deallocate(opts%params_file) - end block sum + call sum_profiles(params, results) else call gray_main(params, data, results, err) end if @@ -192,16 +105,40 @@ program main call log_message(msg, level=INFO, mod='main') end block print_res + write_res: block + integer :: i + + ! Write results to file + do i = 1, size(results%tables%iter) + if (.not. associated(results%tables%iter(i)%ptr)) cycle + + associate (tbl => results%tables%iter(i)%ptr) + call tbl%save(err, make_header=make_header) + if (err /= 0) then + call log_message('failed to save '//trim(tbl%title)//' table', & + level=ERROR, mod='main') + end if + end associate + end do + end block write_res + ! Free memory call deinit_equilibrium(data%equilibrium) call deinit_profiles(data%profiles) call deinit_misc call deinit_cli_options(opts) deallocate(results%dpdv, results%jcd) - call close_units contains + subroutine make_header(file) + ! Writes the gray output files header + use gray_params, only : write_gray_header + integer, intent(in) :: file + call write_gray_header(params, file) + end subroutine + + subroutine init_equilibrium(params, data, err) ! Reads the MHD equilibrium file (either in the G-EQDSK format ! or an analytical description) and initialises the respective @@ -462,83 +399,120 @@ contains end subroutine deinit_misc - subroutine sum_profiles(params, jphi, jcd, dpdv, currins, pins, pabs, icd, & - jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & - drhotjava, drhotpav, ratjamx, ratjbmx) - use const_and_precisions, only : zero, degree - use coreprofiles, only : set_profiles_an, set_profiles_spline, & - temp, fzeff - use dispersion, only : expinit - use gray_params, only : gray_parameters, print_parameters - use beams, only : launchangles2n, xgygcoeff - use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr - use pec, only : pec_init, postproc_profiles, dealloc_pec, & - rhop_tab, rhot_tab - use gray_core, only : print_headers, print_finals, print_pec, & - print_bres, print_prof, print_maps, & - print_surfq + subroutine sum_profiles(params, results) + use gray_params, only : gray_parameters + use magsurf_data, only : compute_flux_averages, dealloc_surfvec + use beamdata, only : init_btr + use pec, only : pec_init, postproc_profiles, dealloc_pec, & + rhot_tab ! subroutine arguments - type(gray_parameters), intent(inout) :: params - real(wp_), intent(in) :: pabs, icd - real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins - real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava, & - rhotp, rhotpav, drhotjava, drhotpav, & - ratjamx, ratjbmx + type(gray_parameters), intent(inout) :: params + type(gray_results), intent(inout) :: results + + real(wp_), dimension(:), allocatable :: jphi + real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin + real(wp_), dimension(:), allocatable :: extracol + integer :: i, j, k, l, n, ngam, nextracol, irt + character(len=255) :: filename, headerline, fmtstr + real(wp_), dimension(5) :: f48v + real(wp_) :: jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & + drhotjava, drhotpav, ratjamx, ratjbmx ! local variables - real(wp_) :: ak0, bres, xgcn - real(wp_) :: chipol, psipol, st real(wp_) :: drhotp, drhotj, dpdvmx, jphimx - ! ======== set environment BEGIN ======== - ! Compute X=(ω_pe/ω)² and Y=ω_ce/ω (with B=1) - call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) - ! Initialise the ray variables (beamtracing) call init_btr(params%raytracing) - ! Initialise the dispersion module - if (params%ecrh_cd%iwarm > 1) call expinit - ! Initialise the magsurf_data module - call flux_average ! requires frhotor for dadrhot,dvdrhot + call compute_flux_averages ! Initialise the output profiles call pec_init(params%output%ipec) - ! ======= set environment END ====== - ! ======== pre-proc prints BEGIN ======== - call print_headers(params) + associate(nrho =>params%output%nrho) + allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho)) + allocate(results%jcd(nrho), results%dpdv(nrho)) + end associate - ! Print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout - call print_surfq([1.5_wp_, 2.0_wp_]) + open(100, file=opts%sum_filelist, action='read', status='old', iostat=err) + if (err /= 0) then + call log_message(level=ERROR, mod='main', & + msg='opening file list ('//opts%sum_filelist//') failed!') + call exit(1) + end if - ! Print ne, Te, q, Jphi versus psi, rhop, rhot - call print_bres(bres) - call print_prof(params%profiles) - call print_maps(bres, xgcn, & - norm2(params%antenna%pos(1:2)) *0.01_wp_ , & - sin(params%antenna%beta*degree)) - ! ========= pre-proc prints END ========= + open(unit=100 -1, file='f48sum.txt', action='write', status='unknown') + open(unit=100 -2, file='f7sum.txt', action='write', status='unknown') - ! Print power and current density profiles - call print_pec(rhop_tab, rhot_tab, jphi, jcd, & - dpdv, currins, pins, index_rt=1) - ! Compute profiles width - call postproc_profiles(pabs, icd, rhot_tab, dpdv, jphi, & - rhotpav, drhotpav, rhotjava, drhotjava, & - dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, & - dpdvmx, jphimx, ratjamx, ratjbmx) - ! Print 0D results - call print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, & - drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, & - drhotj, ratjamx, ratjbmx, st, psipol, chipol, & - 1, params%antenna%power, cpl1=zero, cpl2=zero) + read(100, *) n, ngam, nextracol + allocate(extracol(nextracol)) + do i=1,n + read(100, *) filename + open(100 + i, file=filename, action='read', status='old', iostat=err) + if (err /= 0) then + call log_message(level=ERROR, mod='main', & + msg='opening summand file ('//trim(filename)//') failed!') + call exit(1) + end if + + do j=1,22 + read(100 + i, '(a255)') headerline + if (i == 1) then + write(100 -1, '(a)') trim(headerline) + write(100 -2, '(a)') trim(headerline) + end if + end do + end do + + do k=1,ngam + jphi = 0 + results%jcd = 0 + results%dpdv = 0 + currins = 0 + pins = 0 + do j=1,params%output%nrho + do i=1,n + read(100+i, *) (extracol(l), l=1,nextracol), & + rpin(j), rtin(j), f48v(1:5), irt + jphi(j) = f48v(1) + jphi(j) + results%jcd(j) = f48v(2) + results%jcd(j) + results%dpdv(j) = f48v(3) + results%dpdv(j) + currins(j) = f48v(4) + currins(j) + pins(j) = f48v(5) + pins(j) + end do + write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7 + write(100 -1,trim(fmtstr)) & + (extracol(l), l=1,nextracol), rpin(j), rtin(j), & + jphi(j), results%jcd(j), results%dpdv(j), & + currins(j), pins(j), irt + end do + results%pabs = pins(params%output%nrho) + results%icd = currins(params%output%nrho) + write(100 -1, *) + ! Recompute final results + call postproc_profiles(results%pabs, results%icd, rhot_tab, & + results%dpdv, jphi, rhotpav, drhotpav, & + rhotjava, drhotjava, dpdvp, jphip, & + rhotp, drhotp, rhotj, drhotj, & + dpdvmx, jphimx, ratjamx, ratjbmx) + write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12 + write(100 -2, trim(fmtstr)) & + (extracol(l), l=1,nextracol), results%icd, results%pabs, & + jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, & + drhotjava, drhotpav, ratjamx, ratjbmx + end do + + do i=-2,n + close(100 + i) + end do + deallocate(jphi, currins, pins, rtin, rpin) + deallocate(extracol) + deallocate(opts%params_file) ! Free memory - call dealloc_surfvec ! for fluxval + call dealloc_surfvec call dealloc_pec end subroutine sum_profiles diff --git a/src/types.f90 b/src/types.f90 new file mode 100644 index 0000000..ea53fc9 --- /dev/null +++ b/src/types.f90 @@ -0,0 +1,354 @@ +module types + + use const_and_precisions, only : wp_ + + implicit none + + type item + ! An item is the building block of a queue + class(*), allocatable :: value ! the actual content + type(item), pointer :: next => null() ! pointer to the next item + end type + + type queue + ! A queue is a list of items with O(1) insertion and extraction. + ! The first item inserted (`put` operation) is the first to be + ! extracted (`get` operation). + + ! References to the first and last items in the queue + type(item), pointer :: first => null(), last => null() + + contains + procedure :: put => queue_put ! inserts an item at the end of the queue + procedure :: get => queue_get ! extracts the first item in the queue + procedure :: empty => queue_empty ! checks whether the queue is empty + end type + + type table + ! A type for storing tabular data before serialisation + integer :: id ! the table unique ID + logical :: active ! whether to process the table + character(64) :: title ! title of the table + character(64), allocatable :: labels(:) ! labels for each column + type(queue), allocatable :: columns(:) ! values for each column + + contains + procedure :: init => table_init ! creates an empty table + procedure :: append => table_append ! appends an extra row to the table + procedure :: empty => table_empty ! checks whether the queue is empty + procedure :: save => table_save ! processes and writes the table to file + end type + + type wrap + ! A wrapper type for storing (references of) heterogeneous + ! values in an array or other homogeneuous container type + class(*), pointer :: value + end type + + ! Interface for custom type constructor + interface wrap + procedure :: wrap_init + end interface + + ! Interface for subroutines that write to file + interface + subroutine writer(file) + integer, intent(in) :: file + end subroutine + end interface + +contains + + subroutine queue_put(self, val) + ! Inserts an item of value `val` at the end of the `self` queue + class(queue), intent(inout) :: self + class(*), intent(in) :: val + + type(item), pointer :: i + + ! set up the item and copy the value + allocate(i) + allocate(i%value, source=val) + + if (self%empty()) then + ! first item received, store it separately + self%first => i + else + ! otherwise, link the last item to the new one + self%last%next => i + end if + + ! update the last item + self%last => i + end subroutine queue_put + + + subroutine queue_get(self, val) + ! Extracts the item from the front of the queue. + ! Notes: + ! - this assumes the queue is not empty. + ! - this should be converted to a function once + ! gcc bug 115072 is fixed. + class(queue), intent(inout) :: self + class(*), allocatable, intent(out) :: val + + type(item), pointer :: tmp + + ! extract the value + allocate(val, source=self%first%value) + + tmp => self%first ! get a grip + self%first => tmp%next ! shift the reference + deallocate(tmp) ! clear the item + end subroutine queue_get + + + pure function queue_empty(self) + ! Whether the `self` queue is empty + class(queue), intent(in) :: self + logical :: queue_empty + + queue_empty = .not. associated(self%first) + end function + + + function wrap_init(val) + ! Initialise a `wrap` value + ! + ! Note: this is needed because, apparently, gfortran can't handle + ! unlimited polymorphic arguments in the default constructor. + class(*), target, intent(in) :: val + type(wrap) :: wrap_init + wrap_init%value => val + end function wrap_init + + + subroutine table_init(self, id, title, labels, active) + ! Initialises the table + class(table), intent(out) :: self + integer, intent(in) :: id + character(*), intent(in) :: title + character(*), intent(in) :: labels(:) + + logical, intent(in), optional :: active + + integer :: i + + self%id = id + self%title = title + self%labels = labels + self%columns = [(queue(), i=1,size(labels))] + + if (present(active)) self%active = active + end subroutine table_init + + + pure function table_empty(self) + ! Checks whether the table is empty + class(table), intent(in) :: self + logical :: table_empty + + table_empty = self%columns(1)%empty() + end function table_empty + + + subroutine table_append(self, row) + ! Appends an extra row to the table + ! + ! Heterogenous values are supported, provided they + ! are wrapped in a `wrap` type. + + class(table), intent(inout) :: self + class(*), intent(in) :: row(:) + + integer :: i + + do i = 1, size(self%columns) + select type (row) + type is (wrap) + call self%columns(i)%put(row(i)%value) + class default + call self%columns(i)%put(row(i)) + end select + end do + + end subroutine table_append + + + subroutine table_save(self, error, filepath, make_header) + ! Save the table to a file + ! + ! Note: this operation consumes the table + use utils, only : get_free_unit, digits + + ! suborutine arguments + class(table), intent(inout) :: self + integer, intent(out) :: error + character(*), intent(in), optional :: filepath + procedure(writer), optional :: make_header + + ! local variables + integer :: i, file + class(*), allocatable :: val + character(len=len(self%title) + 4) :: filepath_ + character(len=10) :: fmt + + if (self%empty()) return + + file = get_free_unit() + + ! set default file name + if (present(filepath)) then + filepath_ = filepath + else + write(filepath_, '(a,".",i0,".txt")') trim(self%title), self%id + end if + + ! open output file + open(unit=file, file=filepath_, action='write', iostat=error) + if (error /= 0) return + + ! write a custom file header + if (present(make_header)) call make_header(file) + + ! write table header + write(file, '(a)', advance='no') '#' + do i = 1, size(self%columns) + write(file, '(x,a)', advance='no') trim(self%labels(i)) + end do + write(file, '(x)') + + ! iterate the table + do while (.not. self%empty()) + + ! write a row + do i = 1, size(self%columns) + call self%columns(i)%get(val) + + ! format a record + select type (val) + type is (integer) + write(fmt, '("(i",i0,")")') max(5, 1 + digits(val)) + write(file, fmt, advance='no') val + type is (real(wp_)) + write(file, '(es16.8e3)', advance='no') val + type is (logical) + write(file, '(l1)', advance='no') val + type is (character(len=*)) + write(file, '(a)', advance='no') trim(val) + end select + + ! add space between values + if (i < size(self%columns)) write(file, '(x)', advance='no') + end do + + ! start newline + write(file, '(x)') + end do + + end subroutine table_save + + + subroutine test_queue() + integer, target :: x=1, y=2 + real, target :: z = 1.231 + logical, target :: p=.false. + character(len=4) :: s='ciao' + type(queue) :: q + class(*), allocatable :: c + + call assert(q%empty(), 'queue is empty') + + ! put some items in + call q%put(x) + call q%put(y) + call q%put(z) + call q%put(p) + call q%put(s) + + call assert(.not. q%empty(), 'queue is not empty') + + ! get them out and check them + call q%get(c) + call assert(transfer(c, x) == x, 'got back an integer') + call q%get(c) + call assert(transfer(c, y) == y, 'got back another integer') + call q%get(c) + call assert(transfer(c, z) == z, 'got back a real') + call q%get(c) + call assert(transfer(c, p) .eqv. p, 'got back a logical') + call q%get(c) + call assert(cast_char(c, 4) == s, 'got back a string') + + call assert(q%empty(), 'queue is empty again') + + ! this should be finalised + call q%put(s) + + contains + + subroutine assert(p, msg) + logical, intent(in) :: p + character(len=*), intent(in) :: msg + + if (.not. p) then + print '(a, a)', 'assertion failed: ', msg + stop + end if + end subroutine + + pure function cast_char(x, n) result(y) + class(*), intent(in) :: x + integer, intent(in) :: n + character(len=n) :: y + + select type (x) + type is (character(len=*)) + y = x + class default + error stop 'invalid cast to character' + end select + end function + + end subroutine test_queue + + + subroutine test_table() + type(table) :: tbl + integer :: i + class(*), allocatable :: x, y + + call tbl%init(title='testing table', id=0, labels=['x', 'y']) + + ! fill in some values + call assert(tbl%empty(), 'table is empty') + do i = 1, 4 + call tbl%append([0.2*i, (0.2*i)**2]) + end do + call assert(.not. tbl%empty(), 'table is not empty') + + ! get them back + do i = 1, 4 + call tbl%columns(1)%get(x) + call tbl%columns(2)%get(y) + call assert(transfer(x, 1.0) == (0.2*i), 'got correct x') + call assert(transfer(y, 1.0) == (0.2*i)**2, 'got correct y') + end do + + call assert(tbl%empty(), 'table is empty again') + + contains + + subroutine assert(p, msg) + logical, intent(in) :: p + character(len=*), intent(in) :: msg + + if (.not. p) then + print '(a, a)', 'assertion failed: ', msg + stop + end if + end subroutine + + end subroutine test_table + +end module types diff --git a/src/units.f90 b/src/units.f90 deleted file mode 100644 index 1fffa55..0000000 --- a/src/units.f90 +++ /dev/null @@ -1,70 +0,0 @@ -! This modules contains the definitions of the unit numbers used for -! all the GRAY output ifles and a mechanism to toggle them -module units - - implicit none - - ! Unit numbers -#ifdef JINTRAC - ! JINTRAC - integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644 - integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631 - integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633 - integer, parameter :: udisp =617, upec =648, usumm =607, umaps =632 - integer, parameter :: uhead =638 -#else - ! STANDARD - integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99 - integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71 - integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33 - integer, parameter :: udisp = 17, upec = 48, usumm = 7, umaps = 72 - integer, parameter :: uhead = 78 -#endif - - ! List of active units - logical, save :: all_enabled = .false. - integer, allocatable, save :: active_units(:) - -contains - - subroutine set_active_units(units) - ! Activate the given units - ! - ! All units are inactive by default and no output - ! will be directed to them. - - ! subroutine arguments - integer, intent(in) :: units(:) - - active_units = units - - end subroutine set_active_units - - - subroutine close_units - ! Close all the active units to flush the buffer. - - integer :: i - - if (allocated(active_units)) then - do i = 1, size(active_units) - close(active_units(i)) - end do - deallocate(active_units) - end if - - end subroutine close_units - - - function unit_active(unit) result(on) - ! Checks whether the given `unit` is active - - ! function arguments - integer, intent(in) :: unit - logical :: on - - on = all_enabled .or. any(active_units == unit) - - end function unit_active - -end module units diff --git a/src/utils.f90 b/src/utils.f90 index 7bc75aa..e56d1c8 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -378,4 +378,20 @@ contains isrelative = (filepath(1:1) /= '/') end function isrelative + + pure function digits(n) + ! Returns the number of digits of an integer + + ! function arguments + integer, intent(in) :: n + integer :: digits + + ! How it works: + ! 1. `bit_size` returns the number of bits + 1 (sign) + ! 2. subtracting the number of leading zeros gives the significant bits + ! 3. multiplying by log₁₀(2) gives the approximate number of digits + ! 4. `ceiling` rounds to the next closest integer + digits = ceiling((bit_size(abs(n)) - leadz(abs(n))) * log10(2.0_wp_)) + end function digits + end module utils diff --git a/src/vendor/cniteq.f90 b/src/vendor/cniteq.f90 new file mode 100644 index 0000000..a79b478 --- /dev/null +++ b/src/vendor/cniteq.f90 @@ -0,0 +1,218 @@ +module cniteq + + use const_and_precisions, only : wp_ + + implicit none + +contains + + subroutine cniteq_f(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon) +! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. +! (based on an older code) +! arguments + integer, intent(in) :: nr,nz + real(wp_), dimension(nr), intent(in) :: rqgrid + real(wp_), dimension(nz), intent(in) :: zqgrid + real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid + real(wp_), intent(in) :: h + integer, intent(inout) :: ncon, icount + integer, dimension(ncon), intent(out) :: npts + real(wp_), dimension(icount), intent(out) :: rcon,zcon +! local variables + integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb + integer :: jabs,jnb,kx,ikx,itm,inext,in + integer, dimension(3,2) :: ja + integer, dimension(icount/2-1) :: lx + real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y + real(wp_), dimension(nr*nz) :: a + logical :: flag1 + + px = 0.5_wp_ + + a = reshape(matr2dgrid, [nr*nz]) + + rcon = 0.0_wp_ + zcon = 0.0_wp_ + + nrqmax = nr + drgrd = rqgrid(2) - rqgrid(1) + dzgrd = zqgrid(2) - zqgrid(1) + + ncon = 0 + + npts = 0 + + iclast = 0 + icount = 0 + mpl = 0 + ix = 0 + mxr = nrqmax * (nz - 1) + n1 = nr - 1 + + do jx=2,n1 + do jm=jx,mxr,nrqmax + j = jm + nrqmax + ah=a(j)-h + if (ah <= 0.0_wp_ .and. a(jm) > h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + end do + + do jm=nr,mxr,nrqmax + j = jm + nrqmax + ah=a(j)-h + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + if (ah <= 0.0_wp_ .and. a(jm) > h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + end do + + do jm=1,mxr,nrqmax + j = jm + nrqmax + if (a(j) <= h .and. a(jm) > h .or. & + a(j) > h .and. a(jm) <= h) then + ix=ix+1 + lx(ix) =-j + end if + end do + + do j=2,nr + if (a(j) <= h .and. a(j-1) > h .or. & + a(j) > h .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + + if(ix<=0) return + +bb: do + in=ix + jx=lx(in) + jfor=0 + lda=1 + ldb=2 + + do + if(jx<0) then + jabs=-jx + jnb = jabs - nrqmax + else + jabs=jx + jnb=jabs-1 + end if + + adn=a(jabs)-a(jnb) + if(adn/=0) px=(a(jabs)-h)/adn + kx = (jabs - 1) / nrqmax + ikx = jabs - nrqmax * kx - 1 + + if(jx<0) then + x = drgrd * ikx + y = dzgrd * (kx - px) + else + x = drgrd * (ikx - px) + y = dzgrd * kx + end if + + icount = icount + 1 + rcon(icount) = x + rqgrid(1) + zcon(icount) = y + zqgrid(1) + mpl= icount + itm = 1 + ja(1,1) = jabs + nrqmax + j=1 + + if(jx<=0) then + ja(1,1) = -jabs-1 + j=2 + end if + + ja(2,1) = -ja(1,1) + ja(3,1) = -jx + 1 - nrqmax + ja(3,2) = -jx + ja(j,2) = jabs - nrqmax + k= 3-j + ja(k,2) = 1-jabs + + if (kx<=0 .or. ikx<=0) then + lda=1 + ldb=lda + else if (ikx + 1 - nr >= 0 .and. jx <= 0) then + lda=2 + ldb=lda + else if(jfor/=0) then + lda=2 + do i=1,3 + if(jfor==ja(i,2)) then + lda=1 + exit + end if + end do + ldb=lda + end if + + flag1=.false. + aa: do k=1,3 + do l=lda,ldb + do i=1,ix + if(lx(i)==ja(k,l)) then + itm=itm+1 + inext= i + if(jfor/=0) exit aa + if(itm > 3) then + flag1=.true. + exit aa + end if + end if + end do + end do + end do aa + + if(.not.flag1) then + lx(in)=0 + if(itm == 1) exit + end if + + jfor=jx + jx=lx(inext) + in = inext + end do + + do + if(lx(ix)/=0) then + if(mpl>=4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + exit + end if + ix= ix-1 + if(ix<=0) exit bb + end do + + end do bb + + if(mpl >= 4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + end subroutine cniteq_f + +end module