fix memory errors and use automatic deallocation

This commit is contained in:
Michele Guerini Rocco 2024-05-16 09:19:22 +02:00
parent f9c313323a
commit ba8fc001e5
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
5 changed files with 216 additions and 201 deletions

View File

@ -24,8 +24,7 @@ module gray_cli
private
public :: cli_options, print_cli_options, parse_cli_options, &
deinit_cli_options, parse_param_overrides, &
print_version, get_next_command
parse_param_overrides, print_version, get_next_command
contains
@ -285,18 +284,4 @@ contains
i = i + 1 ! increment counter
end subroutine
subroutine deinit_cli_options(opts)
! Frees all memory allocated by the parse_cli_options subroutine
! subroutine arguments
type(cli_options), intent(inout) :: opts
if (allocated(opts%output_dir)) deallocate(opts%output_dir)
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%tables)) deallocate(opts%tables)
end subroutine deinit_cli_options
end module gray_cli

View File

@ -6,7 +6,6 @@ module gray_params
implicit none
integer, parameter :: lenfnm = 256
integer, parameter :: headw = 132, headl = 21
! Polarisation mode
enum, bind(c)
@ -266,106 +265,149 @@ module gray_params
contains
subroutine write_gray_header(params, file)
function make_gray_header(params) result(header)
! Writes the Gray output header to a file
! subroutine arguments
! function arguments
type(gray_parameters), intent(in) :: params
integer :: file
! local variables
character(len=8) :: date
character(len=10) :: time
character(len=200) :: line
character(len=21*len(line)) :: header
character(len=8) :: date
character(len=10) :: time
#ifndef REVISION
character(len=*), parameter :: REVISION="unknown"
#endif
header = ""
! Date and time
call date_and_time(date, time)
write(file, '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
write(line, '("# 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)
call append(header, line)
! Git revision
write(file, '("# GRAY Git revision: ",a)') REVISION
write(line, '("# GRAY Git revision: ",a)') REVISION
call append(header, line)
! MHD equilibrium parameters
if (params%equilibrium%iequil > 0) then
write(file, '("# EQL input: ",a)') trim(params%equilibrium%filenm)
write(line, '("# EQL input: ",a)') trim(params%equilibrium%filenm)
call append(header, line)
! TODO add missing values
write(file, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
write(line, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
0, 0, 0, 0, 0
call append(header, line)
else
write(file, '("# EQL input: N/A (vacuum)")')
write(file, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
write(line, '("# EQL input: N/A (vacuum)")')
call append(header, line)
write(line, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
call append(header, line)
end if
write(file, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
write(line, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
params%equilibrium%iequil, params%equilibrium%sgnb, &
params%equilibrium%sgni, params%equilibrium%factb
call append(header, line)
if (params%equilibrium%iequil > 1) then
write(file, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
write(line, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
params%equilibrium%icocos, params%equilibrium%ipsinorm, &
params%equilibrium%idesc, params%equilibrium%ifreefmt
write(file, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
call append(header, line)
write(line, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
params%equilibrium%ssplps, params%equilibrium%ssplf, &
params%equilibrium%ixp
call append(header, line)
else
write(file, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
write(file, '("# EQL ssplps ssplf ixp: N/A (analytical)")')
write(line, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
call append(header, line)
write(line, '("# EQL ssplps ssplf ixp: N/A (analytical)")')
call append(header, line)
end if
! Profiles parameters
if (params%equilibrium%iequil > 0) then
write(file, '("# PRF input: ",a)') trim(params%profiles%filenm)
write(file, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
write(line, '("# PRF input: ",a)') trim(params%profiles%filenm)
call append(header, line)
write(line, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
params%profiles%iprof, params%profiles%iscal, &
params%profiles%factne, params%profiles%factte
call append(header, line)
if (params%profiles%iprof > 0) then
write(file, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
write(line, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
params%profiles%irho, params%profiles%psnbnd, params%profiles%sspld
call append(header, line)
else
write(file, '("# PRF irho psnbnd sspld: N/A (analytical)")')
write(line, '("# PRF irho psnbnd sspld: N/A (analytical)")')
call append(header, line)
end if
! TODO: add missing values
write(file, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0
write(line, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0
call append(header, line)
else
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)")')
write(line, '("# PRF input: N/A (vacuum)")')
call append(header, line)
write(line, '("# PRF iprof iscal factne factte: N/A (vacuum)")')
call append(header, line)
write(line, '("# PRF irho psnbnd sspld: N/A (vacuum)")')
call append(header, line)
write(line, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
call append(header, line)
end if
! Antenna parameters
write(file, '("# ANT input: ",a)') trim(params%antenna%filenm)
write(file, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
write(line, '("# ANT input: ",a)') trim(params%antenna%filenm)
call append(header, line)
write(line, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
params%antenna%ibeam, params%antenna%iox, &
params%antenna%psi, params%antenna%chi
write(file, '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
call append(header, line)
write(line, '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
params%antenna%alpha, params%antenna%beta, &
params%antenna%power, params%antenna%fghz
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))') &
call append(header, line)
write(line, '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
call append(header, line)
write(line, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
params%antenna%w, params%antenna%ri, params%antenna%phi
call append(header, line)
! Other parameters
write(file, '("# RFL rwall:",1x,e12.5)') params%misc%rwall
write(line, '("# RFL rwall:",1x,e12.5)') params%misc%rwall
call append(header, line)
! code parameters
write(file, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
write(line, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
params%raytracing%igrad, params%raytracing%idst, &
params%raytracing%ipass, params%raytracing%ipol
write(file, '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
call append(header, line)
write(line, '("# 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(file, '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
call append(header, line)
write(line, '("# 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(file, '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
call append(header, line)
write(line, '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
params%output%ipec, params%output%nrho, &
params%output%istpr, params%output%istpl
end subroutine write_gray_header
header = trim(header) // trim(line)
contains
subroutine append(header, line)
! Appends a line to the header
character(len=*), intent(inout) :: header
character(len=*), intent(in) :: line
header = trim(header) // trim(line) // new_line('a')
end subroutine
end function make_gray_header
function update_parameter(params, name, value) result(err)

View File

@ -3,13 +3,16 @@ program main
use logger, only : INFO, ERROR, set_log_level, log_message
use utils, only : dirname
use gray_cli, only : cli_options, parse_cli_options, &
deinit_cli_options, parse_param_overrides
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, &
make_gray_header, &
params_set_globals => set_globals
implicit none
no_save: block
! CLI options
type(cli_options) :: opts
@ -75,25 +78,17 @@ 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)
! Change the current directory to output files there
if (allocated(opts%output_dir)) then
if (chdir(opts%output_dir) /= 0) then
call log_message(level=ERROR, mod='main', &
msg='chdir to output_dir ('//opts%output_dir//') failed!')
call exit(1)
end if
end if
if (allocated(opts%sum_filelist)) then
call log_message(level=INFO, mod='main', msg='summing profiles')
call sum_profiles(params, results)
call sum_profiles(params, opts%sum_filelist, opts%output_dir, results)
else
! Activate the given output tables
params%misc%active_tables = opts%tables
! Finally run the simulation
call gray_main(params, data, results, err)
end if
@ -105,15 +100,23 @@ program main
call log_message(msg, level=INFO, mod='main')
end block print_res
! Change the current directory to output files there
if (allocated(opts%output_dir)) then
if (chdir(opts%output_dir) /= 0) then
call log_message(level=ERROR, mod='main', &
msg='chdir to output_dir ('//opts%output_dir//') failed!')
call exit(1)
end if
end if
! Write results to file
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)
call tbl%save(err, header=make_gray_header(params))
if (err /= 0) then
call log_message('failed to save '//trim(tbl%title)//' table', &
level=ERROR, mod='main')
@ -123,22 +126,20 @@ program main
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)
free_mem: block
use coreprofiles, only : unset_profiles_spline
use equilibrium, only : unset_equil_spline
use limiter, only : limiter_unset_globals=>unset_globals
call unset_profiles_spline
call unset_equil_spline
call limiter_unset_globals
end block free_mem
end block no_save
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
@ -191,25 +192,6 @@ contains
end subroutine init_equilibrium
subroutine deinit_equilibrium(data)
! Free all memory allocated by the init_equilibrium subroutine.
use gray_params, only : equilibrium_data
use equilibrium, only : unset_equil_spline
! subroutine arguments
type(equilibrium_data), intent(inout) :: data
! Free the MHD equilibrium arrays
if (allocated(data%rv)) deallocate(data%rv, data%zv, data%fpol, data%qpsi)
if (allocated(data%psin)) deallocate(data%psin, data%psinr)
if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd)
if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
! Unset global variables of the `equilibrium` module
call unset_equil_spline
end subroutine deinit_equilibrium
subroutine init_profiles(params, factb, launch_pos, data, err)
! Reads the plasma kinetic profiles file (containing the elecron
! temperature, density and plasma effective charge) and initialises
@ -272,23 +254,6 @@ contains
end subroutine init_profiles
subroutine deinit_profiles(data)
! Free all memory allocated by the init_profiles subroutine.
use gray_params, only : profiles_data
use coreprofiles, only : unset_profiles_spline
! subroutine arguments
type(profiles_data), intent(inout) :: data
! Free the plasma kinetic profiles arrays
if (allocated(data%psrad)) deallocate(data%psrad)
if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
! Unset global variables of the `coreprofiles` module
call unset_profiles_spline
end subroutine deinit_profiles
subroutine init_antenna(params, err)
! Reads the wave launcher file (containing the wave frequency, launcher
! position, direction and beam description) and initialises the respective
@ -390,40 +355,33 @@ contains
end subroutine init_misc
subroutine deinit_misc
! Free all memory allocated by the init_misc subroutine.
use limiter, only : limiter_unset_globals=>unset_globals
! Unset the global variables of the `limiter` module
call limiter_unset_globals
end subroutine deinit_misc
subroutine sum_profiles(params, results)
subroutine sum_profiles(params, filelist, output_dir, results)
! Combines the EC profiles obtained from multiple gray runs
! (of different beams) and recomputes the summary table
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
character(len=*), intent(in) :: filelist, output_dir
type(gray_results), intent(inout) :: results
! local variables
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
integer :: i, j, k, l, nbeams, ncases, nextracols, irt, err
character(len=255) :: filename, line, extralabels, fmtstr
real(wp_), dimension(5) :: sums
real(wp_) :: jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx
! local variables
real(wp_) :: drhotp, drhotj, dpdvmx, jphimx
! Initialise the ray variables (beamtracing)
call init_btr(params%raytracing)
! unit numbers
integer :: list, prof, summ
integer, allocatable :: beams(:)
! Initialise the magsurf_data module
call compute_flux_averages
@ -436,80 +394,113 @@ contains
allocate(results%jcd(nrho), results%dpdv(nrho))
end associate
open(100, file=opts%sum_filelist, action='read', status='old', iostat=err)
! Read the list of files to be summed
open(newunit=list, file=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!')
msg='opening file list ('//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')
! Open combined output files
open(newunit=prof, file=output_dir//'/'//'sum-ec-profiles.txt', action='write')
open(newunit=summ, file=output_dir//'/'//'sum-summary.txt', action='write')
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)
! Move to the filelist directory
! (the filepaths are assumed relative to it)
if (chdir(dirname(filelist)) /= 0) then
call log_message(level=ERROR, mod='main', &
msg='chdir to file list directory failed!')
call exit(1)
end if
! Read the filelist header
read(list, *) nbeams, ncases, nextracols
allocate(beams(nbeams), extracol(nextracols))
! Read the list of files
do i = 1, nbeams
read(list, *) filename
! Open the summand file
open(newunit=beams(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
! Copy the Gray header from the ec-profiles file
do j = 1, 21
read(beams(i), '(a255)') line
if (i == 1) write(prof, '(a)') trim(line)
if (i == 1) write(summ, '(a)') trim(line)
end do
end do
do k=1,ngam
jphi = 0
! Copy the column labels
! Note: for the summary.7 we copy only the first columns
read(beams(i), '(a255)') line
if (i == 1) then
extralabels = line(1:index(line, "ρ_p")-1)
write(prof, '(a)') trim(line)
write(summ, '(a,x,a,x,a)') trim(extralabels), &
"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"
end if
end do
close(list)
do k = 1, ncases
jphi = 0
currins = 0
pins = 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)
! Sum the radial profiles
do j = 1, params%output%nrho
do i = 1, nbeams
read(beams(i), *) &
(extracol(l), l=1,nextracols), rpin(j), rtin(j), sums(1:5), irt
jphi(j) = sums(1) + jphi(j)
results%jcd(j) = sums(2) + results%jcd(j)
results%dpdv(j) = sums(3) + results%dpdv(j)
currins(j) = sums(4) + currins(j)
pins(j) = sums(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
! Store the sums
write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracols + 7
write(prof, trim(fmtstr)) &
(extracol(l), l=1,nextracols), rpin(j), rtin(j), jphi(j), &
results%jcd(j), results%dpdv(j), currins(j), pins(j), irt
end do
! Combine the total power, current
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
! Recompute the summary values
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))")') nextracols+12
write(summ, trim(fmtstr)) &
(extracol(l), l=1,nextracols), results%icd, results%pabs, &
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, drhotjava, &
drhotpav, ratjamx, ratjbmx
! Separate the case
write(prof, *)
end do
do i=-2,n
close(100 + i)
! Close all files
close(prof); close(summ)
do i = 1, nbeams
close(beams(i))
end do
deallocate(jphi, currins, pins, rtin, rpin)
deallocate(extracol)
deallocate(opts%params_file)
! Free memory
call dealloc_surfvec

View File

@ -22,6 +22,8 @@ program gray_convert
integer(kind(config_format)) :: format ! output file format
end type
no_save: block
type(cli_options) :: opts
type(gray_parameters) :: params
integer :: err
@ -47,6 +49,8 @@ program gray_convert
call exit(0)
end block no_save
contains
subroutine print_help()

View File

@ -50,13 +50,6 @@ module types
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)
@ -175,17 +168,17 @@ contains
end subroutine table_append
subroutine table_save(self, error, filepath, make_header)
subroutine table_save(self, error, filepath, 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
class(table), intent(inout) :: self
integer, intent(out) :: error
character(*), intent(in), optional :: filepath
character(*), optional :: header
! local variables
integer :: i, file
@ -209,7 +202,7 @@ contains
if (error /= 0) return
! write a custom file header
if (present(make_header)) call make_header(file)
if (present(header)) write(file, '(a)') trim(header)
! write table header
write(file, '(a)', advance='no') '#'