istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/gray_core.f90: fix possible division by zero

Browse Source 
If the EQDSK grid extends all the way to R=0, evaluating B ~ 1/R on a
regular grid from rmnm to rmxm results in a division by zero.
This commit is contained in:
Michele Guerini Rocco 2024-01-30 12:28:06 +01:00
parent 3bc1efc2a6
commit 2d16617db8
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 13 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/gray_core.f90
View File

@ -2207,6 +2207,7 @@ bb: do
subroutine print_bres(bres) subroutine print_bres(bres)
! Prints the EC resonance surface table (unit 70) ! Prints the EC resonance surface table (unit 70)
use const_and_precisions, only : comp_eps
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield
use units, only : ubres, unit_active use units, only : ubres, unit_active
use magsurf_data, only : npsi use magsurf_data, only : npsi
@ -2228,10 +2229,10 @@ bb: do
if (.not. unit_active(ubres)) return if (.not. unit_active(ubres)) return
! Build a regular (R, z) grid ! Build a regular (R, z) grid
dr = (rmxm - rmnm)/(npsi - 1) dr = (rmxm - rmnm - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1) dz = (zmxm - zmnm)/(npsi - 1)
do j=1,npsi do j=1,npsi
rv(j) = rmnm + dr*(j - 1) rv(j) = comp_eps + rmnm + dr*(j - 1)
zv(j) = zmnm + dz*(j - 1) zv(j) = zmnm + dz*(j - 1)
end do end do
@ -2271,6 +2272,7 @@ bb: do
subroutine print_maps(B_res, xgcn, R0, Npl0) subroutine print_maps(B_res, xgcn, R0, Npl0)
! Prints several input quantities on a regular grid (unit 72) ! 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 equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
use coreprofiles, only : density, temp use coreprofiles, only : density, temp
use units, only : umaps, unit_active use units, only : umaps, unit_active
@ -2291,10 +2293,10 @@ bb: do
if (.not. unit_active(umaps)) return if (.not. unit_active(umaps)) return
! Build a regular (R, z) grid ! Build a regular (R, z) grid
dR = (rmxm - rmnm)/(npsi - 1) dR = (rmxm - rmnm - comp_eps)/(npsi - 1)
dz = (zmxm - zmnm)/(npsi - 1) dz = (zmxm - zmnm)/(npsi - 1)
do j = 1, npsi do j = 1, npsi
R(j) = rmnm + dR*(j - 1) R(j) = comp_eps + rmnm + dR*(j - 1)
z(j) = zmnm + dz*(j - 1) z(j) = zmnm + dz*(j - 1)
end do end do